[cig-commits] r7895 - in mc/3D/CitcomS/trunk: CitcomS/Components
lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Aug 27 12:34:27 PDT 2007
Author: tan2
Date: 2007-08-27 12:34:26 -0700 (Mon, 27 Aug 2007)
New Revision: 7895
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
mc/3D/CitcomS/trunk/CitcomS/Components/Visc.py
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Add pyre input for new parameters introduced by TWB
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2007-08-27 19:34:08 UTC (rev 7894)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2007-08-27 19:34:26 UTC (rev 7895)
@@ -50,10 +50,16 @@
output_format = inv.str("output_format", default="ascii",
validator=inv.choice(["ascii",
+ "ascii-gz",
"hdf5"]))
output_optional = inv.str("output_optional",
default="surf,botm,tracer,comp_el")
+ # experimental vtk output
+ gzdir_vtkio = inv.int("gzdir_vtkio", default=0)
+ # remove net rotation
+ gzdir_rnr = inv.bool("gzdir_rnr", default=False)
+
# max. degree for spherical harmonics output
output_ll_max = inv.int("output_ll_max", default=20)
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-08-27 19:34:08 UTC (rev 7894)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-08-27 19:34:26 UTC (rev 7895)
@@ -82,9 +82,16 @@
# later for many purposes. One of it is to compute composition,
# either using absolute method or ratio method.
tracer_flavors = inv.int("tracer_flavors", default=0)
+
+ # How to initialize tracer flavors
ic_method_for_flavors = inv.int("ic_method_for_flavors", default=0)
z_interface = inv.list("z_interface", default=[0.7])
+ ictracer_grd_file = inv.str("ictracer_grd_file", default="")
+ ictracer_grd_layers = inv.int("ictracer_grd_layers", default=2)
+ # Enriched internal heat production
+ tracer_enriched = inv.bool("tracer_enriched", default=False)
+ Q0_enriched = inv.float("Q0_enriched", default=0.0)
# Regular grid parameters
regular_grid_deltheta = inv.float("regular_grid_deltheta", default=1.0)
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Visc.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Visc.py 2007-08-27 19:34:08 UTC (rev 7894)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Visc.py 2007-08-27 19:34:26 UTC (rev 7895)
@@ -47,6 +47,10 @@
inv.viscT = map(float, inv.viscT)
inv.viscZ = map(float, inv.viscZ)
inv.sdepv_expt = map(float, inv.sdepv_expt)
+ inv.pdepv_a = map(float, inv.pdepv_a)
+ inv.pdepv_b = map(float, inv.pdepv_b)
+ inv.pdepv_y = map(float, inv.pdepv_y)
+ inv.cdepv_ff = map(float, inv.cdepv_ff)
Visc_set_properties(self.all_variables, inv, stream)
@@ -82,8 +86,18 @@
SDEPV = pyre.inventory.bool("SDEPV", default=False)
sdepv_expt = pyre.inventory.list("sdepv_expt", default=[1, 1, 1, 1])
- sdepv_misfit = pyre.inventory.float("sdepv_misfit", default=0.02)
+ sdepv_misfit = pyre.inventory.float("sdepv_misfit", default=0.001)
+ PDEPV = pyre.inventory.bool("PDEPV", default=False)
+ pdepv_eff = pyre.inventory.bool("pdepv_eff", default=True)
+ pdepv_a= pyre.inventory.list("pdepv_a", default=[1e20, 1e20, 1e20, 1e20])
+ pdepv_b= pyre.inventory.list("pdepv_b", default=[0, 0, 0, 0])
+ pdepv_y= pyre.inventory.list("pdepv_y", default=[1e20, 1e20, 1e20, 1e20])
+ pdepv_offset = pyre.inventory.float("pdepv_offset", default=0.0)
+
+ CDEPV = pyre.inventory.bool("CDEPV", default=False)
+ cdepv_ff= pyre.inventory.list("cdepv_ff", default=[1])
+
low_visc_channel = pyre.inventory.bool("low_visc_channel", default=False)
low_visc_wedge = pyre.inventory.bool("low_visc_wedge", default=False)
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-27 19:34:08 UTC (rev 7894)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-27 19:34:26 UTC (rev 7895)
@@ -399,8 +399,9 @@
const int ends=enodes[dims];
const int sphere_key = 1;
- /* adiabatic and dissipative heating*/
- process_heating(E);
+ /* adiabatic, dissipative and latent heating*/
+ if(E->control.disptn_number != 0)
+ process_heating(E);
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-27 19:34:08 UTC (rev 7894)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-27 19:34:26 UTC (rev 7895)
@@ -31,7 +31,7 @@
#include <string.h>
#include <math.h>
#include "global_defs.h"
-
+#include "parallel_related.h"
#include "setProperties.h"
@@ -76,6 +76,9 @@
#define getDoubleVectorProperty(p, a, v, l, o) if (-1 == _getDoubleVectorProperty(p, a, v, l, o)) return NULL
+void myerror(struct All_variables *,char *);
+void report(struct All_variables *,char *);
+
/*==============================================================*/
@@ -309,6 +312,12 @@
getStringProperty(properties, "output_format", E->output.format, fp);
getStringProperty(properties, "output_optional", E->output.optional, fp);
+ getIntProperty(properties, "gzdir_vtkio", E->output.gzdir.vtk_io, fp);
+ getIntProperty(properties, "gzdir_rnr", E->output.gzdir.rnr, fp);
+ E->output.gzdir.vtk_base_init = 0;
+ /* should we save the basis vectors? (memory!) */
+ E->output.gzdir.vtk_base_save = 1;
+
getIntProperty(properties, "output_ll_max", E->output.llmax, fp);
getIntProperty(properties, "cb_block_size", E->output.cb_block_size, fp);
@@ -576,7 +585,7 @@
struct All_variables *E;
FILE *fp;
double tmp;
- void parallel_process_termination();
+ char message[100];
if (!PyArg_ParseTuple(args, "OOO:Tracer_set_properties",
&obj, &properties, &out))
@@ -589,6 +598,17 @@
getIntProperty(properties, "tracer", E->control.tracer, fp);
+ getIntProperty(properties, "tracer_enriched", E->control.tracer_enriched, fp);
+ if(E->control.tracer_enriched) {
+ if(!E->control.tracer)
+ myerror(E,"need to switch on tracers for tracer_enriched");
+
+ getFloatProperty(properties, "Q0_enriched", E->control.Q0ER, fp);
+ snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g (only one composition!)",
+ E->control.Q0,E->control.Q0ER);
+ report(E,message);
+ }
+
getIntProperty(properties, "tracer_ic_method",
E->trace.ic_method, fp);
@@ -611,11 +631,26 @@
getIntProperty(properties, "tracer_flavors", E->trace.nflavors, fp);
getIntProperty(properties, "ic_method_for_flavors", E->trace.ic_method_for_flavors, fp);
- if (E->trace.nflavors > 1 && E->trace.ic_method_for_flavors == 0) {
- E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
- *sizeof(double));
- getDoubleVectorProperty(properties, "z_interface", E->trace.z_interface, E->trace.nflavors-1, fp);
+ if (E->trace.nflavors > 1) {
+ switch(E->trace.ic_method_for_flavors){
+ case 0: /* layer */
+ E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
+ *sizeof(double));
+
+ getDoubleVectorProperty(properties, "z_interface", E->trace.z_interface, E->trace.nflavors-1, fp);
+ break;
+ case 1: /* from grid in top n materials */
+ /* file from which to read */
+ getStringProperty(properties, "ictracer_grd_file", E->trace.ggrd_file, fp);
+ /* which top layers to use */
+ getIntProperty(properties, "ictracer_grd_layers", E->trace.ggrd_layers, fp);
+ break;
+ default:
+ fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
+ }
}
getIntProperty(properties, "chemical_buoyancy",
@@ -672,7 +707,7 @@
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
- int num_mat;
+ int num_mat, i;
if (!PyArg_ParseTuple(args, "OOO:Visc_set_properties",
&obj, &properties, &out))
@@ -715,10 +750,38 @@
E->viscosity.Z, num_mat, fp);
getIntProperty(properties, "SDEPV", E->viscosity.SDEPV, fp);
- getFloatProperty(properties, "sdepv_misfit", E->viscosity.sdepv_misfit, fp);
getFloatVectorProperty(properties, "sdepv_expt",
E->viscosity.sdepv_expt, num_mat, fp);
+ getIntProperty(properties, "PDEPV", E->viscosity.PDEPV, fp);
+ if (E->viscosity.PDEPV) {
+ E->viscosity.pdepv_visited = 0;
+ getIntProperty(properties, "pdepv_eff", E->viscosity.pdepv_eff, fp);
+ getFloatVectorProperty(properties, "pdepv_a",
+ E->viscosity.pdepv_a, num_mat, fp);
+ getFloatVectorProperty(properties, "pdepv_b",
+ E->viscosity.pdepv_b, num_mat, fp);
+ getFloatVectorProperty(properties, "pdepv_y",
+ E->viscosity.pdepv_y, num_mat, fp);
+ getFloatProperty(properties, "pdepv_offset", E->viscosity.pdepv_offset, fp);
+ }
+ if(E->viscosity.PDEPV || E->viscosity.SDEPV)
+ getFloatProperty(properties, "sdepv_misfit", E->viscosity.sdepv_misfit, fp);
+
+ getIntProperty(properties, "CDEPV", E->viscosity.CDEPV, fp);
+ if(E->viscosity.CDEPV){ /* compositional viscosity */
+ if(!E->control.tracer)
+ myerror(E,"error: CDEPV requires tracers, but tracer is off");
+ if(E->trace.nflavors > 10)
+ myerror(E,"error: too many flavors for CDEPV");
+ /* read in flavor factors */
+ getFloatVectorProperty(properties, "cdepv_ff",
+ E->viscosity.cdepv_ff, E->trace.nflavors, fp);
+ /* and take the log because we're using a geometric avg */
+ for(i=0;i<E->trace.nflavors;i++)
+ E->viscosity.cdepv_ff[i] = log(E->viscosity.cdepv_ff[i]);
+ }
+
getIntProperty(properties, "low_visc_channel", E->viscosity.channel, fp);
getIntProperty(properties, "low_visc_wedge", E->viscosity.wedge, fp);
More information about the cig-commits
mailing list