[cig-commits] r7840 - in mc/3D/CitcomS/trunk: CitcomS/Components
lib module tests/checkpoint visual
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Aug 17 11:33:38 PDT 2007
Author: tan2
Date: 2007-08-17 11:33:36 -0700 (Fri, 17 Aug 2007)
New Revision: 7840
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
mc/3D/CitcomS/trunk/lib/Checkpoints.c
mc/3D/CitcomS/trunk/lib/Composition_related.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg
mc/3D/CitcomS/trunk/visual/autocombine.py
mc/3D/CitcomS/trunk/visual/batchcombine.py
mc/3D/CitcomS/trunk/visual/dxgeneral.py
Log:
Multi-component chemical convection.
By setting tracer_flavors > 2, one can enable multi-component chemical convection.
See tests/checkpoint/tracer.cfg for example input.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-08-17 18:33:36 UTC (rev 7840)
@@ -46,6 +46,11 @@
def setProperties(self, stream):
+ # convert lists of strings to lists of ints/floats
+ inv = self.inventory
+ inv.z_interface = map(float, inv.z_interface)
+ inv.buoyancy_ratio = map(float, inv.buoyancy_ratio)
+
from CitcomSLib import Tracer_set_properties
Tracer_set_properties(self.all_variables, self.inventory, stream)
return
@@ -78,7 +83,7 @@
# either using absolute method or ratio method.
tracer_flavors = inv.int("tracer_flavors", default=0)
ic_method_for_flavors = inv.int("ic_method_for_flavors", default=0)
- z_interface = inv.float("z_interface", default=0.7)
+ z_interface = inv.list("z_interface", default=[0.7])
# Regular grid parameters
@@ -93,7 +98,7 @@
# ibuoy_type=0 (absolute method, not implemented)
# ibuoy_type=1 (ratio method)
buoy_type = inv.int("buoy_type", default=1)
- buoyancy_ratio = inv.float("buoyancy_ratio", default=1.0)
+ buoyancy_ratio = inv.list("buoyancy_ratio", default=[1.0])
# This is not used anymore and is left here for backward compatibility
reset_initial_composition = inv.bool("reset_initial_composition",
Modified: mc/3D/CitcomS/trunk/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Checkpoints.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Checkpoints.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -335,20 +335,25 @@
static void composition_checkpoint(struct All_variables *E, FILE *fp)
{
- int m;
+ int i, m;
write_sentinel(fp);
- fwrite(&(E->composition.bulk_composition), sizeof(double), 1, fp);
- fwrite(&(E->composition.initial_bulk_composition), sizeof(double), 1, fp);
+ fwrite(&(E->composition.ncomp), sizeof(int), 1, fp);
+ fwrite(E->composition.bulk_composition, sizeof(double),
+ E->composition.ncomp, fp);
+ fwrite(E->composition.initial_bulk_composition, sizeof(double),
+ E->composition.ncomp, fp);
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- fwrite(E->composition.comp_el[m], sizeof(double),
- E->lmesh.nel+1, fp);
+ for(i=0; i<E->composition.ncomp; i++)
+ fwrite(E->composition.comp_el[m][i], sizeof(double),
+ E->lmesh.nel+1, fp);
}
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- fwrite(E->composition.comp_node[m], sizeof(double),
- E->lmesh.nno+1, fp);
+ for(i=0; i<E->composition.ncomp; i++)
+ fwrite(E->composition.comp_node[m][i], sizeof(double),
+ E->lmesh.nno+1, fp);
}
return;
@@ -358,31 +363,44 @@
static void read_composition_checkpoint(struct All_variables *E, FILE *fp)
{
double tmp;
- int m;
+ int m, i, itmp;
read_sentinel(fp, E->parallel.me);
- fread(&tmp, sizeof(double), 1, fp);
- E->composition.bulk_composition = tmp;
+ fread(&itmp, sizeof(int), 1, fp);
+ if (itmp != E->composition.ncomp) {
+ fprintf(stderr, "Error in reading checkpoint file: ncomp, me=%d\n",
+ E->parallel.me);
+ fprintf(stderr, "%d\n", itmp);
+ exit(-1);
+ }
- fread(&tmp, sizeof(double), 1, fp);
- E->composition.initial_bulk_composition = tmp;
+ fread(E->composition.bulk_composition, sizeof(double),
+ E->composition.ncomp, fp);
+ fread(E->composition.initial_bulk_composition, sizeof(double),
+ E->composition.ncomp, fp);
+
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- fread(E->composition.comp_el[m], sizeof(double),
- E->lmesh.nel+1, fp);
+ for(i=0; i<E->composition.ncomp; i++)
+ fread(E->composition.comp_el[m][i], sizeof(double),
+ E->lmesh.nel+1, fp);
}
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- fread(E->composition.comp_node[m], sizeof(double),
- E->lmesh.nno+1, fp);
+ for(i=0; i<E->composition.ncomp; i++)
+ fread(E->composition.comp_node[m][i], sizeof(double),
+ E->lmesh.nno+1, fp);
}
/* preventing uninitialized access */
E->trace.istat_iempty = 0;
- E->composition.error_fraction = E->composition.bulk_composition
- / E->composition.initial_bulk_composition - 1.0;
+ for (i=0; i<E->composition.ncomp; i++) {
+ E->composition.error_fraction[i] = E->composition.bulk_composition[i]
+ / E->composition.initial_bulk_composition[i] - 1.0;
+ }
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -42,6 +42,7 @@
void composition_input(struct All_variables *E)
{
+ int i;
int m = E->parallel.me;
input_boolean("chemical_buoyancy",
@@ -50,9 +51,6 @@
if (E->composition.ichemical_buoyancy) {
- input_double("buoyancy_ratio",
- &(E->composition.buoyancy_ratio),"1.0",m);
-
/* ibuoy_type=0 (absolute method) */
/* ibuoy_type=1 (ratio method) */
@@ -63,6 +61,21 @@
parallel_process_termination();
}
+ if (E->composition.ibuoy_type==0)
+ E->composition.ncomp = E->trace.nflavors;
+ else if (E->composition.ibuoy_type==1)
+ E->composition.ncomp = E->trace.nflavors - 1;
+
+ E->composition.buoyancy_ratio = (double*) malloc(E->composition.ncomp
+ *sizeof(double));
+
+ /* default values .... */
+ for (i=0; i<E->composition.ncomp; i++)
+ E->composition.buoyancy_ratio[i] = 1.0;
+
+ input_double_vector("buoyancy_ratio", E->composition.ncomp,
+ E->composition.buoyancy_ratio,m);
+
}
@@ -97,6 +110,7 @@
void write_composition_instructions(struct All_variables *E)
{
+ int k;
E->composition.on = 0;
if (E->composition.ichemical_buoyancy ||
@@ -116,12 +130,14 @@
fprintf(E->trace.fpt,"Active Tracers\n");
- if (E->composition.ibuoy_type==1)
+ if (E->composition.ibuoy_type==1)
fprintf(E->trace.fpt,"Ratio Method\n");
- if (E->composition.ibuoy_type==0)
+ if (E->composition.ibuoy_type==0)
fprintf(E->trace.fpt,"Absolute Method\n");
- fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->composition.buoyancy_ratio);
+ for(k=0; k<E->composition.ncomp; k++) {
+ fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->composition.buoyancy_ratio[k]);
+ }
/*
if (E->composition.icompositional_rheology==0) {
@@ -173,22 +189,41 @@
static void allocate_composition_memory(struct All_variables *E)
{
- int j;
+ int i, j;
+ for (i=0; i<E->composition.ncomp; i++) {
+ E->composition.bulk_composition = (double*) malloc(E->composition.ncomp*sizeof(double));
+ E->composition.initial_bulk_composition = (double*) malloc(E->composition.ncomp*sizeof(double));
+ E->composition.error_fraction = (double*) malloc(E->composition.ncomp*sizeof(double));
+ }
+
/* allocat memory for composition fields at the nodes and elements */
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- if ((E->composition.comp_el[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 8989y\n");
+ if ((E->composition.comp_el[j]=(double **)malloc((E->composition.ncomp)*sizeof(double*)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 8987y\n");
fflush(E->trace.fpt);
exit(10);
}
-
- if ((E->composition.comp_node[j]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 983rk\n");
+ if ((E->composition.comp_node[j]=(double **)malloc((E->composition.ncomp)*sizeof(double*)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 8988y\n");
fflush(E->trace.fpt);
exit(10);
}
+
+ for (i=0; i<E->composition.ncomp; i++) {
+ if ((E->composition.comp_el[j][i]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 8989y\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ if ((E->composition.comp_node[j][i]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(allocate_composition_memory)-no memory 983rk\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
}
return;
@@ -197,7 +232,7 @@
void init_composition(struct All_variables *E)
{
- if (E->composition.ichemical_buoyancy &&
+ if (E->composition.ichemical_buoyancy &&
E->composition.ibuoy_type) {
fill_composition(E);
check_initial_composition(E);
@@ -233,18 +268,10 @@
static void compute_elemental_composition_ratio_method(struct All_variables *E)
{
- int j, e, flavor, numtracers;
+ int i, j, e, flavor, numtracers;
int iempty = 0;
- /* XXX: currently only two composition is supported */
- if (E->trace.nflavors != 2) {
- fprintf(E->trace.fpt, "Sorry - Only two flavors of tracers is supported\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
-
for (j=1; j<=E->sphere.caps_per_proc; j++) {
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
@@ -259,10 +286,11 @@
continue;
}
- /* XXX: generalize for more than one composition */
- flavor = 1;
- E->composition.comp_el[j][e] =
- E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
+ for(i=0;i<E->composition.ncomp;i++) {
+ flavor = i + 1;
+ E->composition.comp_el[j][i][e] =
+ E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
+ }
}
@@ -288,8 +316,8 @@
static void map_composition_to_nodes(struct All_variables *E)
{
-
- int kk;
+ double *tmp[NCS];
+ int i, n, kk;
int nelem, nodenum;
int j;
@@ -297,8 +325,10 @@
for (j=1;j<=E->sphere.caps_per_proc;j++) {
/* first, initialize node array */
- for (kk=1;kk<=E->lmesh.nno;kk++)
- E->composition.comp_node[j][kk]=0.0;
+ for(i=0;i<E->composition.ncomp;i++) {
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ E->composition.comp_node[j][i][kk]=0.0;
+ }
/* Loop through all elements */
for (nelem=1;nelem<=E->lmesh.nel;nelem++) {
@@ -308,34 +338,42 @@
/* weight composition */
for (nodenum=1;nodenum<=8;nodenum++) {
+ n = E->ien[j][nelem].node[nodenum];
+ for(i=0;i<E->composition.ncomp;i++) {
- E->composition.comp_node[j][E->ien[j][nelem].node[nodenum]] +=
- E->composition.comp_el[j][nelem]*
- E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
-
+ E->composition.comp_node[j][i][n] +=
+ E->composition.comp_el[j][i][nelem]*
+ E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
+ }
}
} /* end nelem */
} /* end j */
+ for(i=0;i<E->composition.ncomp;i++) {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ tmp[j] = E->composition.comp_node[j][i];
- (E->exchange_node_d)(E,E->composition.comp_node,E->mesh.levmax);
+ (E->exchange_node_d)(E,tmp,E->mesh.levmax);
+ }
-
/* Divide by nodal volume */
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- for (kk=1;kk<=E->lmesh.nno;kk++)
- E->composition.comp_node[j][kk] *= E->MASS[E->mesh.levmax][j][kk];
+ for(i=0;i<E->composition.ncomp;i++)
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ E->composition.comp_node[j][i][kk] *= E->MASS[E->mesh.levmax][j][kk];
/* testing */
/**
- for (kk=1;kk<=E->lmesh.nel;kk++) {
- fprintf(E->trace.fpt,"%d %f\n",kk,E->composition.comp_el[j][kk]);
- }
+ for(i=0;i<E->composition.ncomp;i++)
+ for (kk=1;kk<=E->lmesh.nel;kk++) {
+ fprintf(E->trace.fpt,"%d %f\n",kk,E->composition.comp_el[j][i][kk]);
+ }
- for (kk=1;kk<=E->lmesh.nno;kk++) {
- fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->composition.comp_node[j][kk]);
- }
+ for(i=0;i<E->composition.ncomp;i++)
+ for (kk=1;kk<=E->lmesh.nno;kk++) {
+ fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->composition.comp_node[j][i][kk]);
+ }
fflush(E->trace.fpt);
/**/
@@ -352,15 +390,23 @@
double return_bulk_value_d();
double volume;
- int ival=0;
+ double *tmp[NCS];
+ int i, m;
+ const int ival=0;
- /* ival=0 returns integral not average */
- volume = return_bulk_value_d(E,E->composition.comp_node,ival);
+ for (i=0; i<E->composition.ncomp; i++) {
- E->composition.bulk_composition = volume;
- E->composition.initial_bulk_composition = volume;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ tmp[m] = E->composition.comp_node[m][i];
+ /* ival=0 returns integral not average */
+ volume = return_bulk_value_d(E,tmp,ival);
+
+ E->composition.bulk_composition[i] = volume;
+ E->composition.initial_bulk_composition[i] = volume;
+ }
+
return;
}
@@ -370,15 +416,22 @@
double return_bulk_value_d();
double volume;
+ double *tmp[NCS];
+ int i, m;
const int ival = 0;
- /* ival=0 returns integral not average */
- volume=return_bulk_value_d(E,E->composition.comp_node,ival);
+ for (i=0; i<E->composition.ncomp; i++) {
- E->composition.bulk_composition=volume;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ tmp[m] = E->composition.comp_node[m][i];
- E->composition.error_fraction=((volume-E->composition.initial_bulk_composition)/
- E->composition.initial_bulk_composition);
+ /* ival=0 returns integral not average */
+ volume = return_bulk_value_d(E,tmp,ival);
+ E->composition.bulk_composition[i] = volume;
+
+ E->composition.error_fraction[i] = (volume - E->composition.initial_bulk_composition[i]) / E->composition.initial_bulk_composition[i];
+ }
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -1836,6 +1836,8 @@
/**** WRITE TRACE INSTRUCTIONS ***************/
static void write_trace_instructions(struct All_variables *E)
{
+ int i;
+
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
@@ -1866,7 +1868,8 @@
fprintf(E->trace.fpt,"Initialized tracer flavors by: %d\n", E->trace.ic_method_for_flavors);
if (E->trace.ic_method_for_flavors == 0) {
fprintf(E->trace.fpt,"Layered tracer flavors\n");
- fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
+ for (i=0; i<E->trace.nflavors-1; i++)
+ fprintf(E->trace.fpt,"Interface Height: %d %f\n",i,E->trace.z_interface[i]);
}
#ifdef USE_GGRD
else if(E->trace.ic_method_for_flavors == 1) {
@@ -1882,7 +1885,16 @@
}
}
+ for (i=0; i<E->trace.nflavors-2; i++) {
+ if (E->trace.z_interface[i] < E->trace.z_interface[i+1]) {
+ fprintf(E->trace.fpt,"Sorry - The %d-th z_interface is smaller than the next one.\n", i);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ }
+
+
/* regular grid stuff */
fprintf(E->trace.fpt,"Regular Grid-> deltheta: %f delphi: %f\n",
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -454,7 +454,7 @@
void output_comp_nd(struct All_variables *E, int cycles)
{
- int i, j;
+ int i, j, k;
char output_file[255];
FILE *fp1;
@@ -463,14 +463,21 @@
fp1 = output_open(output_file);
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%3d %7d %.5e %.5e %.5e\n",
+ fprintf(fp1,"%3d %7d %.5e %d\n",
j, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
+ E->monitor.elapsed_time, E->composition.ncomp);
+ for(i=0;i<E->composition.ncomp;i++) {
+ fprintf(fp1,"%.5e %.5e ",
+ E->composition.initial_bulk_composition[i],
+ E->composition.bulk_composition[i]);
+ }
+ fprintf(fp1,"\n");
for(i=1;i<=E->lmesh.nno;i++) {
- fprintf(fp1,"%.6e\n",E->composition.comp_node[j][i]);
+ for(k=0;k<E->composition.ncomp;k++) {
+ fprintf(fp1,"%.6e ",E->composition.comp_node[j][k][i]);
+ }
+ fprintf(fp1,"\n");
}
}
@@ -482,7 +489,7 @@
void output_comp_el(struct All_variables *E, int cycles)
{
- int i, j;
+ int i, j, k;
char output_file[255];
FILE *fp1;
@@ -491,14 +498,21 @@
fp1 = output_open(output_file);
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%3d %7d %.5e %.5e %.5e\n",
+ fprintf(fp1,"%3d %7d %.5e %d\n",
j, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
+ E->monitor.elapsed_time, E->composition.ncomp);
+ for(i=0;i<E->composition.ncomp;i++) {
+ fprintf(fp1,"%.5e %.5e ",
+ E->composition.initial_bulk_composition[i],
+ E->composition.bulk_composition[i]);
+ }
+ fprintf(fp1,"\n");
- for(i=1;i<=E->lmesh.nel;i++) {
- fprintf(fp1,"%.6e\n",E->composition.comp_el[j][i]);
+ for(i=1;i<=E->lmesh.nno;i++) {
+ for(k=0;k<E->composition.ncomp;k++) {
+ fprintf(fp1,"%.6e ",E->composition.comp_el[j][k][i]);
+ }
+ fprintf(fp1,"\n");
}
}
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -156,10 +156,12 @@
/* chemical buoyancy */
if(E->control.tracer &&
(E->composition.ichemical_buoyancy)) {
- temp2 = E->composition.buoyancy_ratio * temp;
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++)
- buoy[m][i] -= temp2 * E->composition.comp_node[m][i];
+ for(j=0;j<E->composition.ncomp;j++) {
+ temp2 = E->composition.buoyancy_ratio[j] * temp;
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.nno;i++)
+ buoy[m][i] -= temp2 * E->composition.comp_node[m][j][i];
+ }
}
phase_change_apply_410(E, buoy);
Modified: mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -146,6 +146,8 @@
/**** WRITE TRACE INSTRUCTIONS ***************/
static void write_trace_instructions(struct All_variables *E)
{
+ int i;
+
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
@@ -166,7 +168,8 @@
fprintf(E->trace.fpt,"Initialized tracer flavors by: %d\n", E->trace.ic_method_for_flavors);
if (E->trace.ic_method_for_flavors == 0) {
fprintf(E->trace.fpt,"Layered tracer flavors\n");
- fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
+ for (i=0; i<E->trace.nflavors-1; i++)
+ fprintf(E->trace.fpt,"Interface Height: %d %f\n",i,E->trace.z_interface[i]);
}
else {
fprintf(E->trace.fpt,"Sorry-This IC methods for Flavors are Unavailable %d\n",E->trace.ic_method_for_flavors);
@@ -175,8 +178,16 @@
}
}
+ for (i=0; i<E->trace.nflavors-2; i++) {
+ if (E->trace.z_interface[i] < E->trace.z_interface[i+1]) {
+ fprintf(E->trace.fpt,"Sorry - The %d-th z_interface is smaller than the next one.\n", i);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ }
+
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -84,6 +84,7 @@
void full_tracer_input();
char message[100];
int m=E->parallel.me;
+ int i;
input_boolean("tracer",&(E->control.tracer),"off",m);
input_boolean("tracer_enriched",
@@ -133,7 +134,13 @@
&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
switch(E->trace.ic_method_for_flavors){
case 0: /* layer */
- input_double("z_interface",&(E->trace.z_interface),"0.7",m);
+ E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
+ *sizeof(double));
+ for(i=0; i<E->trace.nflavors-1; i++)
+ E->trace.z_interface[i] = 0.7;
+
+ input_double_vector("z_interface", E->trace.nflavors-1,
+ E->trace.z_interface, m);
break;
case 1: /* from grid in top n materials */
input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
@@ -227,11 +234,8 @@
void tracer_post_processing(struct All_variables *E)
{
- char output_file[200];
+ int i;
- double convection_time,tracer_time;
- double trace_fraction,total_time;
-
if(E->control.verbose){
fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
@@ -261,9 +265,17 @@
get_bulk_composition(E);
if (E->parallel.me==0) {
- fprintf(E->fp,"composition: %e %e\n",E->monitor.elapsed_time,E->composition.bulk_composition);
- fprintf(E->fp,"composition_error_fraction: %e %e\n",E->monitor.elapsed_time,E->composition.error_fraction);
+ fprintf(E->fp,"composition: %e",E->monitor.elapsed_time);
+ for (i=0; i<E->composition.ncomp; i++)
+ fprintf(E->fp," %e", E->composition.bulk_composition[i]);
+ fprintf(E->fp,"\n");
+
+ fprintf(E->fp,"composition_error_fraction: %e",E->monitor.elapsed_time);
+ for (i=0; i<E->composition.ncomp; i++)
+ fprintf(E->fp," %e", E->composition.error_fraction[i]);
+ fprintf(E->fp,"\n");
+
}
}
@@ -1155,21 +1167,29 @@
static void init_tracer_flavors(struct All_variables *E)
{
int j, kk, number_of_tracers;
+ int i;
+ double flavor;
double rad;
switch(E->trace.ic_method_for_flavors){
case 0:
/* ic_method_for_flavors == 0 (layered structure) */
- /* any tracer above z_interface is of flavor 0 */
- /* any tracer below z_interface is of flavor 1 */
+ /* any tracer above z_interface[i] is of flavor i */
+ /* any tracer below z_interface is of flavor (nflavors-1) */
for (j=1;j<=E->sphere.caps_per_proc;j++) {
number_of_tracers = E->trace.ntracers[j];
for (kk=1;kk<=number_of_tracers;kk++) {
rad = E->trace.basicq[j][2][kk];
- if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
- if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
+ flavor = E->trace.nflavors - 1;
+ for (i=0; i<E->trace.nflavors-1; i++) {
+ if (rad > E->trace.z_interface[i]) {
+ flavor = i;
+ break;
+ }
+ }
+ E->trace.extraq[j][0][kk] = flavor;
}
}
break;
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-17 18:33:36 UTC (rev 7840)
@@ -622,15 +622,16 @@
/* if any of the above flags is true, turn this flag on */
int on;
- double buoyancy_ratio;
int ibuoy_type;
+ int ncomp;
+ double *buoyancy_ratio;
- double *comp_el[13];
- double *comp_node[13];
+ double **comp_el[13];
+ double **comp_node[13];
- double initial_bulk_composition;
- double bulk_composition;
- double error_fraction;
+ double *initial_bulk_composition;
+ double *bulk_composition;
+ double *error_fraction;
double compositional_rheology_prefactor;
};
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-08-17 18:33:36 UTC (rev 7840)
@@ -67,13 +67,13 @@
int **ntracer_flavor[13];
int ic_method_for_flavors;
- double z_interface;
-
- char ggrd_file[255]; /* for grd input */
- int ggrd_layers;
+ double *z_interface;
+ char ggrd_file[255]; /* for grd input */
+ int ggrd_layers;
+
/* statistical parameters */
int istat_ichoice[13][5];
int istat_isend;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-17 18:33:36 UTC (rev 7840)
@@ -71,7 +71,11 @@
float* vector, int len, FILE* fp);
#define getFloatVectorProperty(p, a, v, l, o) if (-1 == _getFloatVectorProperty(p, a, v, l, o)) return NULL
+int _getDoubleVectorProperty(PyObject* properties, char* attribute,
+ double* vector, int len, FILE* fp);
+#define getDoubleVectorProperty(p, a, v, l, o) if (-1 == _getDoubleVectorProperty(p, a, v, l, o)) return NULL
+
/*==============================================================*/
@@ -604,15 +608,28 @@
getIntProperty(properties, "tracer_flavors", E->trace.nflavors, fp);
getIntProperty(properties, "ic_method_for_flavors", E->trace.ic_method_for_flavors, fp);
- if (E->trace.ic_method_for_flavors == 0)
- getDoubleProperty(properties, "z_interface", E->trace.z_interface, fp);
+ if (E->trace.ic_method_for_flavors == 0) {
+ E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
+ *sizeof(double));
+ getDoubleVectorProperty(properties, "z_interface", E->trace.z_interface, E->trace.nflavors-1, fp);
+ }
+
getIntProperty(properties, "chemical_buoyancy",
E->composition.ichemical_buoyancy, fp);
if (E->composition.ichemical_buoyancy==1) {
getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
- getDoubleProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, fp);
+
+ if (E->composition.ibuoy_type==0)
+ E->composition.ncomp = E->trace.nflavors;
+ else if (E->composition.ibuoy_type==1)
+ E->composition.ncomp = E->trace.nflavors - 1;
+
+ E->composition.buoyancy_ratio = (double*) malloc(E->composition.ncomp
+ *sizeof(double));
+
+ getDoubleVectorProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, E->composition.ncomp, fp);
}
Modified: mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg
===================================================================
--- mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/tests/checkpoint/tracer.cfg 2007-08-17 18:33:36 UTC (rev 7840)
@@ -6,6 +6,9 @@
monitoringFrequency = 1
checkpointFrequency = 1
+[CitcomS.solver.output]
+output_optional = comp_el,tracer
+
[CitcomS.solver]
datafile = zzz
@@ -21,12 +24,12 @@
tracer_ic_method = 0
tracers_per_element = 20
-tracer_flavors = 2
+tracer_flavors = 3
ic_method_for_flavors = 0
-z_interface = 0.7
+z_interface = 0.9,0.7
chemical_buoyancy = 1
-buoyancy_ratio = 0.4
+buoyancy_ratio = 0.3,0.4
regular_grid_deltheta = 1.0
regular_grid_delphi = 1.0
Modified: mc/3D/CitcomS/trunk/visual/autocombine.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/autocombine.py 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/visual/autocombine.py 2007-08-17 18:33:36 UTC (rev 7840)
@@ -34,6 +34,8 @@
# default values for CitcomS input
defaults = {'output_format': 'ascii',
'output_optional': 'surf,botm',
+ 'buoy_type': 1,
+ 'tracer_flavors': 2,
'nprocx': 1,
'nprocy': 1,
'nprocz': 1,
@@ -89,6 +91,19 @@
output_optional = parser.getstr('output_optional')
optional_fields = normalize_optional(output_optional)
+ buoy_type = parser.getint('buoy_type')
+ nflavors = parser.getint('tracer_flavors')
+ if buoy_type == 0:
+ ncompositions = nflavors
+ elif buoy_type == 1:
+ ncompositions = nflavors - 1
+ else:
+ print "Error: don't know how to combine the output", \
+ "(buoy_type=%d)" % buoy_type
+ sys.exit(1)
+
+
+
nodex = parser.getint('nodex')
nodey = parser.getint('nodey')
nodez = parser.getint('nodez')
@@ -105,14 +120,14 @@
bc.batchcombine(nodelist, datadir, datafile, timestep,
nodex, nodey, nodez,
ncap, nprocx, nprocy, nprocz,
- 'coord,velo,visc')
+ 'coord,velo,visc', 0)
# combining optional fields, if necessary
if optional_fields:
bc.batchcombine(nodelist, datadir, datafile, timestep,
nodex, nodey, nodez,
ncap, nprocx, nprocy, nprocz,
- optional_fields)
+ optional_fields, ncompositions)
Modified: mc/3D/CitcomS/trunk/visual/batchcombine.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/batchcombine.py 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/visual/batchcombine.py 2007-08-17 18:33:36 UTC (rev 7840)
@@ -29,7 +29,7 @@
'''
Paste and combine Citcom data
-Usage: batchcombine.py <machinefile | node-list> datadir datafile timestep nodex nodey nodez ncap nprocx nprocy nprocz [fields]
+Usage: batchcombine.py <machinefile | node-list> datadir datafile timestep nodex nodey nodez ncap nprocx nprocy nprocz [fields [ncompositions]]
'''
@@ -100,7 +100,8 @@
def batchcombine(nodes, datadir, datafile, timestep, nodex, nodey, nodez,
- ncap, nprocx, nprocy, nprocz, fields):
+ ncap, nprocx, nprocy, nprocz, fields,
+ ncompositions):
# paste
batchpaste(datadir, datafile, fields, timestep, nodes)
@@ -123,7 +124,7 @@
# create .general file
import dxgeneral
- dxgeneral.write(fields, combined_files)
+ dxgeneral.write(fields, ncompositions, combined_files)
return
@@ -133,7 +134,7 @@
import sys
- if not (12 <= len(sys.argv) <= 13):
+ if not (12 <= len(sys.argv) <= 14):
print __doc__
sys.exit(1)
@@ -148,16 +149,20 @@
nprocx = int(sys.argv[9])
nprocy = int(sys.argv[10])
nprocz = int(sys.argv[11])
- if len(sys.argv) == 13:
+
+ fields = 'coord,velo,visc'
+ ncompositions = 0
+ if len(sys.argv) >= 13:
fields = sys.argv[12]
- else:
- fields = 'coord,velo,visc'
+ if len(sys.argv) == 14:
+ ncompositions = int(sys.argv[13])
+
totalnodes = nprocx * nprocy * nprocz * ncap
nodelist = machinefile2nodelist(machinefile, totalnodes)
batchcombine(nodelist, datadir, datafile, timestep, nodex, nodey, nodez,
- ncap, nprocx, nprocy, nprocz, fields)
+ ncap, nprocx, nprocy, nprocz, fields, ncompositions)
Modified: mc/3D/CitcomS/trunk/visual/dxgeneral.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/dxgeneral.py 2007-08-17 02:54:31 UTC (rev 7839)
+++ mc/3D/CitcomS/trunk/visual/dxgeneral.py 2007-08-17 18:33:36 UTC (rev 7840)
@@ -28,13 +28,13 @@
'''Create OpenDX .general file for combined Citcom Data
- Usage: dxgeneral.py combined_fields file1 [file2 [...]]
+ Usage: dxgeneral.py combined_fields ncompositions file1 [file2 [...]]
'''
import os, sys
-def write(opts, filenames):
+def write(opts, ncomp, filenames):
for filename in filenames:
if not os.path.exists(filename):
@@ -47,7 +47,7 @@
f = open(outfile, 'w')
try:
- write_general_file(f, filename, opts, shape)
+ write_general_file(f, filename, opts, shape, ncomp)
finally:
f.close()
@@ -65,7 +65,7 @@
-def write_general_file(f, filename, opts, shape):
+def write_general_file(f, filename, opts, shape, ncomp):
template = '''file = %(filename)s
grid = %(shape_str)s
@@ -96,6 +96,10 @@
'velo': '3-vector, scalar',
'visc': 'scalar'}
+ # 'comp_nd' needs special treatment
+ if ncomp > 1:
+ struct['comp_nd'] = '%d-vector' % ncomp
+
# mapping from opt name to data type
type = {'comp_nd': 'float',
'coord': 'float',
@@ -108,6 +112,7 @@
opt_struct = []
opt_type = []
for opt in opts.split(','):
+ opt = opt.strip()
opt_field.append(field[opt])
opt_struct.append(struct[opt])
opt_type.append(type[opt])
@@ -130,9 +135,10 @@
sys.exit(1)
opts = sys.argv[1]
- filenames = sys.argv[2:]
+ ncomp = int(sys.argv[2])
+ filenames = sys.argv[3:]
- write(opts, filenames)
+ write(opts, ncomp, filenames)
# End of file
More information about the cig-commits
mailing list