[cig-commits] r8101 - in mc/3D/CitcomS/trunk: bin lib
becker at geodynamics.org
becker at geodynamics.org
Wed Oct 10 16:40:39 PDT 2007
Author: becker
Date: 2007-10-10 16:40:38 -0700 (Wed, 10 Oct 2007)
New Revision: 8101
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/convection_variables.h
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Added coor=3 option for radial node spacing a la CitcomCU.
Streamlined ggrd temperature init.
Added file flush to heat flow output.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -161,7 +161,7 @@
E->control.keep_going = 0;
}
- if (E->monitor.T_interior>1.5) {
+ if (E->monitor.T_interior > E->monitor.T_interior_max_for_exit) {
fprintf(E->fp,"quit due to maxT = %.4e sub_iteration%d\n",E->monitor.T_interior,E->advection.last_sub_iterations);
parallel_process_termination();
}
Modified: mc/3D/CitcomS/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -749,31 +749,48 @@
}
/* took this apart to allow call from other subroutines */
+
+/*
+
+
+determine layer number based on radial coordinate r
+
+if E->viscosity.z... set to Earth values, then
+
+1: lithosphere
+2: 100-410
+3: 410-660
+4: lower mantle
+
+*/
int layers_r(E,r)
struct All_variables *E;
float r;
{
- float zlith, z410, zlm;
+ float rlith, r410, rlm;
int llayers = 0;
- zlith = E->viscosity.zlith;
- z410 = E->viscosity.z410;
- zlm = E->viscosity.zlm;
+ /*
+ the z-values, as read in, are non-dimensionalized depth
+ convert to radii
- if (r > (E->sphere.ro-zlith))
+ */
+ rlith = E->sphere.ro - E->viscosity.zlith; /* lithosphere */
+ r410 = E->sphere.ro - E->viscosity.z410;
+ rlm = E->sphere.ro - E->viscosity.zlm;
+
+ if (r > rlith) /* in lithospherre */
llayers = 1;
- else if ((r > (E->sphere.ro-z410))&&
- (r <= (E->sphere.ro-zlith)))
+ else if ((r > r410)&& (r <= rlith)) /* in asthenosphere 100...410 km */
llayers = 2;
- else if ((r > (E->sphere.ro-zlm)) &&
- (r <= (E->sphere.ro-z410)))
+ else if ((r > rlm) && (r <= r410)) /* in transition zone, 410 - 660 km */
llayers = 3;
- else
+ else /* lower mantle */
llayers = 4;
return (llayers);
}
-
+/* determine layer number of node "node" of cap "m" */
int layers(E,m,node)
struct All_variables *E;
int m,node;
@@ -784,6 +801,10 @@
/* ==============================================================
construct array mat
+
+
+
+
============================================================== */
void construct_mat_group(E)
struct All_variables *E;
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -33,7 +33,8 @@
void ggrd_full_temp_init(struct All_variables *);
#endif
-void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
+void get_r_spacing_fine(double *,struct All_variables *);
+void get_r_spacing_at_levels(double *,struct All_variables *);
/* Setup global mesh parameters */
void full_global_derived_values(E)
@@ -149,13 +150,13 @@
switch(E->control.coor){
- case 2:
- /* higher radial spacing in top and bottom fractions */
- get_r_spacing_fine(rr, (double)E->sphere.ri,(double)E->sphere.ro,
- E->mesh.noz,(double)E->control.coor_refine[0] ,
- (double)E->control.coor_refine[1] ,
- (double)E->control.coor_refine[2] ,
- (double)E->control.coor_refine[3], E);
+ case 0:
+ /* generate uniform mesh in radial direction */
+ dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
+
+ for (k=1;k <= E->mesh.noz;k++) {
+ rr[k] = E->sphere.ri + (k-1)*dr;
+ }
break;
case 1: /* read nodal radii from file */
sprintf(output_file,"%s",E->control.coor_file);
@@ -172,21 +173,27 @@
fclose(fp1);
break;
+ case 2:
+ /* higher radial spacing in top and bottom fractions */
+ get_r_spacing_fine(rr,E);
+ break;
+ case 3:
+ /* assign radial spacing CitcomCU style */
+ get_r_spacing_at_levels(rr,E);
+ break;
default:
- /* generate uniform mesh in radial direction */
- dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
-
- for (k=1;k<=E->mesh.noz;k++) {
- rr[k] = E->sphere.ri + (k-1)*dr;
- }
+ myerror(E,"coor flag undefined in Full_version_dependent");
break;
}
-
+
for (i=1;i<=E->lmesh.noz;i++) {
- k = E->lmesh.nzs+i-1;
- RR[i] = rr[k];
- }
+ k = E->lmesh.nzs+i-1;
+ RR[i] = rr[k];
+
+ }
+
+
for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
if (E->control.NMULTIGRID||E->control.EMULTIGRID)
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -47,6 +47,7 @@
void myerror(struct All_variables *,char *);
void ggrd_full_temp_init(struct All_variables *);
void ggrd_reg_temp_init(struct All_variables *);
+void ggrd_temp_init_general(struct All_variables *,char *);
/*
@@ -136,14 +137,24 @@
report(E,"ggrd tracer init done");
}
+void ggrd_full_temp_init(struct All_variables *E)
+{
+ ggrd_temp_init_general(E,"-L");
+}
+void ggrd_reg_temp_init(struct All_variables *E)
+{
+ ggrd_temp_init_general(E,"");
+}
+
+
/*
initialize temperatures from grd files for spherical geometry
*/
-void ggrd_full_temp_init(struct All_variables *E)
+void ggrd_temp_init_general(struct All_variables *E,char *gmtflag)
{
MPI_Status mpi_stat;
@@ -159,7 +170,7 @@
noxnoz = nox * noz;
if(E->parallel.me == 0)
- fprintf(stderr,"ggrd_full_temp_init: using GMT grd files for temperatures\n");
+ fprintf(stderr,"ggrd_temp_init_general: using GMT grd files for temperatures, gmtflag: %s\n",gmtflag);
/*
@@ -190,7 +201,7 @@
*/
E->convection.ggrd_tinit_d[0].init = FALSE;
if(ggrd_grdtrack_init_general(TRUE,E->convection.ggrd_tinit_gfile,
- E->convection.ggrd_tinit_dfile,"-L", /* this could also be -Lg for global */
+ E->convection.ggrd_tinit_dfile,gmtflag,
E->convection.ggrd_tinit_d,(E->parallel.me == 0),
FALSE))
myerror(E,"grd init error");
@@ -198,7 +209,7 @@
/* tell the next processor to go ahead with the init step */
mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
}else{
- fprintf(stderr,"ggrd_full_temp_init: last processor (%i) done with grd init\n",
+ fprintf(stderr,"ggrd_temp_init_general: last processor (%i) done with grd init\n",
E->parallel.me);
}
/*
@@ -236,7 +247,7 @@
(double)E->sx[m][2][node],
E->convection.ggrd_tinit_d,&tadd,
FALSE))
- myerror(E,"ggrd_full_temp_init");
+ myerror(E,"ggrd_temp_init_general");
if(E->convection.ggrd_tinit_scale_with_prem){
/*
get the PREM density at r for additional scaling
@@ -254,22 +265,28 @@
E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale;
}
+ if(E->convection.ggrd_tinit_limit_trange){
+ /* limit to 0 < T < 1 ?*/
+ E->T[m][node] = min(max(E->T[m][node], 0.0),1.0);
+ }
+ //fprintf(stderr,"z: %11g T: %11g\n",E->sx[m][3][node],E->T[m][node]);
if(E->convection.ggrd_tinit_override_tbc){
if((k == 1) && (E->mesh.bottbc == 1)){ /* bottom TBC */
E->sphere.cap[m].TB[1][node] = E->T[m][node];
E->sphere.cap[m].TB[2][node] = E->T[m][node];
E->sphere.cap[m].TB[3][node] = E->T[m][node];
+ //fprintf(stderr,"z: %11g TBB: %11g\n",E->sx[m][3][node],E->T[m][node]);
}
if((k == noz) && (E->mesh.toptbc == 1)){ /* top TBC */
E->sphere.cap[m].TB[1][node] = E->T[m][node];
E->sphere.cap[m].TB[2][node] = E->T[m][node];
E->sphere.cap[m].TB[3][node] = E->T[m][node];
+ //fprintf(stderr,"z: %11g TBT: %11g\n",E->sx[m][3][node],E->T[m][node]);
}
}
+
- //if(E->convection.ggrd_tinit_limit_unity)
- /* limit to >0 */
- //E->T[m][node] = min(max(E->T[m][node], 0.0),1.0);
+
}
/*
free the structure, not needed anymore since T should now
@@ -279,82 +296,11 @@
/*
end temperature/density from GMT grd init
*/
-
+ temperatures_conform_bcs(E);
}
-/*
-initialize temperatures from grd files for spherical geometry, regional version
-(this is identical to global version but for flags for GMT interpolation)
-*/
-void ggrd_reg_temp_init(struct All_variables *E)
-{
-
- MPI_Status mpi_stat;
- int mpi_rc, mpi_inmsg, mpi_success_message = 1;
- double temp1,tbot,tgrad,tmean,tadd,rho_prem;
- int i,j,k,m,ii,node,noxnoz,nox,noz,noy;
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- noxnoz = nox * noz;
- if(E->parallel.me==0)
- fprintf(stderr,"ggrd_reg_temp_init: using GMT grd files for temperatures\n");
- if(E->parallel.me > 0)
- mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
- if(E->convection.ggrd_tinit_scale_with_prem){
- if(prem_read_model(E->convection.prem.model_filename,&E->convection.prem, (E->parallel.me==0)))
- myerror(E,"prem error");
- }
- if(ggrd_grdtrack_init_general(TRUE,E->convection.ggrd_tinit_gfile,
- E->convection.ggrd_tinit_dfile,"", /* this part differnent from global */
- E->convection.ggrd_tinit_d,(E->parallel.me==0),
- FALSE))
- myerror(E,"grdtrack init error");
- if(E->parallel.me < E->parallel.nproc-1){
- mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
- }else{
- fprintf(stderr,"ggrd_reg_temp_init: last processor (%i) done with grd init\n",
- E->parallel.me);
- }
- tmean = (E->control.TBCbotval + E->control.TBCtopval)/2.0 + E->convection.ggrd_tinit_offset;
- for(m=1;m <= E->sphere.caps_per_proc;m++)
- for(i=1;i <= noy;i++)
- for(j=1;j <= nox;j++)
- for(k=1;k <= noz;k++) {
- node=k+(j-1)*noz+(i-1)*noxnoz;
- if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],(double)E->sx[m][1][node],
- (double)E->sx[m][2][node],E->convection.ggrd_tinit_d,&tadd,FALSE))
- myerror(E,"ggrd_reg_temp_init");
- if(E->convection.ggrd_tinit_scale_with_prem){
- prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->convection.prem);
- if(rho_prem < 3200.0)
- rho_prem = 3200.0;
- E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale * rho_prem / E->data.density;
- }else{
- E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale;
- }
- if(E->convection.ggrd_tinit_override_tbc){
- if((k == 1) && (E->mesh.bottbc == 1)){ /* bottom TBC */
- E->sphere.cap[m].TB[1][node] = E->T[m][node];
- E->sphere.cap[m].TB[2][node] = E->T[m][node];
- E->sphere.cap[m].TB[3][node] = E->T[m][node];
- }
- if((k == noz) && (E->mesh.toptbc == 1)){ /* top TBC */
- E->sphere.cap[m].TB[1][node] = E->T[m][node];
- E->sphere.cap[m].TB[2][node] = E->T[m][node];
- E->sphere.cap[m].TB[3][node] = E->T[m][node];
- }
- }
- }
- ggrd_grdtrack_free_gstruc(E->convection.ggrd_tinit_d);
-}
-
-
-
-
-
#endif
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -151,16 +151,15 @@
/* read in some more parameters */
/* scale the anomalies with PREM densities */
input_boolean("ggrd_tinit_scale_with_prem",&(E->convection.ggrd_tinit_scale_with_prem),"off",E->parallel.me);
- /* scaling factor for the grids */
+ /* limit T to 0...1 */
+ input_boolean("ggrd_tinit_limit_trange",&(E->convection.ggrd_tinit_limit_trange),"on",E->parallel.me);
+ /* scaling factor for the grids */
input_double("ggrd_tinit_scale",&(E->convection.ggrd_tinit_scale),"1.0",E->parallel.me); /* scale */
/* temperature offset factor */
input_double("ggrd_tinit_offset",&(E->convection.ggrd_tinit_offset),"0.0",E->parallel.me); /* offset */
/* grid name, without the .i.grd suffix */
input_string("ggrd_tinit_gfile",E->convection.ggrd_tinit_gfile,"",E->parallel.me); /* grids */
- input_string("ggrd_tinit_dfile",E->convection.ggrd_tinit_dfile,"",E->parallel.me); /* depth.dat
- layers
- of
- grids*/
+ input_string("ggrd_tinit_dfile",E->convection.ggrd_tinit_dfile,"",E->parallel.me); /* depth.dat layers of grids*/
/* override temperature boundary condition? */
input_boolean("ggrd_tinit_override_tbc",&(E->convection.ggrd_tinit_override_tbc),"off",E->parallel.me);
input_string("ggrd_tinit_prem_file",E->convection.prem.model_filename,"", E->parallel.me); /* PREM model filename */
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -341,13 +341,35 @@
input_int("levels",&(E->mesh.levels),"0",m);
input_int("coor",&(E->control.coor),"0",m);
- if(E->control.coor == 2){ /* refinement */
+ if(E->control.coor == 2){
+ /*
+ refinement in two layers
+ */
+ /* number of refinement layers */
E->control.coor_refine[0] = 0.10; /* bottom 10% */
E->control.coor_refine[1] = 0.15; /* get 15% of the nodes */
E->control.coor_refine[2] = 0.10; /* top 10% */
E->control.coor_refine[3] = 0.20; /* get 20% of the nodes */
input_float_vector("coor_refine",4,E->control.coor_refine,m);
- }
+ }else if(E->control.coor == 3){
+ /*
+
+ refinement CitcomCU style, by reading in layers, e.g.
+
+ r_grid_layers=3 # minus 1 is number of layers with uniform grid in r
+ rr=0.5,0.75,1.0 # starting and ending r coodinates
+ nr=1,37,97 # starting and ending node in r direction
+
+ */
+ input_int("r_grid_layers", &(E->control.rlayers), "1",m);
+ if(E->control.rlayers > 20)
+ myerror(E,"number of rlayers out of bounds (20) for coor = 3");
+ /* layers radii */
+ input_float_vector("rr", E->control.rlayers, (E->control.rrlayer),m);
+ /* associated node numbers */
+ input_int_vector("nr", E->control.rlayers, (E->control.nrlayer),m);
+ }
+
input_string("coor_file",E->control.coor_file,"",m);
input_int("nprocx",&(E->parallel.nprocx),"1",m);
@@ -369,10 +391,16 @@
input_int("solution_cycles_init",&(E->monitor.solution_cycles_init),"0",m);
/* for layers */
- input_float("z_cmb",&(E->viscosity.zcmb),"0.45",m);
+ /*
+
+ these boundaries are a little wacko
+
+
+ */
+ input_float("z_cmb",&(E->viscosity.zcmb),"0.45",m); /* does this ever get used? */
input_float("z_lmantle",&(E->viscosity.zlm),"0.45",m);
- input_float("z_410",&(E->viscosity.z410),"0.225",m);
- input_float("z_lith",&(E->viscosity.zlith),"0.225",m);
+ input_float("z_410",&(E->viscosity.z410),"0.225",m); /* 0.06434, more like it */
+ input_float("z_lith",&(E->viscosity.zlith),"0.225",m); /* 0.0157, more like it */
/* the start age and initial subduction history */
input_float("start_age",&(E->control.start_age),"0.0",m);
@@ -389,6 +417,9 @@
input_float("topvbyval",&(E->control.VBYtopval),"0.0",m);
input_float("botvbyval",&(E->control.VBYbotval),"0.0",m);
+
+ input_float("T_interior_max_for_exit",&(E->monitor.T_interior_max_for_exit),"1.5",m);
+
input_int("pseudo_free_surf",&(E->control.pseudo_free_surf),"0",m);
input_int("toptbc",&(E->mesh.toptbc),"1",m);
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -56,9 +56,9 @@
void convert_pvec_to_cvec(float ,float , float , float *,float *);
void *safe_malloc (size_t );
void myerror(struct All_variables *,char *);
-void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,
- struct All_variables *);
void xyz2rtp(float ,float ,float ,float *);
+void get_r_spacing_fine(double *,struct All_variables *);
+void get_r_spacing_at_levels(double *,struct All_variables *);
int get_process_identifier()
@@ -443,6 +443,8 @@
parallel_process_termination();
}
+
+
/*
@@ -450,15 +452,21 @@
attempt to space rr[1...nz] such that bfrac*nz nodes will be within the lower
brange fraction of (ro-ri), and similar for the top layer
+function below is more general
+
*/
-void get_r_spacing_fine(double *rr, double ri, double ro,
- int nz,double brange, double bfrac,
- double trange, double tfrac,
- struct All_variables *E)
+void get_r_spacing_fine(double *rr, struct All_variables *E)
{
int k,klim,nb,nt,nm;
- double drb,dr0,drt,dr,drm,range,r,mrange;
- range = ro-ri; /* original range */
+ double drb,dr0,drt,dr,drm,range,r,mrange, brange,bfrac,trange, tfrac;
+
+ brange = (double)E->control.coor_refine[0];
+ bfrac = (double)E->control.coor_refine[1];
+ trange = (double)E->control.coor_refine[2];
+ tfrac = (double)E->control.coor_refine[3];
+
+ range = (double) E->sphere.ro - E->sphere.ri; /* original range */
+
mrange = 1 - brange - trange;
if(mrange <= 0)
myerror(E,"get_r_spacing_fine: bottom and top range too large");
@@ -467,9 +475,9 @@
trange *= range; /* top */
mrange *= range; /* middle */
- nb = nz * bfrac;
- nt = nz * tfrac;
- nm = nz-nb-nt;
+ nb = E->mesh.noz * bfrac;
+ nt = E->mesh.noz * tfrac;
+ nm = E->mesh.noz - nb - nt;
if((nm < 1)||(nt < 2)||(nb < 2))
myerror(E,"get_r_spacing_fine: refinement out of bounds");
@@ -477,14 +485,55 @@
drt = trange/(nt-1);
drm = mrange / (nm + 1);
- for(r=ri,k=1;k<=nb;k++,r+=drb){
+ for(r=E->sphere.ri,k=1;k<=nb;k++,r+=drb){
rr[k] = r;
}
- klim = nz-nt+1;
+ klim = E->mesh.noz - nt + 1;
for(r=r-drb+drm;k < klim;k++,r+=drm){
rr[k] = r;
}
- for(;k <= nz;k++,r+=drt){
+ for(;k <= E->mesh.noz;k++,r+=drt){
rr[k] = r;
}
}
+/*
+
+
+get r spacing at radial locations and node numbers as specified
+CitcomCU style
+
+rr[1...E->mesh.noz]
+
+
+e.g.:
+
+ r_grid_layers=3 # minus 1 is number of layers with uniform grid in r
+ rr=0.5,0.75,1.0 # starting and ending r coodinates
+ nr=1,37,97 # starting and ending node in r direction
+
+*/
+void get_r_spacing_at_levels(double *rr,struct All_variables *E)
+{
+ double ddr;
+ int k,j;
+ /* do some sanity checks */
+ if(E->control.nrlayer[0] != 1)
+ myerror(E,"first node for coor=3 should be unity");
+ if(E->control.nrlayer[E->control.rlayers-1] != E->mesh.noz)
+ myerror(E,"last node for coor = 3 input should match max nr z nodes");
+ if(fabs(E->control.rrlayer[0] -E->sphere.ri) > 1e-6)
+ myerror(E,"inner layer for coor=3 input should be inner radius");
+ if(fabs(E->control.rrlayer[ E->control.rlayers-1] - E->sphere.ro)>1e-6)
+ myerror(E,"outer layer for coor=3 input should be inner radius");
+ if(E->control.rlayers < 2)
+ myerror(E,"number of rlayers needs to be at leats two for coor = 3");
+
+ rr[1] = E->control.rrlayer[0];
+ for(j = 1; j < E->control.rlayers; j++){
+ ddr = (E->control.rrlayer[j] - E->control.rrlayer[j - 1]) /
+ (E->control.nrlayer[j] - E->control.nrlayer[j - 1]);
+ for(k = E->control.nrlayer[j-1]+1;k <= E->control.nrlayer[j];k++)
+ rr[k] = rr[k-1]+ddr;
+ }
+
+}
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -166,9 +166,14 @@
} */
if (E->parallel.me==E->parallel.nprocz-1) {
fprintf(stderr,"surface heat flux= %f\n",sum_h[0]);
- /*fprintf(E->fp,"surface heat flux= %f\n",sum_h[0]);*/
- if(E->output.write_q_files)
- fprintf(E->output.fpqt,"%.5e %.5e\n",E->monitor.elapsed_time,sum_h[0]);
+ /* XXX */
+ fprintf(E->fp,"surface heat flux= %f\n",sum_h[0]); /* commented this back in TWB , why was it out in the first place? */
+
+ if(E->output.write_q_files > 0){
+ /* format: time heat_flow sqrt(v.v) */
+ fprintf(E->output.fpqt,"%13.5e %13.5e %13.5e\n",E->monitor.elapsed_time,sum_h[0],sqrt(E->monitor.vdotv));
+ fflush(E->output.fpqt);
+ }
}
}
@@ -178,8 +183,11 @@
if (E->parallel.me==0) {
fprintf(stderr,"bottom heat flux= %f\n",sum_h[2]);
fprintf(E->fp,"bottom heat flux= %f\n",sum_h[2]);
- if(E->output.write_q_files)
- fprintf(E->output.fpqb,"%.5e %.5e\n",E->monitor.elapsed_time,sum_h[2]);
+ if(E->output.write_q_files > 0){
+ fprintf(E->output.fpqb,"%13.5e %13.5e %13.5e\n",
+ E->monitor.elapsed_time,sum_h[2],sqrt(E->monitor.vdotv));
+ fflush(E->output.fpqb);
+ }
}
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -72,54 +72,57 @@
temp = max(E->mesh.NOY[E->mesh.levmax],E->mesh.NOX[E->mesh.levmax]);
-if(E->control.coor==1) {
- for(i=E->mesh.gridmin;i<=E->mesh.gridmax;i++) {
- theta1[i] = (float *)malloc((temp+1)*sizeof(float));
- fi1[i] = (float *)malloc((temp+1)*sizeof(float));
- }
+ if(E->control.coor==1) {
- temp = E->mesh.NOY[E->mesh.levmax]*E->mesh.NOX[E->mesh.levmax];
+ /* read in node locations from file */
- sprintf(output_file,"%s",E->control.coor_file);
- fp=fopen(output_file,"r");
- if (fp == NULL) {
- fprintf(E->fp,"(Sphere_related #1) Cannot open %s\n",output_file);
- exit(8);
- }
+ for(i=E->mesh.gridmin;i<=E->mesh.gridmax;i++) {
+ theta1[i] = (float *)malloc((temp+1)*sizeof(float));
+ fi1[i] = (float *)malloc((temp+1)*sizeof(float));
+ }
+
+ temp = E->mesh.NOY[E->mesh.levmax]*E->mesh.NOX[E->mesh.levmax];
+
+ sprintf(output_file,"%s",E->control.coor_file);
+ fp=fopen(output_file,"r");
+ if (fp == NULL) {
+ fprintf(E->fp,"(Sphere_related #1) Cannot open %s\n",output_file);
+ exit(8);
+ }
- fscanf(fp,"%s %d",a,&nn);
- for(i=1;i<=gnox;i++) {
- fscanf(fp,"%d %e",&nn,&theta1[E->mesh.gridmax][i]);
- }
- E->control.theta_min = theta1[E->mesh.gridmax][1];
- E->control.theta_max = theta1[E->mesh.gridmax][gnox];
-
- fscanf(fp,"%s %d",a,&nn);
- for(i=1;i<=gnoy;i++) {
- fscanf(fp,"%d %e",&nn,&fi1[E->mesh.gridmax][i]);
+ fscanf(fp,"%s %d",a,&nn);
+ for(i=1;i<=gnox;i++) {
+ fscanf(fp,"%d %e",&nn,&theta1[E->mesh.gridmax][i]);
}
- E->control.fi_min = fi1[E->mesh.gridmax][1];
- E->control.fi_max = fi1[E->mesh.gridmax][gnox];
-
- fclose(fp);
-
-
- for (lev=E->mesh.gridmin;lev<=E->mesh.gridmax;lev++) {
-
- if (E->control.NMULTIGRID||E->control.EMULTIGRID)
+ E->control.theta_min = theta1[E->mesh.gridmax][1];
+ E->control.theta_max = theta1[E->mesh.gridmax][gnox];
+
+ fscanf(fp,"%s %d",a,&nn);
+ for(i=1;i<=gnoy;i++) {
+ fscanf(fp,"%d %e",&nn,&fi1[E->mesh.gridmax][i]);
+ }
+ E->control.fi_min = fi1[E->mesh.gridmax][1];
+ E->control.fi_max = fi1[E->mesh.gridmax][gnox];
+
+ fclose(fp);
+
+
+ for (lev=E->mesh.gridmin;lev<=E->mesh.gridmax;lev++) {
+
+ if (E->control.NMULTIGRID||E->control.EMULTIGRID)
step = (int) pow(2.0,(double)(E->mesh.levmax-lev));
- else
+ else
step = 1;
+
+ for (i=1;i<=E->mesh.NOX[lev];i++)
+ theta1[lev][i] = theta1[E->mesh.gridmax][(i-1)*step+1];
+
+ for (i=1;i<=E->mesh.NOY[lev];i++)
+ fi1[lev][i] = fi1[E->mesh.gridmax][(i-1)*step+1];
+
+ }
- for (i=1;i<=E->mesh.NOX[lev];i++)
- theta1[lev][i] = theta1[E->mesh.gridmax][(i-1)*step+1];
- for (i=1;i<=E->mesh.NOY[lev];i++)
- fi1[lev][i] = fi1[E->mesh.gridmax][(i-1)*step+1];
-
- }
-
-
for (lev=E->mesh.gridmin;lev<=E->mesh.gridmax;lev++) {
elx = E->lmesh.ELX[lev];
ely = E->lmesh.ELY[lev];
@@ -162,10 +165,10 @@
free ((void *)fi1[lev] );
}
-}
+ } /* end of coord = 1 */
+
+ else if((E->control.coor==0) || (E->control.coor==2)|| (E->control.coor==3)) {
- else if((E->control.coor==0) || (E->control.coor==2)) {
-
/*
for(i=1;i<=5;i++) {
x[i] = (double *) malloc((E->parallel.nproc+1)*sizeof(double));
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -29,8 +29,10 @@
#include "global_defs.h"
#include "parallel_related.h"
-void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
+void get_r_spacing_fine(double *,struct All_variables *);
+void get_r_spacing_at_levels(double *,struct All_variables *);
+
#ifdef USE_GGRD
void ggrd_reg_temp_init(struct All_variables *);
#endif
@@ -159,7 +161,16 @@
noz=E->mesh.noz;
- if(E->control.coor==1) { /* get nodal levels from file */
+ switch(E->control.coor) {
+ case 0:
+ /* default: regular node spacing */
+ dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
+ for (k=1;k<=E->mesh.noz;k++) {
+ rr[k] = E->sphere.ri + (k-1)*dr;
+ }
+ break;
+ case 1:
+ /* get nodal levels from file */
sprintf(output_file,"%s",E->control.coor_file);
fp1=fopen(output_file,"r");
if (fp1 == NULL) {
@@ -184,23 +195,18 @@
E->sphere.ro = rr[E->mesh.noz];
fclose(fp1);
-
- } else if(E->control.coor==0) {
- /* default: regular node spacing */
- dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
- for (k=1;k<=E->mesh.noz;k++) {
- rr[k] = E->sphere.ri + (k-1)*dr;
- }
- } else if(E->control.coor==2){
+ break;
+ case 2:
/* higher radial spacing in top and bottom fractions */
- get_r_spacing_fine(rr, (double)E->sphere.ri,(double)E->sphere.ro,
- E->mesh.noz,(double)E->control.coor_refine[0] ,
- (double)E->control.coor_refine[1] ,
- (double)E->control.coor_refine[2] ,
- (double)E->control.coor_refine[3],E);
-
- } else {
+ get_r_spacing_fine(rr, E);
+ break;
+ case 3:
+ /* assign radial spacing CitcomCU style */
+ get_r_spacing_at_levels(rr,E);
+ break;
+ default:
myerror(E,"regional_version_dependent: coor mode not implemented");
+ break;
}
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-10 23:40:38 UTC (rev 8101)
@@ -472,11 +472,13 @@
for(m=1;m <= E->sphere.caps_per_proc;m++)
for(i=1;i <= nel;i++) {
- l = E->mat[m][i];
+
+ l = E->mat[m][i] - 1;
+
if(E->control.mat_control)
- tempa = E->viscosity.N0[l-1] * E->VIP[m][i];
+ tempa = E->viscosity.N0[l] * E->VIP[m][i];
else
- tempa = E->viscosity.N0[l-1];
+ tempa = E->viscosity.N0[l];
j = 0;
for(kk=1;kk<=ends;kk++) {
@@ -491,11 +493,11 @@
temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
}
- //fprintf(stderr,"N0 %11g T %11g T0 %11g E %11g 1-z %11g Z %11g\n",
- // tempa,temp,E->viscosity.T[l-1],E->viscosity.E[l-1], zzz ,E->viscosity.Z[l-1] );
EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( E->viscosity.E[l-1]*(E->viscosity.T[l-1] - temp) +
- zzz * E->viscosity.Z[l-1] );
+ exp( E->viscosity.E[l]*(E->viscosity.T[l] - temp) +
+ zzz * E->viscosity.Z[l]);
+ //fprintf(stderr,"N0 %11g T %11g T0 %11g E %11g z %11g km Z %11g mat: %i log10(eta): %11g\n",
+ // tempa,temp,E->viscosity.T[l],E->viscosity.E[l], zzz *6371 ,E->viscosity.Z[l],l+1,log10(EEta[m][ (i-1)*vpts + jj ]));
}
}
break;
Modified: mc/3D/CitcomS/trunk/lib/convection_variables.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/convection_variables.h 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/convection_variables.h 2007-10-10 23:40:38 UTC (rev 8101)
@@ -54,7 +54,7 @@
#ifdef USE_GGRD
/* for temperature init from grd files */
int ggrd_tinit,ggrd_tinit_scale_with_prem;
- int ggrd_tinit_override_tbc;
+ int ggrd_tinit_override_tbc,ggrd_tinit_limit_trange;
double ggrd_tinit_scale,ggrd_tinit_offset,ggrd_vstage_transition;
char ggrd_tinit_gfile[1000];
char ggrd_tinit_dfile[1000];
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-10-10 23:35:45 UTC (rev 8100)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-10-10 23:40:38 UTC (rev 8101)
@@ -382,6 +382,8 @@
float T_interior;
float T_maxvaried;
+ float T_interior_max_for_exit;
+
};
struct CONTROL {
@@ -472,6 +474,9 @@
int coor;
float coor_refine[4];
+ float rrlayer[20];
+ int nrlayer[20],rlayers;
+
char coor_file[100];
int lith_age;
More information about the cig-commits
mailing list