[cig-commits] r16089 - in mc/3D/CitcomS/trunk: CitcomS/Components/Stokes_solver lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Dec 8 13:44:26 PST 2009
Author: tan2
Date: 2009-12-08 13:44:25 -0800 (Tue, 08 Dec 2009)
New Revision: 16089
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Added two boolean parameters for controlling convergence criterion:
check_continuity_convergence (default=on): include div/v in the criterion.
check_pressure_convergence (default=on): include dp/p in the criterion.
Retired the parameter: only_check_vel_convergence.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2009-12-08 21:43:38 UTC (rev 16088)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2009-12-08 21:44:25 UTC (rev 16089)
@@ -87,6 +87,9 @@
precond = prop.bool("precond", default=True)
accuracy = prop.float("accuracy", default=1.0e-4)
+ check_continuity_convergence = prop.bool("check_continuity_convergence", default=True)
+ check_pressure_convergence = prop.bool("check_pressure_convergence", default=True)
+
mg_cycle = prop.int("mg_cycle", default=1)
down_heavy = prop.int("down_heavy", default=3)
up_heavy = prop.int("up_heavy", default=3)
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2009-12-08 21:43:38 UTC (rev 16088)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2009-12-08 21:44:25 UTC (rev 16089)
@@ -609,9 +609,9 @@
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("check_continuity_convergence",&(E->control.check_continuity_convergence),"on",m);
+ input_boolean("check_pressure_convergence",&(E->control.check_pressure_convergence),"on",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);
input_int("max_mg_cycles",&(E->control.max_mg_cycles),"50,0,nomax",m);
@@ -749,10 +749,6 @@
myerror(E, "number of multigrid levels out of bound\n");
}
- if((E->parallel.me == 0) && (E->control.only_check_vel_convergence)) {
- fprintf(stderr,"solve_Ahat_p_fhat: WARNING: overriding pressure and div check\n");
- }
-
/* remove angular momentum/rigid rotation should only be done in free
convection of global models. */
if(E->sphere.caps == 12 &&
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2009-12-08 21:43:38 UTC (rev 16088)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2009-12-08 21:44:25 UTC (rev 16089)
@@ -148,7 +148,20 @@
}
+static int keep_iterating(struct All_variables *E,
+ double acc, int converging)
+{
+ const int required_converging_loops = 2;
+ if(E->control.check_continuity_convergence)
+ return (E->monitor.incompressibility > acc) ||
+ (converging < required_converging_loops);
+ else
+ return (E->monitor.incompressibility > acc) &&
+ (converging < required_converging_loops);
+}
+
+
static void solve_Ahat_p_fhat(struct All_variables *E,
double **V, double **P, double **F,
double imp, int *steps_max)
@@ -272,9 +285,7 @@
r0dotz0 = 0;
- while( (count < *steps_max) &&
- ((E->monitor.incompressibility > imp) ||
- (converging < 2) )) {
+ while( (count < *steps_max) && keep_iterating(E, imp, converging) ) {
/* require two consecutive converging iterations to quit the while-loop */
/* preconditioner BPI ~= inv(K), z1 = BPI*r1 */
@@ -364,25 +375,25 @@
dvelocity, dpressure,
E->monitor.incompressibility);
}
+
if(!valid){
- converging = 0;
+ /* reset consecutive converging iterations */
+ converging = 0;
}else{
-
-
- if(E->control.only_check_vel_convergence){
- /* disregard pressure and div check */
- if(dvelocity < imp)
- converging++;
- else
- converging = 0;
-
- }else{
- /* how many consecutive converging iterations? */
- if(dvelocity < imp && dpressure < imp)
- converging++;
- else
- converging = 0;
- }
+ /* how many consecutive converging iterations? */
+ if(E->control.check_pressure_convergence) {
+ /* check dv and dp */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+ }else{
+ /* check dv only */
+ if(dvelocity < imp)
+ converging++;
+ else
+ converging = 0;
+ }
}
@@ -528,9 +539,7 @@
valid = 1;
r0dotrt = alpha = omega = 0;
- while( (count < *steps_max) &&
- ((E->monitor.incompressibility > imp) ||
- (converging < 2) )) {
+ while( (count < *steps_max) && keep_iterating(E, imp, converging) ) {
/* require two consecutive converging iterations to quit the while-loop */
/* r1dotrt = <r1, rt> */
@@ -664,27 +673,28 @@
dvelocity, dpressure,
E->monitor.incompressibility);
}
+
if(!valid){
- converging = 0;
+ /* reset consecutive converging iterations */
+ converging = 0;
}else{
- if(E->control.only_check_vel_convergence){
- /*
-
- override pressure and compressibility check
-
- */
- if(dvelocity < imp)
- converging++;
- else
- converging =0;
- }else{
- /* how many consecutive converging iterations? */
- if(dvelocity < imp && dpressure < imp)
- converging++;
- else
- converging = 0;
- }
+ /* how many consecutive converging iterations? */
+ if(E->control.check_pressure_convergence) {
+ /* check dv and dp */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+ }else{
+ /* check dv only */
+ if(dvelocity < imp)
+ converging++;
+ else
+ converging = 0;
+ }
+
}
+
/* shift array pointers */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
shuffle[m] = p1[m];
@@ -794,10 +804,7 @@
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 2009-12-08 21:43:38 UTC (rev 16088)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2009-12-08 21:44:25 UTC (rev 16089)
@@ -530,7 +530,8 @@
ggrd_boolean ggrd_mat_is_3d;
#endif
double accuracy;
- int only_check_vel_convergence;
+ int check_continuity_convergence;
+ int check_pressure_convergence;
char velocity_boundary_file[1000];
char temperature_boundary_file[1000];
char mat_file[1000];
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2009-12-08 21:43:38 UTC (rev 16088)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2009-12-08 21:44:25 UTC (rev 16089)
@@ -881,6 +881,9 @@
getDoubleProperty(properties, "accuracy", E->control.accuracy, fp);
+ getIntProperty(properties, "check_continuity_convergence", E->control.check_continuity_convergence, fp);
+ getIntProperty(properties, "check_pressure_convergence", E->control.check_pressure_convergence, fp);
+
getIntProperty(properties, "mg_cycle", E->control.mg_cycle, fp);
getIntProperty(properties, "down_heavy", E->control.down_heavy, fp);
getIntProperty(properties, "up_heavy", E->control.up_heavy, fp);
More information about the CIG-COMMITS
mailing list