[cig-commits] r16778 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Mon May 24 11:00:39 PDT 2010
Author: becker
Date: 2010-05-24 11:00:39 -0700 (Mon, 24 May 2010)
New Revision: 16778
Modified:
mc/3D/CitcomS/trunk/lib/BC_util.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/ggrd_handling.h
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
toplayerbc > 0 will now assign BCs for all noes with r >= toplayerbc_r, which defaults to
0.984303876942 i.e. 100 km depth. (Still experimental.)
Modified: mc/3D/CitcomS/trunk/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/BC_util.c 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/BC_util.c 2010-05-24 18:00:39 UTC (rev 16778)
@@ -255,8 +255,7 @@
options:
-toplayerbc > 0: assign surface boundary condition down to all nodes within the the toplayerbc
- layer, e.g. layer one for lithosphere
+toplayerbc > 0: assign surface boundary condition down to all nodes with r >= toplayerbc_r
toplayerbc == 0: no action
toplayerbc < 0: assign surface boundary condition within medium at node -toplayerbc depth, ie.
toplayerbc = -1 is one node underneath surface
@@ -279,8 +278,8 @@
ontop = ((k==noz) && (E->parallel.me_loc[3]==E->parallel.nprocz-1))?(1):(0);
onbottom = ((k==1) && (E->parallel.me_loc[3]==0))?(1):(0);
/* node number is k, assuming no dependence on x and y */
- if((lay = layers(E,j,k)) <= E->mesh.toplayerbc){
-
+ if(E->SX[lv][lv][3][k] >= E->mesh.toplayerbc_r){
+ lay = layers(E,j,k);
if((!ontop)&&(!onbottom)&&(lv==E->mesh.gridmax))
ncount++; /* not in top or bottom */
if(E->mesh.topvbc != 1) { /* free slip */
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-05-24 18:00:39 UTC (rev 16778)
@@ -175,7 +175,7 @@
MPI_Status mpi_stat;
int mpi_rc;
int mpi_inmsg, mpi_success_message = 1;
- double temp1,tbot,tgrad,tmean,tadd,rho_prem;
+ double temp1,tbot,tgrad,tmean,tadd,rho_prem,depth,loc_scale;
char gmt_string[10];
int i,j,k,m,node,noxnoz,nox,noy,noz;
static ggrd_boolean shift_to_pos_lon = FALSE;
@@ -248,10 +248,6 @@
*/
tbot = 1.0;
}
- /*
- mean temp is (top+bot)/2 + offset
- */
- tmean = (tbot + E->control.TBCtopval)/2.0 + E->control.ggrd.temp.offset;
for(m=1;m <= E->sphere.caps_per_proc;m++)
@@ -264,31 +260,46 @@
/*
get interpolated velocity anomaly
*/
- if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],(double)E->sx[m][1][node],
+ depth = (1-E->sx[m][3][node])*6371;
+ if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],
+ (double)E->sx[m][1][node],
(double)E->sx[m][2][node],
E->control.ggrd.temp.d,&tadd,
FALSE,shift_to_pos_lon)){
fprintf(stderr,"%g %g %g\n",E->sx[m][2][node]*57.29577951308232087,
- 90-E->sx[m][1][node]*57.29577951308232087,(1-E->sx[m][3][node])*6371);
+ 90-E->sx[m][1][node]*57.29577951308232087,depth);
myerror(E,"ggrd__temp_init_general: interpolation error");
}
+
+ if(depth < E->control.ggrd_lower_depth_km){
+ /*
+ mean temp is (top+bot)/2 + offset
+ */
+ tmean = (tbot + E->control.TBCtopval)/2.0 + E->control.ggrd.temp.offset;
+ loc_scale = E->control.ggrd.temp.scale;
+ }else{ /* lower mantle */
+ tmean = (tbot + E->control.TBCtopval)/2.0 + E->control.ggrd_lower_offset;
+ loc_scale = E->control.ggrd_lower_scale;
+ }
if(E->control.ggrd.temp.scale_with_prem){
/*
get the PREM density at r for additional scaling
*/
prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->control.ggrd.temp.prem);
- if(rho_prem < 3200.0)
- rho_prem = 3200.0; /* we don't want the density of water */
+ if(rho_prem < GGRD_DENS_MIN){
+ fprintf(stderr,"WARNING: restricting minimum density to %g, would have been %g\n",
+ GGRD_DENS_MIN,rho_prem);
+ rho_prem = GGRD_DENS_MIN; /* we don't want the density of water or crust*/
+ }
/*
assign temperature
*/
- E->T[m][node] = tmean + tadd * E->control.ggrd.temp.scale *
- rho_prem / E->data.density;
+ E->T[m][node] = tmean + tadd * loc_scale * rho_prem / E->data.density;
}else{
/* no PREM scaling */
- E->T[m][node] = tmean + tadd * E->control.ggrd.temp.scale;
+ E->T[m][node] = tmean + tadd * loc_scale;
}
if(E->control.ggrd.temp.limit_trange){
@@ -1209,7 +1220,7 @@
if(E->mesh.toplayerbc > 0){
/* check for internal nodes in layers */
for(k=nozl;k >= 1;k--){
- if(layers(E,m,k) > E->mesh.toplayerbc) /* assume regular mesh structure */
+ if(E->SX[level][m][3][k] < E->mesh.toplayerbc_r) /* assume regular mesh structure */
break;
}
if(k == nozl){ /* */
@@ -1227,13 +1238,14 @@
*assign = TRUE;
*topnode = *botnode;
}else{
- fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: toplayerbc %i\n", E->mesh.toplayerbc);
+ fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: toplayerbc %i, r_min: %g\n",
+ E->mesh.toplayerbc,E->mesh.toplayerbc_r);
myerror(E,"ggrd_vtop_helper_decide_on_internal_nodes: logic error");
}
}
if(verbose)
- fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: mixed: internal: %i assign: %i k: %i to %i (%i)\n",
- allow_internal,*assign,*botnode,*topnode,nozl);
+ fprintf(stderr,"ggrd_vtop_helper_decide_on_internal_nodes: mixed: internal: %i assign: %i k: %i to %i (%i), r_min: %g\n",
+ allow_internal,*assign,*botnode,*topnode,nozl,E->mesh.toplayerbc_r);
}
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2010-05-24 18:00:39 UTC (rev 16778)
@@ -164,19 +164,37 @@
*/
/* scale the anomalies with PREM densities */
- input_boolean("ggrd_tinit_scale_with_prem",&(E->control.ggrd.temp.scale_with_prem),"off",E->parallel.me);
+ input_boolean("ggrd_tinit_scale_with_prem",
+ &(E->control.ggrd.temp.scale_with_prem),"off",E->parallel.me);
/* limit T to 0...1 */
- input_boolean("ggrd_tinit_limit_trange",&(E->control.ggrd.temp.limit_trange),"on",E->parallel.me);
+ input_boolean("ggrd_tinit_limit_trange",
+ &(E->control.ggrd.temp.limit_trange),"on",E->parallel.me);
/* scaling factor for the grids */
- input_double("ggrd_tinit_scale",&(E->control.ggrd.temp.scale),"1.0",E->parallel.me); /* scale */
+ input_double("ggrd_tinit_scale",
+ &(E->control.ggrd.temp.scale),"1.0",E->parallel.me); /* scale */
/* temperature offset factor */
- input_double("ggrd_tinit_offset",&(E->control.ggrd.temp.offset),"0.0",E->parallel.me); /* offset */
+ input_double("ggrd_tinit_offset",
+ &(E->control.ggrd.temp.offset),"0.0",E->parallel.me); /* offset */
+ /*
+ do we want a different scaling for the lower mantle?
+ */
+ input_float("ggrd_lower_depth_km",&(E->control.ggrd_lower_depth_km),"7000",
+ E->parallel.me); /* depth, in km, below which
+ different scaling applies */
+ input_float("ggrd_lower_scale",&(E->control.ggrd_lower_scale),"1.0",E->parallel.me);
+ input_float("ggrd_lower_offset",&(E->control.ggrd_lower_offset),"1.0",E->parallel.me);
+
/* grid name, without the .i.grd suffix */
- input_string("ggrd_tinit_gfile",E->control.ggrd.temp.gfile,"",E->parallel.me); /* grids */
- input_string("ggrd_tinit_dfile",E->control.ggrd.temp.dfile,"",E->parallel.me); /* depth.dat layers of grids*/
+ input_string("ggrd_tinit_gfile",
+ E->control.ggrd.temp.gfile,"",E->parallel.me); /* grids */
+ input_string("ggrd_tinit_dfile",
+ E->control.ggrd.temp.dfile,"",E->parallel.me); /* depth.dat layers of grids*/
/* override temperature boundary condition? */
- input_boolean("ggrd_tinit_override_tbc",&(E->control.ggrd.temp.override_tbc),"off",E->parallel.me);
- input_string("ggrd_tinit_prem_file",E->control.ggrd.temp.prem.model_filename,"hc/prem/prem.dat", E->parallel.me); /* PREM model filename */
+ input_boolean("ggrd_tinit_override_tbc",
+ &(E->control.ggrd.temp.override_tbc),"off",E->parallel.me);
+ input_string("ggrd_tinit_prem_file",
+ E->control.ggrd.temp.prem.model_filename,"hc/prem/prem.dat",
+ E->parallel.me); /* PREM model filename */
#else
fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
parallel_process_termination();
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-05-24 18:00:39 UTC (rev 16778)
@@ -434,17 +434,16 @@
input_int("topvbc",&(E->mesh.topvbc),"0",m);
input_int("botvbc",&(E->mesh.botvbc),"0",m);
- input_int("toplayerbc",&(E->mesh.toplayerbc),"0",m); /* > 0: apply
- surface BC
- throughout
- all layer
- nodes of layers toplayerbc (i.e. = 1 for top layer)
+ input_int("toplayerbc",&(E->mesh.toplayerbc),"0",m); /* > 0: apply surface BC
+ throughout all nodes with r > toplayerbc_r
-
< 0: apply to single node layer noz+toplayerbc
*/
+ input_float("toplayerbc_r",&(E->mesh.toplayerbc_r),"0.984303876942",m); /* minimum r to apply BC to,
+ 100 km depth */
+
input_float("topvbxval",&(E->control.VBXtopval),"0.0",m);
input_float("botvbxval",&(E->control.VBXbotval),"0.0",m);
input_float("topvbyval",&(E->control.VBYtopval),"0.0",m);
Modified: mc/3D/CitcomS/trunk/lib/ggrd_handling.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2010-05-24 18:00:39 UTC (rev 16778)
@@ -46,3 +46,4 @@
float find_age_in_MY();
void ggrd_adjust_tbl_rayleigh(struct All_variables *,double **);
+#define GGRD_DENS_MIN 3200 /* minimum density for PREM scaling of input files */
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2010-05-24 15:21:04 UTC (rev 16777)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-05-24 18:00:39 UTC (rev 16778)
@@ -346,7 +346,8 @@
int botvbc;
int sidevbc;
- int toplayerbc; /* apply surface BC throughout top layers, or for a single internal node */
+ int toplayerbc; /* apply surface BC throughout top, or for a single internal node */
+ float toplayerbc_r; /* minimum r to apply BC to */
int periodic_x;
int periodic_y;
@@ -528,6 +529,7 @@
float ggrd_vtop_omega[4];
char ggrd_mat_depth_file[1000];
ggrd_boolean ggrd_mat_is_3d;
+ float ggrd_lower_depth_km,ggrd_lower_scale,ggrd_lower_offset;
#endif
double accuracy,inner_accuracy_scale;
int check_continuity_convergence;
More information about the CIG-COMMITS
mailing list