[cig-commits] r8095 - in mc/3D/CitcomS/trunk:
CitcomS/Components/Stokes_solver examples/Cookbook8
examples/Full examples/Regional lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Oct 10 13:11:01 PDT 2007
Author: tan2
Date: 2007-10-10 13:11:00 -0700 (Wed, 10 Oct 2007)
New Revision: 8095
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
mc/3D/CitcomS/trunk/examples/Full/input.sample
mc/3D/CitcomS/trunk/examples/Regional/input.sample
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Remove rigid body rotation from the velocity solution from global model.
A input parameter 'solver.vsolver.remove_rigid_rotation', default to on, indicate whether to remove
the rigid body rotation. For free-slip model, this parameter should be on. For model with imposed plate
velocity, it is advised to turn off this parameter.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2007-10-10 20:11:00 UTC (rev 8095)
@@ -106,6 +106,7 @@
validator=prop.choice(["cg", "bicg"]))
compress_iter_maxstep = prop.int("compress_iter_maxstep", default=100)
relative_err_accuracy = prop.float("relative_err_accuracy", default=0.001)
+ remove_rigid_rotation = prop.bool("remove_rigid_rotation", default=True)
# version
__id__ = "$Id$"
Modified: mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2007-10-10 20:11:00 UTC (rev 8095)
@@ -95,4 +95,6 @@
compress_iter_maxstep = 100
relative_err_accuracy = 0.001
+remove_rigid_rotation = on
+
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-10 20:11:00 UTC (rev 8095)
@@ -223,6 +223,7 @@
aug_lagr=on
aug_number=2.0e3
+remove_rigid_rotation=on
# Age information
start_age=4.0
Modified: mc/3D/CitcomS/trunk/examples/Regional/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Regional/input.sample 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/examples/Regional/input.sample 2007-10-10 20:11:00 UTC (rev 8095)
@@ -226,6 +226,7 @@
aug_lagr=on
aug_number=2.0e3
+remove_rigid_rotation=on
# Age information
start_age=4.0
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-10 20:11:00 UTC (rev 8095)
@@ -64,6 +64,7 @@
void assemble_forces();
void sphere_harmonics_layer();
void get_system_viscosity();
+ void remove_rigid_rot();
float vmag;
@@ -140,6 +141,10 @@
} /*end if SDEPV or PDEPV */
+ /* remove the rigid rotation component from the velocity solution */
+ if(E->sphere.caps == 12 && E->control.remove_rigid_rotation)
+ remove_rigid_rot(E);
+
return;
}
@@ -151,6 +156,7 @@
void assemble_forces_pseudo_surf();
void get_system_viscosity();
void std_timestep();
+ void remove_rigid_rot();
void get_STD_freesurf(struct All_variables *, float**);
float vmag;
@@ -229,6 +235,11 @@
} /*end if SDEPV or PDEPV */
E->monitor.topo_loop++;
}
+
+ /* remove the rigid rotation component from the velocity solution */
+ if(E->sphere.caps == 12 && E->control.remove_rigid_rotation)
+ remove_rigid_rot(E);
+
get_STD_freesurf(E,E->slice.freesurf);
return;
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-10-10 20:11:00 UTC (rev 8095)
@@ -787,3 +787,103 @@
return;
}
+
+void remove_rigid_rot(struct All_variables *E)
+{
+ void velo_from_element_d();
+ double myatan();
+ double wx, wy, wz, v_theta, v_phi;
+ double vx[9], vy[9], vz[9];
+ double r, t, f;
+
+ double exyz[4], fxyz[4];
+
+ int m, e, i, k, j, node;
+
+ const int nno = E->lmesh.nno;
+ const int ends = ENODES3D;
+ const int ppts = PPOINTS3D;
+ const int vpts = VPOINTS3D;
+ const int sphere_key = 1;
+ double VV[4][9];
+ double rot, fr, tr;
+
+ /* Note: no need to weight in rho(r) here. */
+ double moment_of_inertia = (8.0*M_PI/15.0)*
+ (pow(E->sphere.ro,(double)5.0) - pow(E->sphere.ri,(double)5.0));
+
+ /* compute and add angular momentum components */
+
+ exyz[1] = exyz[2] = exyz[3] = 0.0;
+ fxyz[1] = fxyz[2] = fxyz[3] = 0.0;
+
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+
+ for (e=1;e<=E->lmesh.nel;e++) {
+
+ t = E->eco[m][e].centre[1];
+ f = E->eco[m][e].centre[2];
+ r = E->eco[m][e].centre[3];
+
+ velo_from_element(E,VV,m,e,sphere_key);
+
+ for (j=1;j<=ppts;j++) {
+ vx[j] = 0.0;
+ vy[j] = 0.0;
+ }
+
+ for (j=1;j<=ppts;j++) {
+ for (i=1;i<=ends;i++) {
+ vx[j] += VV[1][i]*E->N.ppt[GNPINDEX(i,j)];
+ vy[j] += VV[2][i]*E->N.ppt[GNPINDEX(i,j)];
+ }
+ }
+
+ wx = -r*vy[1];
+ wy = r*vx[1];
+
+ exyz[1] += (wx*cos(t)*cos(f)-wy*sin(f)) * E->eco[m][e].area;
+ exyz[2] += (wx*cos(t)*sin(f)+wy*cos(f)) * E->eco[m][e].area;
+ exyz[3] -= (wx*sin(t) ) * E->eco[m][e].area;
+ }
+ } /* end cap */
+
+ MPI_Allreduce(exyz,fxyz,4,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+ fxyz[1] = fxyz[1] / moment_of_inertia;
+ fxyz[2] = fxyz[2] / moment_of_inertia;
+ fxyz[3] = fxyz[3] / moment_of_inertia;
+
+
+ rot = sqrt(fxyz[1]*fxyz[1] + fxyz[2]*fxyz[2] + fxyz[3]*fxyz[3]);
+ fr = myatan(fxyz[2], fxyz[1]);
+ tr = acos(fxyz[3] / rot);
+
+ if (E->parallel.me==0) {
+ fprintf(E->fp,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ }
+
+
+ /* remove rigid rotation */
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for (node=1;node<=nno;node++) {
+
+ v_theta = E->sx[m][3][node] * rot * sin(tr)
+ * sin(fr - E->sx[m][2][node]);
+ v_phi = E->sx[m][3][node] * rot
+ * ( sin(E->sx[m][1][node]) * cos(tr)
+ - cos(E->sx[m][1][node]) * sin(tr)
+ * cos(fr-E->sx[m][2][node]) );
+
+ E->sphere.cap[m].V[1][node] -= v_theta;
+ E->sphere.cap[m].V[2][node] -= v_phi;
+
+ }
+ }
+
+
+ return;
+
+}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-10 20:11:00 UTC (rev 8095)
@@ -416,6 +416,8 @@
input_boolean("aug_lagr",&(E->control.augmented_Lagr),"off",m);
input_double("aug_number",&(E->control.augmented),"0.0",m);
+ input_boolean("remove_rigid_rotation",&(E->control.remove_rigid_rotation),"on",m);
+
input_float("tole_compressibility",&(E->control.tole_comp),"0.0",m);
input_int("storage_spacing",&(E->control.record_every),"10",m);
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-10-10 20:11:00 UTC (rev 8095)
@@ -501,6 +501,7 @@
int up_heavy;
int verbose;
+ int remove_rigid_rotation;
int side_sbcs;
int vbcs_file;
int mat_control;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-10-09 12:35:28 UTC (rev 8094)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-10-10 20:11:00 UTC (rev 8095)
@@ -848,6 +848,8 @@
getIntProperty(properties, "aug_lagr", E->control.augmented_Lagr, fp);
getDoubleProperty(properties, "aug_number", E->control.augmented, fp);
+ getIntProperty(properties, "remove_rigid_rotation", E->control.remove_rigid_rotation, fp);
+
if(E->control.inv_gruneisen != 0) {
/* which compressible solver to use: "cg" or "bicg" */
getStringProperty(properties, "uzawa", E->control.uzawa, fp);
More information about the cig-commits
mailing list