[cig-commits] r6796 - in mc/3D/CitcomS/branches/compressible:
CitcomS/Components lib module tests/checkpoint
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon May 7 14:30:14 PDT 2007
Author: tan2
Date: 2007-05-07 14:30:13 -0700 (Mon, 07 May 2007)
New Revision: 6796
Modified:
mc/3D/CitcomS/branches/compressible/CitcomS/Components/Tracer.py
mc/3D/CitcomS/branches/compressible/lib/Checkpoints.c
mc/3D/CitcomS/branches/compressible/lib/Composition_related.c
mc/3D/CitcomS/branches/compressible/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/compressible/lib/Output.c
mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c
mc/3D/CitcomS/branches/compressible/lib/global_defs.h
mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h
mc/3D/CitcomS/branches/compressible/module/setProperties.c
mc/3D/CitcomS/branches/compressible/tests/checkpoint/tracer.cfg
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/branches/compressible/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Components/Tracer.py 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Components/Tracer.py 2007-05-07 21:30:13 UTC (rev 6796)
@@ -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/branches/compressible/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Checkpoints.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Checkpoints.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -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/branches/compressible/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Composition_related.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Composition_related.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -42,6 +42,7 @@
void composition_input(struct All_variables *E)
{
+ int i;
int m = E->parallel.me;
input_int("chemical_buoyancy",&(E->composition.ichemical_buoyancy),
@@ -49,9 +50,6 @@
if (E->composition.ichemical_buoyancy==1) {
- input_double("buoyancy_ratio",
- &(E->composition.buoyancy_ratio),"1.0",m);
-
/* ibuoy_type=0 (absolute method) */
/* ibuoy_type=1 (ratio method) */
@@ -62,6 +60,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);
+
}
@@ -96,6 +109,7 @@
void write_composition_instructions(struct All_variables *E)
{
+ int k;
E->composition.on = 0;
if (E->composition.ichemical_buoyancy==1 ||
@@ -119,7 +133,9 @@
if (E->composition.ibuoy_type==1) fprintf(E->trace.fpt,"Ratio Method\n");
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) {
@@ -171,22 +187,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;
@@ -230,18 +265,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;
@@ -256,10 +283,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;
+ }
}
@@ -285,8 +313,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;
@@ -294,8 +322,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++) {
@@ -305,34 +335,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);
/**/
@@ -349,15 +387,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;
}
@@ -367,15 +413,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/branches/compressible/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Full_tracer_advection.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Full_tracer_advection.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -1829,6 +1829,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");
@@ -1859,7 +1861,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);
@@ -1868,7 +1871,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/branches/compressible/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Output.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Output.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -434,7 +434,7 @@
void output_comp_nd(struct All_variables *E, int cycles)
{
- int i, j;
+ int i, j, k;
char output_file[255];
FILE *fp1;
@@ -443,14 +443,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",
j, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
+ E->monitor.elapsed_time);
+ 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");
}
}
@@ -462,7 +469,7 @@
void output_comp_el(struct All_variables *E, int cycles)
{
- int i, j;
+ int i, j, k;
char output_file[255];
FILE *fp1;
@@ -471,14 +478,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",
j, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
+ E->monitor.elapsed_time);
+ 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/branches/compressible/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -146,10 +146,12 @@
/* chemical buoyancy */
if(E->control.tracer && E->composition.ichemical_buoyancy==1) {
- 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/branches/compressible/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -147,6 +147,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");
@@ -167,7 +169,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);
@@ -176,8 +179,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/branches/compressible/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -76,6 +76,7 @@
{
void full_tracer_input();
int m=E->parallel.me;
+ int i;
input_int("tracer",&(E->control.tracer),"0",m);
if(E->control.tracer) {
@@ -107,9 +108,15 @@
input_int("ic_method_for_flavors",&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
- if (E->trace.ic_method_for_flavors == 0)
- input_double("z_interface",&(E->trace.z_interface),"0.7",m);
+ if (E->trace.ic_method_for_flavors == 0) {
+ 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);
+ }
if(E->parallel.nprocxy == 12)
@@ -190,12 +197,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;
-
-
fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
@@ -222,9 +225,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");
+
}
}
@@ -1095,12 +1106,14 @@
static void init_tracer_flavors(struct All_variables *E)
{
int j, kk, number_of_tracers;
+ int i;
+ double flavor;
double rad;
/* 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 */
+ /* tracers below all z_interface is of flavor (nflavors-1) */
if (E->trace.ic_method_for_flavors == 0) {
for (j=1;j<=E->sphere.caps_per_proc;j++) {
@@ -1108,8 +1121,14 @@
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;
}
}
}
Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-05-07 21:30:13 UTC (rev 6796)
@@ -598,15 +598,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/branches/compressible/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h 2007-05-07 21:30:13 UTC (rev 6796)
@@ -66,7 +66,7 @@
int **ntracer_flavor[13];
int ic_method_for_flavors;
- double z_interface;
+ double *z_interface;
/* statistical parameters */
int istat_ichoice[13][5];
Modified: mc/3D/CitcomS/branches/compressible/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-05-07 21:30:13 UTC (rev 6796)
@@ -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
+
/*==============================================================*/
@@ -601,15 +605,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/branches/compressible/tests/checkpoint/tracer.cfg
===================================================================
--- mc/3D/CitcomS/branches/compressible/tests/checkpoint/tracer.cfg 2007-05-07 21:23:11 UTC (rev 6795)
+++ mc/3D/CitcomS/branches/compressible/tests/checkpoint/tracer.cfg 2007-05-07 21:30:13 UTC (rev 6796)
@@ -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
More information about the cig-commits
mailing list