[cig-commits] r14640 - in mc/3D/CitcomS/trunk: examples/Full examples/Regional lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Apr 8 16:34:44 PDT 2009
Author: tan2
Date: 2009-04-08 16:34:43 -0700 (Wed, 08 Apr 2009)
New Revision: 14640
Modified:
mc/3D/CitcomS/trunk/examples/Full/input.sample
mc/3D/CitcomS/trunk/examples/Regional/input.sample
mc/3D/CitcomS/trunk/lib/Convection.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/advdiffu.c
mc/3D/CitcomS/trunk/module/bindings.c
mc/3D/CitcomS/trunk/module/misc.c
mc/3D/CitcomS/trunk/module/misc.h
Log:
Clean up.
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2009-04-08 23:34:43 UTC (rev 14640)
@@ -98,7 +98,6 @@
# required information
Problem=convection
Geometry=sphere
-Spacing=regular
Solver=cgrad
node_assemble=1
Modified: mc/3D/CitcomS/trunk/examples/Regional/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Regional/input.sample 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/examples/Regional/input.sample 2009-04-08 23:34:43 UTC (rev 14640)
@@ -105,7 +105,6 @@
# required information
Problem=convection
Geometry=sphere
-Spacing=regular
Solver=cgrad
node_assemble=1
Modified: mc/3D/CitcomS/trunk/lib/Convection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Convection.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/lib/Convection.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -48,6 +48,8 @@
void convection_initial_fields();
void twiddle_thumbs();
+ E->control.CONVECTION = 1;
+
E->advection.timestep = 0.0;
E->advection.timesteps = 0;
E->advection.temp_iterations = 2; /* petrov-galerkin iterations: minimum value. */
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -291,68 +291,34 @@
/* first the problem type (defines subsequent behaviour) */
input_string("Problem",E->control.PROBLEM_TYPE,"convection",m);
- if ( strcmp(E->control.PROBLEM_TYPE,"convection") == 0) {
- E->control.CONVECTION = 1;
+ if ( strcmp(E->control.PROBLEM_TYPE,"convection") == 0)
set_convection_defaults(E);
- }
-
- else if ( strcmp(E->control.PROBLEM_TYPE,"convection-chemical") == 0) {
- E->control.CONVECTION = 1;
+ else if ( strcmp(E->control.PROBLEM_TYPE,"convection-chemical") == 0)
set_convection_defaults(E);
- }
-
else {
fprintf(E->fp,"Unable to determine problem type, assuming convection ... \n");
- E->control.CONVECTION = 1;
set_convection_defaults(E);
}
input_string("Geometry",E->control.GEOMETRY,"sphere",m);
- if ( strcmp(E->control.GEOMETRY,"cart2d") == 0)
- { E->control.CART2D = 1;
- (E->solver.set_2dc_defaults)(E);}
- else if ( strcmp(E->control.GEOMETRY,"axi") == 0)
- { E->control.AXI = 1;
- }
- else if ( strcmp(E->control.GEOMETRY,"cart2pt5d") == 0)
- { E->control.CART2pt5D = 1;
- (E->solver.set_2pt5dc_defaults)(E);}
- else if ( strcmp(E->control.GEOMETRY,"cart3d") == 0)
- { E->control.CART3D = 1;
- (E->solver.set_3dc_defaults)(E);}
- else if ( strcmp(E->control.GEOMETRY,"sphere") == 0)
- {
- (E->solver.set_3dsphere_defaults)(E);}
- else
- { fprintf(E->fp,"Unable to determine geometry, assuming cartesian 2d ... \n");
- E->control.CART2D = 1;
- (E->solver.set_2dc_defaults)(E); }
+ if ( strcmp(E->control.GEOMETRY,"sphere") == 0)
+ (E->solver.set_3dsphere_defaults)(E);
+ else {
+ fprintf(E->fp,"Unable to determine geometry, assuming sphere 3d ... \n");
+ (E->solver.set_3dsphere_defaults)(E);
+ }
input_string("Solver",E->control.SOLVER_TYPE,"cgrad",m);
if ( strcmp(E->control.SOLVER_TYPE,"cgrad") == 0)
set_cg_defaults(E);
else if ( strcmp(E->control.SOLVER_TYPE,"multigrid") == 0)
set_mg_defaults(E);
- else
- { if (E->parallel.me==0) fprintf(stderr,"Unable to determine how to solve, specify Solver=VALID_OPTION \n");
+ else {
+ if (E->parallel.me==0) fprintf(stderr,"Unable to determine how to solve, specify Solver=VALID_OPTION \n");
parallel_process_termination();
- }
+ }
- /* admin */
-
- input_string("Spacing",E->control.NODE_SPACING,"regular",m);
- if ( strcmp(E->control.NODE_SPACING,"regular") == 0)
- E->control.GRID_TYPE = 1;
- else if ( strcmp(E->control.NODE_SPACING,"bound_lyr") == 0)
- E->control.GRID_TYPE = 2;
- else if ( strcmp(E->control.NODE_SPACING,"region") == 0)
- E->control.GRID_TYPE = 3;
- else if ( strcmp(E->control.NODE_SPACING,"ortho_files") == 0)
- E->control.GRID_TYPE = 4;
- else
- { E->control.GRID_TYPE = 1; }
-
/* Information on which files to print, which variables of the flow to calculate and print.
Default is no information recorded (apart from special things for given applications.
*/
@@ -1176,8 +1142,6 @@
E->control.augmented_Lagr = 0;
E->control.augmented = 0.0;
- E->control.GRID_TYPE=1;
-
E->trace.fpt = NULL;
E->control.tracer = 0;
E->composition.on = 0;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -280,23 +280,22 @@
float **EEta;
int propogate;
{
- int m,i,j,k,l,z,jj,kk,imark;
- float zero,e_6,one,eta0,Tave,depth,temp,tempa,temp1,TT[9];
+ int m,i,k,l,z,jj,kk;
+ float zero,one,eta0,temp,tempa,TT[9];
float zzz,zz[9],dr;
- float visc1, visc2, tempa_exp;
+ float visc1, visc2;
const int vpts = vpoints[E->mesh.nsd];
const int ends = enodes[E->mesh.nsd];
const int nel = E->lmesh.nel;
- e_6 = 1.e-6;
one = 1.0;
zero = 0.0;
- imark = 0;
- /* consisntent handling : l is material number - 1 to allow
+ /* consistent handling : l is (material number - 1) to allow
addressing viscosity arrays, which are all 0...n-1 */
switch (E->viscosity.RHEOL) {
- case 1: /* eta = N_0 exp( E * (T_0 - T)) */
+ case 1:
+ /* eta = N_0 exp( E * (T_0 - T)) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
@@ -323,7 +322,8 @@
}
break;
- case 2: /* eta = N_0 exp(-T/T_0) */
+ case 2:
+ /* eta = N_0 exp(-T/T_0) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
@@ -350,8 +350,8 @@
}
break;
- case 3: /* eta = N_0 exp(E/(T+T_0) - E/(1+T_0)) */
-
+ case 3:
+ /* eta = N_0 exp(E/(T+T_0) - E/(1+T_0)) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
@@ -359,7 +359,6 @@
tempa = E->viscosity.N0[l] * E->VIP[m][i];
else
tempa = E->viscosity.N0[l];
- j = 0;
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
@@ -391,8 +390,6 @@
else
tempa = E->viscosity.N0[l];
- j = 0;
-
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
@@ -419,26 +416,22 @@
case 5:
- /* same as rheol 3, except alternative margin, VIP, formulation */
+ /* when mat_control=0, same as rheol 3,
+ when mat_control=1, applying viscosity cut-off before mat_control */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
tempa = E->viscosity.N0[l];
/* fprintf(stderr,"\nINSIDE visc_from_T, l=%d, tempa=%g",l+1,tempa);*/
- j = 0;
-
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- /* zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]); */
}
for(jj=1;jj<=vpts;jj++) {
temp=0.0;
- /* zzz=0.0; */
for(kk=1;kk<=ends;kk++) {
TT[kk]=max(TT[kk],zero);
temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- /* zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)]; */
}
if(E->control.mat_control==0)
@@ -466,13 +459,11 @@
break;
- case 6: /*
- like case 1, but allowing for depth-dependence if Z_0 != 0
-
- eta = N_0 exp(E(T_0-T) + (1-z) Z_0 )
+ case 6:
+ /* like case 1, but allowing for depth-dependence if Z_0 != 0
+ eta = N_0 exp(E(T_0-T) + (1-z) Z_0 )
+ */
- */
-
for(m=1;m <= E->sphere.caps_per_proc;m++)
for(i=1;i <= nel;i++) {
@@ -482,7 +473,6 @@
tempa = E->viscosity.N0[l] * E->VIP[m][i];
else
tempa = E->viscosity.N0[l];
- j = 0;
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
@@ -499,17 +489,44 @@
EEta[m][ (i-1)*vpts + jj ] = tempa*
exp( E->viscosity.E[l]*(E->viscosity.T[l] - temp) +
zzz * E->viscosity.Z[l]);
- //if(E->parallel.me == 0)
- // fprintf(stderr,"z %11g km mat %i N0 %11g T %11g T0 %11g E %11g Z %11g mat: %i log10(eta): %11g\n",
- // zzz *E->data.radius_km ,l+1,
- // tempa,temp,E->viscosity.T[l],E->viscosity.E[l], E->viscosity.Z[l],l+1,log10(EEta[m][ (i-1)*vpts + jj ]));
+ /*
+ if(E->parallel.me == 0)
+ fprintf(stderr,"z %11g km mat %i N0 %11g T %11g T0 %11g E %11g Z %11g mat: %i log10(eta): %11g\n",
+ zzz *E->data.radius_km ,l+1,
+ tempa,temp,E->viscosity.T[l],E->viscosity.E[l], E->viscosity.Z[l],l+1,log10(EEta[m][ (i-1)*vpts + jj ]));
+ */
}
}
break;
case 7:
+ /* The viscosity formulation (dimensional) is:
+ visc=visc0*exp[(Ea+p*Va)/(R*T)]
+ Typical values for dry upper mantle are:
+ Ea = 300 KJ/mol ; Va = 1.e-5 m^3/mol
+
+ T=DT*(T0+T');
+ where DT - temperature contrast (from Rayleigh number)
+ T' - nondimensional temperature;
+ T0 - nondimensional surface tempereture;
+
+ =>
+ visc = visc0 * exp{(Ea+p*Va) / [R*DT*(T0 + T')]}
+ = visc0 * exp{[Ea/(R*DT) + p*Va/(R*DT)] / (T0 + T')}
+
+ so:
+ E->viscosity.E = Ea/(R*DT);
+ (1-r) = p/(rho*g);
+ E->viscosity.Z = Va*rho*g/(R*DT);
+ E->viscosity.T = T0;
+
+ after normalizing visc=1 at T'=1 and r=r_CMB:
+ visc = visc0*exp{ [viscE + (1-r)*viscZ] / (viscT+T')
+ - [viscE + (1-r_CMB)*viscZ] / (viscT+1) }
+ */
+
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
@@ -519,8 +536,6 @@
else
tempa = E->viscosity.N0[l];
- j = 0;
-
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
@@ -534,32 +549,7 @@
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
}
- /* The viscosity formulation (dimensional) is:
- visc=visc0*exp[(Ea+p*Va)/(R*T)]
- Typical values for dry upper mantle are:
- Ea = 300 KJ/mol ; Va = 1.e-5 m^3/mol
-
- T=DT*(T0+T');
- where DT - temperature contrast (from Rayleigh number)
- T' - nondimensional temperature;
- T0 - nondimensional surface tempereture;
-
- =>
- visc = visc0 * exp{(Ea+p*Va) / [R*DT*(T0 + T')]}
- = visc0 * exp{[Ea/(R*DT) + p*Va/(R*DT)] / (T0 + T')}
-
- so:
- E->viscosity.E = Ea/(R*DT);
- (1-r) = p/(rho*g);
- E->viscosity.Z = Va*rho*g/(R*DT);
- E->viscosity.T = T0;
-
- after normalizing visc=1 at T'=1 and r=r_CMB:
- visc=visc0*exp{ [viscE + (1-r)*viscZ] / (viscT+T')
- - [viscE + (1-r_CMB)*viscZ] / (viscT+1) }
- */
-
EEta[m][ (i-1)*vpts + jj ] = tempa*
exp( (E->viscosity.E[l] + E->viscosity.Z[l]*zzz )
/ (E->viscosity.T[l] + temp)
@@ -570,20 +560,21 @@
}
break;
- case 8: /*
- eta0 = N_0 exp(E/(T+T_0) - E/(1+T_0))
+ case 8:
+ /*
+ eta0 = N_0 exp(E/(T+T_0) - E/(1+T_0))
- eta = eta0 if T < T_sol0 + 2(1-z)
- eta = ET_red*eta0 if T >= T_sol0 + 2(1-z)
+ eta = eta0 if T < T_sol0 + 2(1-z)
+ eta = ET_red*eta0 if T >= T_sol0 + 2(1-z)
- where z is normalized by layer
- thickness, and T_sol0 is something
- like 0.6, and ET_red = 0.1
+ where z is normalized by layer
+ thickness, and T_sol0 is something
+ like 0.6, and ET_red = 0.1
- (same as case 3, but for viscosity reduction)
+ (same as case 3, but for viscosity reduction)
+ */
- */
- dr = E->sphere.ro - E->sphere.ri;
+ dr = E->sphere.ro - E->sphere.ri;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i] - 1;
@@ -591,7 +582,6 @@
tempa = E->viscosity.N0[l] * E->VIP[m][i];
else
tempa = E->viscosity.N0[l];
- j = 0;
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2009-04-08 23:34:43 UTC (rev 14640)
@@ -410,7 +410,6 @@
int stokes;
int restart;
int post_p;
- int post_topo;
char GEOMETRY[20]; /* one of ... */
int CART2D;
@@ -419,14 +418,9 @@
int AXI;
char SOLVER_TYPE[20]; /* one of ... */
- int DIRECT;
int CONJ_GRAD;
int NMULTIGRID;
- int DIRECTII;
- char NODE_SPACING[20]; /* turns into ... */
- int GRID_TYPE;
-
int pseudo_free_surf;
int tracer;
Modified: mc/3D/CitcomS/trunk/module/advdiffu.c
===================================================================
--- mc/3D/CitcomS/trunk/module/advdiffu.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/module/advdiffu.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -99,7 +99,6 @@
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
- E->control.CONVECTION = 1;
set_convection_defaults(E);
Py_INCREF(Py_None);
Modified: mc/3D/CitcomS/trunk/module/bindings.c
===================================================================
--- mc/3D/CitcomS/trunk/module/bindings.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/module/bindings.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -48,11 +48,8 @@
METH_VARARGS,
pyCitcom_copyright__doc__},
+ /* from misc.h */
- /*////////////////////////////////////////////////////////////////////////
- // This section is for testing or temporatory implementation
- ////////////////////////////////////////////////////////////////////////*/
-
{pyCitcom_return1_test__name__,
pyCitcom_return1_test,
METH_VARARGS,
@@ -63,18 +60,6 @@
METH_VARARGS,
pyCitcom_CPU_time__doc__},
- {pyCitcom_read_instructions__name__,
- pyCitcom_read_instructions,
- METH_VARARGS,
- pyCitcom_read_instructions__doc__},
-
-
- /*////////////////////////////////////////////////////////////////////////
- // This section is for finished implementation
- ////////////////////////////////////////////////////////////////////////*/
-
- /* from misc.h */
-
{pyCitcom_citcom_init__name__,
pyCitcom_citcom_init,
METH_VARARGS,
Modified: mc/3D/CitcomS/trunk/module/misc.c
===================================================================
--- mc/3D/CitcomS/trunk/module/misc.c 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/module/misc.c 2009-04-08 23:34:43 UTC (rev 14640)
@@ -71,12 +71,7 @@
return Py_BuildValue("s", pyCitcom_copyright_note);
}
-/*////////////////////////////////////////////////////////////////////////
-// This section is for testing or temporatory implementation
-////////////////////////////////////////////////////////////////////////*/
-
-
char pyCitcom_return1_test__doc__[] = "";
char pyCitcom_return1_test__name__[] = "return1_test";
@@ -88,30 +83,6 @@
}
-char pyCitcom_read_instructions__doc__[] = "";
-char pyCitcom_read_instructions__name__[] = "read_instructions";
-
-PyObject * pyCitcom_read_instructions(PyObject *self, PyObject *args)
-{
- PyObject *obj;
- struct All_variables* E;
- char *filename;
-
- if (!PyArg_ParseTuple(args, "Os:read_instructions", &obj, &filename))
- return NULL;
-
- E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
-
- read_instructions(E, filename);
-
- /* test */
- fprintf(stderr,"output file prefix: %s\n", E->control.data_file);
-
- Py_INCREF(Py_None);
- return Py_None;
-}
-
-
char pyCitcom_CPU_time__doc__[] = "";
char pyCitcom_CPU_time__name__[] = "CPU_time";
@@ -121,11 +92,6 @@
}
-/*////////////////////////////////////////////////////////////////////////
-// This section is for finished implementation
-////////////////////////////////////////////////////////////////////////*/
-
-
void deleteE(struct All_variables *E)
{
free(E);
Modified: mc/3D/CitcomS/trunk/module/misc.h
===================================================================
--- mc/3D/CitcomS/trunk/module/misc.h 2009-04-08 23:31:11 UTC (rev 14639)
+++ mc/3D/CitcomS/trunk/module/misc.h 2009-04-08 23:34:43 UTC (rev 14640)
@@ -38,11 +38,6 @@
PyObject * pyCitcom_return1_test(PyObject *, PyObject *);
-extern char pyCitcom_read_instructions__name__[];
-extern char pyCitcom_read_instructions__doc__[];
-PyObject * pyCitcom_read_instructions(PyObject *, PyObject *);
-
-
extern char pyCitcom_CPU_time__name__[];
extern char pyCitcom_CPU_time__doc__[];
PyObject * pyCitcom_CPU_time(PyObject *, PyObject *);
More information about the CIG-COMMITS
mailing list