[cig-commits] r12303 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Sat Jun 21 15:28:52 PDT 2008
Author: becker
Date: 2008-06-21 15:28:51 -0700 (Sat, 21 Jun 2008)
New Revision: 12303
Modified:
mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.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/Regional_boundary_conditions.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Assignment of velocity boundary condition only called once at the top multigrid level
during boundary condition assignment.
Experimental implementation of mixed velocity/free slip boundary
condition based on netcdf grids, in Ggrd_handling.
Modified: mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Full_boundary_conditions.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -37,8 +37,8 @@
static void velocity_apply_periodic_bcs();
static void temperature_apply_periodic_bcs();
void read_temperature_boundary_from_file(struct All_variables *);
+void read_velocity_boundary_from_file(struct All_variables *);
-
/* ========================================== */
void full_velocity_boundary_conditions(E)
@@ -46,7 +46,7 @@
{
void velocity_imp_vert_bc();
void velocity_apply_periodicapply_periodic_bcs();
- void read_velocity_boundary_from_file();
+
void apply_side_sbc();
int j,noz,lv;
@@ -54,7 +54,7 @@
for(lv=E->mesh.gridmax;lv>=E->mesh.gridmin;lv--)
for (j=1;j<=E->sphere.caps_per_proc;j++) {
noz = E->mesh.NOZ[lv];
- if(E->mesh.topvbc != 1) {
+ if(E->mesh.topvbc != 1) { /* free slip top */
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,0.0,VBX,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,VBZ,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,0.0,VBY,0,lv,j);
@@ -62,7 +62,7 @@
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,SBZ,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,E->control.VBYtopval,SBY,1,lv,j);
}
- if(E->mesh.botvbc != 1) {
+ if(E->mesh.botvbc != 1) { /* free slip bottom */
horizontal_bc(E,E->sphere.cap[j].VB,1,1,0.0,VBX,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,1,3,0.0,VBZ,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,1,2,0.0,VBY,0,lv,j);
@@ -71,7 +71,7 @@
horizontal_bc(E,E->sphere.cap[j].VB,1,2,E->control.VBYbotval,SBY,1,lv,j);
}
- if(E->mesh.topvbc == 1) {
+ if(E->mesh.topvbc == 1) { /* velocity/no slip BC */
horizontal_bc(E,E->sphere.cap[j].VB,noz,1,E->control.VBXtopval,VBX,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,VBZ,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,E->control.VBYtopval,VBY,1,lv,j);
@@ -79,11 +79,18 @@
horizontal_bc(E,E->sphere.cap[j].VB,noz,3,0.0,SBZ,0,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,0.0,SBY,0,lv,j);
- if(E->control.vbcs_file)
- read_velocity_boundary_from_file(E);
-
- }
- if(E->mesh.botvbc == 1) {
+ if(E->control.vbcs_file){ /* this should either only be called
+ once, or the input routines need
+ to be told what to do for each
+ multigrid level and cap. it might
+ be easiest to call only once and
+ have routines deal with multigrid
+ */
+ if((lv == E->mesh.gridmin) && (j == E->sphere.caps_per_proc))
+ read_velocity_boundary_from_file(E);
+ }
+ }
+ if(E->mesh.botvbc == 1) { /* velocity bottom BC */
horizontal_bc(E,E->sphere.cap[j].VB,1,1,E->control.VBXbotval,VBX,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,1,3,0.0,VBZ,1,lv,j);
horizontal_bc(E,E->sphere.cap[j].VB,1,2,E->control.VBYbotval,VBY,1,lv,j);
Modified: mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -769,7 +769,6 @@
const int dims=E->mesh.nsd;
-
me = E->parallel.me;
nprocz = E->parallel.nprocz;
@@ -797,6 +796,7 @@
else { /* the last FOUR communications are for lines */
E->parallel.NUM_sNODE[lev][m].pass[kkk]=1;
+
for (k=1;k<=E->parallel.NUM_sNODE[lev][m].pass[kkk];k++) {
if (E->parallel.nprocx*E->parallel.nprocy > 1) { /* 4 or more horiz. proc*/
switch(kkk) {
@@ -838,6 +838,7 @@
} /* end for lev */
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -48,7 +48,7 @@
float age, newage1, newage2;
char output_file1[255],output_file2[255];
float *TB1, *TB2, *VB1[4],*VB2[4], inputage1, inputage2;
- int nox,noz,noy,nnn,nox1,noz1,noy1,lev;
+ int nox,noz,noy,nnn,nox1,noz1,noy1;
int i,ii,ll,m,mm,j,k,n,nodeg,nodel,node,cap;
int intage, pos_age;
int nodea;
@@ -69,7 +69,6 @@
nox1=E->lmesh.nox;
noz1=E->lmesh.noz;
noy1=E->lmesh.noy;
- lev=E->mesh.levmax;
elx=E->lmesh.elx;
elz=E->lmesh.elz;
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -49,13 +49,14 @@
noz = E->mesh.mgunitz * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocz + 1;
if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
+ /* multigrid set up */
E->mesh.levmax = E->mesh.levels-1;
E->mesh.gridmax = E->mesh.levmax;
E->mesh.nox = E->mesh.mgunitx * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocx + 1;
E->mesh.noy = E->mesh.mgunity * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocy + 1;
E->mesh.noz = E->mesh.mgunitz * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocz + 1;
- }
- else {
+ } else {
+ /* regular, conjugate gradient setup */
if (nox!=E->mesh.nox || noy!=E->mesh.noy || noz!=E->mesh.noz) {
if (E->parallel.me==0)
fprintf(stderr,"inconsistent mesh for interpolation, quit the run\n");
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -322,7 +322,7 @@
int mpi_rc,timedep,interpolate;
int mpi_inmsg, mpi_success_message = 1;
int m,el,i,j,k,inode,i1,i2,elxlz,elxlylz,ind;
- int llayer,nox,noy,noz,nox1,noz1,noy1,lev,lselect,idim,elx,ely,elz;
+ int llayer,nox,noy,noz,nox1,noz1,noy1,level,lselect,idim,elx,ely,elz;
char gmt_string[10],char_dummy;
double indbl,indbl2,age,f1,f2,vip,rout[3],xloc[4];
char tfilename[1000];
@@ -335,7 +335,7 @@
elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
elxlz = elx * elz;
elxlylz = elxlz * ely;
- lev=E->mesh.levmax;
+
/*
if we have not initialized the time history structure, do it now
*/
@@ -596,29 +596,33 @@
} /* end ray control */
+/*
+read velocity boundary conditions from netcdf grd files
+
+*/
+
void ggrd_read_vtop_from_file(struct All_variables *E, int is_global)
{
MPI_Status mpi_stat;
int mpi_rc,interpolate,timedep;
- int mpi_inmsg, mpi_success_message = 1;
- int m,el,i,k,i1,i2,ind,nodel;
- int nox1,noz1,noy1,lev,lselect,idim,nox1noz1;
+ int mpi_inmsg, mpi_success_message = 1,nsum[2];
+ int m,el,i,k,i1,i2,ind,nodel,j,nsuml[2],level;
+ int noxg,nozg,noyg,noxl,noyl,nozl,lselect,idim,noxgnozg,noxlnozl;
char gmt_string[10],char_dummy;
-
- double vin1[2],vin2[2],age,f1,f2,vp,vt,vscale,rout[3];
+ static int lc =0; /* only for debugging */
+ double vin1[2],vin2[2],age,f1,f2,vp,vt,vscale,rout[3],cutoff;
char tfilename1[1000],tfilename2[1000];
const int dims=E->mesh.nsd;
- nox1 = E->lmesh.nox;noz1=E->lmesh.noz;noy1=E->lmesh.noy;
- nox1noz1 = nox1*noz1;
- lev = E->mesh.levmax;
+ /* global, top level number of nodes */
+ noxg = E->lmesh.nox;nozg=E->lmesh.noz;noyg=E->lmesh.noy;
+ noxgnozg = noxg*nozg;
/* velocity scaling, assuming input is cm/yr */
vscale = E->data.scalev * E->data.timedir;
-
/*
if we have not initialized the time history structure, do it now
if this file is not found, will use constant velocities
@@ -652,6 +656,9 @@
*/
E->control.ggrd.svt = (struct ggrd_gt *)calloc(E->control.ggrd.time_hist.nvtimes,sizeof(struct ggrd_gt));
E->control.ggrd.svp = (struct ggrd_gt *)calloc(E->control.ggrd.time_hist.nvtimes,sizeof(struct ggrd_gt));
+ /* for detecting the actual max */
+ E->control.ggrd.svt->bandlim = E->control.ggrd.svp->bandlim = 1e6;
+
for(i=0;i < E->control.ggrd.time_hist.nvtimes;i++){
if(!timedep){ /* constant */
sprintf(tfilename1,"%s/vt.grd",E->control.ggrd.vtop_dir);
@@ -671,13 +678,13 @@
mpi_rc = MPI_Send(&mpi_success_message, 1,
MPI_INT, (E->parallel.me+1), 0, E->parallel.world);
}else{
- fprintf(stderr,"ggrd_read_vtop_from_file: last processor done with ggrd vtop BC init, %i timesteps\n",
- E->control.ggrd.time_hist.nvtimes);
+ fprintf(stderr,"ggrd_read_vtop_from_file: last processor done with ggrd vtop BC init, %i timesteps, vp band lim max: %g\n",
+ E->control.ggrd.time_hist.nvtimes,E->control.ggrd.svp->fmaxlim[0]);
}
/* end init */
}
- if((E->control.ggrd.time_hist.nvtimes > 1)||
- (!E->control.ggrd.vtop_control_init)){
+ if((E->control.ggrd.time_hist.nvtimes > 1)|| (!E->control.ggrd.vtop_control_init)){
+ /* either first time around, or time-dependent */
age = find_age_in_MY(E);
if(timedep){
/*
@@ -697,6 +704,7 @@
if(E->parallel.me == 0)
fprintf(stderr,"ggrd_read_vtop_from_file: interpolating vtop for age = %g\n",age);
}
+
}else{
interpolate = 0; /* single timestep, use single file */
i1 = 0;
@@ -705,19 +713,96 @@
}
if(E->parallel.me==0)
fprintf(stderr,"ggrd_read_vtop_from_file: assigning velocities BC, timedep: %i time: %g\n",
- timedep,age);
+ timedep,age);
+ /* if mixed BCs are allowed, need to reassign the boundary
+ condition */
+ if(E->control.ggrd_allow_mixed_vbcs){
+ if(E->parallel.me == 0)fprintf(stderr,"WARNING: allowing mixed velocity BCs\n");
+
+ if (E->parallel.me_loc[3] == E->parallel.nprocz-1) { /* top processor onlt */
+
+ cutoff = E->control.ggrd.svp->fmaxlim[0] + 1e-5;
+
+ for(level=E->mesh.gridmax;level>=E->mesh.gridmin;level--){
+ /* multigrid levels */
+ noxl = E->lmesh.NOX[level];noyl = E->lmesh.NOY[level];nozl = E->lmesh.NOZ[level];
+ noxlnozl = noxl*nozl;
+
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+ nsuml[0] = nsuml[1] = 0;
+ /*
+ loop through all horizontal nodes
+ */
+ for(i=1;i<=noyl;i++){
+ for(j=1;j<=noxl;j++) {
+
+ nodel = j * nozl + (i-1)*noxlnozl; /* top node */
+
+ rout[1] = E->SX[level][m][1][nodel]; /* theta,phi */
+ rout[2] = E->SX[level][m][2][nodel];
+ /* find vp */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i1),
+ vin1,FALSE)){
+ fprintf(stderr,"ggrd_read_vtop_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ if(interpolate){ /* second time */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i2),vin2,FALSE)){
+ fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ vp = (f1*vin1[0] + f2*vin2[0]); /* unscaled! */
+ }else{
+ vp = vin1[0];
+ }
+ if(fabs(vp) > cutoff){
+ /* free slip */
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~VBX);
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | SBX;
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~VBY);
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | SBY;
+ nsuml[0]++;
+ }else{
+ /* no slip */
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | VBX;
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~SBX);
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] | VBY;
+ E->NODE[level][m][nodel] = E->NODE[level][m][nodel] & (~SBY);
+ nsuml[1]++;
+ }
+ } /* end x loop */
+ } /* end y loop */
+ //fprintf(stderr,"cpu: %i - %i %i %i level: %i \t sum: %i %i\n",
+ // E->parallel.me,E->parallel.me_loc[1],E->parallel.me_loc[2] ,E->parallel.me_loc[3] ,level,nsuml[0],nsuml[1]);
+ MPI_Allreduce(nsuml,nsum,2,MPI_INT,MPI_SUM,E->parallel.horizontal_comm);
+ if(E->parallel.me == E->parallel.loc2proc_map[0][0][0][E->parallel.nprocz-1])
+ fprintf(stderr,"mixed velocity BC: multi grid level %2i : %7i fixed %7i free\n",level,nsum[1],nsum[0]);
+ } /* cap */
+ } /* level */
+ } /* proc */
+ /* end BC branch */
+ } /* end mixed BC assign */
/*
- loop through all elements and assign
+
+ loop through all nodes and assign velocity boundary condition
+ values
+
*/
- for (m=1;m <= E->sphere.caps_per_proc;m++) {
- if(E->parallel.me_loc[3] == E->parallel.nprocz-1 ) { /* if on top */
- for(k=1;k <= noy1;k++) {/* loop through surface nodes */
- for(i=1;i <= nox1;i++) {
- nodel = (k-1)*nox1noz1 + (i-1)*noz1+noz1;
+ if(E->parallel.me_loc[3] == E->parallel.nprocz-1 ) { /* if on top */
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+ nsuml[0] = nsuml[1] = 0;
+ cutoff = E->control.ggrd.svp->fmaxlim[0] * vscale + 1e-5;
+
+ for(k=1;k <= noyg;k++) {/* loop through surface nodes */
+ for(i=1;i <= noxg;i++) {
+
+ nodel = (k-1)*noxgnozg + i * nozg; /* top node */
+
rout[1] = E->sx[m][1][nodel]; /* theta,phi coordinates */
rout[2] = E->sx[m][2][nodel];
-
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svt+i1),
vin1,FALSE)){
fprintf(stderr,"ggrd_read_vtop_from_ggrd_file: interpolation error at %g, %g\n",
@@ -747,20 +832,35 @@
vt = vin1[0]*vscale;
vp = vin1[1]*vscale;
}
- /* assign as boundary condition */
- E->sphere.cap[m].VB[1][nodel] = vt; /* theta */
- E->sphere.cap[m].VB[2][nodel] = vp; /* phi */
+ if(fabs(vp) > cutoff){
+ /* huge velocitie - free slip */
+ E->sphere.cap[m].VB[1][nodel] = 0; /* theta */
+ E->sphere.cap[m].VB[2][nodel] = 0; /* phi */
+ nsuml[0]++;
+
+ }else{
+ /* regular no slip , assign velocities as BCs */
+ E->sphere.cap[m].VB[1][nodel] = vt; /* theta */
+ E->sphere.cap[m].VB[2][nodel] = vp; /* phi */
+ nsuml[1]++;
+ }
E->sphere.cap[m].VB[3][nodel] = 0.0; /* r */
}
} /* end surface node loop */
} /* end top proc branch */
} /* end cap loop */
+
+ MPI_Allreduce(nsuml,nsum,2,MPI_INT,MPI_SUM,E->parallel.horizontal_comm);
+ if(E->parallel.me == E->parallel.loc2proc_map[0][0][0][E->parallel.nprocz-1])
+ fprintf(stderr,"mixed velocity BC val: %i fixed %i free\n",nsum[1],nsum[0]);
} /* end assignment branch */
+
if((!timedep)&&(!E->control.ggrd.vtop_control_init)){ /* forget the grids */
ggrd_grdtrack_free_gstruc(E->control.ggrd.svt);
ggrd_grdtrack_free_gstruc(E->control.ggrd.svp);
}
E->control.ggrd.vtop_control_init = 1;
+ if(E->parallel.me == 0)fprintf(stderr,"vtop from grd done: %i\n",lc++);
}
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -53,10 +53,15 @@
int m = E->parallel.me;
int noz = E->lmesh.noz;
- int n;
+ int n,tmp;
input_int("tic_method", &(E->convection.tic_method), "0,0,2", m);
+#ifdef USE_GGRD /* for backward capability */
+ input_int("ggrd_tinit", &tmp, "0", m);
+ if(tmp)
+ E->convection.tic_method = 4; /* */
+#endif
/* When tic_method is 0 (default), the temperature is a linear profile +
perturbation at some layers.
@@ -151,7 +156,10 @@
case 4: initial temp from grd files
*/
#ifdef USE_GGRD
- /* read in some more parameters */
+ /*
+ read in some more parameters
+
+ */
/* scale the anomalies with PREM densities */
input_boolean("ggrd_tinit_scale_with_prem",&(E->control.ggrd.temp_init.scale_with_prem),"off",E->parallel.me);
/* limit T to 0...1 */
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -85,7 +85,9 @@
void initial_mesh_solver_setup(struct All_variables *E)
{
-
+ int chatty;
+ chatty = ((E->parallel.me == 0)&&(E->control.verbose))?(1):(0);
+
E->monitor.cpu_time_at_last_cycle =
E->monitor.cpu_time_at_start = CPU_time0();
@@ -99,18 +101,21 @@
allocate_common_vars(E);
(E->problem_allocate_vars)(E);
(E->solver_allocate_vars)(E);
-
+ if(chatty)fprintf(stderr,"memory allocation done\n");
/* logical domain */
construct_ien(E);
construct_surface(E);
(E->solver.construct_boundary)(E);
(E->solver.parallel_domain_boundary_nodes)(E);
+ if(chatty)fprintf(stderr,"parallel setup done\n");
/* physical domain */
(E->solver.node_locations)(E);
allocate_velocity_vars(E);
+ if(chatty)fprintf(stderr,"velocity vars done\n");
+
get_initial_elapsed_time(E); /* Set elapsed time */
set_starting_age(E); /* set the starting age to elapsed time, if desired */
set_elapsed_time(E); /* reset to elapsed time to zero, if desired */
@@ -131,13 +136,17 @@
(E->problem_boundary_conds)(E);
check_bc_consistency(E);
+ if(chatty)fprintf(stderr,"boundary conditions done\n");
construct_masks(E); /* order is important here */
construct_id(E);
construct_lm(E);
+ if(chatty)fprintf(stderr,"id/lm done\n");
(E->solver.parallel_communication_routs_v)(E);
+ if(chatty)fprintf(stderr,"v communications done\n");
(E->solver.parallel_communication_routs_s)(E);
+ if(chatty)fprintf(stderr,"s communications done\n");
reference_state(E);
@@ -156,13 +165,16 @@
construct_surf_det (E);
construct_bdry_det (E);
+ if(chatty)fprintf(stderr,"mass matrix, dets done\n");
+
set_sphere_harmonics (E);
+
if(E->control.tracer) {
tracer_initial_settings(E);
(E->problem_tracer_setup)(E);
}
-
+ if(chatty)fprintf(stderr,"initial_mesh_solver_setup done\n");
}
@@ -189,7 +201,6 @@
global_default_values(E);
read_initial_settings(E);
-
shutdown_parser(E);
return;
@@ -349,6 +360,8 @@
input_int("mgunitz",&(E->mesh.mgunitz),"1",m);
input_int("mgunity",&(E->mesh.mgunity),"1",m);
input_int("levels",&(E->mesh.levels),"0",m);
+ if(E->mesh.levels > MAX_LEVELS)
+ myerror(E,"number of multigrid levels out of bound");
input_int("coor",&(E->control.coor),"0",m);
if(E->control.coor == 2){
@@ -530,6 +543,13 @@
if(E->control.ggrd.vtop_control) /* this will override mat_control setting */
E->control.vbcs_file = 1;
+ /* if set, will check the theta velocities from grid input for
+ scaled (internal non dim) values of > 1e9. if found, those nodes will
+ be set to free slip
+ */
+ input_boolean("allow_mixed_vbcs",&(E->control.ggrd_allow_mixed_vbcs),"off",m);
+
+
#endif
input_int("nodex",&(E->mesh.nox),"essential",m);
@@ -606,6 +626,13 @@
input_float("density_below",&(E->data.density_below),"6600.0",m);
input_float("refvisc",&(E->data.ref_viscosity),"1.0e21",m);
+
+
+ /*
+
+ ellipticity and rotation settings
+
+ */
/* f = (a-c)/a, where c is the short, a=b the long axis
1/298.257 = 0.00335281317789691 for Earth at present day
*/
@@ -614,9 +641,12 @@
/* define ra and rc such that R=1 is the volume equivalanet */
E->data.ra = pow((1.-(double)E->data.ellipticity),(double)-1./3.); /* non dim long axis */
E->data.rc = 1./(E->data.ra * E->data.ra); /* non dim short axis */
- if(E->parallel.me == 0)
+ if(E->parallel.me == 0){
fprintf(stderr,"WARNING: ellipticity: %.5e equivalent radii: r_a: %g r_b: %g\n",
E->data.ellipticity,E->data.ra,E->data.rc);
+ if(E->control.remove_rigid_rotation)
+ fprintf(stderr,"WARNING: remove rigid rotations is switched on for elliptical run !!!\n");
+ }
}else{
E->data.ra = E->data.rc = 1.0;
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Regional_boundary_conditions.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -39,6 +39,7 @@
static void velocity_refl_vert_bc();
static void temperature_refl_vert_bc();
void read_temperature_boundary_from_file(struct All_variables *);
+void read_velocity_boundary_from_file(struct All_variables *);
/* ========================================== */
@@ -46,7 +47,6 @@
struct All_variables *E;
{
void velocity_imp_vert_bc();
- void read_velocity_boundary_from_file();
void renew_top_velocity_boundary();
void apply_side_sbc();
@@ -73,7 +73,8 @@
horizontal_bc(E,E->sphere.cap[j].VB,noz,2,0.0,SBY,0,lv,j);
if(E->control.vbcs_file) {
- read_velocity_boundary_from_file(E); /* read in the velocity boundary condition from file */
+ if((lv == E->mesh.gridmin) && (j == E->sphere.caps_per_proc))
+ read_velocity_boundary_from_file(E); /* read in the velocity boundary condition from file */
}
}
else if(E->mesh.topvbc == 2) {
Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-06-21 22:28:51 UTC (rev 12303)
@@ -48,7 +48,7 @@
float age, newage1, newage2;
char output_file1[255],output_file2[255];
float *TB1, *TB2, *VB1[4],*VB2[4], inputage1, inputage2;
- int nox,noz,noy,nnn,nox1,noz1,noy1,lev;
+ int nox,noz,noy,nnn,nox1,noz1,noy1;
int i,ii,ll,mm,j,k,n,nodeg,nodel,node;
int intage, pos_age;
int nodea;
@@ -73,8 +73,8 @@
nox1=E->lmesh.nox;
noz1=E->lmesh.noz;
noy1=E->lmesh.noy;
- lev=E->mesh.levmax;
+
elx=E->lmesh.elx;
elz=E->lmesh.elz;
ely=E->lmesh.ely;
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-06-20 23:01:28 UTC (rev 12302)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-06-21 22:28:51 UTC (rev 12303)
@@ -516,6 +516,7 @@
#ifdef USE_GGRD
struct ggrd_master ggrd;
float *surface_rayleigh;
+ int ggrd_allow_mixed_vbcs;
#endif
double accuracy;
double vaccuracy;
More information about the cig-commits
mailing list