https://github.com/geodynamics/citcoms
Tip revision: c59234914a80773a2b9b69b63279df4dddff34ad authored by Eh Tan on 26 April 2007, 01:19:04 UTC
Move branches/CitcomS/ to tags/pre-2.0 since no development will be done in this branch.
Move branches/CitcomS/ to tags/pre-2.0 since no development will be done in this branch.
Tip revision: c592349
Drive_solvers.c
#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
void general_stokes_solver(E)
struct All_variables *E;
{
void construct_stiffness_B_matrix();
void velocities_conform_bcs();
void assemble_forces();
void sphere_harmonics_layer();
double global_vdot(),kineticE_radial();
float global_fvdot();
float vnorm_nonnewt();
void get_system_viscosity();
float vmag;
float *delta_U[NCS];
double Udot_mag, dUdot_mag;
double CPU_time0(),time;
int m,count,i,j,k;
static float *oldU[NCS];
static int visits=0;
const int nno = E->lmesh.nno;
const int nel = E->lmesh.nel;
const int nnov = E->lmesh.nnov;
const int neq = E->lmesh.neq;
const int vpts = vpoints[E->mesh.nsd];
const int dims = E->mesh.nsd;
const int addi_dof = additional_dof[dims];
if(visits==0) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
oldU[m] = (float *)malloc((neq+2)*sizeof(float));
for(i=0;i<=neq;i++)
oldU[m][i]=0.0;
}
visits ++;
}
for (m=1;m<=E->sphere.caps_per_proc;m++) {
delta_U[m] = (float *)malloc((neq+2)*sizeof(float));
}
/* FIRST store the old velocity field */
E->monitor.elapsed_time_vsoln = E->monitor.elapsed_time;
if(E->parallel.me==0) time=CPU_time0();
velocities_conform_bcs(E,E->U);
assemble_forces(E,0);
Udot_mag=dUdot_mag=0.0;
count=1;
do {
if(E->viscosity.update_allowed)
get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
solve_constrained_flow_iterative(E);
/* Udot_mag = kineticE_radial(E,E->U,E->mesh.levmax);
if(E->parallel.me==0)
fprintf(E->fp_out,"%g %g \n",E->monitor.elapsed_time,Udot_mag);
fflush(E->fp_out);
if(E->parallel.me==0)
fprintf(stderr,"kinetic energy= %g time4= %g seconds \n",Udot_mag,CPU_time0()-time);
*/
if ( E->viscosity.SDEPV ) {
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=0;i<neq;i++) {
delta_U[m][i] = E->U[m][i] - oldU[m][i];
oldU[m][i] = E->U[m][i];
}
Udot_mag = sqrt(global_fvdot(E,oldU,oldU,E->mesh.levmax));
dUdot_mag = vnorm_nonnewt(E,delta_U,oldU,E->mesh.levmax);
if(E->parallel.me==0){
fprintf(stderr,"Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n",dUdot_mag,Udot_mag,count);
fprintf(E->fp,"Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n",dUdot_mag,Udot_mag,count);
fflush(E->fp);
}
count++;
} /* end for SDEPV */
} while((count < 50) && (dUdot_mag>E->viscosity.sdepv_misfit) && E->viscosity.SDEPV);
for (m=1;m<=E->sphere.caps_per_proc;m++) {
free((void *) delta_U[m]);
}
return;
}