[cig-commits] r11221 - in mc/3D/CitcomS/trunk: CitcomS lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Feb 20 17:03:44 PST 2008
Author: tan2
Date: 2008-02-20 17:03:43 -0800 (Wed, 20 Feb 2008)
New Revision: 11221
Modified:
mc/3D/CitcomS/trunk/CitcomS/BaseApplication.py
mc/3D/CitcomS/trunk/CitcomS/Controller.py
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/module/initial_conditions.c
Log:
Finished interpolating fields onto tracers
Modified: mc/3D/CitcomS/trunk/CitcomS/BaseApplication.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/BaseApplication.py 2008-02-21 01:02:41 UTC (rev 11220)
+++ mc/3D/CitcomS/trunk/CitcomS/BaseApplication.py 2008-02-21 01:03:43 UTC (rev 11221)
@@ -60,7 +60,11 @@
'''
self.initialize()
self.reportConfiguration()
- self.launch()
+
+ try:
+ self.launch()
+ except SystemExit:
+ pass
return
Modified: mc/3D/CitcomS/trunk/CitcomS/Controller.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Controller.py 2008-02-21 01:02:41 UTC (rev 11220)
+++ mc/3D/CitcomS/trunk/CitcomS/Controller.py 2008-02-21 01:03:43 UTC (rev 11221)
@@ -68,6 +68,10 @@
# do io for 0th step
self.save()
+ if self.solver.inventory.ic.inventory.post_p:
+ self.endSimulation()
+ raise SystemExit()
+
### XXX: if stokes: advection tracers and terminate
return
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-02-21 01:02:41 UTC (rev 11220)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-02-21 01:03:43 UTC (rev 11221)
@@ -995,12 +995,15 @@
int iwedge = shp[0];
if (iwedge==1)
- return data[1]*shp[1] + data[2]*shp[2] + shp[3]*data[3]
- + data[6]*shp[4] + data[6]*shp[5] + shp[7]*data[6];
+ return data[1]*shp[1] + data[2]*shp[2] + data[3]*shp[3]
+ + data[5]*shp[4] + data[6]*shp[5] + data[7]*shp[6];
if (iwedge==2)
- return data[1]*shp[1] + data[3]*shp[2] + shp[4]*data[3]
- + data[5]*shp[4] + data[7]*shp[5] + shp[8]*data[6];
+ return data[1]*shp[1] + data[3]*shp[2] + data[4]*shp[3]
+ + data[5]*shp[4] + data[7]*shp[5] + data[8]*shp[6];
+
+ fprintf(stderr, "full_interpolate_data: shouldn't be here\n");
+ exit(2);
}
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2008-02-21 01:02:41 UTC (rev 11220)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2008-02-21 01:03:43 UTC (rev 11221)
@@ -59,7 +59,6 @@
const int m = 1;
int n, ncolumns, ncomp;
double *compositions;
- double (*fields)[9];
snprintf(output_file, 255, "%s.intp_fields.%d",
E->control.data_file, E->parallel.me);
@@ -87,6 +86,13 @@
* [flavor0, flavor1, radius, temperature, composition(s)]
*/
+ if(E->parallel.me == 0) {
+ fprintf(E->fp, "Temperature contrast is %e Kelvin\n",
+ E->data.ref_temperature);
+ fprintf(stderr, "Temperature contrast is %e Kelvin\n",
+ E->data.ref_temperature);
+ }
+
ncolumns = 4;
if(E->composition.on) {
ncolumns += E->composition.ncomp;
@@ -95,14 +101,6 @@
/* get the horizontal average of temperature and composition */
compute_horiz_avg(E);
- /* allocate memory for fields that need to be interpolated,
- * ie. excluding [flavor0, flavor1, radius] */
- fields = malloc((ncolumns-3) * sizeof(fields));
- if(fields == NULL) {
- fprintf(stderr, "output_interpolated_fields(): 2 not enough memory.\n");
- exit(1);
- }
-
fprintf(fp1,"%d %d %d %d\n",
E->trace.ntracers[m], E->trace.itracer_interpolate_fields,
ncolumns, ncomp);
@@ -110,8 +108,9 @@
for(n=1; n<=E->trace.ntracers[m]; n++) {
int i, j, k;
- int nelem;
- double shpfn[9];
+ int nelem, flavor0, flavor1;
+ int node[9], nz[9];
+ double shpfn[9], data[9];
double theta, phi, rad;
double temperature;
@@ -120,32 +119,58 @@
phi = E->trace.basicq[m][1][n];
rad = E->trace.basicq[m][2][n];
+ flavor0 = E->trace.extraq[m][0][n];
+ flavor1 = E->trace.extraq[m][1][n];
+
+ /* get shape functions at the tracer location */
+ if(E->parallel.nprocxy == 12)
+ full_get_shape_functions(E, shpfn, nelem, theta, phi, rad);
+ else
+ regional_get_shape_functions(E, shpfn, nelem, theta, phi, rad);
+
/* fetch element data for interpolation */
for(i=1; i<=ENODES3D; i++) {
- int node = E->ien[m][nelem].node[i];
- int nz = (node - 1) % E->lmesh.noz + 1;
- fields[0][i] = E->T[m][node] - E->Have.T[nz];
- for(j=0, k=1; j<E->composition.ncomp; j++, k++)
- fields[k][i] = E->composition.comp_node[m][j][node]
- - E->Have.C[j][nz];
+ node[i] = E->ien[m][nelem].node[i];
+ nz[i] = (node[i] - 1) % E->lmesh.noz + 1;
}
- if(E->parallel.nprocxy == 12) {
- full_get_shape_functions(E, shpfn, nelem, theta, phi, rad);
- temperature = full_interpolate_data(E, shpfn, fields[0]);
- for(j=0, k=1; j<E->composition.ncomp; j++, k++)
- compositions[j] = full_interpolate_data(E, shpfn, fields[k]);
+ for(i=1; i<=ENODES3D; i++) {
+ data[i] = E->T[m][node[i]] - E->Have.T[nz[i]];
}
- else {
- regional_get_shape_functions(E, shpfn, nelem, theta, phi, rad);
- temperature = regional_interpolate_data(E, shpfn, fields[0]);
- for(j=0, k=1; j<E->composition.ncomp; j++, k++)
- compositions[j] = regional_interpolate_data(E, shpfn, fields[k]);
+
+ if(E->parallel.nprocxy == 12)
+ temperature = full_interpolate_data(E, shpfn, data);
+ else
+ temperature = regional_interpolate_data(E, shpfn, data);
+
+ /** debug **
+ fprintf(E->trace.fpt, "result: %e data: %e %e %e %e %e %e %e %e\n",
+ temperature, data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8]);
+ /**/
+
+ for(j=0; j<E->composition.ncomp; j++) {
+ for(i=1; i<=ENODES3D; i++) {
+ data[i] = E->composition.comp_node[m][j][node[i]]
+ - E->Have.C[j][nz[i]];
+ }
+ if(E->parallel.nprocxy == 12)
+ compositions[j] = full_interpolate_data(E, shpfn, data);
+ else
+ compositions[j] = regional_interpolate_data(E, shpfn, data);
+
+ /** debug **
+ fprintf(E->trace.fpt, "result: %e data: %e %e %e %e %e %e %e %e\n",
+ compositions[j], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8]);
+ /**/
}
+ /* dimensionalize */
+ rad *= 1e3 * E->data.radius_km;
+ temperature *= E->data.ref_temperature;
+
+ /* output */
fprintf(fp1,"%d %d %e %e",
- E->trace.extraq[m][0][n], E->trace.extraq[m][1][n],
- rad, temperature);
+ flavor0, flavor1, rad, temperature);
for(j=0; j<E->composition.ncomp; j++) {
fprintf(fp1," %e", compositions[j]);
@@ -153,7 +178,6 @@
fprintf(fp1, "\n");
}
- free(fields);
break;
case 100:
/* user modification here */
@@ -168,7 +192,7 @@
}
- if(!compositions)
+ if(E->composition.on)
free(compositions);
fclose(fp1);
Modified: mc/3D/CitcomS/trunk/module/initial_conditions.c
===================================================================
--- mc/3D/CitcomS/trunk/module/initial_conditions.c 2008-02-21 01:02:41 UTC (rev 11220)
+++ mc/3D/CitcomS/trunk/module/initial_conditions.c 2008-02-21 01:03:43 UTC (rev 11221)
@@ -211,8 +211,6 @@
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
post_processing(E);
- (E->problem_output)(E, E->monitor.solution_cycles);
- parallel_process_termination();
Py_INCREF(Py_None);
return Py_None;
More information about the cig-commits
mailing list