[cig-commits] r16139 - mc/3D/CitcomS/trunk/lib
danb at geodynamics.org
danb at geodynamics.org
Sat Jan 16 17:19:25 PST 2010
Author: danb
Date: 2010-01-16 17:19:24 -0800 (Sat, 16 Jan 2010)
New Revision: 16139
Modified:
mc/3D/CitcomS/trunk/lib/Composition_related.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/tracer_defs.h
Log:
Preliminary support for the absolute tracer method. Only tested with one flavor, reading in from a file, and the initial volume of the tracers is HARD-CODED! In addition, the code automaticaly truncates the composition to one, although this should be a switch for the user to specify.
Modified: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c 2010-01-16 02:47:04 UTC (rev 16138)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c 2010-01-17 01:19:24 UTC (rev 16139)
@@ -35,6 +35,7 @@
static void allocate_composition_memory(struct All_variables *E);
static void compute_elemental_composition_ratio_method(struct All_variables *E);
+static void compute_elemental_composition_absolute_method(struct All_variables *E);
static void init_bulk_composition(struct All_variables *E);
static void check_initial_composition(struct All_variables *E);
static void fill_composition_from_neighbors(struct All_variables *E);
@@ -55,11 +56,11 @@
/* ibuoy_type=1 (ratio method) */
input_int("buoy_type",&(E->composition.ibuoy_type),"1,0,nomax",m);
- if (E->composition.ibuoy_type!=1) {
+ /* if (E->composition.ibuoy_type!=1) {
fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
fflush(stderr);
parallel_process_termination();
- }
+ } */
if (E->composition.ibuoy_type==0)
E->composition.ncomp = E->trace.nflavors;
@@ -161,9 +162,6 @@
void fill_composition(struct All_variables *E)
{
- /* XXX: Currently, only the ratio method works here. */
- /* Will have to come back here to include the absolute method. */
-
/* ratio method */
if (E->composition.ibuoy_type==1) {
@@ -171,12 +169,15 @@
}
/* absolute method */
+ if (E->composition.ibuoy_type==0) {
+ compute_elemental_composition_absolute_method(E);
+ }
- if (E->composition.ibuoy_type!=1) {
+ /* if (E->composition.ibuoy_type!=1) {
fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
fflush(E->trace.fpt);
exit(10);
- }
+ } */
/* Map elemental composition to nodal points */
@@ -241,20 +242,21 @@
void init_composition(struct All_variables *E)
{
- /* XXX: Currently, only the ratio method works here. */
- /* Will have to come back here to include the absolute method. */
-
/* ratio method */
if (E->composition.ibuoy_type==1) {
compute_elemental_composition_ratio_method(E);
}
/* absolute method */
- if (E->composition.ibuoy_type!=1) {
+ if (E->composition.ibuoy_type==0) {
+ compute_elemental_composition_absolute_method(E);
+ }
+
+ /* if (E->composition.ibuoy_type!=1) {
fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
fflush(E->trace.fpt);
exit(10);
- }
+ } */
/* for empty elements */
check_initial_composition(E);
@@ -336,6 +338,70 @@
return;
}
+
+
+/*********** COMPUTE ELEMENTAL ABSOLUTE METHOD *************************/
+/* */
+/* This function computes the composition per element. */
+/* The concentration of material i in an element is */
+/* defined as: */
+/* tracer i mass * (# of tracers of flavor i) / ( volume of element) */
+/* */
+/* where tracer i mass = tot tracer volume / tot num of tracers */
+/* */
+/* Added by DJB 01/14/10. Known caveats: (will address at some point) */
+/* 1) XXX: needs better error handling */
+/* 2) XXX: only tested with one flavor of tracer (flavor 0) */
+/* 3) XXX: only tested by reading in a file of tracer locations */
+/* 4) XXX: the initial volume of the tracer domain is HARD-CODED! */
+
+static void compute_elemental_composition_absolute_method(struct All_variables *E)
+{
+ int i, j, e, flavor, numtracers;
+ double domain_volume;
+ float comp;
+ float one = 1.0;
+
+ /* This needs to be `manually' changed for the total volume */
+ /* occupied by your tracers */
+ domain_volume = 1e-2;
+
+ for (j=1; j<=E->sphere.caps_per_proc; j++) {
+ for (e=1; e<=E->lmesh.nel; e++) {
+ numtracers = 0;
+ for (flavor=0; flavor<E->trace.nflavors; flavor++)
+ numtracers += E->trace.ntracer_flavor[j][flavor][e];
+
+ /* Check for empty entries */
+ /* If no tracers are in an element, comp = 0.0 (i.e. is ambient) */
+ if (numtracers == 0) {
+ for(i=0;i<E->composition.ncomp;i++) {
+ E->composition.comp_el[j][i][e] = 0.0;
+ }
+ continue;
+ }
+
+ /* Composition is proportional to (local) tracer density */
+ for(i=0;i<E->composition.ncomp;i++) {
+ flavor = i;
+ comp =
+ E->trace.ntracer_flavor[j][flavor][e] / E->eco[j][e].area
+ * domain_volume / E->trace.number_of_tracers;
+
+ /* truncate composition at 1.0 */
+ /* This violates mass conservation but prevents unphysical C */
+ /* XXX: make truncation a switch for the user to specify */
+ E->composition.comp_el[j][i][e] = min(comp,one);
+
+ }
+ }
+ }
+
+ return;
+}
+
+
+
/********** MAP COMPOSITION TO NODES ****************/
/* */
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2010-01-16 02:47:04 UTC (rev 16138)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2010-01-17 01:19:24 UTC (rev 16139)
@@ -938,6 +938,9 @@
icushion=100;
+ /* for absolute tracer method */
+ E->trace.number_of_tracers = number_of_tracers;
+
iestimate=number_of_tracers/E->parallel.nproc + icushion;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2010-01-16 02:47:04 UTC (rev 16138)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2010-01-17 01:19:24 UTC (rev 16139)
@@ -58,6 +58,8 @@
int max_ntracers[13];
int *ielement[13];
+ int number_of_tracers;
+
int ilatersize[13];
int ilater[13];
double *rlater[13][100];
More information about the CIG-COMMITS
mailing list