[cig-commits] r14358 - in mc/3D/CitcomS/trunk: CitcomS/Components/Stokes_solver examples/Cookbook8 lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Mar 16 17:16:11 PDT 2009
Author: tan2
Date: 2009-03-16 17:16:11 -0700 (Mon, 16 Mar 2009)
New Revision: 14358
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
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/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
remove angular momentum from stokes solution by setting solver.vsolver.remove_angular_mementum=1
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2009-03-17 00:16:11 UTC (rev 14358)
@@ -105,6 +105,7 @@
validator=prop.choice(["cg", "bicg"]))
compress_iter_maxstep = prop.int("compress_iter_maxstep", default=100)
remove_rigid_rotation = prop.bool("remove_rigid_rotation", default=True)
+ remove_angular_momentum = prop.bool("remove_angular_momentum", default=True)
# Not used. Retained here for backward compatibility.
tole_compressibility = prop.float("tole_compressibility", default=1.0e-7)
Modified: mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2009-03-17 00:16:11 UTC (rev 14358)
@@ -93,6 +93,6 @@
accuracy = 0.001
compress_iter_maxstep = 100
-remove_rigid_rotation = on
+remove_angular_momentum = on
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -144,7 +144,8 @@
} /*end if SDEPV or PDEPV */
/* remove the rigid rotation component from the velocity solution */
- if((E->sphere.caps == 12) && E->control.remove_rigid_rotation) {
+ if((E->sphere.caps == 12) &&
+ (E->control.remove_rigid_rotation || E->control.remove_angular_momentum)) {
remove_rigid_rot(E);
}
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -881,7 +881,7 @@
/*
- * remove rigid body rotation from the velocity
+ * remove rigid body rotation or angular momentum from the velocity
*/
void remove_rigid_rot(struct All_variables *E)
@@ -903,12 +903,23 @@
const int sphere_key = 1;
double VV[4][9];
double rot, fr, tr;
+ double tmp, moment_of_inertia, rho;
+ if(E->control.remove_angular_momentum) {
+ moment_of_inertia = tmp = 0;
+ for (i=1;i<=E->lmesh.elz;i++)
+ tmp += (8.0*M_PI/15.0)*
+ 0.5*(E->refstate.rho[i] + E->refstate.rho[i+1])*
+ (pow(E->sx[1][3][i+1],5.0) - pow(E->sx[1][3][i],5.0));
- /* 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));
+ MPI_Allreduce(&tmp, &moment_of_inertia, 1, MPI_DOUBLE,
+ MPI_SUM, E->parallel.vertical_comm);
+ } else {
+ /* no need to weight in rho(r) here. */
+ moment_of_inertia = (8.0*M_PI/15.0)*
+ (pow(E->sphere.ro,5.0) - pow(E->sphere.ri,5.0));
+ }
/* compute and add angular momentum components */
@@ -940,11 +951,19 @@
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;
+
+ wx = -r*vy[1];
+ wy = r*vx[1];
+
+ if(E->control.remove_angular_momentum) {
+ int nz = (e-1) % E->lmesh.elz + 1;
+ rho = 0.5 * (E->refstate.rho[nz] + E->refstate.rho[nz+1]);
+ } else {
+ rho = 1;
+ }
+ exyz[1] += (wx*cos_t*cos_f - wy*sin_f) * E->eco[m][e].area * rho;
+ exyz[2] += (wx*cos_t*sin_f + wy*cos_f) * E->eco[m][e].area * rho;
+ exyz[3] -= (wx*sin_t ) * E->eco[m][e].area * rho;
}
} /* end cap */
@@ -959,8 +978,13 @@
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);
- fprintf(stderr,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ if(E->control.remove_angular_momentum) {
+ fprintf(E->fp,"Angular momentum: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ fprintf(stderr,"Angular momentum: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ } else {
+ fprintf(E->fp,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ fprintf(stderr,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ }
}
/*
remove rigid rotation
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -633,6 +633,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("remove_angular_momentum",&(E->control.remove_angular_momentum),"on",m);
input_boolean("self_gravitation",&(E->control.self_gravitation),"off",m);
input_boolean("use_cbf_topo",&(E->control.use_cbf_topo),"off",m); /* make default on later XXX TWB */
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -477,7 +477,7 @@
}
if(E->output.gzdir.rnr){ /* remove the whole model net rotation */
- if((E->control.remove_rigid_rotation)&&
+ if((E->control.remove_rigid_rotation || E->control.remove_angular_momentum) &&
(E->parallel.me == 0)) /* that's not too terrible but wastes time */
fprintf(stderr,"WARNING: both gzdir.rnr and remove_rigid_rotation are switched on!\n");
oamp = determine_model_net_rotation(E,omega);
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -903,7 +903,10 @@
fptracer=fopen(E->trace.tracer_file,"r");
fgets(input_s,200,fptracer);
- sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns);
+ if(sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns) != 2) {
+ fprintf(stderr,"Error while reading file '%s'\n", E->trace.tracer_file);
+ exit(8);
+ }
fprintf(E->trace.fpt,"%d Tracers, %d columns in file \n",
number_of_tracers, ncolumns);
@@ -1074,9 +1077,13 @@
for(j=1;j<=E->sphere.caps_per_proc;j++) {
fgets(input_s,200,fp1);
- sscanf(input_s,"%d %d %d %lf",
- &idum1, &numtracers, &ncolumns, &rdum1);
+ if(sscanf(input_s,"%d %d %d %lf",
+ &idum1, &numtracers, &ncolumns, &rdum1) != 4) {
+ fprintf(stderr,"Error while reading file '%s'\n", output_file);
+ exit(8);
+ }
+
/* some error control */
if (E->trace.number_of_extra_quantities+3 != ncolumns) {
fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-wrong # of columns\n");
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2009-03-17 00:16:11 UTC (rev 14358)
@@ -514,6 +514,8 @@
int verbose;
int remove_rigid_rotation;
+ int remove_angular_momentum;
+
int side_sbcs;
int vbcs_file;
int tbcs_file;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2009-03-16 23:22:25 UTC (rev 14357)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2009-03-17 00:16:11 UTC (rev 14358)
@@ -861,6 +861,7 @@
getDoubleProperty(properties, "aug_number", E->control.augmented, fp);
getIntProperty(properties, "remove_rigid_rotation", E->control.remove_rigid_rotation, fp);
+ getIntProperty(properties, "remove_angular_momentum", E->control.remove_angular_momentum, fp);
if(E->control.inv_gruneisen != 0) {
/* which compressible solver to use: "cg" or "bicg" */
More information about the CIG-COMMITS
mailing list