[cig-commits] r13327 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Mon Nov 17 22:55:47 PST 2008
Author: becker
Date: 2008-11-17 22:55:46 -0800 (Mon, 17 Nov 2008)
New Revision: 13327
Modified:
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Patched up Ggrd_handling for velocity grids close to pole.
Addded flag to suppress check of incompressibility and pressure convergence. Partially,
this is because for kinematic BCs pressure will only be constrained up to a constant, partially
because mixed density/plate flow models show very poor convergence, still to be checked.
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-11-18 06:02:17 UTC (rev 13326)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-11-18 06:55:46 UTC (rev 13327)
@@ -636,9 +636,9 @@
char gmt_string[10],char_dummy;
static int lc =0; /* only for debugging */
double vin1[2],vin2[2],age,f1,f2,vscale,rout[3],cutoff,v[3],sin_theta,vx[4],
- cos_theta,sin_phi,cos_phi;
+ cos_theta,sin_phi,cos_phi,theta_max,theta_min;
char tfilename1[1000],tfilename2[1000];
-
+ static pole_warned = FALSE;
static ggrd_boolean shift_to_pos_lon = FALSE;
const int dims=E->mesh.nsd;
#ifdef USE_GZDIR
@@ -758,7 +758,15 @@
fprintf(stderr,"ggrd_read_vtop_from_file: 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 */
-
+
+ /* geographic bounds */
+ theta_max = (90-E->control.ggrd.svp[i1].south)*M_PI/180-1e-5;
+ theta_min = (90-E->control.ggrd.svp[i1].north)*M_PI/180+1e-5;
+ if((E->parallel.me ==0) && (is_global)){
+ fprintf(stderr,"ggrd_read_vtop_from_file: determined South/North range: %g/%g\n",
+ E->control.ggrd.svp[i1].south,E->control.ggrd.svp[i1].north);
+ }
+
if((E->control.ggrd.time_hist.nvtimes > 1)|| (!E->control.ggrd.vtop_control_init)){
/*
either first time around, or time-dependent assignment
@@ -789,10 +797,12 @@
if(E->parallel.me == 0)
fprintf(stderr,"ggrd_read_vtop_from_file: temporally constant velocity BC \n");
}
- if(E->parallel.me==0)
+
+ if(E->parallel.me==0){
fprintf(stderr,"ggrd_read_vtop_from_file: assigning velocities BC, timedep: %i time: %g\n",
timedep,age);
-
+
+ }
/* if mixed BCs are allowed, need to reassign the boundary
condition */
if(E->control.ggrd_allow_mixed_vbcs){
@@ -825,6 +835,28 @@
nodel = j * nozl + (i-1)*noxlnozl; /* top node = nozl + (j-1) * nozl + (i-1)*noxlnozl; */
/* node location */
rout[1] = E->SX[level][m][1][nodel]; /* theta,phi */
+ /*
+
+ for global grid, shift theta if too close to poles
+
+ */
+ if((is_global)&&(rout[1] > theta_max)){
+ if(!pole_warned){
+ fprintf(stderr,"WARNING: shifting theta from %g (%g) to max theta %g (%g)\n",
+ rout[1],90-180/M_PI*rout[1],theta_max,90-180/M_PI*theta_max);
+ pole_warned = TRUE;
+ }
+ rout[1] = theta_max;
+ }
+ if((is_global)&&(rout[1] < theta_min)){
+ if(!pole_warned){
+ fprintf(stderr,"WARNING: shifting theta from %g (%g) to min theta %g (%g)\n",
+ rout[1],90-180/M_PI*rout[1],theta_min,90-180/M_PI*theta_min);
+ pole_warned = TRUE;
+ }
+ rout[1] = theta_min;
+ }
+ /* */
rout[2] = E->SX[level][m][2][nodel];
/* find vp */
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i1),
@@ -884,6 +916,30 @@
nodel = (k-1)*noxgnozg + i * nozg; /* top node = nozg + (i-1) * nozg + (k-1)*noxgnozg */
/* */
rout[1] = E->sx[m][1][nodel]; /* theta,phi coordinates */
+
+ /*
+
+ for global grid, shift theta if too close to poles
+
+ */
+ if((is_global)&&(rout[1] > theta_max)){
+ if(!pole_warned){
+ fprintf(stderr,"WARNING: shifting theta from %g (%g) to max theta %g (%g)\n",
+ rout[1],90-180/M_PI*rout[1],theta_max,90-180/M_PI*theta_max);
+ pole_warned = TRUE;
+ }
+ rout[1] = theta_max;
+ }
+ if((is_global)&&(rout[1] < theta_min)){
+ if(!pole_warned){
+ fprintf(stderr,"WARNING: shifting theta from %g (%g) to min theta %g (%g)\n",
+ rout[1],90-180/M_PI*rout[1],theta_min,90-180/M_PI*theta_min);
+ pole_warned = TRUE;
+ }
+ rout[1] = theta_min;
+ }
+
+
rout[2] = E->sx[m][2][nodel];
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svt+i1),
vin1,FALSE,shift_to_pos_lon)){
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-11-18 06:02:17 UTC (rev 13326)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-11-18 06:55:46 UTC (rev 13327)
@@ -642,6 +642,7 @@
input_int("down_heavy",&(E->control.down_heavy),"1,0,nomax",m);
input_int("up_heavy",&(E->control.up_heavy),"1,0,nomax",m);
input_double("accuracy",&(E->control.accuracy),"1.0e-4,0.0,1.0",m);
+ input_boolean("only_check_vel_convergence",&(E->control.only_check_vel_convergence),"off",m);
input_int("vhighstep",&(E->control.v_steps_high),"1,0,nomax",m);
input_int("vlowstep",&(E->control.v_steps_low),"250,0,nomax",m);
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2008-11-18 06:02:17 UTC (rev 13326)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2008-11-18 06:55:46 UTC (rev 13327)
@@ -207,7 +207,7 @@
double v_norm, p_norm;
double dvelocity, dpressure;
int converging;
-
+ static int warned = FALSE;
void assemble_c_u();
void assemble_div_u();
void assemble_del2_u();
@@ -216,6 +216,11 @@
int solve_del2_u();
void parallel_process_termination();
+ if((E->parallel.me == 0) && (E->control.only_check_vel_convergence) && (!warned)){
+ fprintf(stderr,"solve_Ahat_p_fhat_BiCG: WARNING: overriding pressure and div check\n");
+ warned = TRUE;
+ }
+
npno = E->lmesh.npno;
neq = E->lmesh.neq;
lev = E->mesh.levmax;
@@ -361,11 +366,7 @@
dvelocity = alpha * sqrt(global_v_norm2(E, E->u1) / (1e-32 + E->monitor.vdotv));
dpressure = alpha * sqrt(global_p_norm2(E, s2) / (1e-32 + E->monitor.pdotp));
- /* how many consecutive converging iterations? */
- if(dvelocity < imp && dpressure < imp)
- converging++;
- else
- converging = 0;
+
assemble_div_u(E, V, z1, lev);
if(E->control.inv_gruneisen != 0)
@@ -385,8 +386,23 @@
dvelocity, dpressure,
E->monitor.incompressibility);
}
+ if(E->control.only_check_vel_convergence){
+ /* disregard pressure and div check */
+ if(dvelocity < imp)
+ converging++;
+ else
+ converging = 0;
+ E->monitor.incompressibility = dvelocity;
+ }else{
+ /* how many consecutive converging iterations? */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+ }
+
/* shift array pointers */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
shuffle[m] = s1[m];
@@ -445,6 +461,8 @@
double global_v_norm2(), global_p_norm2(), global_div_norm2();
double CPU_time0();
+ static int warned = FALSE;
+
int npno, neq;
int m, j, count, lev;
int valid;
@@ -463,6 +481,11 @@
double time0, v_res;
+ if((E->parallel.me == 0) && (E->control.only_check_vel_convergence) && (!warned)){
+ fprintf(stderr,"solve_Ahat_p_fhat_BiCG: WARNING: overriding pressure and div check\n");
+ warned = TRUE;
+ }
+
npno = E->lmesh.npno;
neq = E->lmesh.neq;
lev = E->mesh.levmax;
@@ -647,17 +670,16 @@
dvelocity = sqrt(global_v_norm2(E, F) / (1e-32 + E->monitor.vdotv));
dpressure = sqrt(global_p_norm2(E, s0) / (1e-32 + E->monitor.pdotp));
- /* how many consecutive converging iterations? */
- if(dvelocity < imp && dpressure < imp)
- converging++;
- else
- converging = 0;
+
+
assemble_div_rho_u(E, V, t0, lev);
E->monitor.incompressibility = sqrt(global_div_norm2(E, t0)
/ (1e-32 + E->monitor.vdotv));
+
+
count++;
if(E->control.print_convergence && E->parallel.me==0) {
@@ -667,7 +689,24 @@
E->monitor.incompressibility);
}
+ if(E->control.only_check_vel_convergence){
+ /*
+ override pressure and compressibility check
+
+ */
+ if(dvelocity < imp)
+ converging++;
+ else
+ converging =0;
+ E->monitor.incompressibility = dvelocity;
+ }else{
+ /* how many consecutive converging iterations? */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+ }
/* shift array pointers */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
shuffle[m] = p1[m];
@@ -721,7 +760,7 @@
double relative_err_v, relative_err_p;
double *old_v[NCS], *old_p[NCS],*diff_v[NCS],*diff_p[NCS];
double div_res;
-
+ static int warned = FALSE;
const int npno = E->lmesh.npno;
const int neq = E->lmesh.neq;
const int lev = E->mesh.levmax;
@@ -729,6 +768,10 @@
double global_v_norm2(),global_p_norm2();
double global_div_norm2();
void assemble_div_rho_u();
+ if((E->parallel.me == 0) && (E->control.only_check_vel_convergence) && (!warned)){
+ fprintf(stderr,"solve_Ahat_p_fhat_iterCG: WARNING: overriding pressure and div check\n");
+ warned = TRUE;
+ }
for (m=1;m<=E->sphere.caps_per_proc;m++) {
old_v[m] = (double *)malloc(neq*sizeof(double));
@@ -778,7 +821,10 @@
fprintf(stderr, "itercg -- div(rho*v)/v=%.2e dv/v=%.2e and dp/p=%.2e loop %d\n\n", div_res, relative_err_v, relative_err_p, num_of_loop);
fprintf(E->fp, "itercg -- div(rho*v)/v=%.2e dv/v=%.2e and dp/p=%.2e loop %d\n\n", div_res, relative_err_v, relative_err_p, num_of_loop);
}
-
+ if(E->control.only_check_vel_convergence){
+ /* override pressure and compressibility check */
+ relative_err_p = div_res = relative_err_v;
+ }
num_of_loop++;
} /* end of while */
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-11-18 06:02:17 UTC (rev 13326)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-11-18 06:55:46 UTC (rev 13327)
@@ -525,6 +525,7 @@
float ggrd_vtop_omega[4];
#endif
double accuracy;
+ int only_check_vel_convergence;
char velocity_boundary_file[1000];
char temperature_boundary_file[1000];
char mat_file[1000];
More information about the CIG-COMMITS
mailing list