[cig-commits] r19414 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Sun Jan 22 22:11:00 PST 2012
Author: becker
Date: 2012-01-22 22:11:00 -0800 (Sun, 22 Jan 2012)
New Revision: 19414
Modified:
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
Log:
Further improved restart (one hopes), now relying on T, V (and P)
solutions.
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -60,6 +60,11 @@
iterate = need_to_iterate(E);
if(visits == 0){ /* initialization stage */
+ if(E->control.restart){ /* we have read in a velocity
+ solution, use this for the force
+ vector initial solution guess */
+ vector_from_v(E, E->U,E->V);
+ }
if(iterate){
/* damping factors */
alpha = E->viscosity.sdepv_iter_damp;
@@ -73,8 +78,13 @@
}
/* allocate oldU only if iterations are needed */
oldU = (double *)malloc(neq * sizeof(double));
- for(i = 0; i < neq; i++)
- oldU[i] = 0.0;
+ if(!E->control.restart){
+ for(i = 0; i < neq; i++)
+ oldU[i] = 0.0;
+ }else{
+ for(i = 0; i < neq; i++)
+ oldU[i] = E->U[i];
+ }
} /* end iterate */
visits++;
}
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -455,6 +455,8 @@
E->control.composition = 0;
+ E->sphere.vtk_base_init = 0;
+
E->control.GRID_TYPE = 1;
E->mesh.hwidth[1] = E->mesh.hwidth[2] = E->mesh.hwidth[3] = 1.0; /* divide by this one ! */
E->mesh.magnitude[1] = E->mesh.magnitude[2] = E->mesh.magnitude[3] = 0.0;
@@ -1171,14 +1173,19 @@
void common_initial_fields(struct All_variables *E)
{
- report(E, "Initialize pressure field");
- initial_pressure(E);
- report(E, "Initialize velocity field");
- initial_velocity(E);
- report(E, "Initialize viscosity field");
- get_viscosity_option(E);
- return;
+ if(E->control.restart){
+ if(E->parallel.me == 0)fprintf(stderr,"skipping regular P V init because of restart\n");
+ }else{
+ report(E, "Initialize pressure field");
+ initial_pressure(E);
+ report(E, "Initialize velocity field");
+ initial_velocity(E);
+ }
+ report(E, "Initialize viscosity field");
+ get_viscosity_option(E);
+
+ return;
}
/* ========================================== */
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -64,8 +64,8 @@
FILE *safe_fopen(char *,char *);
void *safe_malloc (size_t );
void calc_cbase_at_tp(float , float , float *);
-void convert_pvec_to_cvec(float ,float , float , float *,
- float *);
+void convert_pvec_to_cvec(float ,float , float , float *,float *);
+void convert_cvec_to_pvec(float ,float , float , float *,float *);
void rtp2xyz(float r, float , float , float *);
@@ -86,7 +86,6 @@
static float *SV, *EV;
static int been_here = 0;
float vs;
- static int vtk_base_init = 0; /* for spherical */
int vtk_comp_out;
const int dims = E->mesh.nsd;
@@ -288,17 +287,12 @@
gzprintf(gzout,"%.6e\n",E->T[i]);
gzclose(gzout);
/* velocities */
- if(E->control.Rsphere && (!vtk_base_init)){
- /*
-
- STILL NEED TO CHECK THE CONVENTION OF ORDERING FOR SPHERICAL HERE!
- */
-
+ if(E->control.Rsphere && (!E->sphere.vtk_base_init)){
/* need to init spherical/cartesian conversion vectors */
E->sphere.vtk_base = (float *)safe_malloc(sizeof(float)*E->lmesh.nno*9);
for(k=0,i=1;i<=E->lmesh.nno;i++,k+=9)
calc_cbase_at_tp(E->SX[1][i],E->SX[2][i],(E->sphere.vtk_base+k));
- vtk_base_init = 1;
+ E->sphere.vtk_base_init = 1;
}
/* write velocities to file */
sprintf(output_file,"%s/%d/vtk_v.%d.%d.gz",
@@ -691,15 +685,28 @@
int node, i, j, k, p;
FILE *fp;
- float temp1, temp2, *temp;
+ float temp1, temp2, *temp,cvec[3],pvec[3];
char input_s[200], input_file[255];
gzFile gzin;
+ if(E->control.Rsphere && (!E->sphere.vtk_base_init)){
+ /* need to init spherical/cartesian conversion vectors */
+ E->sphere.vtk_base = (float *)safe_malloc(sizeof(float)*E->lmesh.nno * 9);
+ for(k=0,i=1;i <= E->lmesh.nno;i++,k+=9)
+ calc_cbase_at_tp(E->SX[1][i],E->SX[2][i],(E->sphere.vtk_base+k));
+ E->sphere.vtk_base_init = 1;
+ }
+
temp = (float *)malloc((E->mesh.noz + 1) * sizeof(float));
if(E->control.restart == 1 || E->control.restart == 3)
{
-
+ if(E->parallel.me == 0)
+ fprintf(stderr,"restarting from %s timestep %i time %g\n",E->convection.old_T_file,E->monitor.solution_cycles,E->monitor.elapsed_time);
+
+ /*
+ temperatures
+ */
sprintf(input_file,"%s/%d/t.%d.%d.gz",
E->convection.old_T_file,
E->monitor.solution_cycles,E->parallel.me, E->monitor.solution_cycles);
@@ -718,15 +725,61 @@
{
gzgets (gzin,input_s, 200);
sscanf(input_s, "%f", &E->T[node]);
- //if(E->SX[3][node] == 0)fprintf(stderr,"%g %g\n",E->SX[3][node],E->T[node]);
+ //if(E->parallel.me==0)fprintf(stderr,"rt: %6i %11g\n",node,E->T[node]);
E->C[node] = 0;
}
gzclose(gzin);
- if(E->parallel.me == 0)
- fprintf(stderr,"restarting from %s timestep %i time %g\n",
- E->convection.old_T_file,E->monitor.solution_cycles,
- E->monitor.elapsed_time);
-
+ if(E->parallel.me == 0)fprintf(stderr,"restart T OK\n");
+ /*
+
+ velocities
+
+ */
+ sprintf(input_file,"%s/%d/vtk_v.%d.%d.gz",
+ E->convection.old_T_file,
+ E->monitor.solution_cycles,E->parallel.me, E->monitor.solution_cycles);
+ gzin = safe_gzopen(input_file, "r");
+ if(E->control.Rsphere){
+ for(node = 1,k=0; node <= E->lmesh.nno; node++,k+=9)
+ {
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%f %f %f", (cvec),(cvec+1),(cvec+2));
+ convert_cvec_to_pvec(cvec[0],cvec[1],cvec[2],(E->sphere.vtk_base+k),pvec);
+ E->V[3][i] = pvec[0];E->V[1][i]=pvec[1];E->V[2][i]=pvec[2];
+ }
+ }else{ /* cartesian */
+ for(i=1;i<=E->lmesh.nno;i++) {
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%f %f %f",&(E->V[1][i]),&(E->V[2][i]),&(E->V[3][i]));
+ //if(E->parallel.me==0)fprintf(stderr,"rv: %6i %11g %11g %11g\n",node,E->V[1][i],E->V[2][i],E->V[3][i]);
+ }
+ }
+ gzclose(gzin);
+ if(E->parallel.me == 0)fprintf(stderr,"restart V OK\n");
+ if(E->control.vtk_pressure_out){
+ /* pressures */
+ sprintf(input_file,"%s/%d/p.%d.%d.gz",
+ E->convection.old_T_file,
+ E->monitor.solution_cycles,E->parallel.me, E->monitor.solution_cycles);
+ gzin = safe_gzopen(input_file, "r");
+ if(gzgets (gzin,input_s, 200) == Z_NULL){
+ fprintf(stderr,"read error\n");
+ parallel_process_termination();
+ }
+ sscanf(input_s, "%i %i %f", &i, &j, &E->monitor.elapsed_time);
+ if(j != E->lmesh.nno)
+ myerror("mismatch of total node number upon restart",E);
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%i %i", &i, &j);
+ /* */
+ for(node = 1; node <= E->lmesh.nno; node++)
+ {
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%f", &E->NP[node]);
+ }
+ gzclose(gzin);
+ if(E->parallel.me == 0)fprintf(stderr,"restart P OK\n");
+ }
}
else
{
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -868,10 +868,7 @@
*/
void calc_cbase_at_tp(float theta, float phi, float *base)
{
-
-
- double ct,cp,st,sp;
-
+ double ct,cp,st,sp;
ct=cos(theta);
cp=cos(phi);
st=sin(theta);
@@ -889,6 +886,9 @@
base[7]= cp;
base[8]= 0.0;
}
+
+
+
/* given a base from calc_cbase_at_tp, convert a polar vector to
cartesian */
void convert_pvec_to_cvec(float vr,float vt, float vp, float *base,
@@ -901,6 +901,17 @@
cvec[i] += base[6+i]* vp;
}
}
+/* convert cartesian to polae */
+void convert_cvec_to_pvec(float vx,float vy, float vz, float *base,
+ float *pvec)
+{
+ int i,i3;
+ for(i=i3=0;i<3;i++,i3+=3){
+ pvec[i] = base[i3] * vx;
+ pvec[i] += base[i3+1]* vy;
+ pvec[i] += base[i3+2]* vz;
+ }
+}
/* convert r,theta,phi system to cartesian, xout[3] */
void rtp2xyz(float r, float theta, float phi, float *xout)
{
Modified: mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -550,3 +550,20 @@
}
return;
}
+/* for a restart, we might have a velocity solution, so assign this to
+ the F vector */
+void vector_from_v(struct All_variables *E, double *F, float **V)
+{
+ int node;
+
+ const int nno = E->lmesh.nno;
+
+ for(node = 1; node <= nno; node++)
+ {
+ F[E->id[node].doff[1]] = V[1][node];
+ F[E->id[node].doff[2]] = V[2][node];
+ F[E->id[node].doff[3]] = V[3][node];
+
+ }
+ return;
+}
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-01-23 06:11:00 UTC (rev 19414)
@@ -707,14 +707,14 @@
one = 1.0;
two = 2.0;
- if(visits == 0)
+ if((!E->control.restart) && (visits == 0))
{
for(e = 1; e <= nel; e++)
eedot[e] = one;
}
else
strain_rate_2_inv(E, eedot, 1);
- if((!E->viscosity.sdepv_start_from_newtonian)||(visits)){
+ if((!E->viscosity.sdepv_start_from_newtonian)||(E->control.restart)||(visits)){
switch(E->viscosity.sdepv_rheology){
case 1: /* old default, i think the factors
don't make sense, leave in for
@@ -1291,13 +1291,13 @@
myerror("Byerlee is to apply only to certain flavor, but no flavor is set",E);
}
}
- /*
- get strain rates for all elements
- */
+ }
+ if((!E->control.restart)&&(!visited)){
+ /* start with uniform strain-rates */
for(i=1;i<=nel;i++)
eedot[i] = 1.0;
-
}else{
+ /* either second time or restart with velocities */
if(E->viscosity.psrw)
strain_rate_2_inv(E,eedot,0);
else
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2012-01-23 06:11:00 UTC (rev 19414)
@@ -451,7 +451,8 @@
double dircos[4][4];
float *vtk_base;
-
+ int vtk_base_init;
+
int output_llmax;
int lnox;
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2012-01-21 17:43:49 UTC (rev 19413)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2012-01-23 06:11:00 UTC (rev 19414)
@@ -265,6 +265,7 @@
void *safe_malloc(size_t);
void calc_cbase_at_tp(float, float, float *);
void convert_pvec_to_cvec(float, float, float, float *, float *);
+void convert_cvec_to_pvec(float, float, float, float *, float *);
void rtp2xyz(float, float, float, float *);
void xyz2rtp(float, float, float, float *);
void myerror(char *, struct All_variables *);
@@ -372,6 +373,7 @@
float solve_Ahat_p_fhat(struct All_variables *, double *, double *, double *, double, int *);
void generate_log_message(int, double, double, double, double, struct All_variables *);
void v_from_vector(struct All_variables *, float **, double *);
+void vector_from_v(struct All_variables *, double *, float **);
/* Topo_gravity.c */
void get_CBF_topo(struct All_variables *, float *, float *);
void get_STD_topo(struct All_variables *, float *, float *, int);
More information about the CIG-COMMITS
mailing list