https://github.com/geodynamics/citcoms
Revision 25a7493514f4d4b222ed1562b3566a057b52e548 authored by Eh Tan on 18 March 2009, 19:39:54 UTC, committed by Eh Tan on 18 March 2009, 19:39:54 UTC
In pyre version, the mesh size is always specified by nodex etc. In C version, the mesh size is specified by nodex if Solver=cgrad, and by mgunitx and levels if Solver=multigrid.
1 parent 04c4125
Tip revision: 25a7493514f4d4b222ed1562b3566a057b52e548 authored by Eh Tan on 18 March 2009, 19:39:54 UTC
Remove mgunitx etc from pyre input. This restores the behavior in v3.0 and earlier version.
Remove mgunitx etc from pyre input. This restores the behavior in v3.0 and earlier version.
Tip revision: 25a7493
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);
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, "mantle_temp", E->control.lith_age_mantle_temp, 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->mesh.mgunity = E->mesh.noy - 1;
E->mesh.mgunitz = E->mesh.noz - 1;
}
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, "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 */
Computing file changes ...