setProperties.c
/*
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//<LicenseText>
//
// CitcomS.py by Eh Tan, Eun-seo Choi, and Pururav Thoutireddy.
// Copyright (C) 2002-2005, California Institute of Technology.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
//</LicenseText>
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <Python.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "global_defs.h"
#include "parallel_related.h"
#include "setProperties.h"
/* See PEP 353. */
#if PY_VERSION_HEX < 0x02050000 && !defined(PY_SSIZE_T_MIN)
typedef int Py_ssize_t;
#define PY_SSIZE_T_MAX INT_MAX
#define PY_SSIZE_T_MIN INT_MIN
#endif
/*==============================================================*/
/* functions and macros which fetch properties from 'inventory' */
FILE *get_output_stream(PyObject *out, struct All_variables*E);
#define PUTS(s) if (fp) fprintf(fp, s)
int _getStringProperty(PyObject* properties, char* attribute,
char* value, size_t valueSize, FILE* fp);
#define getStringProperty(p, a, v, o) if (-1 == _getStringProperty(p, a, v, sizeof(v), o)) return NULL
int _getIntProperty(PyObject* properties, char* attribute, int *value, FILE* fp);
#define getIntProperty(p, a, v, o) if (-1 == _getIntProperty(p, a, &(v), o)) return NULL
int _getFloatProperty(PyObject* properties, char* attribute, float *value, FILE* fp);
#define getFloatProperty(p, a, v, o) if (-1 == _getFloatProperty(p, a, &(v), o)) return NULL
int _getDoubleProperty(PyObject* properties, char* attribute, double *value, FILE* fp);
#define getDoubleProperty(p, a, v, o) if (-1 == _getDoubleProperty(p, a, &(v), o)) return NULL
int _getIntVectorProperty(PyObject* properties, char* attribute,
int* vector, int len, FILE* fp);
#define getIntVectorProperty(p, a, v, l, o) if (-1 == _getIntVectorProperty(p, a, v, l, o)) return NULL
int _getFloatVectorProperty(PyObject* properties, char* attribute,
float* vector, int len, FILE* fp);
#define getFloatVectorProperty(p, a, v, l, o) if (-1 == _getFloatVectorProperty(p, a, v, l, o)) return NULL
int _getDoubleVectorProperty(PyObject* properties, char* attribute,
double* vector, int len, FILE* fp);
#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 *);
/*==============================================================*/
char pyCitcom_Advection_diffusion_set_properties__doc__[] = "";
char pyCitcom_Advection_diffusion_set_properties__name__[] = "Advection_diffusion_set_properties";
PyObject * pyCitcom_Advection_diffusion_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:Advection_diffusion_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.tsolver]\n"));
getIntProperty(properties, "ADV", E->advection.ADVECTION, fp);
getIntProperty(properties, "filter_temp", E->advection.filter_temperature, fp);
getIntProperty(properties, "monitor_max_T", E->advection.monitor_max_T, fp);
getFloatProperty(properties, "finetunedt", E->advection.fine_tune_dt, fp);
getFloatProperty(properties, "fixed_timestep", E->advection.fixed_timestep, fp);
getFloatProperty(properties, "adv_gamma", E->advection.gamma, fp);
getIntProperty(properties, "adv_sub_iterations", E->advection.temp_iterations, fp);
getFloatProperty(properties, "inputdiffusivity", E->control.inputdiff, fp);
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_BC_set_properties__doc__[] = "";
char pyCitcom_BC_set_properties__name__[] = "BC_set_properties";
PyObject * pyCitcom_BC_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:BC_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.bc]\n"));
getIntProperty(properties, "side_sbcs", E->control.side_sbcs, fp);
getIntProperty(properties, "pseudo_free_surf", E->control.pseudo_free_surf, fp);
getIntProperty(properties, "topvbc", E->mesh.topvbc, fp);
getFloatProperty(properties, "topvbxval", E->control.VBXtopval, fp);
getFloatProperty(properties, "topvbyval", E->control.VBYtopval, fp);
getIntProperty(properties, "botvbc", E->mesh.botvbc, fp);
getFloatProperty(properties, "botvbxval", E->control.VBXbotval, fp);
getFloatProperty(properties, "botvbyval", E->control.VBYbotval, fp);
getIntProperty(properties, "toptbc", E->mesh.toptbc, fp);
getFloatProperty(properties, "toptbcval", E->control.TBCtopval, fp);
getIntProperty(properties, "bottbc", E->mesh.bottbc, fp);
getFloatProperty(properties, "bottbcval", E->control.TBCbotval, fp);
getIntProperty(properties, "temperature_bound_adj", E->control.temperature_bound_adj, fp);
getFloatProperty(properties, "depth_bound_adj", E->control.depth_bound_adj, fp);
getFloatProperty(properties, "width_bound_adj", E->control.width_bound_adj, fp);
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Const_set_properties__doc__[] = "";
char pyCitcom_Const_set_properties__name__[] = "Const_set_properties";
PyObject * pyCitcom_Const_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
float radius;
if (!PyArg_ParseTuple(args, "OOO:Const_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.const]\n"));
getFloatProperty(properties, "radius", radius, fp);
getFloatProperty(properties, "density", E->data.density, fp);
getFloatProperty(properties, "thermdiff", E->data.therm_diff, fp);
getFloatProperty(properties, "gravacc", E->data.grav_acc, fp);
getFloatProperty(properties, "thermexp", E->data.therm_exp, fp);
getFloatProperty(properties, "refvisc", E->data.ref_viscosity, fp);
getFloatProperty(properties, "cp", E->data.Cp, fp);
getFloatProperty(properties, "density_above", E->data.density_above, fp);
getFloatProperty(properties, "density_below", E->data.density_below, fp);
E->data.therm_cond = E->data.therm_diff * E->data.density * E->data.Cp;
E->data.ref_temperature = E->control.Atemp * E->data.therm_diff
* E->data.ref_viscosity / (radius * radius * radius)
/ (E->data.density * E->data.grav_acc * E->data.therm_exp);
getFloatProperty(properties, "z_lith", E->viscosity.zlith, fp);
getFloatProperty(properties, "z_410", E->viscosity.z410, fp);
getFloatProperty(properties, "z_lmantle", E->viscosity.zlm, fp);
getFloatProperty(properties, "z_cmb", E->viscosity.zcmb, fp); /* this is used as the D" phase change depth */
E->viscosity.zbase_layer[0] = E->viscosity.zlith;
E->viscosity.zbase_layer[1] = E->viscosity.z410;
E->viscosity.zbase_layer[2] = E->viscosity.zlm;
E->viscosity.zbase_layer[3] = E->viscosity.zcmb;
/* convert meter to kilometer */
E->data.radius_km = radius / 1e3;
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_IC_set_properties__doc__[] = "";
char pyCitcom_IC_set_properties__name__[] = "IC_set_properties";
PyObject * pyCitcom_IC_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
int num_perturb;
if (!PyArg_ParseTuple(args, "OOO:IC_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.ic]\n"));
getIntProperty(properties, "restart", E->control.restart, fp);
getIntProperty(properties, "post_p", E->control.post_p, fp);
getIntProperty(properties, "solution_cycles_init", E->monitor.solution_cycles_init, fp);
getIntProperty(properties, "zero_elapsed_time", E->control.zero_elapsed_time, fp);
getIntProperty(properties, "tic_method", E->convection.tic_method, fp);
getIntProperty(properties, "num_perturbations", num_perturb, fp);
if(num_perturb > PERTURB_MAX_LAYERS) {
fprintf(stderr, "'num_perturb' greater than allowed value, set to %d\n", PERTURB_MAX_LAYERS);
num_perturb = PERTURB_MAX_LAYERS;
}
E->convection.number_of_perturbations = num_perturb;
getIntVectorProperty(properties, "perturbl", E->convection.perturb_ll,
num_perturb, fp);
getIntVectorProperty(properties, "perturbm", E->convection.perturb_mm,
num_perturb, fp);
getIntVectorProperty(properties, "perturblayer", E->convection.load_depth,
num_perturb, fp);
getFloatVectorProperty(properties, "perturbmag", E->convection.perturb_mag,
num_perturb, fp);
getFloatProperty(properties, "half_space_age", E->convection.half_space_age, fp);
getFloatProperty(properties, "mantle_temp", E->control.mantle_temp, fp);
getFloatVectorProperty(properties, "blob_center", E->convection.blob_center, 3, fp);
if( E->convection.blob_center[0] == -999.0 && E->convection.blob_center[1] == -999.0 && E->convection.blob_center[2] == -999.0 ) {
E->convection.blob_center[0] = 0.5*(E->control.theta_min+E->control.theta_max);
E->convection.blob_center[1] = 0.5*(E->control.fi_min+E->control.fi_max);
E->convection.blob_center[2] = 0.5*(E->sphere.ri+E->sphere.ro);
}
getFloatProperty(properties, "blob_radius", E->convection.blob_radius, fp);
getFloatProperty(properties, "blob_dT", E->convection.blob_dT, fp);
PUTS(("\n"));
if (PyErr_Occurred())
return NULL;
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Output_set_properties__doc__[] = "";
char pyCitcom_Output_set_properties__name__[] = "Output_set_properties";
PyObject * pyCitcom_Output_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:Output_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.output]\n"));
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, "self_gravitation", E->control.self_gravitation, fp);
getIntProperty(properties, "use_cbf_topo", E->control.use_cbf_topo, fp);
getIntProperty(properties, "cb_block_size", E->output.cb_block_size, fp);
getIntProperty(properties, "cb_buffer_size", E->output.cb_buffer_size, fp);
getIntProperty(properties, "sieve_buf_size", E->output.sieve_buf_size, fp);
getIntProperty(properties, "output_alignment", E->output.alignment, fp);
getIntProperty(properties, "output_alignment_threshold", E->output.alignment_threshold, fp);
getIntProperty(properties, "cache_mdc_nelmts", E->output.cache_mdc_nelmts, fp);
getIntProperty(properties, "cache_rdcc_nelmts", E->output.cache_rdcc_nelmts, fp);
getIntProperty(properties, "cache_rdcc_nbytes", E->output.cache_rdcc_nbytes, fp);
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Param_set_properties__doc__[] = "";
char pyCitcom_Param_set_properties__name__[] = "Param_set_properties";
PyObject * pyCitcom_Param_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:Param_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.param]\n"));
getIntProperty(properties, "reference_state", E->refstate.choice, fp);
if(E->refstate.choice == 0) {
getStringProperty(properties, "refstate_file", E->refstate.filename, fp);
}
getIntProperty(properties, "mineral_physics_model", E->control.mineral_physics_model, fp);
getIntProperty(properties, "file_vbcs", E->control.vbcs_file, fp);
getStringProperty(properties, "vel_bound_file", E->control.velocity_boundary_file, fp);
getIntProperty(properties, "file_tbcs", E->control.tbcs_file, fp);
getStringProperty(properties, "temp_bound_file", E->control.temperature_boundary_file, fp);
getIntProperty(properties, "mat_control", E->control.mat_control, fp);
getStringProperty(properties, "mat_file", E->control.mat_file, fp);
getIntProperty(properties, "lith_age", E->control.lith_age, fp);
getStringProperty(properties, "lith_age_file", E->control.lith_age_file, fp);
getIntProperty(properties, "lith_age_time", E->control.lith_age_time, fp);
getFloatProperty(properties, "lith_age_depth", E->control.lith_age_depth, fp);
getFloatProperty(properties, "start_age", E->control.start_age, fp);
getIntProperty(properties, "reset_startage", E->control.reset_startage, fp);
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Phase_set_properties__doc__[] = "";
char pyCitcom_Phase_set_properties__name__[] = "Phase_set_properties";
PyObject * pyCitcom_Phase_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
float width;
if (!PyArg_ParseTuple(args, "OOO:Phase_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.phase]\n"));
getFloatProperty(properties, "Ra_410", E->control.Ra_410, fp);
getFloatProperty(properties, "clapeyron410", E->control.clapeyron410, fp);
getFloatProperty(properties, "transT410", E->control.transT410, fp);
getFloatProperty(properties, "width410", width, fp);
if (width!=0.0)
E->control.inv_width410 = 1.0 / width;
getFloatProperty(properties, "Ra_670", E->control.Ra_670 , fp);
getFloatProperty(properties, "clapeyron670", E->control.clapeyron670, fp);
getFloatProperty(properties, "transT670", E->control.transT670, fp);
getFloatProperty(properties, "width670", width, fp);
if (width!=0.0)
E->control.inv_width670 = 1.0 / width;
getFloatProperty(properties, "Ra_cmb", E->control.Ra_cmb, fp);
getFloatProperty(properties, "clapeyroncmb", E->control.clapeyroncmb, fp);
getFloatProperty(properties, "transTcmb", E->control.transTcmb, fp);
getFloatProperty(properties, "widthcmb", width, fp);
if (width!=0.0)
E->control.inv_widthcmb = 1.0 / width;
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Solver_set_properties__doc__[] = "";
char pyCitcom_Solver_set_properties__name__[] = "Solver_set_properties";
PyObject * pyCitcom_Solver_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
float tmp;
if (!PyArg_ParseTuple(args, "OOO:Solver_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver]\n"));
getStringProperty(properties, "datadir", E->control.data_dir, fp);
getStringProperty(properties, "datafile", E->control.data_prefix, fp);
getStringProperty(properties, "datadir_old", E->control.data_dir_old, fp);
getStringProperty(properties, "datafile_old", E->control.data_prefix_old, fp);
getFloatProperty(properties, "rayleigh", E->control.Atemp, fp);
getFloatProperty(properties, "dissipation_number", E->control.disptn_number, fp);
getFloatProperty(properties, "gruneisen", tmp, fp);
/* special case: if tmp==0, set gruneisen as inf */
if(tmp != 0)
E->control.inv_gruneisen = 1/tmp;
else
E->control.inv_gruneisen = 0;
getFloatProperty(properties, "surfaceT", E->control.surface_temp, fp);
/*getFloatProperty(properties, "adiabaticT0", E->control.adiabaticT0, fp);*/
getFloatProperty(properties, "Q0", E->control.Q0, fp);
getIntProperty(properties, "stokes_flow_only", E->control.stokes, fp);
getIntProperty(properties, "verbose", E->control.verbose, fp);
getIntProperty(properties, "see_convergence", E->control.print_convergence, fp);
/* parameters not used in pyre version,
assigned value here to prevent uninitialized access */
E->advection.min_timesteps = 1;
E->advection.max_timesteps = 1;
E->advection.max_total_timesteps = 1;
E->control.checkpoint_frequency = 1;
E->control.record_every = 1;
E->control.record_all_until = 1;
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Sphere_set_properties__doc__[] = "";
char pyCitcom_Sphere_set_properties__name__[] = "Sphere_set_properties";
PyObject * pyCitcom_Sphere_set_properties(PyObject *self, PyObject *args)
{
void full_set_3dsphere_defaults2(struct All_variables *);
void regional_set_3dsphere_defaults2(struct All_variables *);
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:Sphere_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.mesher]\n"));
getIntProperty(properties, "nproc_surf", E->parallel.nprocxy, fp);
getIntProperty(properties, "nprocx", E->parallel.nprocx, fp);
getIntProperty(properties, "nprocy", E->parallel.nprocy, fp);
getIntProperty(properties, "nprocz", E->parallel.nprocz, fp);
if (E->parallel.nprocxy == 12)
if (E->parallel.nprocx != E->parallel.nprocy) {
char errmsg[] = "!!!! nprocx must equal to nprocy";
PyErr_SetString(PyExc_SyntaxError, errmsg);
return NULL;
}
getIntProperty(properties, "coor", E->control.coor, fp);
getFloatVectorProperty(properties, "coor_refine", E->control.coor_refine, 4, fp);
getStringProperty(properties, "coor_file", E->control.coor_file, fp);
getIntProperty(properties, "nodex", E->mesh.nox, fp);
getIntProperty(properties, "nodey", E->mesh.noy, fp);
getIntProperty(properties, "nodez", E->mesh.noz, fp);
getIntProperty(properties, "levels", E->mesh.levels, fp);
if (E->control.CONJ_GRAD) {
if (E->mesh.levels != 1) {
fprintf(stderr, "!!! cgrad solver must have levels=1\n");
abort();
}
E->mesh.mgunitx = (E->mesh.nox - 1) / E->parallel.nprocx;
E->mesh.mgunity = (E->mesh.noy - 1) / E->parallel.nprocy;
E->mesh.mgunitz = (E->mesh.noz - 1) / E->parallel.nprocz;
}
else {
double levmax;
int nox, noy, noz;
if (E->mesh.levels <= 1) {
fprintf(stderr, "!!! multigrid solver must have levels>1\n");
abort();
}
levmax = E->mesh.levels - 1;
E->mesh.mgunitx = (E->mesh.nox - 1) / E->parallel.nprocx / pow(2.0, levmax);
E->mesh.mgunity = (E->mesh.noy - 1) / E->parallel.nprocy / pow(2.0, levmax);
E->mesh.mgunitz = (E->mesh.noz - 1) / E->parallel.nprocz / pow(2.0, levmax);
nox = E->mesh.mgunitx * (int) pow(2.0,levmax) * E->parallel.nprocx + 1;
noy = E->mesh.mgunity * (int) pow(2.0,levmax) * E->parallel.nprocy + 1;
noz = E->mesh.mgunitz * (int) pow(2.0,levmax) * E->parallel.nprocz + 1;
if ((nox != E->mesh.nox) || (noy != E->mesh.noy) || (noz != E->mesh.noz)) {
fprintf(stderr, "!!! inconsistent mesh size and levels.\n");
abort();
}
}
if (E->parallel.nprocxy == 12) {
if (E->mesh.nox != E->mesh.noy) {
char errmsg[] = "!!!! nodex must equal to nodey";
PyErr_SetString(PyExc_SyntaxError, errmsg);
return NULL;
}
}
getDoubleProperty(properties, "radius_outer", E->sphere.ro, fp);
getDoubleProperty(properties, "radius_inner", E->sphere.ri, fp);
if (E->sphere.caps == 1) {
getDoubleProperty(properties, "theta_min", E->control.theta_min, fp);
getDoubleProperty(properties, "theta_max", E->control.theta_max, fp);
getDoubleProperty(properties, "fi_min", E->control.fi_min, fp);
getDoubleProperty(properties, "fi_max", E->control.fi_max, fp);
}
E->mesh.layer[1] = 1;
E->mesh.layer[2] = 1;
E->mesh.layer[3] = 1;
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Tracer_set_properties__doc__[] = "";
char pyCitcom_Tracer_set_properties__name__[] = "Tracer_set_properties";
PyObject * pyCitcom_Tracer_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
double tmp;
char message[100];
if (!PyArg_ParseTuple(args, "OOO:Tracer_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.tracer]\n"));
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);
if (E->trace.ic_method==0) {
getIntProperty(properties, "tracers_per_element",
E->trace.itperel, fp);
}
else if (E->trace.ic_method==1) {
getStringProperty(properties, "tracer_file",
E->trace.tracer_file, fp);
}
else if (E->trace.ic_method==2) {
}
else {
fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
fflush(stderr);
parallel_process_termination();
}
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) {
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, "itracer_warnings", E->trace.itracer_warnings, fp);
getIntProperty(properties, "chemical_buoyancy",
E->composition.ichemical_buoyancy, fp);
if (E->control.tracer && E->composition.ichemical_buoyancy==1) {
getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
if (E->composition.ibuoy_type==0)
E->composition.ncomp = E->trace.nflavors;
else if (E->composition.ibuoy_type==1)
E->composition.ncomp = E->trace.nflavors - 1;
E->composition.buoyancy_ratio = (double*) malloc(E->composition.ncomp
*sizeof(double));
getDoubleVectorProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, E->composition.ncomp, fp);
}
if(E->parallel.nprocxy == 12) {
getDoubleProperty(properties, "regular_grid_deltheta", tmp, fp);
E->trace.deltheta[0] = tmp;
getDoubleProperty(properties, "regular_grid_delphi", tmp, fp);
E->trace.delphi[0] = tmp;
E->trace.ianalytical_tracer_test = 0;
//getIntProperty(properties, "analytical_tracer_test", E->trace.ianalytical_tracer_test, fp);
E->composition.icompositional_rheology = 0;
/*
getIntProperty(properties, "compositional_rheology", E->composition.icompositional_rheology, fp);
if (E->composition.icompositional_rheology==1) {
getDoubleProperty(properties, "compositional_prefactor", E->composition.compositional_rheology_prefactor, fp);
}
*/
}
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Visc_set_properties__doc__[] = "";
char pyCitcom_Visc_set_properties__name__[] = "Visc_set_properties";
PyObject * pyCitcom_Visc_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
int num_mat, i;
if (!PyArg_ParseTuple(args, "OOO:Visc_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.visc]\n"));
getStringProperty(properties, "Viscosity", E->viscosity.STRUCTURE, fp);
if (strcmp(E->viscosity.STRUCTURE,"system") == 0)
E->viscosity.FROM_SYSTEM = 1;
else
E->viscosity.FROM_SYSTEM = 0;
getIntProperty(properties, "visc_smooth_method", E->viscosity.smooth_cycles, fp);
getIntProperty(properties, "VISC_UPDATE", E->viscosity.update_allowed, fp);
#define MAX_MAT 40
getIntProperty(properties, "num_mat", num_mat, fp);
if(num_mat > MAX_MAT) {
/* max. allowed material types = 40 */
fprintf(stderr, "'num_mat' greater than allowed value, set to %d\n", MAX_MAT);
num_mat = MAX_MAT;
}
E->viscosity.num_mat = num_mat;
getFloatVectorProperty(properties, "visc0",
E->viscosity.N0, num_mat, fp);
getIntProperty(properties, "TDEPV", E->viscosity.TDEPV, fp);
getIntProperty(properties, "rheol", E->viscosity.RHEOL, fp);
getFloatVectorProperty(properties, "viscE",
E->viscosity.E, num_mat, fp);
getFloatVectorProperty(properties, "viscT",
E->viscosity.T, num_mat, fp);
getFloatVectorProperty(properties, "viscZ",
E->viscosity.Z, num_mat, fp);
getIntProperty(properties, "SDEPV", E->viscosity.SDEPV, 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);
getFloatProperty(properties, "lv_min_radius", E->viscosity.lv_min_radius, fp);
getFloatProperty(properties, "lv_max_radius", E->viscosity.lv_max_radius, fp);
getFloatProperty(properties, "lv_channel_thickness", E->viscosity.lv_channel_thickness, fp);
getFloatProperty(properties, "lv_reduction", E->viscosity.lv_reduction, fp);
getIntProperty(properties, "VMIN", E->viscosity.MIN, fp);
getFloatProperty(properties, "visc_min", E->viscosity.min_value, fp);
getIntProperty(properties, "VMAX", E->viscosity.MAX, fp);
getFloatProperty(properties, "visc_max", E->viscosity.max_value, fp);
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
char pyCitcom_Incompressible_set_properties__doc__[] = "";
char pyCitcom_Incompressible_set_properties__name__[] = "Incompressible_set_properties";
PyObject * pyCitcom_Incompressible_set_properties(PyObject *self, PyObject *args)
{
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
if (!PyArg_ParseTuple(args, "OOO:Incompressible_set_properties",
&obj, &properties, &out))
return NULL;
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
fp = get_output_stream(out, E);
PUTS(("[CitcomS.solver.vsolver]\n"));
getStringProperty(properties, "Solver", E->control.SOLVER_TYPE, fp);
getIntProperty(properties, "node_assemble", E->control.NASSEMBLE, fp);
getIntProperty(properties, "precond", E->control.precondition, fp);
getDoubleProperty(properties, "accuracy", E->control.accuracy, fp);
getIntProperty(properties, "mg_cycle", E->control.mg_cycle, fp);
getIntProperty(properties, "down_heavy", E->control.down_heavy, fp);
getIntProperty(properties, "up_heavy", E->control.up_heavy, fp);
getIntProperty(properties, "vlowstep", E->control.v_steps_low, fp);
getIntProperty(properties, "vhighstep", E->control.v_steps_high, fp);
getIntProperty(properties, "max_mg_cycles", E->control.max_mg_cycles, fp);
getIntProperty(properties, "piterations", E->control.p_iterations, fp);
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);
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" */
getStringProperty(properties, "uzawa", E->control.uzawa, fp);
if(strcmp(E->control.uzawa, "cg") == 0) {
/* more convergence parameters for "cg" */
getIntProperty(properties, "compress_iter_maxstep", E->control.compress_iter_maxstep, fp);
}
}
PUTS(("\n"));
Py_INCREF(Py_None);
return Py_None;
}
/*==========================================================*/
/* helper functions */
FILE *get_output_stream(PyObject *out, struct All_variables*E)
{
if (PyFile_Check(out)) {
return PyFile_AsFile(out);
}
return NULL;
}
int _getStringProperty(PyObject* properties, char* attribute,
char* value, size_t valueSize, FILE* fp)
{
PyObject *prop;
char *buffer;
Py_ssize_t length;
if (!(prop = PyObject_GetAttrString(properties, attribute)))
return -1;
if (-1 == PyString_AsStringAndSize(prop, &buffer, &length))
return -1;
if (length >= (Py_ssize_t)valueSize) {
PyErr_Format(PyExc_ValueError,
"value of '%s' cannot exceed %zu characters in length",
attribute, valueSize);
return -1;
}
strcpy(value, buffer);
if (fp)
fprintf(fp, "%s=%s\n", attribute, value);
return 0;
}
#define getTYPEProperty _getIntProperty
#define getTYPEVectorProperty _getIntVectorProperty
#define PyTYPE_Check PyInt_Check
#define CTYPE int
#define PyTYPE_AsCTYPE PyInt_AsLong
#define MESSAGE "an integer is required"
#define FORMAT "%d"
#include "getProperty.h"
#undef getTYPEProperty
#undef getTYPEVectorProperty
#undef PyTYPE_Check
#undef CTYPE
#undef PyTYPE_AsCTYPE
#undef MESSAGE
#undef FORMAT
#define getTYPEProperty _getFloatProperty
#define getTYPEVectorProperty _getFloatVectorProperty
#define PyTYPE_Check PyFloat_Check
#define CTYPE float
#define PyTYPE_AsCTYPE PyFloat_AsDouble
#define MESSAGE "a float is required"
#define FORMAT "%g"
#include "getProperty.h"
#undef getTYPEProperty
#undef getTYPEVectorProperty
#undef PyTYPE_Check
#undef CTYPE
#undef PyTYPE_AsCTYPE
#undef MESSAGE
#undef FORMAT
#define getTYPEProperty _getDoubleProperty
#define getTYPEVectorProperty _getDoubleVectorProperty
#define PyTYPE_Check PyFloat_Check
#define CTYPE double
#define PyTYPE_AsCTYPE PyFloat_AsDouble
#define MESSAGE "a float is required"
#define FORMAT "%g"
#include "getProperty.h"
/* $Id$ */
/* End of file */
