[cig-commits] r16138 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Fri Jan 15 18:47:04 PST 2010
Author: becker
Date: 2010-01-15 18:47:04 -0800 (Fri, 15 Jan 2010)
New Revision: 16138
Modified:
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Implementation of two, perhaps entirely unnecessary features, in
search of why certain prescribed velocity models do not converge like
they should.
- allow for each inner loop solution to have a higher accuracy than
the outer, Uzawa loop. This is set by new parameter
"inner_accuracy_scale", which pre-multiplies the inner accuracy (imp)
setting (defaults to unity)
- remove any rigid body rotation at each Uzawa iteration
step. Parmaeter "inner_remove_rigid_rotation", by default off. For
this, I also added an experimental routine assign_v_to_vector() in
order to reasign the V velocities to the U solution vector after
removing the net rotations from V (how is supposed to work for the
regular operational mode?)
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-01-15 22:08:42 UTC (rev 16137)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-01-16 02:47:04 UTC (rev 16138)
@@ -585,6 +585,7 @@
input_double("aug_number",&(E->control.augmented),"0.0",m);
input_boolean("remove_rigid_rotation",&(E->control.remove_rigid_rotation),"on",m);
+ input_boolean("inner_remove_rigid_rotation",&(E->control.inner_remove_rigid_rotation),"off",m);
input_boolean("remove_angular_momentum",&(E->control.remove_angular_momentum),"on",m);
input_boolean("self_gravitation",&(E->control.self_gravitation),"off",m);
@@ -608,6 +609,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_double("inner_accuracy_scale",&(E->control.inner_accuracy_scale),"1.0,0.000001,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);
Modified: mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-01-15 22:08:42 UTC (rev 16137)
+++ mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-01-16 02:47:04 UTC (rev 16138)
@@ -58,6 +58,22 @@
return;
}
+void assign_v_to_vector(E)
+ struct All_variables *E;
+{
+ int m,node;
+ const int nno = E->lmesh.nno;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(node=1;node<=nno;node++) {
+ E->U[m][E->id[m][node].doff[1]] = E->sphere.cap[m].V[1][node];
+ E->U[m][E->id[m][node].doff[2]] = E->sphere.cap[m].V[2][node];
+ E->U[m][E->id[m][node].doff[3]] = E->sphere.cap[m].V[3][node];
+ }
+ }
+ return;
+}
+
void v_from_vector_pseudo_surf(E)
struct All_variables *E;
{
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2010-01-15 22:08:42 UTC (rev 16137)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2010-01-16 02:47:04 UTC (rev 16138)
@@ -196,7 +196,7 @@
double *shuffle[NCS];
double alpha, delta, r0dotz0, r1dotz1;
double v_res;
-
+ double inner_imp;
double global_pdot();
double global_v_norm2(), global_p_norm2(), global_div_norm2();
@@ -211,7 +211,12 @@
void strip_bcs_from_residual();
int solve_del2_u();
void parallel_process_termination();
+ void v_from_vector();
+ void v_from_vector_pseudo_surf();
+ void assign_v_to_vector();
+ inner_imp = imp * E->control.inner_accuracy_scale; /* allow for different innner loop accuracy */
+
npno = E->lmesh.npno;
neq = E->lmesh.neq;
lev = E->mesh.levmax;
@@ -251,7 +256,7 @@
/* In the compressible case, the initial guess of P might be bad.
* Do not correct V with it. */
if(E->control.inv_gruneisen == 0)
- initial_vel_residual(E, V, P, F, imp*v_res);
+ initial_vel_residual(E, V, P, F, inner_imp*v_res);
/* initial residual r1 = div(V) */
@@ -313,7 +318,7 @@
/* solve K*u1 = grad(s2) for u1 */
assemble_grad_p(E, s2, F, lev);
- valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ valid = solve_del2_u(E, E->u1, F, inner_imp*v_res, lev);
if(!valid && (E->parallel.me==0)) {
fputs("Warning: solver not converging! 1\n", stderr);
fputs("Warning: solver not converging! 1\n", E->fp);
@@ -410,7 +415,16 @@
/* shift <r0, z0> = <r1, z1> */
r0dotz0 = r1dotz1;
-
+ if((E->sphere.caps == 12) && (E->control.inner_remove_rigid_rotation)){
+ /* allow for removal of net rotation at each iterative step
+ (expensive) */
+ if(E->control.pseudo_free_surf) /* move from U to V */
+ v_from_vector_pseudo_surf(E);
+ else
+ v_from_vector(E);
+ remove_rigid_rot(E); /* correct V */
+ assign_v_to_vector(E); /* assign V to U */
+ }
} /* end loop for conjugate gradient */
assemble_div_u(E, V, z1, lev);
@@ -459,7 +473,7 @@
int m, j, count, lev;
int valid;
- double alpha, beta, omega;
+ double alpha, beta, omega,inner_imp;
double r0dotrt, r1dotrt;
double v_norm, p_norm;
double dvelocity, dpressure;
@@ -472,6 +486,8 @@
double *shuffle[NCS];
double time0, v_res;
+
+ inner_imp = imp * E->control.inner_accuracy_scale; /* allow for different innner loop accuracy */
npno = E->lmesh.npno;
neq = E->lmesh.neq;
@@ -505,7 +521,7 @@
/* calculate the initial velocity residual */
- initial_vel_residual(E, V, P, F, imp*v_res);
+ initial_vel_residual(E, V, P, F, inner_imp*v_res);
/* initial residual r1 = div(rho_ref*V) */
@@ -575,7 +591,7 @@
/* solve K*u0 = grad(pt) for u1 */
assemble_grad_p(E, pt, F, lev);
- valid = solve_del2_u(E, u0, F, imp*v_res, lev);
+ valid = solve_del2_u(E, u0, F, inner_imp*v_res, lev);
if(!valid && (E->parallel.me==0)) {
fputs("Warning: solver not converging! 1\n", stderr);
fputs("Warning: solver not converging! 1\n", E->fp);
@@ -605,7 +621,7 @@
/* solve K*u1 = grad(st) for u1 */
assemble_grad_p(E, st, F, lev);
- valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ valid = solve_del2_u(E, E->u1, F, inner_imp*v_res, lev);
if(!valid && (E->parallel.me==0)) {
fputs("Warning: solver not converging! 2\n", stderr);
fputs("Warning: solver not converging! 2\n", E->fp);
@@ -755,7 +771,7 @@
double global_v_norm2(),global_p_norm2();
double global_div_norm2();
void assemble_div_rho_u();
-
+
for (m=1;m<=E->sphere.caps_per_proc;m++) {
old_v[m] = (double *)malloc(neq*sizeof(double));
diff_v[m] = (double *)malloc(neq*sizeof(double));
@@ -766,7 +782,7 @@
cycles = E->control.p_iterations;
initial_vel_residual(E, V, P, F,
- imp * E->monitor.fdotf);
+ imp * E->control.inner_accuracy_scale * E->monitor.fdotf);
div_res = 1.0;
relative_err_v = 1.0;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-01-15 22:08:42 UTC (rev 16137)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-01-16 02:47:04 UTC (rev 16138)
@@ -876,7 +876,7 @@
where \eta_0 is the regular viscosity
- psrw = TRUE
+ psrw = 1
a strain-rate weakening rheology applies
based on a steady state stress-strain
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2010-01-15 22:08:42 UTC (rev 16137)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-01-16 02:47:04 UTC (rev 16138)
@@ -513,7 +513,7 @@
int up_heavy;
int verbose;
- int remove_rigid_rotation;
+ int remove_rigid_rotation,inner_remove_rigid_rotation;
int remove_angular_momentum;
int side_sbcs;
@@ -529,7 +529,7 @@
char ggrd_mat_depth_file[1000];
ggrd_boolean ggrd_mat_is_3d;
#endif
- double accuracy;
+ double accuracy,inner_accuracy_scale;
int check_continuity_convergence;
int check_pressure_convergence;
char velocity_boundary_file[1000];
More information about the CIG-COMMITS
mailing list