[cig-commits] r13268 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Thu Nov 6 19:29:17 PST 2008
Author: becker
Date: 2008-11-06 19:29:17 -0800 (Thu, 06 Nov 2008)
New Revision: 13268
Added:
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Makefile.gzdir
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
Modified:
mc/3D/CitcomCU/trunk/src/Advection_diffusion.c
mc/3D/CitcomCU/trunk/src/Construct_arrays.c
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Convection_variables.h
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Makefile
mc/3D/CitcomCU/trunk/src/Nodal_mesh.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
mc/3D/CitcomCU/trunk/src/Parsing.c
mc/3D/CitcomCU/trunk/src/Process_velocity.c
mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
Log:
Replaced Parsing.c with the CIG CitcomS version to ensure that
defaults settings are actually used. Before, the respective function
always returned NULL to avoid problems with strtok on Linux.
Modified layers() such that the full range of 1...4 layers gets
assigned to E->mat[], otherwise not all four layers would be used in
Viscosity_structures.c
Added two new temp dep rheology options, rheol==3 is very similar to
0, and case has a non-dim depth dependence in it
Added ggrd netcdf support for temperature and composition using the HC
package.
Added gzdir output, in which case (-DUSE_GZDIR), get_system_viscosity
will compute the visocosities at the nodes.
Modified: mc/3D/CitcomCU/trunk/src/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Advection_diffusion.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Advection_diffusion.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -64,6 +64,7 @@
void advection_diffusion_parameters(struct All_variables *E)
{
/* Set intial values, defaults & read parameters */
+ int m = E->parallel.me;
E->advection.fixed_timestep = 0.0;
E->advection.temp_iterations = 2; /* petrov-galerkin iterations: minimum value. */
@@ -72,32 +73,32 @@
E->advection.last_sub_iterations = 1;
E->advection.gamma = 0.5;
E->advection.ADVECTION = 1;
+
+ input_boolean("ADV", &(E->advection.ADVECTION), "on",m);
- input_boolean("ADV", &(E->advection.ADVECTION), "on");
+ input_int("minstep", &(E->advection.min_timesteps), "1",m);
+ input_int("maxstep", &(E->advection.max_timesteps), "1000",m);
+ input_int("maxtotstep", &(E->advection.max_total_timesteps), "1000000",m);
+ input_float("finetunedt", &(E->advection.fine_tune_dt), "0.9",m);
+ input_float("fixed_timestep", &(E->advection.fixed_timestep), "0.0",m);
+ input_int("adv_sub_iterations", &(E->advection.temp_iterations), "2,2,nomax",m);
+ input_float("maxadvtime", &(E->advection.max_dimensionless_time), "10.0",m);
- input_int("minstep", &(E->advection.min_timesteps), "1");
- input_int("maxstep", &(E->advection.max_timesteps), "1000");
- input_int("maxtotstep", &(E->advection.max_total_timesteps), "1000000");
- input_float("finetunedt", &(E->advection.fine_tune_dt), "0.9");
- input_float("fixed_timestep", &(E->advection.fixed_timestep), "0.0");
- input_int("adv_sub_iterations", &(E->advection.temp_iterations), "2,2,nomax");
- input_float("maxadvtime", &(E->advection.max_dimensionless_time), "10.0");
+ input_float("sub_tolerance", &(E->advection.vel_substep_aggression), "0.005",m);
+ input_int("maxsub", &(E->advection.max_substeps), "25",m);
- input_float("sub_tolerance", &(E->advection.vel_substep_aggression), "0.005");
- input_int("maxsub", &(E->advection.max_substeps), "25");
-
- input_float("liddefvel", &(E->advection.lid_defining_velocity), "0.01");
- input_float("sublayerfrac", &(E->advection.sub_layer_sample_level), "0.5");
+ input_float("liddefvel", &(E->advection.lid_defining_velocity), "0.01",m);
+ input_float("sublayerfrac", &(E->advection.sub_layer_sample_level), "0.5",m);
E->control.adi_heating = 0;
E->control.visc_heating = 0;
- input_int("adi_heating", &(E->control.adi_heating), "1");
- input_int("visc_heating", &(E->control.visc_heating), "1");
+ input_int("adi_heating", &(E->control.adi_heating), "1",m);
+ input_int("visc_heating", &(E->control.visc_heating), "1",m);
E->viscosity.zcomp = E->control.Q0ER = 0;
- input_int("markers_per_ele", &(E->advection.markers_per_ele), "0");
- input_float("comp_depth", &(E->viscosity.zcomp), "0.0");
- input_float("Q0_enriched", &(E->control.Q0ER), "0.0");
+ input_int("markers_per_ele", &(E->advection.markers_per_ele), "0",m);
+ input_float("comp_depth", &(E->viscosity.zcomp), "0.0",m);
+ input_float("Q0_enriched", &(E->control.Q0ER), "0.0",m);
E->viscosity.zcomp = 1.0 - E->viscosity.zcomp;
Modified: mc/3D/CitcomCU/trunk/src/Construct_arrays.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -47,6 +47,7 @@
NOTE: this is not really general enough for new elements:
it should be done through a pre-calculated lookup table.
======================================================== */
+int layers(struct All_variables *, float);
void construct_ien(struct All_variables *E)
{
@@ -784,12 +785,11 @@
x3 = x3 / ends;
llayer = layers(E, x3);
- if(llayer)
- { /* for layers:1-lith,2-upper and 3-lower mantle */
- E->mat[el] = llayer;
+ if(llayer){
+ /* for layers:1-lith,2-upper and 3-lower mantle */
+ E->mat[el] = llayer;
}
-
/* if (E->mat[el]==1 || E->mat[el]==2) {
* for (a=1;a<=ends;a++) {
* nodea = E->ien[el].node[a];
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -45,10 +45,13 @@
#include "global_defs.h"
#include <stdlib.h> /* for "system" command */
#include <string.h>
+void convection_initial_markers(struct All_variables *,int );
+
void set_convection_defaults(struct All_variables *E)
{
- input_int("composition", &(E->control.composition), "0");
+
+ input_int("composition", &(E->control.composition), "0", E->parallel.me);
if(E->control.composition)
E->next_buoyancy_field = PG_timestep_particle;
@@ -79,52 +82,56 @@
char tmp_string[100], tmp1_string[100];
/* parameters */
+ int m;
- input_float("rayleigh", &(E->control.Atemp), "essential");
+ m = E->parallel.me;
- input_float("rayleigh_comp", &(E->control.Acomp), "essential");
+ input_float("rayleigh", &(E->control.Atemp), "essential", m);
- input_boolean("halfspace", &(E->convection.half_space_cooling), "off");
- input_float("halfspage", &(E->convection.half_space_age), "nodefault");
+ input_float("rayleigh_comp", &(E->control.Acomp), "essential", m);
+ input_boolean("halfspace", &(E->convection.half_space_cooling), "off", m);
+ input_float("halfspage", &(E->convection.half_space_age), "nodefault", m);
+
/*
- input_int("temperature_blobs", &(E->convection.temp_blobs), "0");
- input_float_vector("temperature_blobx", E->convection.temp_blobs, E->convection.temp_blob_x);
- input_float_vector("temperature_bloby", E->convection.temp_blobs, E->convection.temp_blob_y);
- input_float_vector("temperature_blobz", E->convection.temp_blobs, E->convection.temp_blob_z);
- input_float_vector("temperature_blobsize", E->convection.temp_blobs, E->convection.temp_blob_radius);
- input_float_vector("temperature_blobDT", E->convection.temp_blobs, E->convection.temp_blob_T);
- input_float_vector("temperature_blobbg", E->convection.temp_blobs, E->convection.temp_blob_bg);
- input_int_vector("temperature_blobsticky", E->convection.temp_blobs, E->convection.temp_blob_sticky);
+ input_int("temperature_blobs", &(E->convection.temp_blobs), "0", m);
+ input_float_vector("temperature_blobx", E->convection.temp_blobs, E->convection.temp_blob_x, m);
+ input_float_vector("temperature_bloby", E->convection.temp_blobs, E->convection.temp_blob_y, m);
+ input_float_vector("temperature_blobz", E->convection.temp_blobs, E->convection.temp_blob_z, m);
+ input_float_vector("temperature_blobsize", E->convection.temp_blobs, E->convection.temp_blob_radius, m);
+ input_float_vector("temperature_blobDT", E->convection.temp_blobs, E->convection.temp_blob_T, m);
+ input_float_vector("temperature_blobbg", E->convection.temp_blobs, E->convection.temp_blob_bg, m);
+ input_int_vector("temperature_blobsticky", E->convection.temp_blobs, E->convection.temp_blob_sticky, m);
- input_int("temperature_zones", &(E->convection.temp_zones), "0");
- input_float_vector("temperature_zonex1", E->convection.temp_zones, E->convection.temp_zonex1);
- input_float_vector("temperature_zonex2", E->convection.temp_zones, E->convection.temp_zonex2);
- input_float_vector("temperature_zonez1", E->convection.temp_zones, E->convection.temp_zonez1);
- input_float_vector("temperature_zonez2", E->convection.temp_zones, E->convection.temp_zonez2);
- input_float_vector("temperature_zoney1", E->convection.temp_zones, E->convection.temp_zoney1);
- input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2);
- input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2);
- input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2);
- input_float_vector("temperature_zonehw", E->convection.temp_zones, E->convection.temp_zonehw);
- input_float_vector("temperature_zonemag", E->convection.temp_zones, E->convection.temp_zonemag);
- input_int_vector("temperature_zonesticky", E->convection.temp_zones, E->convection.temp_zone_sticky);
+ input_int("temperature_zones", &(E->convection.temp_zones), "0", m);
+ input_float_vector("temperature_zonex1", E->convection.temp_zones, E->convection.temp_zonex1, m);
+ input_float_vector("temperature_zonex2", E->convection.temp_zones, E->convection.temp_zonex2, m);
+ input_float_vector("temperature_zonez1", E->convection.temp_zones, E->convection.temp_zonez1, m);
+ input_float_vector("temperature_zonez2", E->convection.temp_zones, E->convection.temp_zonez2, m);
+ input_float_vector("temperature_zoney1", E->convection.temp_zones, E->convection.temp_zoney1, m);
+ input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2, m);
+ input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2, m);
+ input_float_vector("temperature_zoney2", E->convection.temp_zones, E->convection.temp_zoney2, m);
+ input_float_vector("temperature_zonehw", E->convection.temp_zones, E->convection.temp_zonehw, m);
+ input_float_vector("temperature_zonemag", E->convection.temp_zones, E->convection.temp_zonemag, m);
+ input_int_vector("temperature_zonesticky", E->convection.temp_zones, E->convection.temp_zone_sticky, m);
*/
- input_int("num_perturbations", &(E->convection.number_of_perturbations), "0,0,32");
- input_float_vector("perturbmag", E->convection.number_of_perturbations, E->convection.perturb_mag);
- input_float_vector("perturbk", E->convection.number_of_perturbations, E->convection.perturb_k);
- input_float_vector("perturbm", E->convection.number_of_perturbations, E->convection.perturb_mm);
- input_float_vector("perturbl", E->convection.number_of_perturbations, E->convection.perturb_ll);
+ input_int("num_perturbations", &(E->convection.number_of_perturbations), "0,0,32", m);
+ input_boolean("random_t_init",&E->convection.random_t_init,"off", m);
+ input_float_vector("perturbmag", E->convection.number_of_perturbations, E->convection.perturb_mag, m);
+ input_float_vector("perturbk", E->convection.number_of_perturbations, E->convection.perturb_k, m);
+ input_float_vector("perturbm", E->convection.number_of_perturbations, E->convection.perturb_mm, m);
+ input_float_vector("perturbl", E->convection.number_of_perturbations, E->convection.perturb_ll, m);
- input_string("prevT", E->convection.old_T_file, "initialize");
+ input_string("prevT", E->convection.old_T_file, "initialize", m);
if(E->control.restart)
{
- input_int("restart_timesteps", &(E->monitor.solution_cycles), "0");
+ input_int("restart_timesteps", &(E->monitor.solution_cycles), "0", m);
- input_string("oldfile", tmp1_string, "initialize");
- input_string("use_scratch", tmp_string, "local");
+ input_string("oldfile", tmp1_string, "initialize", m);
+ input_string("use_scratch", tmp_string, "local", m);
if(strcmp(tmp_string, "local") == 0)
strcpy(E->convection.old_T_file, tmp1_string);
else
@@ -164,7 +171,7 @@
E->advection.markers = E->advection.markers_per_ele * E->mesh.nel;
E->advection.markers = E->advection.markers * E->lmesh.volume / E->mesh.volume;
E->advection.markers_uplimit = E->advection.markers * 2;
- fprintf(stderr, "aaaa %d %g %g\n", E->advection.markers, E->lmesh.volume, E->mesh.volume);
+ //fprintf(stderr, "aaaa %d %g %g\n", E->advection.markers, E->lmesh.volume, E->mesh.volume);
for(i = 1; i <= E->mesh.nsd; i++)
{
E->VO[i] = (float *)malloc((E->advection.markers_uplimit + 1) * sizeof(float));
@@ -178,8 +185,11 @@
}
report(E, "convection, initial temperature");
+#ifdef USE_GGRD
+ convection_initial_temperature_ggrd(E);
+#else
convection_initial_temperature(E);
-
+#endif
return;
}
@@ -197,6 +207,11 @@
/* ===============================
Initialization of fields .....
+
+ there's a different routine for GGRD handling in Ggrd_handling
+ which has differnet logic for marker init etc.
+
+
=============================== */
void convection_initial_temperature(struct All_variables *E)
@@ -208,7 +223,7 @@
//int noz2, nfz, in1, in2, in3, instance, nox, noy, noz;
int noz2, nox, noy, noz;
- //char input_s[200], output_file[255];
+ //char input_s[200], output_file[255, m);
//float weight, para1, plate_velocity, delta_temp, age;
//const int dims = E->mesh.nsd;
@@ -288,7 +303,7 @@
if(E->control.composition)
- convection_initial_markers(E);
+ convection_initial_markers(E,0);
} // end for restart==0
@@ -403,7 +418,7 @@
if(E->control.composition && (E->control.restart == 1 || E->control.restart == 2))
{
- convection_initial_markers(E);
+ convection_initial_markers(E,0);
}
else if(E->control.composition && (E->control.restart == 3))
@@ -565,15 +580,19 @@
return;
}
-void convection_initial_markers(struct All_variables *E)
+void convection_initial_markers(struct All_variables *E,int use_element_nodes_for_init_c)
{
//int el, i, j, k, p, node, ii, jj;
- int el, node;
+ int el, node,j;
//double x, y, z, r, t, f, dX[4], dx, dr;
double x, y, z, r, t, f, dX[4];
//char input_s[100], output_file[255];
//FILE *fp;
+ float temp;
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+
if(E->control.CART3D)
{
node = 0;
@@ -592,10 +611,20 @@
el = get_element(E, E->XMC[1][node], E->XMC[2][node], E->XMC[3][node], dX);
E->CElement[node] = el;
- if(E->XMC[3][node] > E->viscosity.zcomp)
- E->C12[node] = 0;
- else
- E->C12[node] = 1;
+ if(use_element_nodes_for_init_c){
+ for(temp=0.0,j = 1; j <= ends; j++)
+ temp += E->C[E->ien[el].node[j]];
+ temp /= ends;
+ if(temp >0.5)
+ E->C12[node] = 1;
+ else
+ E->C12[node] = 0;
+ }else{ /* use depth */
+ if(E->XMC[3][node] > E->viscosity.zcomp)
+ E->C12[node] = 0;
+ else
+ E->C12[node] = 1;
+ }
}
} while(node < E->advection.markers);
}
@@ -619,10 +648,20 @@
E->XMC[3][node] = r;
el = get_element(E, E->XMC[1][node], E->XMC[2][node], E->XMC[3][node], dX);
E->CElement[node] = el;
- if(r > E->viscosity.zcomp)
- E->C12[node] = 0;
- else
- E->C12[node] = 1;
+ if(use_element_nodes_for_init_c){
+ for(temp=0.0,j = 1; j <= ends; j++)
+ temp += E->C[E->ien[el].node[j]];
+ temp /= ends;
+ if(temp >0.5)
+ E->C12[node] = 1;
+ else
+ E->C12[node] = 0;
+ }else{ /* use depth */
+ if(r > E->viscosity.zcomp)
+ E->C12[node] = 0;
+ else
+ E->C12[node] = 1;
+ }
}
} while(node < E->advection.markers);
}
@@ -640,24 +679,21 @@
//char output_file[255];
int l, noz;
- double temp, temp1, time_scale, velo_scale;
+ double temp, temp1;
noz = E->lmesh.noz;
E->control.plate_vel = E->control.plate_vel * 0.01 / 3.1536e7;
- /* now with a unit of m/sec */
- velo_scale = E->data.therm_diff / (E->monitor.length_scale);
/* in m/sec */
- E->control.plate_vel = E->control.plate_vel / velo_scale;
+ E->control.plate_vel = E->control.plate_vel / E->monitor.velo_scale;
/* dimensionless */
- time_scale = E->monitor.length_scale * E->monitor.length_scale / E->data.therm_diff;
/* in sec */
E->control.plate_age = E->control.plate_age * 1.0e6 * 3.1536e7;
/* in sec */
- E->control.plate_age = E->control.plate_age / time_scale;
+ E->control.plate_age = E->control.plate_age / E->monitor.time_scale;
/* dimensionless */
Modified: mc/3D/CitcomCU/trunk/src/Convection_variables.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection_variables.h 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Convection_variables.h 2008-11-07 03:29:17 UTC (rev 13268)
@@ -58,6 +58,8 @@
int temp_zone_sticky[40];
int temp_zones;
+ int random_t_init;
+
float half_space_age;
int half_space_cooling;
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -50,16 +50,33 @@
//int count, i, j, k;
int count, i;
+ static double alpha,alpha1;
static double *oldU;
- static int visits = 0;
+ static int damp=0,visits = 0;
//const int nno = E->lmesh.nno;
const int neq = E->lmesh.neq;
//const int dims = E->mesh.nsd;
//const int vpts = vpoints[E->mesh.nsd];
+ static int powerlaw;
+
if(visits == 0)
{
+
+ powerlaw = (E->viscosity.SDEPV || E->viscosity.BDEPV)?(1):(0);
+ if(powerlaw){
+ /* damping factors */
+ alpha = E->viscosity.sdepv_iter_damp;
+ alpha1 = 1 - alpha;
+ if(fabs(alpha-1) > 1e-7){
+ if(E->parallel.me == 0)
+ fprintf(stderr,"damping stress dependent iteration velocities by %g\n",alpha);
+ damp = 1;
+ } else{
+ damp = 0;
+ }
+ } /* end powerlaw */
oldU = (double *)malloc((neq + 2) * sizeof(double));
for(i = 1; i <= neq; i++)
oldU[i] = 0.0;
@@ -102,8 +119,14 @@
solve_constrained_flow_iterative(E);
- if(E->viscosity.SDEPV)
+ if(powerlaw)
{
+
+ if(damp){
+ /* add some of the old solution */
+ for(i = 1; i <= neq; i++)
+ E->U[i] = alpha * E->U[i] + alpha1 * oldU[i];
+ }
for(i = 1; i <= neq; i++)
{
delta_U[i] = E->U[i] - oldU[i];
@@ -115,17 +138,17 @@
if(Udot_mag != 0.0)
dUdot_mag /= Udot_mag;
- if(E->control.sdepv_print_convergence < E->monitor.solution_cycles && E->parallel.me == 0)
+ if(E->control.sdepv_print_convergence && (E->parallel.me == 0))
{
fprintf(stderr, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, count);
fprintf(E->fp, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, count);
fflush(E->fp);
}
count++;
- } /* end for SDEPV */
+ } /* end for SDEPV / BDEPV */
+
+ } while((count < 50) && (dUdot_mag > E->viscosity.sdepv_misfit) && powerlaw);
- } while((count < 50) && (dUdot_mag > E->viscosity.sdepv_misfit) && E->viscosity.SDEPV);
-
free((void *)delta_U);
return;
Added: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c (rev 0)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -0,0 +1,484 @@
+/*
+ * CitcomCU is a Finite Element Code that solves for thermochemical
+ * convection within a three dimensional domain appropriate for convection
+ * within the Earth's mantle. Cartesian and regional-spherical geometries
+ * are implemented. See the file README contained with this distribution
+ * for further details.
+ *
+ * Copyright (C) 1994-2005 California Institute of Technology
+ * Copyright (C) 2000-2005 The University of Colorado
+ *
+ * Authors: Louis Moresi, Shijie Zhong, and Michael Gurnis
+ *
+ * For questions or comments regarding this software, you may contact
+ *
+ * Luis Armendariz <luis at geodynamics.org>
+ * http://geodynamics.org
+ * Computational Infrastructure for Geodynamics (CIG)
+ * California Institute of Technology
+ * 2750 East Washington Blvd, Suite 210
+ * Pasadena, CA 91007
+ *
+ * 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 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
+ */
+
+#include <math.h>
+#include <malloc.h>
+#include <sys/types.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+#include <stdlib.h> /* for "system" command */
+#include <string.h>
+#include <mpi.h>
+#include "hc.h"
+
+
+/* ===============================
+ Initialization of fields .....
+ =============================== */
+
+void convection_initial_temperature_ggrd(struct All_variables *E)
+{
+ int ll, mm, i, j, k, p, node, ii;
+ double temp, temp1, temp2, temp3, base, radius, radius2;
+ unsigned int type;
+
+ FILE *fp;
+ double x1, y1, z1, con;
+
+ int noz2, nfz, in1, in2, in3, instance, nox, noy, noz,nxs,setflag;
+ char input_s[200], output_file[255];
+ float weight, para1, plate_velocity, delta_temp, age;
+
+ char *char_dummy="";
+
+ /* twb additions */
+ double rho_prem;
+ char pfile[1000];
+ double t1,f1,r1,tgrad,tbot,tadd,tz,tmean;
+
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc, mpi_tag=1;
+ int mpi_inmsg, mpi_success_message = 1;
+
+
+ const int dims = E->mesh.nsd;
+
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+ nox = E->lmesh.nox;
+
+
+ setflag = 0;
+
+
+ for(i=1;i<=E->lmesh.nno;i++) {
+ /* init as zeros */
+ E->T[i] = E->C[i] = 0.0;
+ }
+ if(!E->control.restart)
+ { /*
+
+ regular init
+
+ */
+
+ if(E->control.ggrd.use_temp){ /* read T init from grid */
+ if(E->parallel.me==0)
+ fprintf(stderr,"convection_initial_temperature: using GMT grd files for temperatures\n");
+ /*
+
+ read in tempeatures/density from GMT grd files
+
+
+ */
+ if(E->parallel.me > 0){
+ /*
+ wait for the previous processor
+ */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
+ MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(E->parallel.me == 0)
+ fprintf(stderr,"Init: initializing PREM and ggrd files\n");
+ if(E->control.ggrd.temp.scale_with_prem){
+ if(!E->control.Rsphere)
+ myerror("cannot use PREM with Cartesian setting",E);
+ /* initialize PREM */
+ if(prem_read_model(PREM_MODEL_FILE,&E->control.ggrd.temp.prem,
+ (E->parallel.me==0)))
+ myerror("prem initialization",E);
+ }
+ /*
+ initialize the GMT grid files
+ */
+ if(E->control.slab_slice){ /* only slice of x - depth */
+
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.temp.gfile,
+ char_dummy,"",E->control.ggrd.temp.d,
+ (E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* 3-D */
+ if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.temp.gfile,
+ E->control.ggrd.temp.dfile,"",
+ E->control.ggrd.temp.d,(E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack 3-D init error",E);
+ }
+
+
+
+ if(E->parallel.me < E->parallel.nproc-1){
+ /*
+ tell the next processor to go ahead
+ */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
+ (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"convection_initial_temperature: last processor done with grd init\n");
+ }
+ /*
+
+ end MPI synchornization part
+
+ */
+ /*
+
+ interpolate densities to temperature given PREM variations
+
+ */
+ if(E->mesh.bottbc != 0){
+ /* bottom has specified temperature */
+ tbot = E->control.TBCbotval;
+ }else{
+ /*
+ bottom has specified heat flux
+ start with unity bottom temperature
+ */
+ tbot = 1.0;
+ }
+ tmean = (tbot + E->control.TBCtopval)/2.0 + E->control.ggrd.temp.offset;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"tinit: top: %g bot: %g mean: %g scale: %g PREM %i\n",
+ E->control.TBCtopval,tbot,tmean,E->control.ggrd.temp.scale,
+ E->control.ggrd.temp.scale_with_prem);
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
+ if(E->control.slab_slice){
+ if(E->control.Rsphere) {
+ if(E->SX[1][node] <= E->control.slab_theta_bound)
+ /* spherical interpolation */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE);
+ else{
+ if(E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1])
+ tadd = E->control.TBCtopval;
+ else
+ tadd = 1.0;
+ }
+ }else{ /* cartesian interpolation */
+ if(E->X[2][node] <= E->control.slab_theta_bound)
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
+ (double)E->X[3][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE);
+ else{
+ if(E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1])
+ tadd = E->control.TBCtopval;
+ else
+ tadd = 1.0;
+ }
+ }
+ }else{
+ /*
+
+ 3-D
+
+ */
+ if(E->control.Rsphere) /* spherical interpolation */
+ ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],
+ (double)E->SX[1][node],
+ (double)E->SX[2][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE,FALSE);
+ else{ /* cartesian interpolation */
+ ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],
+ (double)E->X[2][node],
+ (double)E->X[3][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE);
+ }
+ }
+ if(E->control.ggrd.temp.scale_with_prem){ /* only works for spherical! */
+ /*
+ get the PREM density at r for additional scaling
+ */
+ prem_get_rho(&rho_prem,(double)E->SX[3][node],&E->control.ggrd.temp.prem);
+ if(rho_prem < 3200.0)
+ rho_prem = 3200.0; /* we don't want the viscosity of water */
+ /*
+ assign temperature
+ */
+ E->T[node] = tmean + tadd * E->control.ggrd.temp.scale *
+ rho_prem / E->data.density;
+ }else{
+ /* no PREM scaling */
+ E->T[node] = tmean + tadd * E->control.ggrd.temp.scale;
+
+ //fprintf(stderr,"z: %11g mean: %11g tadd: %11g scale: %11g T: %11g\n", E->X[3][node],tmean, tadd ,E->control.ggrd_tinit_scale, E->T[node] );
+
+ }
+ /*
+ if we're on the surface or bottom, reassign T and
+ temperature boundary condition if toptbc or bottbsc == 2
+ */
+ if(E->control.Rsphere){
+ if((E->mesh.toptbc == 2) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1])){
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ if((E->mesh.bottbc == 2) && (E->SX[3][node] == E->segment.zzlayer[0])){
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ }else{
+ if((E->mesh.toptbc == 2) && (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ if((E->mesh.bottbc == 2) && (E->X[3][node] == E->segment.zzlayer[0]))
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ /*
+ boundary condition flags
+ */
+ if(!setflag)
+ E->node[node] = E->node[node] | (INTX | INTZ | INTY);
+ } /* end node loop */
+ setflag=1;
+ /* free the structure, not needed anymore */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.temp.d);
+ /*
+ end temperature from GMT grd init
+ */
+ /* end grid init branch */
+ } else {
+
+
+
+ /*
+
+ add a linear profile and potentially perturbations to the temperature
+
+ */
+ if(E->control.CART3D)
+ {
+
+ for(i = 1; i <= noy; i++)
+ for(j = 1; j <= nox; j++)
+ for(k = 1; k <= noz; k++)
+ {
+ node = k + (j - 1) * noz + (i - 1) * noz * nox;
+ x1 = E->X[1][node];
+ z1 = E->X[3][node];
+ y1 = E->X[2][node];
+ /* linear gradient */
+ tz = z1* (E->control.TBCtopval-E->control.TBCbotval) +
+ E->control.TBCbotval;
+ E->T[node] += tz;
+
+ if(E->convection.random_t_init){
+ /* random init */
+ E->T[node] += (-0.5 + drand48()) * E->convection.perturb_mag[0];
+ }else{
+
+ for(p=0;p<E->convection.number_of_perturbations;p++){ /* perturbations */
+
+ E->T[node] += E->convection.perturb_mag[p] *
+ sin(M_PI * (1.0 - z1)) * cos(E->convection.perturb_k[p] * M_PI * x1) *
+ ((E->mesh.nsd != 3) ? 1.0 : cos(E->convection.perturb_k[p] * M_PI * y1));
+ }
+ }
+ if(!setflag)
+ E->node[node] = E->node[node] | (INTX | INTZ | INTY);
+ }
+ setflag=1;
+ }
+ else if(E->control.Rsphere)
+ {
+
+ con = (E->mesh.noz - 1) / (E->sphere.ro - E->sphere.ri);
+ noz2 = (E->mesh.noz - 1) / 2 + 1;
+
+ for(i = 1; i <= noy; i++)
+ for(j = 1; j <= nox; j++)
+ for(k = 1; k <= noz; k++)
+ {
+ ii = k + E->lmesh.nzs - 1;
+ node = k + (j - 1) * noz + (i - 1) * noz * nox;
+ x1 = E->SX[1][node];
+ z1 = E->SX[3][node];
+ y1 = E->SX[2][node];
+ /* linear */
+ tz = (z1 - E->sphere.ri)/(E->sphere.ro - E->sphere.ri) *
+ (E->control.TBCtopval-E->control.TBCbotval) + E->control.TBCbotval;
+ E->T[node] += tz;
+ if(E->convection.random_t_init){
+ /* random init */
+ E->T[node] += (-0.5 + drand48())*E->convection.perturb_mag[0];
+ }else{
+ for(p=0;p < E->convection.number_of_perturbations;p++){ /* perturbations */
+ mm = E->convection.perturb_mm[p];
+ ll = E->convection.perturb_ll[p];
+ con = E->convection.perturb_mag[p];
+
+ E->T[node] += E->convection.perturb_mag[p] * modified_plgndr_a(ll, mm, x1) * cos(mm * y1) *
+ sin(M_PI * (z1 - E->sphere.ri) / (E->sphere.ro - E->sphere.ri));
+ }
+ }
+ if(!setflag)
+ E->node[node] = E->node[node] | (INTX | INTZ | INTY);
+ }
+ setflag =1;
+ } // end for if else if of geometry
+
+ } /* end perturnbation branch */
+ /*
+
+ now deal with composition
+
+ */
+ if(E->control.ggrd.use_comp){ /* read composition init from grid */
+ if(!E->control.composition)
+ myerror("comp grd init but no composition control set!?",E);
+ if(E->parallel.me==0)
+ fprintf(stderr,"convection_initial_temperature: using GMT grd files for composition\n");
+ if(E->parallel.me > 0){
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
+ MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(E->control.slab_slice){
+ /* x - z slice */
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.comp.gfile,
+ char_dummy,"",
+ E->control.ggrd.comp.d,
+ /* (E->parallel.me==0)*/
+ FALSE,FALSE))
+ myerror("grdtrack init error",E);
+ }else{
+ /* 3-D */
+ if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.comp.gfile,
+ E->control.ggrd.comp.dfile,"",
+ E->control.ggrd.comp.d,
+ /* (E->parallel.me==0)*/
+ FALSE,FALSE))
+ myerror("grdtrack init error",E);
+ }
+ if(E->parallel.me < E->parallel.nproc-1){
+ /*
+ tell the next processor to go ahead
+ */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
+ (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"convection_initial_temperature: last processor done with grd init\n");
+ }
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
+ if(E->control.slab_slice){
+ /* slab */
+ if(E->control.Rsphere) {
+ if(E->SX[1][node] <= E->control.slab_theta_bound)
+ /* spherical interpolation */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE);
+ else
+ tadd = 0.0;
+ }else{ /* cartesian interpolation */
+ if(E->X[2][node] <= E->control.slab_theta_bound){
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
+ (double)E->X[3][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE);
+ }else{
+ tadd = 0.0;
+ }
+ }
+ }else{
+ /* 3-D */
+ if(E->control.Rsphere) /* spherical interpolation */
+ ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],
+ (double)E->SX[1][node],
+ (double)E->SX[2][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE,FALSE);
+ else /* cartesian interpolation */
+ ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],
+ (double)E->X[2][node],
+ (double)E->X[3][node],
+ E->control.ggrd.comp.d,&tadd,
+ FALSE);
+ }
+
+ E->C[node] = E->control.ggrd.comp.offset + tadd *
+ E->control.ggrd.comp.scale;
+ }
+ /* free the structure, not needed anymore */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.comp.d);
+ } /* end grid cinit */
+ /*
+
+ check T and C range
+
+ */
+ if(E->control.check_t_irange)
+ for(i=1;i<=E->lmesh.nno;i++){
+ if(E->T[i]>1)E->T[i]=1;
+ if(E->T[i]<0)E->T[i]=0;
+ }
+ if(E->control.check_c_irange)
+ for(i=1;i<=E->lmesh.nno;i++){
+ if(E->C[i]>1)E->C[i]=1;
+ if(E->C[i]<0)E->C[i]=0;
+ }
+
+ if(E->control.composition)
+ convection_initial_markers(E,1);
+ } // end for restart==0
+
+ else if(E->control.restart)
+ {
+ /* restart part */
+ process_restart_tc(E);
+
+ }
+
+ /* make sure temps are assigned as TBC */
+ //assign_T_as_TBC(E,1e-4);
+
+
+ temperatures_conform_bcs(E);
+
+
+ thermal_buoyancy(E);
+
+ return;
+}
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -46,9 +46,12 @@
#include <string.h>
#include "element_definitions.h"
#include "global_defs.h"
-
+#ifdef USE_GGRD
+#include "hc.h"
+#endif
int Emergency_stop;
+
void read_instructions(struct All_variables *E, int argc, char **argv)
{
double start_time;
@@ -73,7 +76,7 @@
* from startup files. (See Parsing.c).
* ================================================== */
- setup_parser(E, argc, argv);
+ setup_parser(E, argv[1]);
if(E->parallel.me == 0)
fprintf(stderr, "ok1\n");
@@ -524,15 +527,18 @@
* record information about the progress of the
* program as it runs
*/
-
+
+#ifdef USE_GZDIR
+ sprintf(logfile, "%s/log%d", E->control.data_file, E->parallel.me);
+#else
sprintf(logfile, "%s.log%d", E->control.data_file, E->parallel.me);
-
+#endif
if((fp = fopen(logfile, "w")) == NULL)
E->fp = stdout;
else
E->fp = fp;
- fprintf(E->fp, "visc_f %g %g\n", E->data.visc_factor, E->control.Q0);
+ //fprintf(E->fp, "visc_f %g %g\n", E->data.visc_factor, E->control.Q0);
if(E->control.NMULTIGRID || E->control.EMULTIGRID)
@@ -637,8 +643,11 @@
char tmp_string[100];
/* first the problem type (defines subsequent behaviour) */
+ int m;
- input_string("Problem", E->control.PROBLEM_TYPE, NULL);
+ m = E->parallel.me;
+
+ input_string("Problem", E->control.PROBLEM_TYPE, NULL, m);
if(strcmp(E->control.PROBLEM_TYPE, "convection") == 0)
{
E->control.CONVECTION = 1;
@@ -659,7 +668,7 @@
set_convection_defaults(E);
}
- input_string("Geometry", E->control.GEOMETRY, NULL);
+ input_string("Geometry", E->control.GEOMETRY, NULL, m);
if(strcmp(E->control.GEOMETRY, "cart2d") == 0)
{
E->control.CART2D = 1;
@@ -691,7 +700,7 @@
set_2dc_defaults(E);
}
- input_string("Solver", E->control.SOLVER_TYPE, NULL);
+ input_string("Solver", E->control.SOLVER_TYPE, NULL, m);
if(strcmp(E->control.SOLVER_TYPE, "cgrad") == 0)
{
E->control.CONJ_GRAD = 1;
@@ -721,57 +730,61 @@
* Default is no information recorded (apart from special things for given applications.
*/
- input_string("datatypes", E->control.which_data_files, "");
- input_string("averages", E->control.which_horiz_averages, "");
- input_string("timelog", E->control.which_running_data, "");
- input_string("observables", E->control.which_observable_data, "");
+ input_string("datatypes", E->control.which_data_files, "", m);
+ input_string("averages", E->control.which_horiz_averages, "", m);
+ input_string("timelog", E->control.which_running_data, "", m);
+ input_string("observables", E->control.which_observable_data, "", m);
- input_string("datafile", E->control.data_file, "initialize");
- input_string("process_command", E->control.output_written_external_command, "");
- input_string("use_scratch", tmp_string, "local");
+ input_string("datafile", E->control.data_file, "initialize", m);
+ input_string("process_command", E->control.output_written_external_command, "", m);
+
+#ifdef USE_GZDIR
+ input_boolean("gzdir",&(E->control.gzdir),"on",m);
+#endif
+ input_string("use_scratch", tmp_string, "local", m);
if(strcmp(tmp_string, "local") == 0)
strcpy(E->control.data_file2, E->control.data_file);
else
sprintf(E->control.data_file2, "/scratch_%s/%s/%s", E->parallel.machinename, tmp_string, E->control.data_file);
- input_boolean("AVS", &(E->control.AVS), "off");
- input_boolean("CONMAN", &(E->control.CONMAN), "off");
+ input_boolean("AVS", &(E->control.AVS), "off", m);
+ input_boolean("CONMAN", &(E->control.CONMAN), "off", m);
if(E->control.NMULTIGRID || E->control.EMULTIGRID)
{
- input_int("mgunitx", &(E->mesh.mgunitx), "1");
- input_int("mgunitz", &(E->mesh.mgunitz), "1");
- input_int("mgunity", &(E->mesh.mgunity), "1");
- input_int("mgunitxl", &(E->lmesh.mgunitx), "1");
- input_int("mgunitzl", &(E->lmesh.mgunitz), "1");
- input_int("mgunityl", &(E->lmesh.mgunity), "1");
- input_int("levels", &(E->mesh.levels), "0");
+ input_int("mgunitx", &(E->mesh.mgunitx), "1", m);
+ input_int("mgunitz", &(E->mesh.mgunitz), "1", m);
+ input_int("mgunity", &(E->mesh.mgunity), "1", m);
+ input_int("mgunitxl", &(E->lmesh.mgunitx), "1", m);
+ input_int("mgunitzl", &(E->lmesh.mgunitz), "1", m);
+ input_int("mgunityl", &(E->lmesh.mgunity), "1", m);
+ input_int("levels", &(E->mesh.levels), "0", m);
}
- input_int("restart", &(E->control.restart), "0");
+ input_int("restart", &(E->control.restart), "0", m);
- input_boolean("regular_grid", &(E->control.ORTHOZ), "off");
+ input_boolean("regular_grid", &(E->control.ORTHOZ), "off", m);
- input_int("stokes_flow_only", &(E->control.stokes), "0");
+ input_int("stokes_flow_only", &(E->control.stokes), "0", m);
- input_boolean("node_assemble", &(E->control.NASSEMBLE), "off");
+ input_boolean("node_assemble", &(E->control.NASSEMBLE), "off", m);
/* general mesh structure */
- input_boolean("parallel_auto", &(E->parallel.automa), "off");
+ input_boolean("parallel_auto", &(E->parallel.automa), "off", m);
if(E->parallel.automa == 0)
{
- input_int("nprocx", &(E->parallel.nprocx), "1");
- input_int("nprocz", &(E->parallel.nprocz), "1");
- input_int("nprocy", &(E->parallel.nprocy), "1");
+ input_int("nprocx", &(E->parallel.nprocx), "1", m);
+ input_int("nprocz", &(E->parallel.nprocz), "1", m);
+ input_int("nprocy", &(E->parallel.nprocy), "1", m);
}
- input_boolean("verbose", &(E->control.verbose), "off");
- input_boolean("see_convergence", &(E->control.print_convergence), "off");
- input_boolean("COMPRESS", &(E->control.COMPRESS), "on");
- input_float("sobtol", &(E->control.sob_tolerance), "0.0001");
+ input_boolean("verbose", &(E->control.verbose), "off", m);
+ input_boolean("see_convergence", &(E->control.print_convergence), "off", m);
+ input_boolean("COMPRESS", &(E->control.COMPRESS), "on", m);
+ input_float("sobtol", &(E->control.sob_tolerance), "0.0001", m);
- input_int("obs_maxlongk", &(E->slice.maxlong), "100,1");
- input_int("obs_minlongk", &(E->slice.minlong), "1,1");
+ input_int("obs_maxlongk", &(E->slice.maxlong), "100,1", m);
+ input_int("obs_minlongk", &(E->slice.minlong), "1,1", m);
/* for layers */
E->viscosity.zlm = 1.0;
@@ -780,79 +793,104 @@
if(E->control.CART3D)
{
- input_float("z_lmantle", &(E->viscosity.zlm), "1.0");
- input_float("z_410", &(E->viscosity.z410), "1.0");
- input_float("z_lith", &(E->viscosity.zlith), "0.0");
+ input_float("z_lmantle", &(E->viscosity.zlm), "1.0", m);
+ input_float("z_410", &(E->viscosity.z410), "1.0", m);
+ input_float("z_lith", &(E->viscosity.zlith), "0.0", m);
}
else if(E->control.Rsphere)
{
- input_float("r_lmantle", &(E->viscosity.zlm), "1.0");
- input_float("r_410", &(E->viscosity.z410), "1.0");
- input_float("r_lith", &(E->viscosity.zlith), "0.0");
+ input_float("r_lmantle", &(E->viscosity.zlm), "1.0", m);
+ input_float("r_410", &(E->viscosity.z410), "1.0", m);
+ input_float("r_lith", &(E->viscosity.zlith), "0.0", m);
}
+
+#ifdef USE_GGRD
+ /* ggrd control */
+ ggrd_init_master(&(E->control.ggrd));
+ input_boolean("ggrd_tinit",&(E->control.ggrd.use_temp),"off", m);
+ input_double("ggrd_tinit_scale",&(E->control.ggrd.temp.scale),"1.0", m);
+ input_boolean("ggrd_scale_with_prem",&(E->control.ggrd.temp.scale_with_prem),"off", m);
+ input_string("ggrd_tinit_gfile",E->control.ggrd.temp.gfile,"", m);
+ input_string("ggrd_tinit_dfile",E->control.ggrd.temp.dfile,"", m);
+ input_double("ggrd_tinit_offset",&(E->control.ggrd.temp.offset),"0.0", m);
+ /* comp */
+ input_boolean("ggrd_cinit",&(E->control.ggrd.use_comp),"off", m);
+ input_double("ggrd_cinit_scale",&(E->control.ggrd.comp.scale),"1.0", m);
+ input_string("ggrd_cinit_gfile",E->control.ggrd.comp.gfile,"", m);
+ input_string("ggrd_cinit_dfile",E->control.ggrd.comp.dfile,"", m);
+ input_double("ggrd_cinit_offset",&(E->control.ggrd.comp.offset),"0.0", m);
+ /* slab slice handling */
+ input_boolean("slab_slice",&(E->control.slab_slice),"off", m);
+ input_float("slab_theta_bound",&(E->control.slab_theta_bound),"1.0", m);
+
+#endif
+
E->control.transT670 = 1500;
E->control.transT410 = 1500;
- input_float("Ra_670", &(E->control.Ra_670), "0.0");
- input_float("clapeyron670", &(E->control.clapeyron670), "0.0");
- input_float("transT670", &(E->control.transT670), "0.0");
- input_float("width670", &(E->control.width670), "0.0");
+ input_float("Ra_670", &(E->control.Ra_670), "0.0", m);
+ input_float("clapeyron670", &(E->control.clapeyron670), "0.0", m);
+ input_float("transT670", &(E->control.transT670), "0.0", m);
+ input_float("width670", &(E->control.width670), "0.0", m);
- input_float("Ra_410", &(E->control.Ra_410), "0.0");
- input_float("clapeyron410", &(E->control.clapeyron410), "0.0");
- input_float("transT410", &(E->control.transT410), "0.0");
- input_float("width410", &(E->control.width410), "0.0");
+ input_float("Ra_410", &(E->control.Ra_410), "0.0", m);
+ input_float("clapeyron410", &(E->control.clapeyron410), "0.0", m);
+ input_float("transT410", &(E->control.transT410), "0.0", m);
+ input_float("width410", &(E->control.width410), "0.0", m);
- input_int("ll_max", &(E->sphere.llmax), "1");
- input_int("nlong", &(E->sphere.noy), "1");
- input_int("nlati", &(E->sphere.nox), "1");
+ input_int("ll_max", &(E->sphere.llmax), "1", m);
+ input_int("nlong", &(E->sphere.noy), "1", m);
+ input_int("nlati", &(E->sphere.nox), "1", m);
- input_int("topvbc", &(E->mesh.topvbc), "0");
- input_int("botvbc", &(E->mesh.botvbc), "0");
- input_int("sidevbc", &(E->mesh.sidevbc), "0");
+ input_int("topvbc", &(E->mesh.topvbc), "0", m);
+ input_int("botvbc", &(E->mesh.botvbc), "0", m);
+ input_int("sidevbc", &(E->mesh.sidevbc), "0", m);
- input_boolean("periodicx", &(E->mesh.periodic_x), "off");
- input_boolean("periodicy", &(E->mesh.periodic_y), "off");
- input_boolean("depthdominated", &(E->control.depth_dominated), "off");
- input_boolean("eqnzigzag", &(E->control.eqn_zigzag), "off");
- input_boolean("eqnviscosity", &(E->control.eqn_viscosity), "off");
+ input_boolean("periodicx", &(E->mesh.periodic_x), "off", m);
+ input_boolean("periodicy", &(E->mesh.periodic_y), "off", m);
+ input_boolean("depthdominated", &(E->control.depth_dominated), "off", m);
+ input_boolean("eqnzigzag", &(E->control.eqn_zigzag), "off", m);
+ input_boolean("eqnviscosity", &(E->control.eqn_viscosity), "off", m);
- input_float("topvbxval", &(E->control.VBXtopval), "0.0");
- input_float("botvbxval", &(E->control.VBXbotval), "0.0");
- input_float("topvbyval", &(E->control.VBYtopval), "0.0");
- input_float("botvbyval", &(E->control.VBYbotval), "0.0");
+ input_float("topvbxval", &(E->control.VBXtopval), "0.0", m);
+ input_float("botvbxval", &(E->control.VBXbotval), "0.0", m);
+ input_float("topvbyval", &(E->control.VBYtopval), "0.0", m);
+ input_float("botvbyval", &(E->control.VBYbotval), "0.0", m);
- input_int("toptbc", &(E->mesh.toptbc), "1");
- input_int("bottbc", &(E->mesh.bottbc), "1");
- input_float("toptbcval", &(E->control.TBCtopval), "0.0");
- input_float("bottbcval", &(E->control.TBCbotval), "1.0");
+ input_int("toptbc", &(E->mesh.toptbc), "1", m);
+ input_int("bottbc", &(E->mesh.bottbc), "1", m);
+ input_float("toptbcval", &(E->control.TBCtopval), "0.0", m);
+ input_float("bottbcval", &(E->control.TBCbotval), "1.0", m);
- input_float("plate_velocity", &(E->control.plate_vel), "0.0");
- input_float("plate_age", &(E->control.plate_age), "0.0");
- input_float("plume_radius", &(E->segment.plume_radius), "0.0");
- input_float("plume_DT", &(E->segment.plume_DT), "0.0");
- input_float("plume_x", &(E->segment.plume_coord[1]), "0.0");
- input_float("plume_y", &(E->segment.plume_coord[2]), "0.0");
- input_float("plume_z", &(E->segment.plume_coord[3]), "0.0");
+ input_float("plate_velocity", &(E->control.plate_vel), "0.0", m);
+ input_float("plate_age", &(E->control.plate_age), "0.0", m);
+ input_float("plume_radius", &(E->segment.plume_radius), "0.0", m);
+ input_float("plume_DT", &(E->segment.plume_DT), "0.0", m);
+ input_float("plume_x", &(E->segment.plume_coord[1]), "0.0", m);
+ input_float("plume_y", &(E->segment.plume_coord[2]), "0.0", m);
+ input_float("plume_z", &(E->segment.plume_coord[3]), "0.0", m);
+ /* check input temp / comp ranges */
+ input_boolean("check_t_irange", &(E->control.check_t_irange), "on", m);
+ input_boolean("check_c_irange", &(E->control.check_c_irange), "on", m);
- input_string("gridxfile", E->mesh.gridfile[1], " ");
- input_string("gridyfile", E->mesh.gridfile[2], " ");
- input_string("gridzfile", E->mesh.gridfile[3], " ");
+ input_string("gridxfile", E->mesh.gridfile[1], " ", m);
+ input_string("gridyfile", E->mesh.gridfile[2], " ", m);
+ input_string("gridzfile", E->mesh.gridfile[3], " ", m);
- input_float("dimenx", &(E->mesh.layer[1]), "nodefault");
- input_float("dimeny", &(E->mesh.layer[2]), "nodefault");
- input_float("dimenz", &(E->mesh.layer[3]), "nodefault");
+ input_float("dimenx", &(E->mesh.layer[1]), "nodefault", m);
+ input_float("dimeny", &(E->mesh.layer[2]), "nodefault", m);
+ input_float("dimenz", &(E->mesh.layer[3]), "nodefault", m);
- input_float("radius_inner", &(E->sphere.ri), "nodefault");
- input_float("radius_outer", &(E->sphere.ro), "nodefault");
- input_float("theta_north", &(E->sphere.ti), "nodefault");
- input_float("theta_south", &(E->sphere.to), "nodefault");
- input_float("fi_west", &(E->sphere.fi), "nodefault");
- input_float("fi_east", &(E->sphere.fo), "nodefault");
+ input_float("radius_inner", &(E->sphere.ri), "nodefault", m);
+ input_float("radius_outer", &(E->sphere.ro), "nodefault", m);
+ input_float("theta_north", &(E->sphere.ti), "nodefault", m);
+ input_float("theta_south", &(E->sphere.to), "nodefault", m);
+ input_float("fi_west", &(E->sphere.fi), "nodefault", m);
+ input_float("fi_east", &(E->sphere.fo), "nodefault", m);
E->sphere.corner[1][1] = E->sphere.ti;
E->sphere.corner[2][1] = E->sphere.to;
E->sphere.corner[1][2] = E->sphere.fi;
@@ -860,79 +898,91 @@
E->sphere.corner[1][3] = E->sphere.ri;
E->sphere.corner[2][3] = E->sphere.ro;
- input_int("nodex", &(E->mesh.nox), "nodefault,1,nomax");
- input_int("nodez", &(E->mesh.noz), "nodefault,1,nomax");
- input_int("nodey", &(E->mesh.noy), "1,1,nomax");
- input_boolean("aug_lagr", &(E->control.augmented_Lagr), "off");
- input_double("aug_number", &(E->control.augmented), "0.0");
+ input_int("nodex", &(E->mesh.nox), "nodefault,1,nomax", m);
+ input_int("nodez", &(E->mesh.noz), "nodefault,1,nomax", m);
+ input_int("nodey", &(E->mesh.noy), "1,1,nomax", m);
+ input_boolean("aug_lagr", &(E->control.augmented_Lagr), "off", m);
+ input_double("aug_number", &(E->control.augmented), "0.0", m);
- input_float("tole_compressibility", &(E->control.tole_comp), "0.0");
- input_boolean("orthogonal", &(E->control.ORTHO), "on");
- input_boolean("crust", &(E->control.crust), "off");
- input_float("crust_width", &(E->crust.width), "0.0");
+ input_boolean("sdepv_print_convergence",&(E->control.sdepv_print_convergence),"off", m);
+ //input_float("tole_compressibility", &(E->control.tole_comp), "0.0", m);
+ input_boolean("orthogonal", &(E->control.ORTHO), "on", m);
+ input_boolean("crust", &(E->control.crust), "off", m);
+ input_float("crust_width", &(E->crust.width), "0.0", m);
- input_int("storage_spacing", &(E->control.record_every), "10");
- input_int("storage_always_before", &(E->control.record_all_until), "5");
+ input_int("storage_spacing", &(E->control.record_every), "10", m);
+ input_int("storage_always_before", &(E->control.record_all_until), "5", m);
- input_boolean("precond", &(E->control.precondition), "off");
- input_boolean("vprecond", &(E->control.vprecondition), "on");
- input_int("mg_cycle", &(E->control.mg_cycle), "2,0,nomax");
- input_int("down_heavy", &(E->control.down_heavy), "1,0,nomax");
- input_int("up_heavy", &(E->control.up_heavy), "1,0,nomax");
- input_double("accuracy", &(E->control.accuracy), "1.0e-4,0.0,1.0");
- input_int("viterations", &(E->control.max_vel_iterations), "250,0,nomax");
+ input_boolean("precond", &(E->control.precondition), "off", m);
+ input_boolean("vprecond", &(E->control.vprecondition), "on", m);
+ input_int("mg_cycle", &(E->control.mg_cycle), "2,0,nomax", m);
+ input_int("down_heavy", &(E->control.down_heavy), "1,0,nomax", m);
+ input_int("up_heavy", &(E->control.up_heavy), "1,0,nomax", m);
+ input_double("accuracy", &(E->control.accuracy), "1.0e-4,0.0,1.0", m);
+ input_int("viterations", &(E->control.max_vel_iterations), "250,0,nomax", m);
- input_int("vhighstep", &(E->control.v_steps_high), "1,0,nomax");
- input_int("vlowstep", &(E->control.v_steps_low), "250,0,nomax");
- input_int("vupperstep", &(E->control.v_steps_upper), "1,0,nomax");
- input_int("piterations", &(E->control.p_iterations), "100,0,nomax");
- input_int("maxsamevisc", &(E->control.max_same_visc), "25,0,nomax");
+ input_int("vhighstep", &(E->control.v_steps_high), "1,0,nomax", m);
+ input_int("vlowstep", &(E->control.v_steps_low), "250,0,nomax", m);
+ input_int("vupperstep", &(E->control.v_steps_upper), "1,0,nomax", m);
+ input_int("piterations", &(E->control.p_iterations), "100,0,nomax", m);
+ input_int("maxsamevisc", &(E->control.max_same_visc), "25,0,nomax", m);
/* data section */
E->data.therm_exp_factor = 1.0;
E->data.therm_diff_factor = 1.0;
E->data.visc_factor = 1.0;
- input_float("ReferenceT", &(E->data.ref_temperature), "2600.0");
- input_float("Q0", &(E->control.Q0), "0.0");
+ input_float("ReferenceT", &(E->data.ref_temperature), "2600.0", m);
+ input_float("Q0", &(E->control.Q0), "0.0", m);
-// input_float("layerd",&(E->data.layer_km),"2870000.0");
+// input_float("layerd",&(E->data.layer_km),"2870000.0", m);
if(E->control.CART3D)
- input_float("layerd", &(E->monitor.length_scale), "2870000.0");
+ input_float("layerd", &(E->monitor.length_scale), "2870000.0", m);
else if(E->control.Rsphere)
- input_float("radius", &(E->monitor.length_scale), "6370000.0");
+ input_float("radius", &(E->monitor.length_scale), "6370000.0", m);
- input_float("gravacc", &(E->data.grav_acc), "9.81");
- input_float("thermexp", &(E->data.therm_exp), "3.28e-5");
- input_float("thermexp_factor", &(E->data.therm_exp_factor), "1.0");
- input_float("visc_factor", &(E->data.visc_factor), "1.0");
- input_float("thermdiff_factor", &(E->data.therm_diff_factor), "1.0");
- input_float("cp", &(E->data.Cp), "1200.0");
- input_float("thermdiff", &(E->data.therm_diff), "8.0e-7");
- input_float("thermcond", &(E->data.therm_cond), "3.168");
- input_float("density", &(E->data.density), "3340.0");
- input_float("mdensity", &(E->data.melt_density), "2800.0");
- input_float("wdensity", &(E->data.density_above), "1030.0");
- input_float("rdensity", &(E->data.res_density), "3295.0");
- input_float("heatflux", &(E->data.surf_heat_flux), "4.4e-2");
- input_float("refvisc", &(E->data.ref_viscosity), "nodefault");
- input_float("meltvisc", &(E->data.melt_viscosity), "nodefault");
- input_float("surf_temp", &(E->data.surf_temp), "0.0");
- input_float("dissipation_number", &(E->data.disptn_number), "0.0");
- input_float("youngs", &(E->data.youngs_mod), "1.0e11");
- input_float("Te", &(E->data.Te), "0.0");
- input_float("Tsol0", &(E->data.T_sol0), "1373.0");
- input_float("dTsoldz", &(E->data.dTsol_dz), "3.4e-3");
- input_float("dTsoldF", &(E->data.dTsol_dF), "440.0");
- input_float("dTdz", &(E->data.dT_dz), "0.48e-3");
- input_float("deltaS", &(E->data.delta_S), "250.0");
- input_float("gasconst", &(E->data.gas_const), "8.3"); /* not much cause to change these ! */
- input_float("gravconst", &(E->data.grav_const), "6.673e-11");
- input_float("permeability", &(E->data.permeability), "3.0e-10");
+
+ input_float("gravacc", &(E->data.grav_acc), "9.81", m);
+ input_float("thermexp", &(E->data.therm_exp), "3.28e-5", m);
+ input_float("thermexp_factor", &(E->data.therm_exp_factor), "1.0", m);
+ input_float("visc_factor", &(E->data.visc_factor), "1.0", m);
+ input_float("thermdiff_factor", &(E->data.therm_diff_factor), "1.0", m);
+ input_float("cp", &(E->data.Cp), "1200.0", m);
+ input_float("thermdiff", &(E->data.therm_diff), "8.0e-7", m);
+ input_float("thermcond", &(E->data.therm_cond), "3.168", m);
+ input_float("density", &(E->data.density), "3340.0", m);
+ input_float("mdensity", &(E->data.melt_density), "2800.0", m);
+ input_float("wdensity", &(E->data.density_above), "1030.0", m);
+ input_float("rdensity", &(E->data.res_density), "3295.0", m);
+ input_float("heatflux", &(E->data.surf_heat_flux), "4.4e-2", m);
+ input_float("refvisc", &(E->data.ref_viscosity), "nodefault", m);
+ input_float("meltvisc", &(E->data.melt_viscosity), "nodefault", m);
+ input_float("surf_temp", &(E->data.surf_temp), "0.0", m);
+ input_float("dissipation_number", &(E->data.disptn_number), "0.0", m);
+ input_float("youngs", &(E->data.youngs_mod), "1.0e11", m);
+ input_float("Te", &(E->data.Te), "0.0", m);
+ input_float("Tsol0", &(E->data.T_sol0), "1373.0", m);
+ input_float("dTsoldz", &(E->data.dTsol_dz), "3.4e-3", m);
+ input_float("dTsoldF", &(E->data.dTsol_dF), "440.0", m);
+ input_float("dTdz", &(E->data.dT_dz), "0.48e-3", m);
+ input_float("deltaS", &(E->data.delta_S), "250.0", m);
+ input_float("gasconst", &(E->data.gas_const), "8.3",m); /* not much cause to change these !
+ input_float("gravconst", &(E->data.grav_const), "6.673e-11", m);
+ input_float("permeability", &(E->data.permeability), "3.0e-10", m);
+
+ /* scaling */
+ E->monitor.time_scale = pow(E->monitor.length_scale, 2.0) /E->data.therm_diff;
+ /* million years */
+ E->monitor.time_scale_ma = E->monitor.time_scale/(3600.0 * 24.0 * 365.25 * 1.0e6);
+
+ E->monitor.velo_scale = E->data.therm_diff / (E->monitor.length_scale);
+ E->monitor.tau_scale = (E->data.ref_viscosity/E->monitor.time_scale); /* scaling stress */
+
+
(E->problem_settings) (E);
viscosity_parameters(E);
Modified: mc/3D/CitcomCU/trunk/src/Makefile
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Makefile 2008-11-07 03:29:17 UTC (rev 13268)
@@ -16,7 +16,8 @@
# Individual Machine Variations
##############################################################################
-COMPRESS=/usr/bin/compress
+#COMPRESS=/usr/bin/compress
+COMPRESS=/bin/gzip
LIB_PATH=
#LIB_LIST= -lmpi
LIB_LIST=
@@ -293,4 +294,7 @@
$(CPP) $(OPTIM) -P Blas_lapack_interfaces.F Blas_lapack_interfaces.f
$(F77) $(OPTIM) $(FOPTIM) -c Blas_lapack_interfaces.f
+Output_gzdir.o: $(HEADER) Output_gzdir.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output_gzdir.c
+
Added: mc/3D/CitcomCU/trunk/src/Makefile.gzdir
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile.gzdir (rev 0)
+++ mc/3D/CitcomCU/trunk/src/Makefile.gzdir 2008-11-07 03:29:17 UTC (rev 13268)
@@ -0,0 +1,314 @@
+#
+# Makefile for "CitcomCU" allowing considerable variations
+#
+
+#
+# NOTES:
+# 1) Macro substitution is different for different version of make, for SGI,
+# CRAY, HP, etc. It seems necessary to go through and replace the
+# CEXT ... FOBJFLAG macros manually.
+#
+# 2) Put FLAGS=$(YOUR_MACHINE_FLAGS) chosen from the list that follows, add
+# machines, where necessary, try and find a better way than just defining
+# "const" to be null.
+#
+##############################################################################
+# Individual Machine Variations
+##############################################################################
+
+#COMPRESS=/usr/bin/compress
+COMPRESS=/bin/gzip
+GMT_DIR = ../../../GMTdev/GMT4.2.1/x86_64/
+NETCDF_DIR = ../../../netcdf-3.6.2/
+HC_DIR = ../../../hc/
+LIB_PATH= -L$(HC_DIR)/objects/x86_64/ -L$(GMT_DIR) -L$(NETCDF_DIR)/lib/
+#LIB_LIST= -lmpi
+LIB_LIST= -lz -lggrd -lhc -lgmt -lnetcdf
+LIB= $(LIB_PATH) $(LIB_LIST) -lm
+
+DEFINES = -DUSE_GZDIR -DUSE_GGRD -DUSE_GMT4 -I$(HC_DIR)/\
+ -I$(GMT_DIR)/include/ -I$(NETCDF_DIR)/include/\
+ -DPREM_MODEL_FILE=\"$(PWD)/prem/prem.dat\"
+
+###################################################################
+# Operating System Variables
+###################################################################
+
+#CC=/usr/local/mpich/bin/mpicc
+#CC=/opt/mpich/bin/mpicc
+#CC=/usr/lib/cmplrs/cc/gemc_cc
+#CC=$(HOME)/usr/local/mpich/bin/mpicc
+CC=mpicc
+F77=f77
+CPP=
+
+CEXT=c
+FEXT=F # which implies further action of cpp
+OBJEXT=o
+FOBJEXT=o
+OBJFLAG=-c
+FOBJFLAG=-c
+
+
+###################################
+# Choose your machine from here.
+###################################
+
+####################################
+# Dec Alpha, OSF 1
+AXPFLAGS= -unsigned -non_shared \
+ -math_library=fast -prefix=all -reentrancy=none -assume=noaccuracy_sensitive \
+ -unsigned_char -extern=strict_refdef -trapuv # -D_INTRINSICS
+AXPLDFLAGS= -unsigned -assume=noaccuracy_sensitive -non_shared -D_INTRINSICS
+AXPOPTIM= -O -Ublas -float -Olimit 1000 # -cord -feedback citcom.feedback
+####################################
+
+####################################
+# CRAY Unicos systems
+CRAYFLAGS=
+CRAYLDFLAGS=
+CRAYOPTIM=
+####################################
+
+####################################
+#IBM AIX systems
+AIXFLAGS= -D__aix__
+AIXLDFLAGS=
+AIXOPTIM= -O2 -qarch=pwr -qopt=pwr -Ublas
+####################################
+
+####################################
+# SUNOS systems
+SUNFLAGS= -D__sunos__ -Dconst=""
+SUNLDFLAGS=
+SUNOPTIM=-O -fsingle
+####################################
+
+####################################
+# Solaris systems
+SOLARISFLAGS= -D__solaris -Dconst="" -I/opt/mpi/include
+SOLARISLDFLAGS=-fast -lsocket -lnsl -lthread
+SOLARISOPTIM=-fast -xO4 -dalign -xtarget=ultra -xarch=v8plus -fsingle
+####################################
+
+####################################
+#HP running HPUX
+HPUXFLAGS=-Dconst=""
+HPUXLDFLAGS=
+HPUXOPTIM=+O3 +Onolimit +Odataprefetch
+#HPUXOPTIM=+O4 +Onolimit +Odataprefetch +Ofastaccess
+#HPUXOPTIM3=+O3 +Onolimit +Odataprefetch
+####################################
+
+####################################
+# SGI with IRIX
+SGIFLAGS=
+SGILDFLAGS=
+SGIOPTIM=-O -fsingle
+####################################
+
+#LinuxFLAGS=-I/usr/local/mpi/include
+#LinuxFLAGS=-I/opt/mpi/include
+#LinuxFLAGS=-I$(HOME)/usr/include
+LinuxFLAGS=
+LinuxLDFLAGS=
+#LinuxOPTIM=-g
+LinuxOPTIM=-O2 $(DEFINES)
+
+####################################
+#PARAGON
+PARAGONFLAGS=-nx
+PARAGONLDFLAGS=
+PARAGONOPTIM=-O4
+####################################
+
+####################################################################
+# Choose the actual flags for your machine here ...
+####################################################################
+
+FLAGS= $(LinuxFLAGS) -DCOMPRESS_BINARY=\"$(COMPRESS)\"
+LDFLAGS= $(LinuxLDFLAGS)
+OPTIM= $(LinuxOPTIM)
+
+#FLAGS= $(SOLARISFLAGS) -DCOMPRESS_BINARY=\"$(COMPRESS)\"
+#LDFLAGS= $(SOLARISLDFLAGS)
+#OPTIM= $(SOLARISOPTIM)
+
+####################################################################
+# CitcomCU's files .....
+####################################################################
+
+CFILES= Advection_diffusion.c\
+ Boundary_conditions.c\
+ Citcom.c\
+ Composition_adv.c\
+ Construct_arrays.c\
+ Convection.c\
+ Drive_solvers.c\
+ Element_calculations.c\
+ Geometry_cartesian.c\
+ General_matrix_functions.c\
+ Global_operations.c\
+ Stokes_flow_Incomp.c\
+ Instructions.c\
+ Nodal_mesh.c\
+ Output.c\
+ Pan_problem_misc_functions.c\
+ Parallel_related.c\
+ Parsing.c\
+ Phase_change.c\
+ Process_buoyancy.c\
+ Process_velocity.c\
+ Profiling.c\
+ Shape_functions.c\
+ Size_does_matter.c\
+ Solver_multigrid.c\
+ Solver_conj_grad.c\
+ Sphere_harmonics.c\
+ Topo_gravity.c\
+ Viscosity_structures.c\
+ Output_gzdir.c\
+ Ggrd_handling.c
+
+
+
+HEADER = element_definitions.h\
+ global_defs.h\
+ viscosity_descriptions.h\
+ advection.h
+
+
+FFILES=#Blas_lapack_interfaces.F
+
+OBJFILES=$(CFILES:.c=.o) $(FFILES:.f=.o)
+
+
+####################################################################
+# Makefile rules follow
+####################################################################
+
+default: citcom.mpi
+
+citcom.mpi: $(OBJFILES) $(HEADER) Makefile
+ $(CC) $(OPTIM) $(FLAGS) $(LDFLAGS) -o citcom.mpi $(OBJFILES) $(FFTLIB) $(LIB)
+
+clean:
+ rm -f *.o
+
+clean-all:
+ rm -f *.o citcom.mpi
+
+smaller:
+ compress $(CFILES)
+
+larger:
+ uncompress $(CFILES)
+
+
+####################################################################
+# File dependencies follow
+####################################################################
+
+global_defs.h: viscosity_descriptions.h advection.h Convection_variables.h
+
+# The following entries can probably be automated from $CFILES etc
+
+Advection_diffusion.o: $(HEADER) Advection_diffusion.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Advection_diffusion.c
+#
+Boundary_conditions.o: $(HEADER) Boundary_conditions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Boundary_conditions.c
+#
+Citcom.o: $(HEADER) Citcom.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Citcom.c
+#
+Composition_adv.o: $(HEADER) Composition_adv.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Composition_adv.c
+#
+Construct_arrays.o: $(HEADER) Construct_arrays.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Construct_arrays.c
+#
+Convection.o: $(HEADER) Convection.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Convection.c
+#
+Drive_solvers.o: $(HEADER) Drive_solvers.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Drive_solvers.c
+#
+Element_calculations.o: $(HEADER) Element_calculations.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Element_calculations.c
+
+General_matrix_functions.o: $(HEADER) General_matrix_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) General_matrix_functions.c
+
+Geometry_cartesian.o: $(HEADER) Geometry_cartesian.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Geometry_cartesian.c
+
+Global_operations.o: $(HEADER) Global_operations.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Global_operations.c
+
+Instructions.o: $(HEADER) Instructions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Instructions.c
+
+Nodal_mesh.o: $(HEADER) Nodal_mesh.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Nodal_mesh.c
+
+Output.o: $(HEADER) Output.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output.c
+
+Pan_problem_misc_functions.o: $(HEADER) Pan_problem_misc_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Pan_problem_misc_functions.c
+
+Parallel_related.o: $(HEADER) Parallel_related.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Parallel_related.c
+
+Parsing.o: $(HEADER) Parsing.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Parsing.c
+
+Phase_change.o: $(HEADER) Phase_change.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Phase_change.c
+
+Process_velocity.o: $(HEADER) Process_velocity.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Process_velocity.c
+
+Process_buoyancy.o: $(HEADER) Process_buoyancy.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Process_buoyancy.c
+
+Profiling.o: $(HEADER) Profiling.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Profiling.c
+
+Shape_functions.o: $(HEADER) Shape_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Shape_functions.c
+
+Size_does_matter.o: $(HEADER) Size_does_matter.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Size_does_matter.c
+
+Solver_conj_grad.o: $(HEADER) Solver_conj_grad.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Solver_conj_grad.c
+
+Solver_multigrid.o: $(HEADER) Solver_multigrid.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Solver_multigrid.c
+
+Sphere_harmonics.o: $(HEADER) Sphere_harmonics.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Sphere_harmonics.c
+
+Stokes_flow_Incomp.o: $(HEADER) Stokes_flow_Incomp.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Stokes_flow_Incomp.c
+
+Topo_gravity.o: $(HEADER) Topo_gravity.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Topo_gravity.c
+
+Viscosity_structures.o: $(HEADER) Viscosity_structures.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Viscosity_structures.c
+
+Blas_lapack_interfaces.f: Blas_lapack_interfaces.F
+ $(CPP) $(OPTIM) -P Blas_lapack_interfaces.F Blas_lapack_interfaces.f
+ $(F77) $(OPTIM) $(FOPTIM) -c Blas_lapack_interfaces.f
+
+Output_gzdir.o: $(HEADER) Output_gzdir.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output_gzdir.c
+
+
+Ggrd_handling.o: $(HEADER) Ggrd_handling.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Ggrd_handling.c
+
+
Modified: mc/3D/CitcomCU/trunk/src/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Nodal_mesh.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Nodal_mesh.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -59,9 +59,11 @@
//int n00, nox, noz, noy, fn;
int nox, noz, noy;
double rad_conv;
-
+ int m;
//const int dims = E->mesh.nsd;
+ m = E->parallel.me;
+
rad_conv = M_PI / 180;
nox = E->mesh.nox;
@@ -96,31 +98,31 @@
if(E->control.CART3D)
{
- input_int("z_grid_layers", &(E->segment.zlayers), "1");
- input_float_vector("zz", E->segment.zlayers, (E->segment.zzlayer));
- input_int_vector("nz", E->segment.zlayers, (E->segment.nzlayer));
+ input_int("z_grid_layers", &(E->segment.zlayers), "1", m);
+ input_float_vector("zz", E->segment.zlayers, (E->segment.zzlayer), m);
+ input_int_vector("nz", E->segment.zlayers, (E->segment.nzlayer), m);
- input_int("x_grid_layers", &(E->segment.xlayers), "1");
- input_float_vector("xx", E->segment.xlayers, (E->segment.xxlayer));
- input_int_vector("nx", E->segment.xlayers, (E->segment.nxlayer));
+ input_int("x_grid_layers", &(E->segment.xlayers), "1", m);
+ input_float_vector("xx", E->segment.xlayers, (E->segment.xxlayer), m);
+ input_int_vector("nx", E->segment.xlayers, (E->segment.nxlayer), m);
- input_int("y_grid_layers", &(E->segment.ylayers), "1");
- input_float_vector("yy", E->segment.ylayers, (E->segment.yylayer));
- input_int_vector("ny", E->segment.ylayers, (E->segment.nylayer));
+ input_int("y_grid_layers", &(E->segment.ylayers), "1", m);
+ input_float_vector("yy", E->segment.ylayers, (E->segment.yylayer), m);
+ input_int_vector("ny", E->segment.ylayers, (E->segment.nylayer), m);
}
else if(E->control.Rsphere)
{
- input_int("r_grid_layers", &(E->segment.zlayers), "1");
- input_float_vector("rr", E->segment.zlayers, (E->segment.zzlayer));
- input_int_vector("nr", E->segment.zlayers, (E->segment.nzlayer));
+ input_int("r_grid_layers", &(E->segment.zlayers), "1", m);
+ input_float_vector("rr", E->segment.zlayers, (E->segment.zzlayer), m);
+ input_int_vector("nr", E->segment.zlayers, (E->segment.nzlayer), m);
- input_int("t_grid_layers", &(E->segment.xlayers), "1");
- input_float_vector("tt", E->segment.xlayers, (E->segment.xxlayer));
- input_int_vector("nt", E->segment.xlayers, (E->segment.nxlayer));
+ input_int("t_grid_layers", &(E->segment.xlayers), "1", m);
+ input_float_vector("tt", E->segment.xlayers, (E->segment.xxlayer), m);
+ input_int_vector("nt", E->segment.xlayers, (E->segment.nxlayer), m);
- input_int("f_grid_layers", &(E->segment.ylayers), "1");
- input_float_vector("ff", E->segment.ylayers, (E->segment.yylayer));
- input_int_vector("nf", E->segment.ylayers, (E->segment.nylayer));
+ input_int("f_grid_layers", &(E->segment.ylayers), "1", m);
+ input_float_vector("ff", E->segment.ylayers, (E->segment.yylayer), m);
+ input_int_vector("nf", E->segment.ylayers, (E->segment.nylayer), m);
}
Added: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c (rev 0)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -0,0 +1,557 @@
+/*
+ * CitcomCU is a Finite Element Code that solves for thermochemical
+ * convection within a three dimensional domain appropriate for convection
+
+* within the Earth's mantle. Cartesian and regional-spherical geometries
+ * are implemented. See the file README contained with this distribution
+ * for further details.
+ *
+ * Copyright (C) 1994-2005 California Institute of Technology
+ * Copyright (C) 2000-2005 The University of Colorado
+ *
+ * Authors: Louis Moresi, Shijie Zhong, and Michael Gurnis
+ *
+ * For questions or comments regarding this software, you may contact
+ *
+ * Luis Armendariz <luis at geodynamics.org>
+ * http://geodynamics.org
+ * Computational Infrastructure for Geodynamics (CIG)
+ * California Institute of Technology
+ * 2750 East Washington Blvd, Suite 210
+ * Pasadena, CA 91007
+ *
+ * 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 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
+
+
+
+
+ minor changes by twb to change output format
+
+*/
+
+/* Routine to process the output of the finite element cycles
+ and to turn them into a coherent suite files */
+
+
+#include <fcntl.h>
+#include <math.h>
+#include <malloc.h>
+#include <stdlib.h> /* for "system" command */
+#ifndef __sunos__ /* string manipulations */
+#include <strings.h>
+#else
+#include <string.h>
+#endif
+
+
+#include <zlib.h>
+
+#include "element_definitions.h"
+#include "global_defs.h"
+
+gzFile *safe_gzopen(char *,char *);
+FILE *safe_fopen(char *,char *);
+void *safe_malloc (size_t );
+void calc_cbase_at_tp(float , float , float *);
+void convert_pvec_to_cvec(float ,float , float , float *,
+ float *);
+void rtp2xyz(float r, float , float , float *);
+
+
+void output_velo_related_gzdir(E, file_number)
+ struct All_variables *E;
+ int file_number;
+{
+ int el, els, i, j, k, ii, m, node, fd,snode;
+ int nox, noz, noy, nfx, nfz, nfy1, nfy2, size1, size2,offset;
+ static int vtkout = 1;
+
+ char output_file[1000];
+ float cvec[3],locx[3];
+ static float *SV, *EV;
+ static int been_here = 0;
+ float vs;
+ static int vtk_base_init = 0; /* for spherical */
+ static int vtk_pressure_out = 1,
+ vtk_viscosity_out = 1;
+ int vtk_comp_out;
+
+
+
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+
+ const int nno = E->mesh.nno;
+
+ gzFile gzout;
+
+
+
+
+
+
+ if(E->control.composition)
+ vtk_comp_out = 1;
+ else
+ vtk_comp_out = 0;
+
+ /* make a directory */
+ if(E->parallel.me == 0){
+ fprintf(stderr,"Output_gzdir: processing output\n");
+
+ sprintf(output_file,"mkdir %s/%d 2> /dev/null",
+ E->control.data_file2,file_number);
+ system(output_file);
+ fprintf(stderr,"making directory: %s\n",output_file);
+ }
+ /* and wait for the other jobs */
+ parallel_process_sync();
+
+
+ if(been_here == 0){
+ /*
+ switch off composition, if not needed
+ */
+ if(!E->control.composition)
+ vtk_comp_out = 0;
+
+ /*
+
+ initial nodeal coordinate output
+
+ */
+ if(E->control.COMPRESS){
+ /* compressed I/O */
+ /*
+ all nodes
+ */
+ sprintf(output_file, "%s/coord.%d.gz", E->control.data_file2, E->parallel.me);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d\n", E->lmesh.nno);
+ if(!E->control.Rsphere)
+ for(i = 1; i <= E->lmesh.nno; i++){
+ gzprintf(gzout, "%.5e %.5e %.5e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
+ }
+ else
+ for(i = 1; i <= E->lmesh.nno; i++)
+ gzprintf(gzout, "%.5e %.5e %.5e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
+ gzclose(gzout);
+ /*
+ surface nodes
+ */
+ sprintf(output_file, "%s/scoord.%d.gz", E->control.data_file2, E->parallel.me);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d\n", E->lmesh.nsf);
+ for(i = 1; i <= E->lmesh.nsf; i++){
+ node = E->surf_node[i];
+ if(!E->control.Rsphere) /* leave in redundancy for now */
+ gzprintf(gzout, "%.5e %.5e %.5e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
+ else
+ gzprintf(gzout, "%.5e %.5e %.5e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
+ }
+ gzclose(gzout);
+ }else{ /* regular I/O */
+ /*
+
+ all nodes
+
+ */
+ sprintf(output_file, "%s/coord.%d", E->control.data_file2, E->parallel.me);
+ E->filed[13] = safe_fopen(output_file, "w");
+ fprintf(E->filed[13], "%6d\n", E->lmesh.nno);
+ if(!E->control.Rsphere)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
+ else
+ for(i = 1; i <= E->lmesh.nno; i++)
+ fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
+ fclose(E->filed[13]);
+ if((E->parallel.me_loc[3] == E->parallel.nprocz - 1) || /* top */
+ (E->parallel.me_loc[3] == 0)){ /* bottom */
+ /*
+ surface nodes
+ */
+ sprintf(output_file, "%s/scoord.%d", E->control.data_file2, E->parallel.me);
+ E->filed[13] = safe_fopen(output_file, "w");
+ fprintf(E->filed[13], "%6d\n", E->lmesh.nsf);
+ for(i = 1; i <= E->lmesh.nsf; i++){
+ node = E->surf_node[snode];
+ if(!E->control.Rsphere) /* leave in redundancy for now */
+ fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
+ else
+ fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
+ }
+ fclose(E->filed[13]);
+ }
+ }
+ if(vtkout){
+ /*
+ vtk output, always compressed
+ */
+ /*
+ node coordinates
+ */
+
+
+ sprintf(output_file,"%s/vtk_ecor.%d.gz",E->control.data_file2,E->parallel.me);
+ gzout = safe_gzopen(output_file,"w");
+ if(E->control.Rsphere){ /* spherical */
+ if(E->parallel.me == 0)
+ fprintf(stderr, "converting spherical to cartesian vtk coords to %s \n",output_file);
+ for(i=1;i <= E->lmesh.nno;i++) {
+ rtp2xyz((float)E->SX[3][i],(float)E->SX[1][i],(float)E->SX[2][i],locx);
+ gzprintf(gzout,"%9.6f %9.6f %9.6f\n",locx[0],locx[1],locx[2]);
+ }
+ }else{ /* cartesian */
+ if(E->parallel.me == 0)
+ fprintf(stderr, "writing cartesian vtk coords to %s \n",output_file);
+ for(i=1;i <= E->lmesh.nno;i++){
+ gzprintf(gzout,"%9.6f %9.6f %9.6f\n",E->X[1][i],E->X[2][i],E->X[3][i]);
+ }
+ }
+ gzclose(gzout);
+
+
+ /*
+ connectivity
+ */
+ offset = E->lmesh.nno * E->parallel.me - 1;
+ sprintf(output_file,"%s/vtk_econ.%d.gz",E->control.data_file2,E->parallel.me);
+ gzout = safe_gzopen(output_file,"w");
+ for(i=1;i<= E->lmesh.nel;i++) {
+ gzprintf(gzout,"%2i\t",ends);
+ /*
+ need to add offset according to the processor for global
+ node numbers
+ */
+ gzprintf(gzout,"%6i %6i %6i %6i %6i %6i %6i %6i\n",
+ E->ien[i].node[1]+offset,E->ien[i].node[2]+offset,
+ E->ien[i].node[3]+offset,E->ien[i].node[4]+offset,
+ E->ien[i].node[5]+offset,E->ien[i].node[6]+offset,
+ E->ien[i].node[7]+offset,E->ien[i].node[8]+offset);
+ }
+ gzclose(gzout);
+ if(E->parallel.me == 0)
+ fprintf(stderr, "writing element connectivity %s \n",output_file);
+ /* end vtk branch */
+ }
+ /* end init loop */
+ been_here++;
+ }
+
+ if(E->parallel.me < E->parallel.nprocz)
+ {
+ sprintf(output_file, "%s/%d/ave_r.%d.%d", E->control.data_file2, file_number,
+ E->parallel.me, file_number);
+ E->filed[10] = safe_fopen(output_file, "w");
+ /* header line */
+ fprintf(E->filed[10], "# %6d %6d %.5e %.5e %.5e %.4e %.4e %.5e %.5e\n",
+ E->lmesh.noz, E->advection.timesteps, E->monitor.elapsed_time,
+ E->slice.Nut, E->slice.Nub, E->data.T_adi0, E->data.T_adi1,
+ E->monitor.Sigma_interior, E->monitor.Sigma_max);
+ /*
+ added first column of radial coordinate TWB
+ */
+ if(!E->control.Rsphere){ /* cartesian */
+ for(i = 1; i <= E->lmesh.noz; i++){
+ fprintf(E->filed[10], "%.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e\n",
+ E->X[3][i],E->Have.T[i], E->Have.vrms[i], E->Have.Vi[i], E->Have.Rho[i],
+ E->Have.F[i], E->Have.f[i], E->Have.C[i], E->Have.Tadi[i]);
+ }
+ }else{ /* spherical */
+ for(i = 1; i <= E->lmesh.noz; i++){
+ fprintf(E->filed[10], "%.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e\n",
+ E->SX[3][i],E->Have.T[i], E->Have.vrms[i], E->Have.Vi[i], E->Have.Rho[i],
+ E->Have.F[i], E->Have.f[i], E->Have.C[i], E->Have.Tadi[i]);
+ }
+ }
+ fclose(E->filed[10]);
+ }
+
+
+ if((file_number % E->control.record_every) == 0)
+ {
+ if(vtkout){
+ /*
+ new vtk output
+
+ */
+ /* temperature at nodes */
+ sprintf(output_file,"%s/%d/t.%d.%d.gz",
+ E->control.data_file2,file_number,E->parallel.me,file_number);
+ gzout=safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%d %d %.5e\n",
+ file_number,E->lmesh.nno,
+ E->monitor.elapsed_time);
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno); /* for backward compatibility */
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e\n",E->T[i]);
+ gzclose(gzout);
+ /* velocities */
+ if(E->control.Rsphere && (!vtk_base_init)){
+ /*
+
+ STILL NEED TO CHECK THE CONVENTION OF ORDERING FOR SPHERICAL HERE!
+ */
+
+ /* need to init spherical/cartesian conversion vectors */
+ E->sphere.vtk_base = (float *)safe_malloc(sizeof(float)*E->lmesh.nno*9);
+ for(k=0,i=1;i<=E->lmesh.nno;i++,k+=9)
+ calc_cbase_at_tp(E->SX[1][i],E->SX[2][i],(E->sphere.vtk_base+k));
+ vtk_base_init = 1;
+ }
+ /* write velocities to file */
+ sprintf(output_file,"%s/%d/vtk_v.%d.%d.gz",
+ E->control.data_file2,file_number, E->parallel.me,file_number);
+
+ /*
+ */
+ gzout=safe_gzopen(output_file,"w");
+ if(E->control.Rsphere){
+ for(k=0,i=1;i<=E->lmesh.nno;i++,k+=9) {
+ /*
+ convert r,theta,phi vector to x,y,z at base location, output is in
+ cm/yr
+ */
+ convert_pvec_to_cvec(E->V[3][i],E->V[1][i],E->V[2][i],(E->sphere.vtk_base+k),cvec);
+ /* output of cartesian vector */
+ gzprintf(gzout,"%.6e %.6e %.6e\n",cvec[0],cvec[1],cvec[2]);
+ }
+ }else{ /* cartesian output in cm/yr */
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e %.6e %.6e\n",E->V[1][i],E->V[2][i],E->V[3][i]);
+ }
+ gzclose(gzout);
+ if(vtk_pressure_out){
+ /*
+ pressure at nodes
+ */
+ sprintf(output_file,"%s/%d/p.%d.%d.gz",E->control.data_file2,
+ file_number, E->parallel.me,file_number);
+ gzout = safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%d %d %.5e\n",file_number,E->lmesh.nno,E->monitor.elapsed_time);
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e\n",E->NP[i]);
+ gzclose(gzout);
+ }
+ if(vtk_viscosity_out){
+ /*
+
+ viscosity output
+
+ */
+ /* write to file */
+ sprintf(output_file,"%s/%d/visc.%d.%d.gz",
+ E->control.data_file2,file_number, E->parallel.me,file_number);
+ gzout=safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++){
+ gzprintf(gzout,"%.4e\n",E->VI[E->mesh.levmax][i]);
+ }
+ gzclose(gzout);
+ }
+ if(vtk_comp_out){
+ /*
+
+ composition output
+
+ */
+
+ /* write to file */
+ sprintf(output_file,"%s/%d/c.%d.%d.gz",
+ E->control.data_file2,file_number, E->parallel.me,file_number);
+ gzout=safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.4e\n",E->C[i]);
+ gzclose(gzout);
+ }
+ /*
+ end vtk out
+ */
+ }else{
+ /*
+
+
+ old output format
+
+ */
+ if(E->control.COMPRESS){
+ sprintf(output_file, "%s/temp.%d.%d.gz", E->control.data_file2, E->parallel.me, file_number);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d %6d %.5e\n", E->lmesh.nno, E->advection.timesteps,
+ E->monitor.elapsed_time);
+ if((file_number % 2* E->control.record_every) == 0)
+ {
+ if(E->control.composition == 0)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // gzprintf(gzout,"%.5e %.4e %.4e %.5e %.5e %.5e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
+ // gzprintf(gzout,"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ gzprintf(gzout, "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ else if(E->control.composition)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // gzprintf(gzout,"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ gzprintf(gzout, "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ }
+ else
+ {
+ if(E->control.composition == 0)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // gzprintf(gzout,"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ gzprintf(gzout, "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ else if(E->control.composition)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // gzprintf(gzout,"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ gzprintf(gzout, "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ }
+ gzclose(gzout);
+
+ }else{
+ sprintf(output_file, "%s/temp.%d.%d", E->control.data_file2, E->parallel.me, file_number);
+ E->filed[10] = safe_fopen(output_file, "w");
+ fprintf(E->filed[10], "%6d %6d %.5e\n", E->lmesh.nno, E->advection.timesteps,
+ E->monitor.elapsed_time);
+ if(file_number % (2 * E->control.record_every) == 0)
+ {
+ if(E->control.composition == 0)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // fprintf(E->filed[10],"%.5e %.4e %.4e %.5e %.5e %.5e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
+ // fprintf(E->filed[10],"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ fprintf(E->filed[10], "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ else if(E->control.composition)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // fprintf(E->filed[10],"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ fprintf(E->filed[10], "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ }
+ else
+ {
+ if(E->control.composition == 0)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // fprintf(E->filed[10],"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ fprintf(E->filed[10], "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ else if(E->control.composition)
+ for(i = 1; i <= E->lmesh.nno; i++)
+ // fprintf(E->filed[10],"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ fprintf(E->filed[10], "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ }
+ fclose(E->filed[10]);
+ }
+
+ /*
+
+ end old field output
+
+ */
+ }
+ /*
+
+ surface and bottom properties
+
+ */
+ if(E->parallel.me_loc[3] == E->parallel.nprocz - 1) /* top */
+ {
+ if(!E->control.COMPRESS){
+ sprintf(output_file, "%s/%d/th_t.%d.%d",
+ E->control.data_file2, file_number, E->parallel.me, file_number);
+ E->filed[11] = safe_fopen(output_file, "w");
+ fprintf(E->filed[11], "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
+ for(i = 1; i <= E->lmesh.nsf; i++)
+ fprintf(E->filed[11], "%.5e %.5e %.5e %.5e\n",
+ E->slice.tpg[i], E->slice.shflux[i], E->Fas410_b[i], E->Fas670_b[i]);
+ fclose(E->filed[11]);
+ }else{
+ sprintf(output_file, "%s/%d/th_t.%d.%d.gz",
+ E->control.data_file2, file_number, E->parallel.me, file_number);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
+ for(i = 1; i <= E->lmesh.nsf; i++)
+ gzprintf(gzout, "%.5e %.5e %.5e %.5e\n",
+ E->slice.tpg[i], E->slice.shflux[i], E->Fas410_b[i], E->Fas670_b[i]);
+ gzclose(gzout);
+ }
+ }
+ if(E->parallel.me_loc[3] == 0) /* bottom */
+ {
+ if(!E->control.COMPRESS){
+ sprintf(output_file, "%s/%d/th_b.%d.%d",
+ E->control.data_file2, file_number, E->parallel.me, file_number);
+ E->filed[11] = safe_fopen(output_file, "w");
+ fprintf(E->filed[11], "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nub);
+ for(i = 1; i <= E->lmesh.nsf; i++)
+ fprintf(E->filed[11], "%.5e %.5e %.5e %.5e\n",
+ E->slice.tpgb[i], E->slice.bhflux[i], E->Fas410_b[i], E->Fas670_b[i]);
+ fclose(E->filed[11]);
+ }else{
+ sprintf(output_file, "%s/%d/th_b.%d.%d.gz",
+ E->control.data_file2, file_number, E->parallel.me, file_number);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d %6d %.5e %.5e\n",
+ E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nub);
+ for(i = 1; i <= E->lmesh.nsf; i++)
+ gzprintf(gzout, "%.5e %.5e %.5e %.5e\n", E->slice.tpgb[i], E->slice.bhflux[i], E->Fas410_b[i], E->Fas670_b[i]);
+ gzclose(gzout);
+ }
+ }
+
+ /* end output number control */
+ }
+
+ /*
+
+ tracer output only every 10th step
+
+ */
+ if((0) && (E->control.composition) && ((file_number % (10*E->control.record_every)) == 0))
+ {
+ /* output of markers */
+
+ if(E->control.COMPRESS){ /* compressed output */
+ sprintf(output_file, "%s/%d/traces.%d.gz", E->control.data_file2, file_number,E->parallel.me);
+ gzout = safe_gzopen(output_file, "w");
+ gzprintf(gzout, "%6d %6d %.5e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
+ for(i = 1; i <= E->advection.markers; i++)
+ gzprintf(gzout, "%g %g %g %d %d\n", E->XMC[1][i], E->XMC[2][i], E->XMC[3][i], E->CElement[i], E->C12[i]);
+ for(i = 1; i <= E->lmesh.nel; i++)
+ fprintf(gzout, "%g\n", E->CE[i]);
+ gzclose(gzout);
+ }else{
+ sprintf(output_file, "%s/%d/traces.%d", E->control.data_file2, file_number, E->parallel.me);
+ E->filed[10] = safe_fopen(output_file, "w");
+ fprintf(E->filed[10], "%6d %6d %.5e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
+ for(i = 1; i <= E->advection.markers; i++)
+ fprintf(E->filed[10], "%g %g %g %d %d\n", E->XMC[1][i], E->XMC[2][i], E->XMC[3][i], E->CElement[i], E->C12[i]);
+ for(i = 1; i <= E->lmesh.nel; i++)
+ fprintf(E->filed[10], "%g\n", E->CE[i]);
+ fclose(E->filed[10]);
+ }
+ }
+
+if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
+ return;
+}
+
+/* safe gzopen function */
+gzFile *safe_gzopen(char *name,char *mode)
+{
+ gzFile *tmp;
+ char m2[300];
+ if((tmp=(gzFile *)gzopen(name,mode))==NULL){
+ fprintf(stderr,"error: cannot gzopen file %s, exiting\n",
+ name);
+ parallel_process_termination();
+ }
+ return ((gzFile *)tmp);
+}
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -180,7 +180,7 @@
/* Define field name, read parameter file to determine file name and column number */
sprintf(input_token, "previous_%s_file", name);
- if(!input_string(input_token, filename, "initialize"))
+ if(!input_string(input_token, filename, "initialize",E->parallel.me))
{
fprintf(E->fp, "No previous %s information found in input file\n", name);
fflush(E->fp);
@@ -484,27 +484,29 @@
//int in1, in2, in3;
//int number, node;
//int combine_option;
-
+ int m;
+ m=E->parallel.me;
+
sprintf(read_string, "%s_rect", name);
- input_int(read_string, &(RECT->numb), "0");
+ input_int(read_string, &(RECT->numb), "0",m);
sprintf(read_string, "%s_rectx1", name);
- input_float_vector(read_string, RECT->numb, RECT->x1);
+ input_float_vector(read_string, RECT->numb, RECT->x1,m);
sprintf(read_string, "%s_rectx2", name);
- input_float_vector(read_string, RECT->numb, RECT->x2);
+ input_float_vector(read_string, RECT->numb, RECT->x2,m);
sprintf(read_string, "%s_rectz1", name);
- input_float_vector(read_string, RECT->numb, RECT->z1);
+ input_float_vector(read_string, RECT->numb, RECT->z1,m);
sprintf(read_string, "%s_rectz2", name);
- input_float_vector(read_string, RECT->numb, RECT->z2);
+ input_float_vector(read_string, RECT->numb, RECT->z2,m);
sprintf(read_string, "%s_recty1", name);
- input_float_vector(read_string, RECT->numb, RECT->y1);
+ input_float_vector(read_string, RECT->numb, RECT->y1,m);
sprintf(read_string, "%s_recty2", name);
- input_float_vector(read_string, RECT->numb, RECT->y2);
+ input_float_vector(read_string, RECT->numb, RECT->y2,m);
sprintf(read_string, "%s_recthw", name);
- input_float_vector(read_string, RECT->numb, RECT->halfw);
+ input_float_vector(read_string, RECT->numb, RECT->halfw,m);
sprintf(read_string, "%s_rectmag", name);
- input_float_vector(read_string, RECT->numb, RECT->mag);
+ input_float_vector(read_string, RECT->numb, RECT->mag,m);
sprintf(read_string, "%s_rectovl", name);
- input_char_vector(read_string, RECT->numb, RECT->overlay);
+ input_char_vector(read_string, RECT->numb, RECT->overlay,m);
if(parse_and_apply)
field_arbitrary_rectangle(E, RECT, field, BC, bcbitf, bcmask_on, bcmask_off);
@@ -599,23 +601,24 @@
//int in1;
//int number, node;
//int combine_option;
+ int m = E->parallel.me;
sprintf(read_string, "%s_circ", name);
- input_int(read_string, &(CIRC->numb), "0");
+ input_int(read_string, &(CIRC->numb), "0",m);
sprintf(read_string, "%s_circx", name);
- input_float_vector(read_string, CIRC->numb, CIRC->x);
+ input_float_vector(read_string, CIRC->numb, CIRC->x,m);
sprintf(read_string, "%s_circz", name);
- input_float_vector(read_string, CIRC->numb, CIRC->z);
+ input_float_vector(read_string, CIRC->numb, CIRC->z,m);
sprintf(read_string, "%s_circy", name);
- input_float_vector(read_string, CIRC->numb, CIRC->y);
+ input_float_vector(read_string, CIRC->numb, CIRC->y,m);
sprintf(read_string, "%s_circrad", name);
- input_float_vector(read_string, CIRC->numb, CIRC->rad);
+ input_float_vector(read_string, CIRC->numb, CIRC->rad,m);
sprintf(read_string, "%s_circmag", name);
- input_float_vector(read_string, CIRC->numb, CIRC->mag);
+ input_float_vector(read_string, CIRC->numb, CIRC->mag,m);
sprintf(read_string, "%s_circhw", name);
- input_float_vector(read_string, CIRC->numb, CIRC->halfw);
+ input_float_vector(read_string, CIRC->numb, CIRC->halfw,m);
sprintf(read_string, "%s_circovl", name);
- input_char_vector(read_string, CIRC->numb, CIRC->overlay);
+ input_char_vector(read_string, CIRC->numb, CIRC->overlay,m);
if(parse_and_apply)
field_arbitrary_circle(E, CIRC, field, BC, bcbitf, bcmask_on, bcmask_off);
@@ -705,44 +708,45 @@
{
char read_string[500];
int i;
+ int m = E->parallel.me;
sprintf(read_string, "%s_harm", name);
- input_int(read_string, &(HARM->numb), "0");
+ input_int(read_string, &(HARM->numb), "0",m);
sprintf(read_string, "%s_harms", name);
- input_int(read_string, &(HARM->harms), "0,0,19");
+ input_int(read_string, &(HARM->harms), "0,0,19",m);
sprintf(read_string, "%s_harmoff", name);
- input_float_vector(read_string, HARM->numb, HARM->off);
+ input_float_vector(read_string, HARM->numb, HARM->off,m);
sprintf(read_string, "%s_harmx1", name);
- input_float_vector(read_string, HARM->numb, HARM->x1);
+ input_float_vector(read_string, HARM->numb, HARM->x1,m);
sprintf(read_string, "%s_harmx2", name);
- input_float_vector(read_string, HARM->numb, HARM->x2);
+ input_float_vector(read_string, HARM->numb, HARM->x2,m);
sprintf(read_string, "%s_harmz1", name);
- input_float_vector(read_string, HARM->numb, HARM->z1);
+ input_float_vector(read_string, HARM->numb, HARM->z1,m);
sprintf(read_string, "%s_harmz2", name);
- input_float_vector(read_string, HARM->numb, HARM->z2);
+ input_float_vector(read_string, HARM->numb, HARM->z2,m);
sprintf(read_string, "%s_harmy1", name);
- input_float_vector(read_string, HARM->numb, HARM->z1);
+ input_float_vector(read_string, HARM->numb, HARM->z1,m);
sprintf(read_string, "%s_harmy2", name);
- input_float_vector(read_string, HARM->numb, HARM->z2);
+ input_float_vector(read_string, HARM->numb, HARM->z2,m);
sprintf(read_string, "%s_harmovl", name);
- input_char_vector(read_string, HARM->numb, HARM->overlay);
+ input_char_vector(read_string, HARM->numb, HARM->overlay,m);
for(i = 0; i < HARM->harms; i++)
{
sprintf(read_string, "%s_harmkx%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->kx[i]);
+ input_float_vector(read_string, HARM->numb, HARM->kx[i],m);
sprintf(read_string, "%s_harmkz%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->kz[i]);
+ input_float_vector(read_string, HARM->numb, HARM->kz[i],m);
sprintf(read_string, "%s_harmky%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->ky[i]);
+ input_float_vector(read_string, HARM->numb, HARM->ky[i],m);
sprintf(read_string, "%s_harmka%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->ka[i]);
+ input_float_vector(read_string, HARM->numb, HARM->ka[i],m);
sprintf(read_string, "%s_harmphx%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->phx[i]);
+ input_float_vector(read_string, HARM->numb, HARM->phx[i],m);
sprintf(read_string, "%s_harmphz%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->phz[i]);
+ input_float_vector(read_string, HARM->numb, HARM->phz[i],m);
sprintf(read_string, "%s_harmphy%02d", name, i + 1);
- input_float_vector(read_string, HARM->numb, HARM->phy[i]);
+ input_float_vector(read_string, HARM->numb, HARM->phy[i],m);
}
if(parse_and_apply)
@@ -837,3 +841,85 @@
return fi;
}
+/* safe file open function */
+FILE *safe_fopen(char *name,char *mode)
+{
+ FILE *tmp;
+ char m2[300];
+ if((tmp=(FILE *)fopen(name,mode))==NULL){
+ fprintf(stderr,"error: cannot fopen file %s, exiting\n",
+ name);
+ parallel_process_termination();
+ }
+ return ((FILE *)tmp);
+}
+/*
+ like malloc, but with test
+*/
+void *safe_malloc (size_t size)
+{
+ void *tmp;
+
+ if ((tmp = malloc(size)) == NULL) {
+ fprintf(stderr, "safe_malloc: could not allocate memory, n = %i, exiting\n",
+ (int)size);
+ parallel_process_termination();
+ }
+ return (tmp);
+}
+/* compute base vectors for conversion of polar to cartesian vectors
+ base[9]
+*/
+void calc_cbase_at_tp(float theta, float phi, float *base)
+{
+
+
+ double ct,cp,st,sp;
+
+ ct=cos(theta);
+ cp=cos(phi);
+ st=sin(theta);
+ sp=sin(phi);
+ /* r */
+ base[0]= st * cp;
+ base[1]= st * sp;
+ base[2]= ct;
+ /* theta */
+ base[3]= ct * cp;
+ base[4]= ct * sp;
+ base[5]= -st;
+ /* phi */
+ base[6]= -sp;
+ base[7]= cp;
+ base[8]= 0.0;
+}
+/* given a base from calc_cbase_at_tp, convert a polar vector to
+ cartesian */
+void convert_pvec_to_cvec(float vr,float vt, float vp, float *base,
+ float *cvec)
+{
+ int i;
+ for(i=0;i<3;i++){
+ cvec[i] = base[i] * vr;
+ cvec[i] += base[3+i]* vt;
+ cvec[i] += base[6+i]* vp;
+ }
+}
+/* convert r,theta,phi system to cartesian, xout[3] */
+void rtp2xyz(float r, float theta, float phi, float *xout)
+{
+ double rst;
+ rst = (double)r * sin((double)theta);
+ xout[0] = rst * cos((double)phi); /* x */
+ xout[1] = rst * sin((double)phi); /* y */
+ xout[2] = (double)r * cos((double)theta);
+}
+void myerror(char *message, struct All_variables *E)
+{
+ E->control.verbose = 1;
+ record(E,message);
+ fprintf(stderr,"node %3i: error: %s\n",
+ E->parallel.me,message);
+ parallel_process_termination();
+}
+
Modified: mc/3D/CitcomCU/trunk/src/Parsing.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Parsing.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Parsing.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -1,54 +1,63 @@
/*
- * CitcomCU is a Finite Element Code that solves for thermochemical
- * convection within a three dimensional domain appropriate for convection
- * within the Earth's mantle. Cartesian and regional-spherical geometries
- * are implemented. See the file README contained with this distribution
- * for further details.
- *
- * Copyright (C) 1994-2005 California Institute of Technology
- * Copyright (C) 2000-2005 The University of Colorado
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
- * Authors: Louis Moresi, Shijie Zhong, and Michael Gurnis
+ *<LicenseText>
*
- * For questions or comments regarding this software, you may contact
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
*
- * Luis Armendariz <luis at geodynamics.org>
- * http://geodynamics.org
- * Computational Infrastructure for Geodynamics (CIG)
- * California Institute of Technology
- * 2750 East Washington Blvd, Suite 210
- * Pasadena, CA 91007
+ * 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 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 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.
*
- * 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
*
- * 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>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
-
/* Routines which read filenames from the command line and
then parse the contents as parameters for citcom */
#include <stdio.h>
-#include <malloc.h>
#include <sys/types.h>
+#ifndef __sunos__
+#include <strings.h>
+#else
#include <string.h>
+#endif
#include "global_defs.h"
-#define MAXLINE 1024 /* max length of line in input file */
-#define MAXNAME 64 /* max length of name */
-#define MAXVALUE 1024 /* max length of value */
-#define MAXFILENAME 64 /* max length of par file name */
-#define MAXVECTOR 10 /* max # of elements for unspecified vectors */
+void setup_parser(struct All_variables *E, char *filename);
+void shutdown_parser(struct All_variables *E);
+int input_string(char *name, char *value, char *Default, int m);
+int input_boolean(char *name, int *value, char *interpret, int m);
+int input_int(char *name, int *value, char *interpret, int m);
+int input_float(char *name, float *value, char *interpret, int m);
+int input_double(char *name, double *value, char *interpret, int m);
+int input_int_vector(char *name, int number, int *value, int m);
+int input_char_vector(char *name, int number, char *value, int m);
+int input_float_vector(char *name,int number, float *value, int m);
+int input_double_vector(char *name, int number, double *value, int m);
+
+
+#define MAXLINE 1024 /* max length of line in input file */
+#define MAXNAME 64 /* max length of name */
+#define MAXVALUE 1024 /* max length of value */
+#define MAXFILENAME 64 /* max length of par file name */
+#define MAXVECTOR 10 /* max # of elements for unspecified vectors */
+#define STRANGE_NUM -98765.4321
+
/* abbreviations: */
#define AL struct arglist
#define PROGNAME ext_par.progname
@@ -62,785 +71,783 @@
#define BUFMAX ext_par.bufmax
#define LISTFILE ext_par.listout
-#define LISTINC 32 /* increment size for arglist */
-#define BUFINC 1024 /* increment size for argbuf */
+#define LISTINC 32 /* increment size for arglist */
+#define BUFINC 1024 /* increment size for argbuf */
-struct ext_par /* global variables for getpar */
+struct ext_par /* global variables for getpar */
{
- char *progname;
- int argflags;
- struct arglist *arglist;
- struct arglist *arghead;
- char *argbuf;
- int nlist;
- int nbuf;
- int listmax;
- int bufmax;
- FILE *listout;
-} ext_par;
+ char *progname;
+ int argflags;
+ struct arglist *arglist;
+ struct arglist *arghead;
+ char *argbuf;
+ int nlist;
+ int nbuf;
+ int listmax;
+ int bufmax;
+ FILE *listout;
+} ext_par;
-struct arglist /* structure of list set up by setpar */
+struct arglist /* structure of list set up by setpar */
{
- int argname_offset;
- int argval_offset;
- int hash;
+ int argname_offset;
+ int argval_offset;
+ int hash;
};
int VERBOSE = 0;
int DESCRIBE = 0;
int BEGINNER = 0;
+int interpret_control_string();
-void setup_parser(struct All_variables *E, int ac, char **av)
+
+void setup_parser(E,filename)
+ struct All_variables *E;
+ char *filename;
{
- FILE *fp;
- char *pl, *pn, *pv;
- char t1, t2, line[MAXLINE], name[MAXNAME], value[MAXVALUE];
- int i, j, k;
+ void unique_copy_file();
+ void add_to_parameter_list();
- /* should get file length & cpp &c before any further parsing */
+ FILE * fp;
+ char *pl,*pn,*pv;
+ char t1, t2, line[MAXLINE], name[MAXNAME], value[MAXVALUE];
+ int i,j,k;
+ int m=E->parallel.me;
- /* for now, read one filename from the command line, we'll parse that ! */
+ /* should get file length & cpp &c before any further parsing */
- if(ac < 2)
- {
- fprintf(stderr, "Usage: citcom PARAMETERFILE\n");
- exit(10);
- }
+ /* for now, read one filename from the command line, we'll parse that ! */
- if((fp = fopen(av[1], "r")) == NULL)
- {
- fprintf(stderr, "File: %s is unreadable\n", av[1]);
- exit(11);
- }
+ /* this section moved to main() */
+/* if (ac < 2) { */
+/* fprintf(stderr,"Usage: citcom PARAMETERFILE\n"); */
+/* exit(10); */
+/* } */
- unique_copy_file(E, av[1], "copy");
+ if ((fp = fopen(filename,"r")) == NULL) {
+ fprintf(stderr,"(Parsing #1) File: %s is unreadable\n",filename);
+ exit(11);
+ }
- /* now the parameter file is open, read into memory */
- while(fgets(line, MAXLINE, fp) != NULL)
- {
- pl = line;
- /* loop over entries on each line */
- loop:
- while(*pl == ' ' || *pl == '\t')
- pl++;
- if(*pl == '\0' || *pl == '\n')
- continue; /* end of line */
- if(*pl == '#')
- continue; /* end of interpretable part of line */
- /* get name */
- pn = name;
- while(*pl != '=' && *pl != '\0' && *pl != ' ' && *pl != '\n' /* FIX by Glenn Nelson */
- && *pl != '\t')
- *pn++ = *pl++;
- *pn = '\0';
- if(*pl == '=')
- pl++;
+ /* now the parameter file is open, read into memory */
- /* get value */
- *value = '\0';
- pv = value;
- if(*pl == '"' || *pl == '\'')
- t1 = t2 = *pl++;
- else
- {
- t1 = ' ';
- t2 = '\t';
- }
- while(*pl != t1 && *pl != t2 && *pl != '\0' && *pl != '\n')
- *pv++ = *pl++;
- *pv = '\0';
- if(*pl == '"' || *pl == '\'')
- pl++;
- add_to_parameter_list(name, value);
+ while( fgets(line,MAXLINE,fp) != NULL )
+ { pl= line;
+ /* loop over entries on each line */
+ loop:
+ while(*pl==' ' || *pl=='\t') pl++;
+ if(*pl=='\0'|| *pl=='\n') continue; /* end of line */
+ if(*pl=='#') continue; /* end of interpretable part of line */
- goto loop;
+ /* get name */
+ pn= name;
+ while(*pl != '=' && *pl != '\0' && *pl != ' '
+ && *pl != '\n' /* FIX by Glenn Nelson */
+ && *pl != '\t')
+ *pn++ = *pl++;
+ *pn = '\0';
+ if(*pl == '=') pl++;
+
+ /* get value */
+ *value= '\0';
+ pv= value;
+ if(*pl=='"' || *pl=='\'')
+ t1= t2= *pl++;
+ else
+ { t1= ' ';
+ t2= '\t';
}
+ while(*pl!=t1 && *pl!=t2 &&
+ *pl!='\0' && *pl!='\n') *pv++= *pl++;
+ *pv= '\0';
+ if(*pl=='"' || *pl=='\'')
+ pl++;
+ add_to_parameter_list(name,value);
- fclose(fp);
+ goto loop;
+ }
- ARGHEAD = ARGLIST;
+ fclose(fp);
- /* Now we can use our routines to check & set their own flags ! */
+ ARGHEAD= ARGLIST;
- printf("%s %d\n", av[1], E->parallel.me);
+ /* Now we can use our routines to check & set their own flags ! */
- input_boolean("VERBOSE", &i, "off");
- input_boolean("DESCRIBE", &j, "off");
- input_boolean("BEGINNER", &k, "off");
- VERBOSE = i;
- DESCRIBE = j;
- BEGINNER = k;
+ input_boolean("VERBOSE",&i,"off",m);
+ input_boolean("DESCRIBE",&j,"off",m);
+ input_boolean("BEGINNER",&k,"off",m);
+ VERBOSE=i;
+ DESCRIBE=j;
+ BEGINNER=k;
- return;
+ /* make this an optional behavior */
+ input_boolean("copy_input_file",&k,"on",m);
+ if(k)
+ unique_copy_file(E,filename,"copy");
+
+
}
-void shutdown_parser(struct All_variables *E)
+void shutdown_parser(E)
+ struct All_variables *E;
+
{
- if(ARGLIST != NULL)
- free(ARGLIST);
- if(ARGBUF != NULL)
- free(ARGBUF);
- ARGBUF = NULL;
- ARGLIST = NULL;
+ if(ARGLIST != NULL) free(ARGLIST);
+ if(ARGBUF != NULL) free(ARGBUF);
+ ARGBUF= NULL;
+ ARGLIST= NULL;
+
}
-/* add an entry to arglist, expanding memory if necessary */
-/* FIXME: return type should be void */
-int add_to_parameter_list(register char *name, register char *value)
+
+/* add an entry to arglist, expanding memory */
+void add_to_parameter_list(name,value)
+ char *name, *value; /* if necessary */
{
- struct arglist *alptr;
- int len;
- register char *ptr;
+ struct arglist *alptr;
+ int len;
+ char *ptr;
+ int compute_parameter_hash_table();
- /* check arglist memory */
- if(NLIST >= LISTMAX)
- {
- LISTMAX += LISTINC;
- if(ARGLIST == NULL)
- ARGLIST = (AL *) malloc(LISTMAX * sizeof(AL));
- else
- ARGLIST = (AL *) realloc(ARGLIST, LISTMAX * sizeof(AL));
- }
- /* check argbuf memory */
- len = strlen(name) + strlen(value) + 2; /* +2 for terminating nulls */
- if(NBUF + len >= BUFMAX)
- {
- BUFMAX += BUFINC;
- if(ARGBUF == NULL)
- ARGBUF = (char *)malloc(BUFMAX);
- else
- ARGBUF = (char *)realloc(ARGBUF, BUFMAX);
- }
- if(ARGBUF == NULL || ARGLIST == NULL)
- fprintf(stderr, "cannot allocate memory\n");
+ /* check arglist memory */
+ if(NLIST >= LISTMAX)
+ { LISTMAX += LISTINC;
+ if(ARGLIST == NULL)
+ ARGLIST= (AL *)malloc(LISTMAX * sizeof(AL));
+ else
+ ARGLIST= (AL *)realloc(ARGLIST,LISTMAX * sizeof(AL));
+ }
+ /* check argbuf memory */
+ len= strlen(name) + strlen(value) + 2; /* +2 for terminating nulls */
+ if(NBUF+len >= BUFMAX)
+ { BUFMAX += BUFINC;
+ if(ARGBUF == NULL)
+ ARGBUF= (char *)malloc(BUFMAX);
+ else ARGBUF= (char *)realloc(ARGBUF,BUFMAX);
+ }
+ if(ARGBUF == NULL || ARGLIST == NULL)
+ fprintf(stderr,"cannot allocate memory\n");
- /* add name */
- alptr = ARGLIST + NLIST;
- alptr->hash = compute_parameter_hash_table(name);
- alptr->argname_offset = NBUF;
- ptr = ARGBUF + NBUF;
- do
- *ptr++ = *name;
- while(*name++);
+ /* add name */
+ alptr= ARGLIST + NLIST;
+ alptr->hash= compute_parameter_hash_table(name);
+ alptr->argname_offset = NBUF;
+ ptr= ARGBUF + NBUF;
+ do
+ *ptr++ = *name;
+ while(*name++);
- /* add value */
- NBUF += len;
- alptr->argval_offset = ptr - ARGBUF;
- do
- *ptr++ = *value;
- while(*value++);
- NLIST++;
-
- return 0; /* successfully added value to list */
+ /* add value */
+ NBUF += len;
+ alptr->argval_offset= ptr - ARGBUF;
+ do
+ *ptr++ = *value;
+ while(*value++);
+ NLIST++;
}
-int compute_parameter_hash_table(register char *s)
-{
- register int h;
+int compute_parameter_hash_table(s)
+ char *s;
+{ int h;
- h = s[0];
- if(s[1])
- h |= (s[1]) << 8;
- else
- return (h);
- if(s[2])
- h |= (s[2]) << 16;
- else
- return (h);
- if(s[3])
- h |= (s[3]) << 24;
- return (h);
+ h= s[0];
+ if(s[1])
+ h |= (s[1])<<8;
+ else
+ return(h);
+ if(s[2])
+ h |= (s[2])<<16;
+ else
+ return(h);
+ if(s[3])
+ h |= (s[3])<<24;
+ return(h);
}
-int input_int(char *name, int *value, char *interpret)
+int input_int(name,value,interpret,m)
+ char *name;
+ int *value;
+ char *interpret;
+ int m;
+
{
- struct arglist *alptr;
- int h, found;
- char *str;
+ struct arglist *alptr;
+ int h, found;
+ char *str;
- int exists, essential;
- double *Default, *minvalue, *maxvalue;
+ int exists,essential;
+ double Default,minvalue,maxvalue;
- if(DESCRIBE)
- fprintf(stderr, "input_int: searching for '%s' with default/range '%s'\n", name, (interpret == NULL) ? "**EMPTY**" : interpret);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_int: searching for '%s' with default/range '%s'\n",
+ name,(interpret == NULL) ? "**EMPTY**" : interpret);
- exists = interpret_control_string(interpret, &essential, &Default, &minvalue, &maxvalue);
+ exists = interpret_control_string(interpret,&essential,&Default,&minvalue,&maxvalue);
- if(Default != NULL)
- *value = (int)(*Default);
+ *value = (int)(Default);
- h = compute_parameter_hash_table(name);
- found = 0;
+ h=compute_parameter_hash_table(name);
+ found=0;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
- str = ARGBUF + alptr->argval_offset;
- sscanf(str, "%d", value);
- found = 1;
- break;
- }
+ str= ARGBUF + alptr->argval_offset;
+ sscanf(str,"%d",value);
+ found=1;
+ break;
+ }
- if(essential && !found)
- {
- fprintf(stderr, "There MUST be an entry for the parameter %s\n", name);
- exit(12);
- }
- if((minvalue != NULL) && (*value < (int)*minvalue))
- {
- *value = (int)*minvalue;
- }
- if((maxvalue != NULL) && (*value > (int)*maxvalue))
- {
- *value = (int)*maxvalue;
- }
+ if(essential && !found)
+ { fprintf(stderr,"There MUST be an entry for the parameter %s\n",name);
+ exit(12);
+ }
+ if((minvalue!=STRANGE_NUM) && (*value < (int) minvalue))
+ { *value = (int) minvalue;
+ }
+ if((maxvalue!=STRANGE_NUM) && (*value > (int) maxvalue))
+ { *value = (int) maxvalue;
+ }
- if(VERBOSE)
- {
- if(found)
- fprintf(stderr, "%25s: (int) = %d \n", name, *value);
- else if(Default != NULL)
- fprintf(stderr, "%25s: (int) = not found (%d) \n", name, (int)(*Default));
- else
- {
- fprintf(stderr, "%25s: (int) = not found (no default) \n", name);
- if(BEGINNER)
- {
- fprintf(stderr, "\t\t Previously set value gives ...");
- fprintf(stderr, "%d\n", *value);
- }
- }
- }
+ if(m==0)
+ if(VERBOSE)
+ { if (found)
+ fprintf(stderr,"%25s: (int) = %d \n",name,*value);
+ else
+ if (Default != STRANGE_NUM)
+ fprintf(stderr,"%25s: (int) = not found (%d) \n",name,(int)(Default));
+ else
+ { fprintf(stderr,"%25s: (int) = not found (no default) \n",name);
+ if(BEGINNER)
+ { fprintf(stderr,"\t\t Previously set value gives ...");
+ fprintf(stderr,"%d\n",*value);
+ }
+ }
+ }
- return found;
+ return(found);
}
-/* in the case of a string default=NULL forces input */
-int input_string(char *name, char *value, char *Default)
+int input_string(name,value,Default,m) /* in the case of a string default=NULL forces input */
+ char *name;
+ char *value;
+ char *Default;
+ int m;
{
- //char *sptr;
- struct arglist *alptr;
- //int h, hno, hyes, found;
- int h, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
- int essential;
+ char *sptr;
+ struct arglist *alptr;
+ int h, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
+ int essential;
- if(DESCRIBE)
- fprintf(stderr, "input_string: searching for '%s' with default '%s'\n", name, (Default == NULL) ? "no default" : Default);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_string: searching for '%s' with default '%s'\n",
+ name,(Default == NULL) ? "no default" : Default);
- h = compute_parameter_hash_table(name);
- essential = found = 0;
+ h=compute_parameter_hash_table(name);
+ essential=found=0;
- if(Default != NULL) /* Cannot use "Essential" as this is a valid input */
- strcpy(value, Default);
- else
- essential = 1;
+ if (Default != NULL) /* Cannot use "Essential" as this is a valid input */
+ strcpy(value,Default);
+ else
+ essential=1;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
- str = ARGBUF + alptr->argval_offset;
- strcpy(value, str);
- found = 1;
- break;
- }
+ str= ARGBUF + alptr->argval_offset;
+ strcpy(value,str);
+ found=1;
+ break;
+ }
- if(essential && !found)
- {
- fprintf(stderr, "There MUST be an entry for the parameter %s\n", name);
- exit(12);
- }
+ if(essential && !found)
+ { fprintf(stderr,"There MUST be an entry for the parameter %s\n",name);
+ exit(12);
+ }
- if(VERBOSE)
- fprintf(stderr, "%25s: (string) = %s (%s)\n", name, (found ? value : "not found"), (Default != NULL ? Default : "no default"));
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (string) = %s (%s)\n",name,
+ (found ? value : "not found"),
+ (Default != NULL ? Default : "no default"));
- return found;
+ return(found);
}
-/* supports name=on/off too */
-int input_boolean (char *name, int *value, char *interpret)
+int input_boolean(name,value,interpret,m) /* supports name=on/off too */
+ char *name;
+ int *value;
+ char *interpret;
+ int m;
{
- //char *sptr;
- struct arglist *alptr;
- //int h, hno, hyes, found;
- int h, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+ char *sptr;
+ struct arglist *alptr;
+ int h, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
- int essential;
- double *Default, *minvalue, *maxvalue;
+ int essential;
+ double Default,minvalue,maxvalue;
- if(DESCRIBE)
- fprintf(stderr, "input_boolean: searching for '%s' with default/range '%s'\n", name, (interpret == NULL) ? "**EMPTY**" : interpret);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_boolean: searching for '%s' with default/range '%s'\n",
+ name,(interpret == NULL) ? "**EMPTY**" : interpret);
- interpret_control_string(interpret, &essential, &Default, &minvalue, &maxvalue);
- if(Default != NULL)
- *value = (int)(*Default);
+ interpret_control_string(interpret,&essential,&Default,&minvalue,&maxvalue);
- h = compute_parameter_hash_table(name);
- found = 0;
+ *value = (int)(Default);
+ h=compute_parameter_hash_table(name);
+ found=0;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
- str = ARGBUF + alptr->argval_offset;
- found = 1;
- break;
- }
+ str= ARGBUF + alptr->argval_offset;
+ found=1;
+ break;
+ }
+ if(!found)
+ {if(m==0)
+ if(VERBOSE)
+ if (Default != STRANGE_NUM)
+ fprintf(stderr,"%25s: (boolean int) = not found (%d) \n",name,(int)(Default));
+ else
+ { fprintf(stderr,"%25s: (boolean int) = not found (no default) \n",name);
+ if(BEGINNER)
+ { fprintf(stderr,"\t\t Previously set value gives ...");
+ fprintf(stderr,"%d\n",*value);
+ }
+ }
- if(!found)
- {
- if(VERBOSE)
- {
- if(Default != NULL)
- fprintf(stderr, "%25s: (boolean int) = not found (%d) \n", name, (int)(*Default));
- else
- {
- fprintf(stderr, "%25s: (boolean int) = not found (no default) \n", name);
- if(BEGINNER)
- {
- fprintf(stderr, "\t\t Previously set value gives ...");
- fprintf(stderr, "%d\n", *value);
- }
- }
- }
- return 0;
- }
+ return(0);
+ }
+ if((strstr(str,"on")!=NULL) || (strstr(str,"ON")!=NULL))
+ *value=1;
+ else if ((strstr(str,"off") != NULL) || (strstr(str,"OFF")!=NULL))
+ *value=0;
+ else /* assume some numerical value */
+ *value=atoi(str);
- if((strstr(str, "on") != NULL) || (strstr(str, "ON") != NULL))
- *value = 1;
- else if((strstr(str, "off") != NULL) || (strstr(str, "OFF") != NULL))
- *value = 0;
- else
- *value = atoi(str);
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (boolean int) = %d \n",name,*value);
- if(VERBOSE)
- fprintf(stderr, "%25s: (boolean int) = %d \n", name, *value);
-
- return (found);
+ return(found);
}
-int input_float(char *name, float *value, char *interpret)
-{
- //char *sptr;
- struct arglist *alptr;
+int input_float(name,value,interpret,m)
+ char *name;
+ float *value;
+ char *interpret;
+ int m;
- //int h, hno, hyes, found;
- int h, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
- int exists, essential;
- double *Default, *minvalue, *maxvalue;
+{ char *sptr;
+ struct arglist *alptr;
+ int h, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
+ int exists,essential;
+ double Default,minvalue,maxvalue;
- if(DESCRIBE)
- fprintf(stderr, "input_float: searching for '%s' with default/range '%s'\n", name, (interpret == NULL) ? "**EMPTY**" : interpret);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_float: searching for '%s' with default/range '%s'\n",
+ name,(interpret == NULL) ? "**EMPTY**" : interpret);
- exists = interpret_control_string(interpret, &essential, &Default, &minvalue, &maxvalue);
+ exists=interpret_control_string(interpret,&essential,&Default,&minvalue,&maxvalue);
- if(Default != NULL)
- *value = (float)*Default;
+ *value = (float) Default;
- h = compute_parameter_hash_table(name);
- found = 0;
+ h=compute_parameter_hash_table(name);
+ found=0;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
- sscanf(str, "%f", value);
- found = 1;
- break;
- }
+ sscanf(str,"%f",value);
+ found=1;
+ break;
+ }
- if(essential && !found)
- {
- fprintf(stderr, "There MUST be an entry for the parameter %s\n", name);
- exit(12);
- }
+ if(essential && !found)
+ { fprintf(stderr,"There MUST be an entry for the parameter %s\n",name);
+ exit(12);
+ }
- if((minvalue != NULL) && (*value < (float)*minvalue))
- *value = (float)*minvalue;
- if((maxvalue != NULL) && (*value > (float)*maxvalue))
- *value = (float)*maxvalue;
+ if((minvalue!=STRANGE_NUM) && (*value < (float) minvalue))
+ *value = (float) minvalue;
+ if((maxvalue!=STRANGE_NUM) && (*value > (float) maxvalue))
+ *value = (float) maxvalue;
- if(VERBOSE)
- {
- if(found)
- fprintf(stderr, "%25s: (float) = %f \n", name, *value);
- else if(Default != NULL)
- fprintf(stderr, "%25s: (float) = not found (%f) \n", name, *Default);
- else
- {
- fprintf(stderr, "%25s: (float) = not found (no default) \n", name);
- if(BEGINNER)
- {
- fprintf(stderr, "\t\t Previously set value gives ...");
- fprintf(stderr, "%g\n", *value);
- }
- }
- }
- return (found);
+ if(m==0)
+ if(VERBOSE)
+ { if (found)
+ fprintf(stderr,"%25s: (float) = %g \n",name,*value);
+ else
+ if (Default != STRANGE_NUM)
+ fprintf(stderr,"%25s: (float) = not found (%g) \n",name,Default);
+ else
+ { fprintf(stderr,"%25s: (float) = not found (no default) \n",name);
+ if(BEGINNER)
+ { fprintf(stderr,"\t\t Previously set value gives ...");
+ fprintf(stderr,"%g\n",*value);
+ }
+ }
+ }
+ return(found);
}
-int input_double(char *name, double *value, char *interpret)
-{
- //char *sptr;
- struct arglist *alptr;
+int input_double(name,value,interpret,m)
+ char *name;
+ double *value;
+ char *interpret;
+ int m;
- //int h, hno, hyes, found;
- int h, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+{ char *sptr;
+ struct arglist *alptr;
- int exists, essential;
- double *Default, *minvalue, *maxvalue;
+ int h, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
+ int exists,essential;
+ double Default,minvalue,maxvalue;
- if(DESCRIBE)
- fprintf(stderr, "input_double: searching for '%s' with default/range '%s'\n", name, (interpret == NULL) ? "**EMPTY**" : interpret);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_double: searching for '%s' with default/range '%s'\n",
+ name,(interpret == NULL) ? "**EMPTY**" : interpret);
- exists = interpret_control_string(interpret, &essential, &Default, &minvalue, &maxvalue);
- if(Default != NULL)
- *value = *Default;
+ exists=interpret_control_string(interpret,&essential,&Default,&minvalue,&maxvalue);
- h = compute_parameter_hash_table(name);
- found = 0;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
- sscanf(str, "%lf", value);
- found = 1;
- break;
- }
+ *value = Default;
- if(essential && !found)
- {
- fprintf(stderr, "There MUST be an entry for the parameter %s\n", name);
- exit(12);
- }
- if((minvalue != NULL) && (*value < *minvalue))
- *value = *minvalue;
- if((maxvalue != NULL) && (*value > *maxvalue))
- *value = *maxvalue;
+ h=compute_parameter_hash_table(name);
+ found=0;
- if(VERBOSE)
- {
- if(found)
- fprintf(stderr, "%25s: (double) = %g \n", name, *value);
- else if(Default != NULL)
- fprintf(stderr, "%25s: (double) = not found (%g) \n", name, *Default);
- else
- {
- fprintf(stderr, "%25s: (double) = not found (no default)\n", name);
- if(BEGINNER)
- {
- fprintf(stderr, "\t\t Previously set value gives ...");
- fprintf(stderr, "%g\n", *value);
- }
- }
- }
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
+ sscanf(str,"%lf",value);
+ found=1;
+ break;
+ }
+ if(essential && !found)
+ { fprintf(stderr,"There MUST be an entry for the parameter %s\n",name);
+ exit(12);
+ }
+ if((minvalue!=STRANGE_NUM) && (*value < minvalue))
+ *value = minvalue;
+ if((maxvalue!=STRANGE_NUM) && (*value > maxvalue))
+ *value = maxvalue;
- return (found);
+ if(m==0)
+ if(VERBOSE)
+ { if (found)
+ fprintf(stderr,"%25s: (double) = %g \n",name,*value);
+ else
+ if (Default != STRANGE_NUM)
+ fprintf(stderr,"%25s: (double) = not found (%g) \n",name,Default);
+ else
+ { fprintf(stderr,"%25s: (double) = not found (no default) \n",name);
+ if(BEGINNER)
+ { fprintf(stderr,"\t\t Previously set value gives ...");
+ fprintf(stderr,"%g\n",*value);
+ }
+ }
+ }
+
+
+ return(found);
}
-/* value is a comma-separated list of ints */
-int input_int_vector(char *name, int number, int *value)
+int input_int_vector(char *name, int number,int *value,int m)
{
- //char *sptr;
- struct arglist *alptr;
- char control_string[500];
+ char *sptr;
+ struct arglist *alptr;
+ char control_string[500];
- //int h, i, hno, hyes, found;
- int h, i, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+ int h,i, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
- if(DESCRIBE)
- fprintf(stderr, "input_int_vector: searching for %s (%d times)\n", name, number);
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_int_vector: searching for %s (%d times)\n",name,number);
- h = compute_parameter_hash_table(name);
- found = 0;
+ h=compute_parameter_hash_table(name);
+ found=0;
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
- found = 1;
- break;
- }
- /* now interpret vector */
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
+ found=1;
+ break;
+ }
+ /* now interpret vector */
- if(!found)
- return 0;
+ if(!found) return(0);
- for(h = 0; h < number; h++)
- {
- sprintf(control_string, "");
- for(i = 0; i < h; i++)
- strcat(control_string, "%*f,");
- strcat(control_string, "%d");
- sscanf(str, control_string, &(value[h]));
- }
+ for(h=0;h<number;h++)
+ { sprintf(control_string,"");
+ for(i=0;i<h;i++)
+ strcat(control_string,"%*f,");
+ strcat(control_string,"%d");
+ sscanf(str,control_string,&(value[h]));
+ }
- if(VERBOSE)
- fprintf(stderr, "%25s: (vector) = %s\n", name, str);
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (vector) = %s\n",name,str);
- return (found);
+ return(found);
}
-/* value is a comma-separated list of ints */
-int input_char_vector(char *name, int number, char *value)
-{
- //char *sptr;
- struct arglist *alptr;
- char control_string[500];
+int input_char_vector(name,number,value,m)
+ char *name;
+ int number;
+ char *value; /* comma-separated list of ints */
+ int m;
- //int h, i, hno, hyes, found;
- int h, i, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+{ char *sptr;
+ struct arglist *alptr;
+ char control_string[500];
- if(DESCRIBE)
- fprintf(stderr, "input_char_vector: searching for %s (%d times)\n", name, number);
+ int h,i, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
- h = compute_parameter_hash_table(name);
- found = 0;
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_char_vector: searching for %s (%d times)\n",name,number);
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
- found = 1;
- break;
- }
- /* now interpret vector */
+ h=compute_parameter_hash_table(name);
+ found=0;
- if(!found)
- return (0);
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
+ found=1;
+ break;
+ }
+ /* now interpret vector */
- for(h = 0; h < number; h++)
- {
- sprintf(control_string, "");
- for(i = 0; i < h; i++)
- strcat(control_string, "%*c,");
- strcat(control_string, "%c");
- sscanf(str, control_string, &(value[h]));
- }
+ if(!found) return(0);
- if(VERBOSE)
- fprintf(stderr, "%25s: (vector) = %s\n", name, str);
+ for(h=0;h<number;h++)
+ { sprintf(control_string,"");
+ for(i=0;i<h;i++)
+ strcat(control_string,"%*c,");
+ strcat(control_string,"%c");
+ sscanf(str,control_string,&(value[h]));
+ }
- return (found);
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (vector) = %s\n",name,str);
+
+ return(found);
}
-/* value is a comma-separated list of floats */
-int input_float_vector(char *name, int number, float *value)
-{
- //char *sptr;
- struct arglist *alptr;
- char control_string[500];
+int input_float_vector(name,number,value,m)
+ char *name;
+ int number;
+ float *value; /* comma-separated list of floats */
+ int m;
- //int h, i, hno, hyes, found;
- int h, i, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+{ char *sptr;
+ struct arglist *alptr;
+ char control_string[500];
- if(0 == number)
- return (0);
+ int h,i, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
- if(DESCRIBE)
- fprintf(stderr, "input_float_vector: searching for %s (%d times)\n", name, number);
+ if(0==number)
+ return(0);
- h = compute_parameter_hash_table(name);
- found = 0;
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_float_vector: searching for %s (%d times)\n",name,number);
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
- found = 1;
- break;
- }
- /* now interpret vector */
+ h=compute_parameter_hash_table(name);
+ found=0;
- if(!found)
- return (0);
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
+ found=1;
+ break;
+ }
+ /* now interpret vector */
- for(h = 0; h < number; h++)
- {
- sprintf(control_string, "");
- for(i = 0; i < h; i++)
- strcat(control_string, "%*f,");
- strcat(control_string, "%f");
- sscanf(str, control_string, &(value[h]));
- }
+ if(!found) return(0);
- if(VERBOSE)
- fprintf(stderr, "%25s: (float vector) = %s\n", name, str);
+ for(h=0;h<number;h++)
+ { sprintf(control_string,"");
+ for(i=0;i<h;i++)
+ strcat(control_string,"%*f,");
+ strcat(control_string,"%f");
+ sscanf(str,control_string,&(value[h]));
+ }
- return found;
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (float vector) = %s\n",name,str);
+
+ return(found);
}
-/* value is a comma-separated list of doubles */
-int input_double_vector(char *name, int number, double *value)
-{
- //char *sptr;
- struct arglist *alptr;
- char control_string[500];
+int input_double_vector(name,number,value,m)
+ char *name;
+ int number;
+ double *value; /* comma-separated list of floats */
+ int m;
- //int h, i, hno, hyes, found;
- int h, i, found;
- //char line[MAXLINE], *str, *noname;
- char *str;
+{ char *sptr;
+ struct arglist *alptr;
+ char control_string[500];
- if(DESCRIBE)
- fprintf(stderr, "input_double_vector: searching for %s (%d times)\n", name, number);
+ int h,i, hno, hyes, found;
+ char line[MAXLINE], *str, *noname;
- h = compute_parameter_hash_table(name);
- found = 0;
+ if(m==0)
+ if(DESCRIBE)
+ fprintf(stderr,"input_double_vector: searching for %s (%d times)\n",name,number);
- /* search list backwards, stopping at first find */
- for(alptr = ARGLIST + (NLIST - 1); alptr >= ARGHEAD; alptr--)
- {
- if(alptr->hash != h)
- continue;
- if(strcmp(ARGBUF + alptr->argname_offset, name) != 0)
- continue;
- str = ARGBUF + alptr->argval_offset;
- found = 1;
- break;
- }
+ h=compute_parameter_hash_table(name);
+ found=0;
- if(!found)
- return (0);
+ /* search list backwards, stopping at first find */
+ for(alptr= ARGLIST +(NLIST-1); alptr >= ARGHEAD; alptr--)
+ { if(alptr->hash != h)
+ continue;
+ if(strcmp(ARGBUF+alptr->argname_offset,name) != 0)
+ continue;
+ str= ARGBUF + alptr->argval_offset;
+ found=1;
+ break;
+ }
- /* now interpret vector */
+ if(!found) return(0);
- for(h = 0; h < number; h++)
- {
- sprintf(control_string, "");
- for(i = 0; i < h; i++)
- strcat(control_string, "%*f,");
- strcat(control_string, "%lf");
- sscanf(str, control_string, &(value[h]));
- }
+ /* now interpret vector */
- if(VERBOSE)
- fprintf(stderr, "%25s: (double vector) = %s\n", name, str);
+ for(h=0;h<number;h++)
+ { sprintf(control_string,"");
+ for(i=0;i<h;i++)
+ strcat(control_string,"%*f,");
+ strcat(control_string,"%lf");
+ sscanf(str,control_string,&(value[h]));
+ }
- return found;
+ if(m==0)
+ if(VERBOSE)
+ fprintf(stderr,"%25s: (double vector) = %s\n",name,str);
+
+ return(found);
}
/* =================================================== */
+/* This is needed to be fixed on Linux machine
+ The function strtok does not work on linux machine
+*/
-int interpret_control_string(char *interpret, int *essential, double **Default, double **minvalue, double **maxvalue)
-{
- char *substring;
+int interpret_control_string(interpret,essential,Default,minvalue,maxvalue)
+ char *interpret;
+ int *essential;
+ double *Default,*minvalue,*maxvalue;
- *Default = *maxvalue = *minvalue = NULL;
- *essential = 0;
+{ char *substring;
- return (0); /* nothing to interpret */
+ *Default=*maxvalue=*minvalue=STRANGE_NUM;
+ *essential=0;
- if((substring = (char *)strtok(interpret, ",")) == NULL)
- return (0); /* nothing to interpret */
+ if (strstr(interpret,"essential")!=NULL)
+ { *essential=1; /* no default possible, must read a value */
+ return(0);
+ }
+ if (strstr(interpret,"nodefault")==NULL)
+ { if((strstr(interpret,"on")!=NULL) || (strstr(interpret,"ON")!=NULL))
+ *Default = 1.0;
+ else
+ if ((strstr(interpret,"off") != NULL) || (strstr(interpret,"OFF")!=NULL))
+ *Default = 0.0;
+ else
+ sscanf(interpret,"%lf",Default); /* read number as a default value */
+ }
- if(strstr(substring, "essential") != NULL)
- *essential = 1; /* no default possible, must read a value */
- else if(strstr(substring, "nodefault") == NULL)
- {
- *Default = (double *)malloc(sizeof(double));
- if((strstr(substring, "on") != NULL) || (strstr(substring, "ON") != NULL))
- **Default = 1.0;
- else if((strstr(substring, "off") != NULL) || (strstr(substring, "OFF") != NULL))
- **Default = 0.0;
- else
- sscanf(substring, "%lf", *Default); /* read number as a default value */
+ if ((substring=strstr(interpret,",")) == NULL) /* minvalue */
+ { /* no minimum, no maximum */
+ return(1);
+ }
- }
+ if (strstr(substring,"nomin")==NULL)
+ sscanf(substring,"%lf",minvalue);
- if((substring = (char *)strtok(NULL, ",")) == NULL) /* minvalue */
- { /* no minimum, no maximum */
- return (1);
- }
+ if ((substring=strstr(substring,",")) == NULL) /* maxvalue */
+ { /* no maximum */
+/* if (DESCRIBE) */
+/* fprintf(stderr,"minimum but no maximum\n"); */
+ return(2);
+ }
- if(strstr(substring, "nomin") == NULL)
- {
- *minvalue = (double *)malloc(sizeof(double));
- sscanf(substring, "%lf", *minvalue);
- }
+ if (strstr(substring,"nomax")==NULL)
+ sscanf(substring,"%lf",maxvalue);
- if((substring = (char *)strtok(NULL, ",")) == NULL) /* maxvalue */
- { /* no maximum */
- if(DESCRIBE)
- fprintf(stderr, "minimum but no maximum\n");
- return (1);
- }
- if(strstr(substring, "nomax") == NULL)
- {
- *maxvalue = (double *)malloc(sizeof(double));
- sscanf(substring, "%lf", *maxvalue);
- }
+ return(0);
- return 1;
}
Modified: mc/3D/CitcomCU/trunk/src/Process_velocity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Process_velocity.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Process_velocity.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -50,28 +50,26 @@
void process_new_velocity(struct All_variables *E, int ii)
{
- static int been_here = 0;
- if(been_here == 0)
- {
- E->monitor.time_scale = pow(E->monitor.length_scale, 2.0) / /* Million years */
- (E->data.therm_diff * 3600.0 * 24.0 * 365.25 * 1.0e6);
- been_here++;
- }
+ if(E->control.stokes || ((ii % E->control.record_every) == 0))
+ {
+ /* get_CBF_topo(E,E->slice.tpg,E->slice.tpgb); */
+
+ get_STD_topo(E, E->slice.tpg, E->slice.tpgb, ii);
+ averages(E);
- if(((ii % E->control.record_every) == 0))
- {
- /* get_CBF_topo(E,E->slice.tpg,E->slice.tpgb); */
+#ifdef USE_GZDIR
+ if(E->control.gzdir)
+ output_velo_related_gzdir(E, ii); /* also topo */
+ else
+ output_velo_related(E, ii); /* also topo */
+#else
+ output_velo_related(E, ii); /* also topo */
+#endif
+ }
- get_STD_topo(E, E->slice.tpg, E->slice.tpgb, ii);
-
- averages(E);
-
- output_velo_related(E, ii); /* also topo */
- }
-
- return;
+ return;
}
/* =============================================== */
Modified: mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -184,13 +184,16 @@
if(E->control.print_convergence && E->parallel.me == 0)
{
- fprintf(E->fp, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n", count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
- /**/ fprintf(stderr, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n", count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
+ fprintf(E->fp, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n",
+ count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
+ /**/ fprintf(stderr, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n",
+ count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
/**/
}
/* while( (count < *steps_max) && (E->monitor.incompressibility >= E->control.tole_comp || dvelocity >= imp) ) {
-*/ while((count < *steps_max) && (dpressure >= imp || dvelocity >= imp))
+*/
+ while((count < *steps_max) && (dpressure >= imp || dvelocity >= imp))
{
for(j = 1; j <= npno; j++)
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2008-11-07 03:29:17 UTC (rev 13268)
@@ -46,17 +46,23 @@
#include "element_definitions.h"
#include "global_defs.h"
+void visc_from_B(struct All_variables *, float *, float *, int );
+void visc_from_C(struct All_variables *, float *, float *, int );
+void *safe_malloc (size_t );
+
void viscosity_parameters(struct All_variables *E)
{
int i, l;
float temp;
-
+ int m = E->parallel.me ;
/* default values .... */
E->viscosity.update_allowed = 0;
- E->viscosity.SDEPV = E->viscosity.TDEPV = 0;
+ E->viscosity.SDEPV = E->viscosity.TDEPV = E->viscosity.BDEPV =
+ E->viscosity.CDEPV = 0;
E->viscosity.EXPX = 0;
E->viscosity.SMOOTH = 0;
+ input_boolean("plasticity_dimensional",&(E->viscosity.plasticity_dimensional),"on",m);
for(i = 0; i < 40; i++)
{
@@ -65,59 +71,116 @@
E->viscosity.Z[i] = 0.0;
E->viscosity.E[i] = 0.0;
E->viscosity.T0[i] = 0.0;
- }
+
+ /*
+ byerlee or plasticity law
+ */
+ if(E->viscosity.plasticity_dimensional){
+ /* for byerlee, dimensional values are used */
+ E->viscosity.abyerlee[i]=2.1e4; /* t_y = (a * z[m] + b) *p
+ whereas a is rho * g * 0.6(from friction)
+ and b is 60MPa cohesion
+ */
+ E->viscosity.bbyerlee[i]=6e7;
+
+ E->viscosity.lbyerlee[i]=.1;
+ }else{
+ /* for plasticity: t_y = min (a + (1-x) * b, l) */
+ E->viscosity.abyerlee[i]=0.0;
+ E->viscosity.bbyerlee[i]=1.0;
+
+ E->viscosity.lbyerlee[i]=1e20;
+ }
+
+ /* comp dep visc */
+ E->viscosity.pre_comp[i] = 1.0;
+
+ } /* end layer loop */
+
/* read in information */
- input_int("rheol", &(E->viscosity.RHEOL), "essential");
- input_int("num_mat", &(E->viscosity.num_mat), "1");
+ input_int("rheol", &(E->viscosity.RHEOL), "essential",m);
+ input_int("num_mat", &(E->viscosity.num_mat), "1",m);
- input_float_vector("viscT", E->viscosity.num_mat, (E->viscosity.T)); /* redundant */
- input_float_vector("viscT1", E->viscosity.num_mat, (E->viscosity.T));
- input_float_vector("viscZ", E->viscosity.num_mat, (E->viscosity.Z));
- input_float_vector("viscE", E->viscosity.num_mat, (E->viscosity.E));
- input_float_vector("viscT0", E->viscosity.num_mat, (E->viscosity.T0));
- input_float_vector("visc0", E->viscosity.num_mat, (E->viscosity.N0)); /* redundant */
- input_float_vector("viscN0", E->viscosity.num_mat, (E->viscosity.N0));
- input_boolean("TDEPV", &(E->viscosity.TDEPV), "on");
- input_boolean("SDEPV", &(E->viscosity.SDEPV), "off");
+ input_float_vector("viscT", E->viscosity.num_mat, (E->viscosity.T),m); /* redundant */
+ input_float_vector("viscT1", E->viscosity.num_mat, (E->viscosity.T),m);
+ input_float_vector("viscZ", E->viscosity.num_mat, (E->viscosity.Z),m);
+ input_float_vector("viscE", E->viscosity.num_mat, (E->viscosity.E),m);
+ input_float_vector("viscT0", E->viscosity.num_mat, (E->viscosity.T0),m);
+ input_float_vector("visc0", E->viscosity.num_mat, (E->viscosity.N0),m); /* redundant */
+ input_float_vector("viscN0", E->viscosity.num_mat, (E->viscosity.N0),m);
- input_float("sdepv_misfit", &(E->viscosity.sdepv_misfit), "0.001");
- input_float_vector("sdepv_expt", E->viscosity.num_mat, (E->viscosity.sdepv_expt));
- input_float_vector("sdepv_trns", E->viscosity.num_mat, (E->viscosity.sdepv_trns));
+ input_boolean("TDEPV", &(E->viscosity.TDEPV), "on",m);
+ input_boolean("SDEPV", &(E->viscosity.SDEPV), "off",m);
+ input_boolean("BDEPV",&(E->viscosity.BDEPV),"off",m);
+ input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
- input_boolean("TDEPV_AVE", &(E->viscosity.TDEPV_AVE), "off");
- input_boolean("VFREEZE", &(E->viscosity.FREEZE), "off");
- input_boolean("VMAX", &(E->viscosity.MAX), "off");
- input_boolean("VMIN", &(E->viscosity.MIN), "off");
- input_boolean("VISC_UPDATE", &(E->viscosity.update_allowed), "on");
+ /* plasticity offset viscosity */
+ input_float("plasticity_viscosity_offset", &(E->viscosity.plasticity_viscosity_offset),"0.0",m);
- input_float("freeze_thresh", &(E->viscosity.freeze_thresh), "0.0");
- input_float("freeze_value", &(E->viscosity.freeze_value), "1.0");
- input_float("visc_max", &(E->viscosity.max_value), "nodefault");
- input_float("visc_min", &(E->viscosity.min_value), "nodefault");
+ /*
+ byerlee
+ */
+ input_float_vector("abyerlee",E->viscosity.num_mat,
+ (E->viscosity.abyerlee),m);
+ input_float_vector("bbyerlee",E->viscosity.num_mat,
+ (E->viscosity.bbyerlee),m);
+ input_float_vector("lbyerlee",E->viscosity.num_mat,
+ (E->viscosity.lbyerlee),m);
+
+
+ /* 1: transition 0: min/max transition */
+ input_boolean("plasticity_trans",&(E->viscosity.plasticity_trans),"on",m);
+ /*
+
+ */
+
+ /* composition factors */
+ input_float_vector("pre_comp",2,(E->viscosity.pre_comp),m);
- input_boolean("VISC_GUESS", &(E->viscosity.guess), "off");
- input_string("visc_old_file", E->viscosity.old_file, " ");
+
+ input_float("sdepv_misfit", &(E->viscosity.sdepv_misfit), "0.001",m);
+ input_float_vector("sdepv_expt", E->viscosity.num_mat, (E->viscosity.sdepv_expt),m);
+ input_float_vector("sdepv_trns", E->viscosity.num_mat, (E->viscosity.sdepv_trns),m);
+
+ /* iteration damping for alpha < 1 */
+ input_float("sdepv_iter_damp", &(E->viscosity.sdepv_iter_damp), "1.0",m);
+
+ input_boolean("TDEPV_AVE", &(E->viscosity.TDEPV_AVE), "off",m);
+ input_boolean("VFREEZE", &(E->viscosity.FREEZE), "off",m);
+ input_boolean("VMAX", &(E->viscosity.MAX), "off",m);
+ input_boolean("VMIN", &(E->viscosity.MIN), "off",m);
+ input_boolean("VISC_UPDATE", &(E->viscosity.update_allowed), "on",m);
+
+ input_float("freeze_thresh", &(E->viscosity.freeze_thresh), "0.0",m);
+ input_float("freeze_value", &(E->viscosity.freeze_value), "1.0",m);
+ input_float("visc_max", &(E->viscosity.max_value), "nodefault",m);
+ input_float("visc_min", &(E->viscosity.min_value), "nodefault",m);
+
+ input_boolean("VISC_GUESS", &(E->viscosity.guess), "off",m);
+ input_string("visc_old_file", E->viscosity.old_file, " ",m);
+
return;
}
void get_viscosity_option(struct All_variables *E)
{
/* general, essential default */
+ int m = E->parallel.me ;
- input_string("Viscosity", E->viscosity.STRUCTURE, NULL); /* Which form of viscosity */
+ input_string("Viscosity", E->viscosity.STRUCTURE, NULL,m); /* Which form of viscosity */
- input_boolean("VISC_EQUIVDD", &(E->viscosity.EQUIVDD), "off"); /* Whether to average it */
- input_int("equivdd_opt", &(E->viscosity.equivddopt), "1");
- input_int("equivdd_x", &(E->viscosity.proflocx), "1");
- input_int("equivdd_y", &(E->viscosity.proflocy), "1");
+ input_boolean("VISC_EQUIVDD", &(E->viscosity.EQUIVDD), "off",m); /* Whether to average it */
+ input_int("equivdd_opt", &(E->viscosity.equivddopt), "1",m);
+ input_int("equivdd_x", &(E->viscosity.proflocx), "1",m);
+ input_int("equivdd_y", &(E->viscosity.proflocy), "1",m);
- input_int("update_every_steps", &(E->control.KERNEL), "1");
+ input_int("update_every_steps", &(E->control.KERNEL), "1",m);
- input_boolean("VISC_SMOOTH", &(E->viscosity.SMOOTH), "off");
- input_int("visc_smooth_cycles", &(E->viscosity.smooth_cycles), "0");
+ input_boolean("VISC_SMOOTH", &(E->viscosity.SMOOTH), "off",m);
+ input_int("visc_smooth_cycles", &(E->viscosity.smooth_cycles), "0",m);
if(strcmp(E->viscosity.STRUCTURE, "system") == 0) /* Interpret */
{
@@ -156,9 +219,16 @@
else
visc_from_mat(E, visc, evisc);
+ if(E->viscosity.CDEPV)
+ visc_from_C(E, visc, evisc, propogate);
+
if(E->viscosity.SDEPV)
visc_from_S(E, visc, evisc, propogate);
+ if(E->viscosity.BDEPV)
+ visc_from_B(E, visc, evisc, propogate);
+
+
if(E->viscosity.SMOOTH)
apply_viscosity_smoother(E, visc, evisc);
@@ -180,8 +250,12 @@
evisc[(i - 1) * vpts + j] = E->viscosity.min_value;
}
- /* v_to_nodes(E,evisc,visc,E->mesh.levmax); */
+#ifdef USE_GZDIR
+ /* this is much preferred over v_to_nodes */
+ visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);
+#endif
+ /* v_to_nodes(E,evisc,visc,E->mesh.levmax); */
return;
}
@@ -208,14 +282,14 @@
void visc_from_mat(struct All_variables *E, float *Eta, float *EEta)
{
//int i, j, k, l, z, jj, kk;
- int i, jj;
-
- for(i = 1; i <= E->lmesh.nel; i++)
- for(jj = 1; jj <= vpoints[E->mesh.nsd]; jj++)
- {
- EEta[(i - 1) * vpoints[E->mesh.nsd] + jj] = E->viscosity.N0[E->mat[i] - 1];
- }
-
+ int i, jj,l;
+ for(i = 1; i <= E->lmesh.nel; i++){
+ l = E->mat[i] - 1;
+ for(jj = 1; jj <= vpoints[E->mesh.nsd]; jj++){
+ EEta[(i - 1) * vpoints[E->mesh.nsd] + jj] = E->viscosity.N0[l];
+ }
+ }
+
return;
}
@@ -271,79 +345,140 @@
if(visits == 0 || E->monitor.solution_cycles % E->control.KERNEL == 0)
{
- if(E->viscosity.RHEOL == 0)
- {
- for(i = 1; i <= nel; i++)
- {
- l = E->mat[i];
- e = (i - 1) % E->lmesh.elz + 1;
+ switch(E->viscosity.RHEOL){
+ case 0:
+ /* eta = eta0 exp(E * (1-T)) */
+ for(i = 1; i <= nel; i++)
+ {
+ l = E->mat[i] - 1; /* moved -1 up here */
+ e = (i - 1) % E->lmesh.elz + 1;
+
+ tempa = E->viscosity.N0[l];
+
+ for(kk = 1; kk <= ends; kk++)
+ TT[kk] = E->T[E->ien[i].node[kk]];
+
+ for(jj = 1; jj <= vpts; jj++)
+ {
+ temp = 1.0e-32;
+ for(kk = 1; kk <= ends; kk++)
+ {
+ temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
+ }
+ EEta[(i - 1) * vpts + jj] = tempa * exp(E->viscosity.E[l] * (one - temp ));
+ }
+ }
+ break;
+ case 1:
+ /*
+ eta = eta0 * exp(E/(T+T0))
+ */
+ for(i = 1; i <= nel; i++)
+ {
+ l = E->mat[i] - 1 ;
+ tempa = E->viscosity.N0[l];
+
+ for(kk = 1; kk <= ends; kk++)
+ TT[kk] = E->T[E->ien[i].node[kk]];
+
+ for(jj = 1; jj <= vpts; jj++)
+ {
+ temp = 1.0e-32;
+ for(kk = 1; kk <= ends; kk++)
+ {
+ temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
+ }
+ EEta[(i - 1) * vpts + jj] = tempa * exp((E->viscosity.E[l]) / (temp + E->viscosity.T[l]));
+ }
+ }
+ case 2:
+ /*
+ eta = eta0 * exp(E + (1-z)*Z0/(T+T0))
+ */
+ for(i = 1; i <= nel; i++)
+ {
+ l = E->mat[i] - 1 ;
- tempa = E->viscosity.N0[l - 1];
+ tempa = E->viscosity.N0[l];
+
+ for(kk = 1; kk <= ends; kk++)
+ TT[kk] = E->T[E->ien[i].node[kk]];
+
+ for(jj = 1; jj <= vpts; jj++)
+ {
+ temp = 1.0e-32;
+ zz = 0;
+ for(kk = 1; kk <= ends; kk++)
+ {
+ temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
+ zz += Xtmp[3][E->ien[i].node[kk]] * E->N.vpt[GNVINDEX(kk, jj)];
+ }
+ EEta[(i - 1) * vpts + jj] = tempa * exp((E->viscosity.E[l] + (1 - zz) * E->viscosity.Z[l]) /(temp + E->viscosity.T[l]) );
+ }
+ }
+ break;
+ case 3:
+ /*
+ eta = eta0 exp(E * (T0 - T))
- for(kk = 1; kk <= ends; kk++)
- TT[kk] = E->T[E->ien[i].node[kk]];
+ (this is like rheol=0 with T0=1, but I left rheol==0
+ unchanged for backward compatibility
- for(jj = 1; jj <= vpts; jj++)
- {
- temp = 1.0e-32;
- for(kk = 1; kk <= ends; kk++)
- {
- temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
- }
- EEta[(i - 1) * vpts + jj] = tempa * exp(E->viscosity.E[l - 1] * (one - temp ));
- }
- }
- }
- else if(E->viscosity.RHEOL == 1)
- {
- for(i = 1; i <= nel; i++)
- {
- l = E->mat[i];
- tempa = E->viscosity.N0[l - 1];
+ */
+ for(i = 1; i <= nel; i++)
+ {
+ l = E->mat[i] - 1;
+ e = (i - 1) % E->lmesh.elz + 1;
+
+ tempa = E->viscosity.N0[l];
+
+ for(kk = 1; kk <= ends; kk++)
+ TT[kk] = E->T[E->ien[i].node[kk]];
+
+ for(jj = 1; jj <= vpts; jj++)
+ {
+ temp = 1.0e-32;
+ for(kk = 1; kk <= ends; kk++)
+ {
+ temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
+ }
+ EEta[(i - 1) * vpts + jj] = tempa * exp(E->viscosity.E[l] * (E->viscosity.T[l] - temp ));
+ }
+ }
+ break;
+ case 4:
+ /*
+ eta = eta0 * exp(E(T_c - T) + (1-z) * Z )
+ */
+ for(i = 1; i <= nel; i++) /* loop over all elements */
+ {
+ l = E->mat[i] - 1; /* material element, determines [,,,,] */
+ tempa = E->viscosity.N0[l]; /* prefactor */
+
+ for(kk = 1; kk <= ends; kk++) /* temperature at nodes of element */
+ TT[kk] = E->T[E->ien[i].node[kk]];
+
+ for(jj = 1; jj <= vpts; jj++) /* loop through integration points of element */
+ {
+ zz = 0;
+ temp = 1.0e-32;
+ for(kk = 1; kk <= ends; kk++) /* loop through points in element */
+ {
+ temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
+ zz += Xtmp[3][E->ien[i].node[kk]] * E->N.vpt[GNVINDEX(kk, jj)]; /* depth */
+ }
+ /* now we have the temperature at this integration point */
+ EEta[(i - 1) * vpts + jj] =
+ tempa * exp(E->viscosity.E[l] * (E->viscosity.T[l] - temp) + (1 - zz) * E->viscosity.Z[l] );
+ }
+ }
+ default:
+ myerror(E,"RHEOL option undefined");
+ break;
+ } /* end swith */
- for(kk = 1; kk <= ends; kk++)
- TT[kk] = E->T[E->ien[i].node[kk]];
-
- for(jj = 1; jj <= vpts; jj++)
- {
- temp = 1.0e-32;
- for(kk = 1; kk <= ends; kk++)
- {
- temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
- }
- EEta[(i - 1) * vpts + jj] = tempa * exp((E->viscosity.E[l - 1]) / (temp + E->viscosity.T[l - 1]));
- }
- }
- }
- else if(E->viscosity.RHEOL == 2)
- {
- for(i = 1; i <= nel; i++)
- {
- l = E->mat[i];
- tempa = E->viscosity.N0[l - 1];
-
- for(kk = 1; kk <= ends; kk++)
- TT[kk] = E->T[E->ien[i].node[kk]];
-
- for(jj = 1; jj <= vpts; jj++)
- {
- temp = 1.0e-32;
- zz = 0;
- for(kk = 1; kk <= ends; kk++)
- {
- temp += max(zero, TT[kk]) * E->N.vpt[GNVINDEX(kk, jj)];;
- zz += Xtmp[3][E->ien[i].node[kk]] * E->N.vpt[GNVINDEX(kk, jj)];
- }
- EEta[(i - 1) * vpts + jj] = tempa * exp((E->viscosity.E[l - 1] + (1 - zz) * E->viscosity.Z[l - 1]) /(temp + E->viscosity.T[l - 1]) );
- }
- }
- }
- else if(E->viscosity.RHEOL == 3)
- {
- }
-
- visits++;
-
+ visits++;
+
}
@@ -522,13 +657,23 @@
{
int llayers = 0;
- if(x3 >= E->viscosity.zlith)
- llayers = 1;
- else if(x3 < E->viscosity.zlith && x3 >= E->viscosity.zlm)
- llayers = 2;
- else if(x3 < E->viscosity.zlm)
- llayers = 3;
-
+ /* this is the old logic, llayer = 4 would never get assigned */
+ /* if(x3 >= E->viscosity.zlith) */
+ /* llayers = 1; */
+ /* else if((x3 < E->viscosity.zlith) && (x3 >= E->viscosity.zlm)) */
+ /* llayers = 2; */
+ /* else if(x3 < E->viscosity.zlm) */
+ /* llayers = 3; */
+
+ if(x3 >= E->viscosity.zlith) /* above 410 */
+ llayers = 1;
+ else if(x3 >= E->viscosity.z410) /* above 410 */
+ llayers = 2;
+ else if(x3 >= E->viscosity.zlm) /* above 660 */
+ llayers = 3;
+ else /* lower mantle */
+ llayers = 4;
+
return (llayers);
}
@@ -572,3 +717,201 @@
return (thickness);
}
+
+/*
+
+limit the viscosity by a plastic-type yielding with depth dependent
+yield stress (aka byerlee law)
+
+*/
+//#define DEBUG
+void visc_from_B(struct All_variables *E, float *Eta, float *EEta, int propogate)
+{
+ static int visited = 0;
+ float scale,stress_magnitude,depth,exponent1;
+ float *eedot;
+ float zzz,zz[9];
+ float tau,ettby,ettnew;
+ int m,l,z,jj,kk,i;
+ static float ndz_to_m;
+#ifdef DEBUG
+ FILE *out;
+#endif
+ const int vpts = vpoints[E->mesh.nsd];
+ const int nel = E->lmesh.nel;
+ const int ends = enodes[E->mesh.nsd];
+
+ eedot = (float *) safe_malloc((2+nel)*sizeof(float));
+
+#ifdef DEBUG
+ out=fopen("tmp.visc","w");
+#endif
+ if (!visited) {
+ /*
+ scaling from nod-dim radius (0...1) to m
+
+ (only used for dimensional version)
+
+ */
+ ndz_to_m = E->monitor.length_scale;
+
+ /* */
+ if(E->parallel.me==0){ /* control output */
+ for(l=1;l <= E->viscosity.num_mat;l++) {
+ fprintf(stderr,"Plasticity: %d/%d: a=%g b=%g p=%g\n",
+ l,E->viscosity.num_mat,
+ E->viscosity.abyerlee[l-1],
+ E->viscosity.bbyerlee[l-1],
+ E->viscosity.lbyerlee[l-1]);
+ }
+ fprintf(stderr,"\tdim: %i trans: %i offset: %g\n",
+ E->viscosity.plasticity_dimensional,
+ E->viscosity.plasticity_trans,
+ E->viscosity.plasticity_viscosity_offset);
+ }
+ /*
+ get strain rates for all elements
+ */
+ for(i=1;i<=nel;i++)
+ eedot[i] = 1.0;
+
+ }else{
+ strain_rate_2_inv(E,eedot,1);
+ }
+
+ for(i=1;i <= nel;i++) {
+ /*
+ loop through all elements
+ */
+ l = E->mat[i] - 1; /* material of element */
+ for(kk=1;kk <= ends;kk++){/* loop through integration points*/
+ /* depth in meters */
+ if(E->control.Rsphere)
+ zz[kk] = (1.0 - E->SX[3][E->ien[i].node[kk]]);
+ else
+ zz[kk] = (1.0 - E->X[3][E->ien[i].node[kk]]);
+ if(E->viscosity.plasticity_dimensional)
+ zz[kk] *= ndz_to_m; /* scale to meters */
+ }
+ for(jj=1;jj <= vpts;jj++) {
+ /* loop over nodes in element */
+ zzz = 0.0;
+ for(kk=1;kk <= ends;kk++) {
+ /*
+ depth [m]
+ */
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
+ if(E->viscosity.plasticity_dimensional){
+ /* byerlee type */
+
+ /*
+ yield stress in [Pa]
+ */
+ tau=(E->viscosity.abyerlee[l] * zzz +
+ E->viscosity.bbyerlee[l]) * E->viscosity.lbyerlee[l];
+ /*
+ scaled stress
+ */
+ tau /= E->monitor.tau_scale;
+ }else{
+
+ tau = E->viscosity.abyerlee[l] * zzz + E->viscosity.bbyerlee[l];
+
+ tau = min(tau, E->viscosity.lbyerlee[l]);
+ }
+ /*
+
+ `byerlee viscosity' : tau = 2 eps eta, this is non-dim
+ plus some offset as in Stein et al.
+
+ */
+ ettby = tau/(2.0 * (eedot[i]+1e-7)) + E->viscosity.plasticity_viscosity_offset;
+ /*
+
+ decide on the plasticity transition
+
+
+ */
+ if(E->viscosity.plasticity_trans){
+ /*
+ eta = 1./(1./eta(k)+1./ettby)
+ */
+
+ ettnew = 1.0/(1.0/EEta[ (i-1)*vpts + jj ] + 1.0/ettby);
+ //fprintf(stderr,"a: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
+ }else{
+ /*
+ min(\eta_p, \eta_visc )
+ */
+ ettnew = min(EEta[ (i-1)*vpts + jj ],ettby);
+ //fprintf(stderr,"m: %g %g %g\n",EEta[ (i-1)*vpts + jj ],ettby,ettnew);
+ }
+#ifdef DEBUG
+ /* output format
+
+ z[km] tau[MPa] eps[s^-1] eta_b[Pas] eta_T[Pas] eta_c[Pas]
+
+ */
+ if(visited)
+ fprintf(out,"%10.2f %17.4e %17.4e %17.4e %17.4e %17.4e\n",
+ zzz/1e3,tau*E->monitor.tau_scale/1e6,
+ eedot[i]/E->monitor.time_scale,
+ ettby*E->data.ref_viscosity,
+ EEta[ (i-1)*vpts + jj ]*E->data.ref_viscosity,
+ ettnew*E->data.ref_viscosity);
+#endif
+ // if(visited)
+ // fprintf(stderr,"%11g %11g %11g %11g\n",ettnew,EEta[ (i-1)*vpts + jj ] ,ettby,eedot[i]);
+ EEta[ (i-1)*vpts + jj ] = ettnew;
+ }
+ }
+#ifdef DEBUG
+ fclose(out);
+#endif
+ visited = 1;
+ free ((void *)eedot);
+ return;
+}
+/*
+
+multiply with composition factor
+
+*/
+void visc_from_C(struct All_variables *E, float *Eta, float *EEta, int propogate)
+{
+ float comp,comp_fac,CC[9],tcomp;
+ double vmean,cc_loc;
+ int m,l,z,jj,kk,i;
+ static int visited = 0;
+ static double logv[2];
+ const int vpts = vpoints[E->mesh.nsd];
+ const int nel = E->lmesh.nel;
+ const int ends = enodes[E->mesh.nsd];
+ if(!visited){
+ /* log of the material viscosities */
+ logv[0] = log(E->viscosity.pre_comp[0]);
+ logv[1] = log(E->viscosity.pre_comp[1]);
+ }
+ for(i = 1; i <= nel; i++)
+ {
+ /* determine composition of each of the nodes of the element */
+ for(kk = 1; kk <= ends; kk++){
+ CC[kk] = E->C[E->ien[i].node[kk]];
+ if(CC[kk] < 0)CC[kk]=0.0;
+ if(CC[kk] > 1)CC[kk]=1.0;
+ }
+ for(jj = 1; jj <= vpts; jj++)
+ {
+ /* compute mean composition */
+ cc_loc = 0.0;
+ for(kk = 1; kk <= ends; kk++){/* the vpt takes care of averaging */
+ cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
+ }
+ /* geometric mean of viscosity */
+ vmean = exp(cc_loc * logv[1] + (1.0-cc_loc) * logv[0]);
+ EEta[ (i-1)*vpts + jj ] *= vmean;
+ } /* end jj loop */
+ }
+ visited++;
+}
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2008-11-07 03:29:17 UTC (rev 13268)
@@ -41,6 +41,10 @@
#include <stdio.h>
#include <stdlib.h>
+#ifdef USE_GGRD
+#include "hc.h"
+#endif
+
#if defined(__osf__)
void *Malloc1();
#endif
@@ -434,6 +438,7 @@
struct SIEN *sien;
double dircos[4][4];
+ float *vtk_base;
int output_llmax;
@@ -643,9 +648,12 @@
int solution_cycles;
- float time_scale;
+ float time_scale,time_scale_ma;
float length_scale;
- float viscosity_scale;
+ float viscosity_scale;
+ float velo_scale;
+ float tau_scale;
+
float geoscale;
float tpgscale;
float grvscale;
@@ -738,9 +746,19 @@
int CONMAN;
int stokes;
+ int check_t_irange,check_c_irange;
+
float Ra_670, clapeyron670, transT670, width670;
float Ra_410, clapeyron410, transT410, width410;
+#ifdef USE_GGRD
+ struct ggrd_master ggrd;
+ int slab_slice;
+ float slab_theta_bound;
+
+#endif
+
+
int adi_heating, visc_heating;
int composition;
float z_comp;
@@ -756,7 +774,7 @@
float plate_vel;
float plate_age;
- float tole_comp;
+ //float tole_comp;
float sob_tolerance;
@@ -797,6 +815,7 @@
int total_iteration_cycles;
int total_v_solver_calls;
+ int gzdir;
int record_every;
int record_all_until;
@@ -945,6 +964,7 @@
float *Fas670, *Fas410, *Fas670_b, *Fas410_b;
+
float *Vi, *EVi;
float *diffusivity, *expansivity;
float *T, *C, *CE, *buoyancy;
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2008-11-07 03:29:17 UTC (rev 13268)
@@ -66,7 +66,7 @@
void convection_initial_temperature(struct All_variables *E);
void process_restart_tc(struct All_variables *E);
void convection_initial_markers1(struct All_variables *E);
-void convection_initial_markers(struct All_variables *E);
+void convection_initial_markers(struct All_variables *,int );
void setup_plume_problem(struct All_variables *E);
void PG_process(struct All_variables *E, int ii);
/* Drive_solvers.c */
@@ -220,20 +220,20 @@
double CPU_time0(void);
void parallel_process_sync(void);
/* Parsing.c */
-void setup_parser(struct All_variables *E, int ac, char **av);
-void shutdown_parser(struct All_variables *E);
-int add_to_parameter_list(register char *name, register char *value);
-int compute_parameter_hash_table(register char *s);
-int input_int(char *name, int *value, char *interpret);
-int input_string(char *name, char *value, char *Default);
-int input_boolean(char *name, int *value, char *interpret);
-int input_float(char *name, float *value, char *interpret);
-int input_double(char *name, double *value, char *interpret);
-int input_int_vector(char *name, int number, int *value);
-int input_char_vector(char *name, int number, char *value);
-int input_float_vector(char *name, int number, float *value);
-int input_double_vector(char *name, int number, double *value);
-int interpret_control_string(char *interpret, int *essential, double **Default, double **minvalue, double **maxvalue);
+void setup_parser(struct All_variables *, char *);
+void shutdown_parser(struct All_variables *);
+void add_to_parameter_list(char *, char *);
+int compute_parameter_hash_table(char *);
+int input_int(char *, int *, char *, int);
+int input_string(char *, char *, char *, int);
+int input_boolean(char *, int *, char *, int);
+int input_float(char *, float *, char *, int);
+int input_double(char *, double *, char *, int);
+int input_int_vector(char *, int, int *, int);
+int input_char_vector(char *, int, char *, int);
+int input_float_vector(char *, int, float *, int);
+int input_double_vector(char *, int, double *, int);
+int interpret_control_string(char *, int *, double *, double *, double *);
/* Phase_change.c */
void phase_change(struct All_variables *E, float *B6, float *B_b6, float *B4, float *B_b4);
/* Process_buoyancy.c */
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2008-11-07 00:02:47 UTC (rev 13267)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2008-11-07 03:29:17 UTC (rev 13268)
@@ -80,6 +80,7 @@
int SDEPV;
float sdepv_misfit;
+ float sdepv_iter_damp;
int sdepv_normalize;
float sdepv_expt[40];
float sdepv_trns[40];
@@ -90,6 +91,46 @@
float E[40], T0[40];
float T[40], Z[40];
+
+ /* byerlee :
+ abyerlee, bbyerlee, pbyerlee
+ Byerlees law, dimensional:
+ Yield strength is exceeded if the
+ second invariant of the stress reaches:
+
+ (abyerlee * (h-z_dim) + bbyerlee) * lbyerlee
+
+ abyerlee: gradient, in Pa/m
+ bbyerlee: strength at surface (Pa)
+ lbyerlee: factor due to pore pressure, reduching the strength by the
+ amount of Plam. 1: dry, 0.7: hydrostatic pressure, 0:
+ lithostatic pressure. Caution: 0 might produce crap.
+ */
+ int BDEPV;
+ float abyerlee[40],bbyerlee[40],
+ lbyerlee[40];
+
+
+ float plasticity_viscosity_offset;
+
+ int plasticity_trans; /* 1: effective viscosity approach
+ 0: min viscosity approach
+
+ */
+ int plasticity_dimensional; /* 1: use Byerlee type setting with
+ dimensional values
+ 0: use non-dimensional values for yield stress
+
+ */
+
+
+ int CDEPV; /* composition dependent viscosity */
+ float pre_comp[40]; /* prefactors */
+
+
+
+
+
int weak_blobs;
float weak_blobx[40];
float weak_bloby[40];
More information about the CIG-COMMITS
mailing list