Skip to main content
  • Home
  • Development
  • Documentation
  • Donate
  • Operational login
  • Browse the archive

swh logo
SoftwareHeritage
Software
Heritage
Archive
Features
  • Search

  • Downloads

  • Save code now

  • Add forge now

  • Help

Revision 0461d2e117ce88704a56dd8bcbf6bf7787991b15 authored by Eh Tan on 08 November 2007, 23:28:46 UTC, committed by Eh Tan on 08 November 2007, 23:28:46 UTC
Merged revisions 8194-8245 via svnmerge from
svn+ssh://svn@geodynamics.org/cig/mc/3D/CitcomS/trunk

........
  r8194 | tan2 | 2007-10-30 14:49:58 -0700 (Tue, 30 Oct 2007) | 1 line
  
  Compute d(rho)/dr/rho from rho(r)
........
  r8195 | tan2 | 2007-10-30 14:50:52 -0700 (Tue, 30 Oct 2007) | 1 line
  
  Fixed a bug in dimensionalizing density. Provided the formula of geoid calculation in the comments. Rearranged the order of functions.
........
  r8196 | tan2 | 2007-10-30 14:53:50 -0700 (Tue, 30 Oct 2007) | 1 line
  
  A post-processing program to project geoid coefficents onto a regular (longitude, latitude) mesh
........
  r8197 | tan2 | 2007-10-30 14:54:14 -0700 (Tue, 30 Oct 2007) | 1 line
  
  Added the C program project_geoid to the makefile
........
  r8199 | tan2 | 2007-10-30 15:29:44 -0700 (Tue, 30 Oct 2007) | 1 line
  
  Minor modification
........
  r8201 | tan2 | 2007-11-01 16:33:30 -0700 (Thu, 01 Nov 2007) | 1 line
  
  Print dv/v=dp/p=1.0 for the 1st Uzawa iteraion
........
  r8202 | tan2 | 2007-11-01 16:33:50 -0700 (Thu, 01 Nov 2007) | 1 line
  
  Fixed an error in comment
........
  r8204 | tan2 | 2007-11-05 17:03:35 -0800 (Mon, 05 Nov 2007) | 1 line
  
  Scaled topo with variable gravity. Fixed an error in comment. Rearranged computation.
........
  r8205 | tan2 | 2007-11-05 17:03:55 -0800 (Mon, 05 Nov 2007) | 1 line
  
  Removed functions related sph. harm in lib/Regional_obsolete.c
........
  r8206 | tan2 | 2007-11-05 17:04:20 -0800 (Mon, 05 Nov 2007) | 1 line
  
  Shrank the size of sph. harm arrays
........
  r8207 | tan2 | 2007-11-05 17:04:43 -0800 (Mon, 05 Nov 2007) | 1 line
  
  Init'd some variables about vtk_io, which might be accessed with uninit'd values in output_finalize()
........
  r8212 | tan2 | 2007-11-06 15:17:54 -0800 (Tue, 06 Nov 2007) | 1 line
  
  Fixed a few memory errors
........
  r8213 | tan2 | 2007-11-06 15:18:12 -0800 (Tue, 06 Nov 2007) | 1 line
  
  Increase vlowstep to match the default value in pyre
........
  r8214 | tan2 | 2007-11-06 15:18:35 -0800 (Tue, 06 Nov 2007) | 1 line
  
  Removed unused multigrid parameters
........
  r8215 | tan2 | 2007-11-06 15:18:54 -0800 (Tue, 06 Nov 2007) | 1 line
  
  Added cgrad solver convergence parameters, increased buoyancy_ratio and lower the # of steps
........
  r8226 | tan2 | 2007-11-07 11:51:56 -0800 (Wed, 07 Nov 2007) | 1 line
  
  Print a warning when matrix eqn solver not converging
........
  r8227 | tan2 | 2007-11-07 11:52:17 -0800 (Wed, 07 Nov 2007) | 1 line
  
  Removed comp_el from default output, since it is not required for restart anymore.
........
  r8228 | tan2 | 2007-11-07 11:52:39 -0800 (Wed, 07 Nov 2007) | 1 line
  
  Decreased the # of processors. This is the only way I can reproduce single-cell convection as in the manual.
........
  r8235 | tan2 | 2007-11-08 11:18:26 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Dereased the timestep size to reduce artifacts in advection
........
  r8236 | tan2 | 2007-11-08 11:18:52 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Update NEWS
........
  r8237 | tan2 | 2007-11-08 11:19:12 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Update the version number
........
  r8241 | tan2 | 2007-11-08 13:17:14 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Updated file ChangeLog to r8240
........
  r8242 | tan2 | 2007-11-08 13:36:55 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Removed binary checkpoint files from makefile, as the file size is too big for distribution.
........
  r8243 | tan2 | 2007-11-08 13:38:09 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Updated file ChangeLog to r8242
........
  r8244 | tan2 | 2007-11-08 14:31:21 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Replaced a system call by std C library remove() and disabled another system call (backup input file). Partially fixed issue130. All remaining system calls are in lib/Output_gzdir.c.
........
  r8245 | tan2 | 2007-11-08 14:41:31 -0800 (Thu, 08 Nov 2007) | 1 line
  
  Updated file ChangeLog to r8244
........

1 parent a828fa9
  • Files
  • Changes
  • dd8a4ed
  • /
  • module
  • /
  • setProperties.c
Raw File Download
Permalinks

To reference or cite the objects present in the Software Heritage archive, permalinks based on SoftWare Hash IDentifiers (SWHIDs) must be used.
Select below a type of object currently browsed in order to display its associated SWHID and permalink.

  • revision
  • directory
  • content
revision badge
swh:1:rev:0461d2e117ce88704a56dd8bcbf6bf7787991b15
directory badge Iframe embedding
swh:1:dir:2e1cacc666e3080231c4c6fdc2695da4b1eaa0a9
content badge Iframe embedding
swh:1:cnt:1939df3ed70a99b6ab09412b635dc7cd9321ce3e
Citations

This interface enables to generate software citations, provided that the root directory of browsed objects contains a citation.cff or codemeta.json file.
Select below a type of object currently browsed in order to generate citations for them.

  • revision
  • directory
  • content
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
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 */

    /* 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;

    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);

    if (E->convection.tic_method == 0 || E->convection.tic_method == 3) {
	int num_perturb;

	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);
    }
    else if (E->convection.tic_method == 1) {
        getFloatProperty(properties, "half_space_age", E->convection.half_space_age, fp);
    }
    else if (E->convection.tic_method == 2) {
        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, "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, "file_vbcs", E->control.vbcs_file, fp);
    getStringProperty(properties, "vel_bound_file", E->control.velocity_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);

    E->mesh.mgunitx = (E->mesh.nox - 1) / E->parallel.nprocx /
	(int) pow(2.0, E->mesh.levels - 1);
    E->mesh.mgunity = (E->mesh.noy - 1) / E->parallel.nprocy /
	(int) pow(2.0, E->mesh.levels - 1);
    E->mesh.mgunitz = (E->mesh.noz - 1) / E->parallel.nprocz /
	(int) pow(2.0, E->mesh.levels - 1);

    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->parallel.nprocxy == 12) {
        full_set_3dsphere_defaults2(E);
    }
    else {
	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);

        regional_set_3dsphere_defaults2(E);
    }

    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);
    getFloatProperty(properties, "tole_compressibility", E->control.tole_comp, 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);

    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);
            getFloatProperty(properties, "relative_err_accuracy", E->control.relative_err_accuracy, 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 */
The diff you're trying to view is too large. Only the first 1000 changed files have been loaded.
Showing with 0 additions and 0 deletions (0 / 0 diffs computed)
swh spinner

Computing file changes ...

back to top

Software Heritage — Copyright (C) 2015–2025, The Software Heritage developers. License: GNU AGPLv3+.
The source code of Software Heritage itself is available on our development forge.
The source code files archived by Software Heritage are available under their own copyright and licenses.
Terms of use: Archive access, API— Contact— JavaScript license information— Web API