[cig-commits] r6222 - in mc/3D/CitcomS/trunk: CitcomS/Components
lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Sat Mar 10 21:20:44 PST 2007
Author: tan2
Date: 2007-03-10 21:20:43 -0800 (Sat, 10 Mar 2007)
New Revision: 6222
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
mc/3D/CitcomS/trunk/lib/Composition_related.c
mc/3D/CitcomS/trunk/lib/Full_obsolete.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/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
* Unifying output for extra tracer quantities
* Renamed rtrac to basicq, etrac to extraq, to match the basic/extra quantities for advection
* Renamed itrac to ielement, itracsize to ielementsize, ielementsize to max_ntracers
* A new member for the # of tracers
* Removed some unused members in trace struct
* Untabify and indent the code
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Tracer.py 2007-03-11 05:20:43 UTC (rev 6222)
@@ -71,9 +71,9 @@
# (tracer_ic_method == 1)
tracer_file = inv.str("tracer_file", default="tracer.dat")
+ # How many types of tracers, must be >= 1
+ tracer_types = inv.int("tracer_types", default=1)
- # TODO: remove
- z_interface = inv.float("z_interface", default=0.5)
# Advection Scheme
@@ -104,7 +104,10 @@
buoyancy_ratio = inv.float("buoyancy_ratio", default=1.0)
reset_initial_composition = inv.bool("reset_initial_composition",
default=False)
+ # TODO: remove
+ z_interface = inv.float("z_interface", default=0.5)
+
# compositional_rheology=1 (not implemented in this version, TODO:remove)
compositional_rheology = inv.bool("compositional_rheology",
default=False)
Modified: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -33,7 +33,11 @@
#include "composition_related.h"
+static void allocate_composition_memory(struct All_variables *E);
+static void compute_elemental_composition_ratio_method(struct All_variables *E);
+static void init_composition(struct All_variables *E);
static void initialize_old_composition(struct All_variables *E);
+static void map_composition_to_nodes(struct All_variables *E);
void composition_input(struct All_variables *E)
@@ -45,23 +49,21 @@
if (E->composition.ichemical_buoyancy==1) {
- input_double("buoyancy_ratio",
- &(E->composition.buoyancy_ratio),"1.0",m);
+ input_double("buoyancy_ratio",
+ &(E->composition.buoyancy_ratio),"1.0",m);
- /* ibuoy_type=0 (absolute method) */
- /* ibuoy_type=1 (ratio method) */
+ /* ibuoy_type=0 (absolute method) */
+ /* ibuoy_type=1 (ratio method) */
- input_int("buoy_type",&(E->composition.ibuoy_type),"1,0,nomax",m);
- if (E->composition.ibuoy_type!=1) {
- fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
- fflush(stderr);
- parallel_process_termination();
- }
+ input_int("buoy_type",&(E->composition.ibuoy_type),"1,0,nomax",m);
+ if (E->composition.ibuoy_type!=1) {
+ fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
- if (E->trace.ic_method==2) {
- input_int("reset_initial_composition",
- &(E->composition.ireset_initial_composition),"0",m);
- }
+ input_int("reset_initial_composition",
+ &(E->composition.ireset_initial_composition),"0",m);
input_double("z_interface",&(E->composition.z_interface),"0.5",m);
@@ -74,35 +76,37 @@
/* icompositional_rheology=1 (on) */
input_int("compositional_rheology",
- &(E->composition.icompositional_rheology),"1,0,nomax",m);
+ &(E->composition.icompositional_rheology),"1,0,nomax",m);
if (E->composition.icompositional_rheology==1) {
- input_double("compositional_prefactor",
- &(E->composition.compositional_rheology_prefactor),
- "1.0",m);
+ input_double("compositional_prefactor",
+ &(E->composition.compositional_rheology_prefactor),
+ "1.0",m);
}
+ return;
}
void composition_setup(struct All_variables *E)
{
+ int j;
E->composition.on = 0;
if (E->composition.ichemical_buoyancy==1 ||
E->composition.icompositional_rheology)
E->composition.on = 1;
- if (E->composition.ichemical_buoyancy==1 && E->composition.ibuoy_type==1) {
- initialize_old_composition(E);
- fill_composition(E);
+ if (E->composition.on) {
+ allocate_composition_memory(E);
+ init_composition(E);
}
+ return;
}
-
void init_tracer_composition(struct All_variables *E)
{
int j, kk, number_of_tracers;
@@ -110,12 +114,12 @@
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- number_of_tracers = E->trace.itrac[j][0];
+ number_of_tracers = E->trace.ntracers[j];
for (kk=1;kk<=number_of_tracers;kk++) {
- rad = E->trace.rtrac[j][2][kk];
+ rad = E->trace.basicq[j][2][kk];
- if (rad<=E->composition.z_interface) E->trace.etrac[j][0][kk]=1.0;
- if (rad>E->composition.z_interface) E->trace.etrac[j][0][kk]=0.0;
+ if (rad<=E->composition.z_interface) E->trace.extraq[j][0][kk]=1.0;
+ if (rad>E->composition.z_interface) E->trace.extraq[j][0][kk]=0.0;
}
}
return;
@@ -126,7 +130,7 @@
{
if (E->composition.ichemical_buoyancy==0)
- fprintf(E->trace.fpt,"Passive Tracers\n");
+ fprintf(E->trace.fpt,"Passive Tracers\n");
if (E->composition.ichemical_buoyancy==1)
fprintf(E->trace.fpt,"Active Tracers\n");
@@ -139,25 +143,26 @@
if (E->composition.ireset_initial_composition==0)
- fprintf(E->trace.fpt,"Using old initial composition from tracer files\n");
+ fprintf(E->trace.fpt,"Using old initial composition from tracer files\n");
else
- fprintf(E->trace.fpt,"Resetting initial composition\n");
+ fprintf(E->trace.fpt,"Resetting initial composition\n");
if (E->composition.icompositional_rheology==0) {
- fprintf(E->trace.fpt,"Compositional Rheology - OFF\n");
+ fprintf(E->trace.fpt,"Compositional Rheology - OFF\n");
}
else if (E->composition.icompositional_rheology>0) {
- fprintf(E->trace.fpt,"Compositional Rheology - ON\n");
- fprintf(E->trace.fpt,"Compositional Prefactor: %f\n",
- E->composition.compositional_rheology_prefactor);
+ fprintf(E->trace.fpt,"Compositional Rheology - ON\n");
+ fprintf(E->trace.fpt,"Compositional Prefactor: %f\n",
+ E->composition.compositional_rheology_prefactor);
}
fflush(E->trace.fpt);
fflush(stderr);
+ return;
}
@@ -165,27 +170,25 @@
void fill_composition(struct All_variables *E)
{
- void compute_elemental_composition_ratio_method();
- void map_composition_to_nodes();
-
- /* Currently, only the ratio method works here. */
+ /* XXX: Currently, only the ratio method works here. */
/* Will have to come back here to include the absolute method. */
+ accumulate_tracers_in_element(E);
+
+
/* ratio method */
- if (E->composition.ibuoy_type==1)
- {
- compute_elemental_composition_ratio_method(E);
- }
+ if (E->composition.ibuoy_type==1) {
+ compute_elemental_composition_ratio_method(E);
+ }
/* absolute method */
- if (E->composition.ibuoy_type!=1)
- {
- fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ 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 */
@@ -196,6 +199,58 @@
+static void allocate_composition_memory(struct All_variables *E)
+{
+ int j;
+
+ /* 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(compute_elemental_composition)-no memory 8989y\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(map_compostion_to_nodes)-no memory 983rk\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ if (E->composition.ibuoy_type==1) {
+ /* allocat memory for tracer ratio method */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ if ((E->composition.ieltrac[j]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 5u83a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((E->composition.celtrac[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 58hy8\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+ return;
+}
+
+
+static void init_composition(struct All_variables *E)
+{
+ if (E->composition.ichemical_buoyancy==1 && E->composition.ibuoy_type==1) {
+ initialize_old_composition(E);
+ fill_composition(E);
+ }
+ return;
+}
+
+
/************ INITIALIZE OLD COMPOSITION ************************/
static void initialize_old_composition(struct All_variables *E)
{
@@ -213,46 +268,46 @@
FILE *fp;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- if ((E->trace.oldel[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(fill old composition)-no memory 324c\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ if ((E->composition.oldel[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(fill old composition)-no memory 324c\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
if ((E->trace.ic_method==0)||(E->trace.ic_method==1))
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
- ibottom_node=E->ien[j][kk].node[1];
- zbottom=E->sx[j][3][ibottom_node];
+ ibottom_node=E->ien[j][kk].node[1];
+ zbottom=E->sx[j][3][ibottom_node];
- if (zbottom<E->composition.z_interface) E->trace.oldel[j][kk]=1.0;
- if (zbottom>=E->composition.z_interface) E->trace.oldel[j][kk]=0.0;
+ if (zbottom<E->composition.z_interface) E->composition.oldel[j][kk]=1.0;
+ if (zbottom>=E->composition.z_interface) E->composition.oldel[j][kk]=0.0;
- } /* end kk */
- } /* end j */
- }
+ } /* end kk */
+ } /* end j */
+ }
/* Else read from file */
else if (E->trace.ic_method==2)
- {
+ {
- /* first look for backing file */
+ /* first look for backing file */
- sprintf(output_file,"%s.comp_el.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
- if ( (fp=fopen(output_file,"r"))==NULL)
- {
+ sprintf(output_file,"%s.comp_el.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
+ if ( (fp=fopen(output_file,"r"))==NULL)
+ {
fprintf(E->trace.fpt,"AKMerror(Initialize Old Composition)-FILE NOT EXIST: %s\n", output_file);
fflush(E->trace.fpt);
exit(10);
@@ -260,36 +315,32 @@
fgets(input_s,200,fp);
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- fgets(input_s,200,fp);
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- fgets(input_s,200,fp);
- sscanf(input_s,"%lf",&E->trace.oldel[j][kk]);
- }
- }
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fgets(input_s,200,fp);
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%lf",&E->composition.oldel[j][kk]);
+ }
+ }
- fclose(fp);
+ fclose(fp);
- } /* endif */
+ } /* endif */
return;
}
-
-
/*********** COMPUTE ELEMENTAL COMPOSITION RATIO METHOD ***/
/* */
/* This function computes the composition per element. */
-/* This function computes the composition per element. */
/* Integer array ieltrac stores tracers per element. */
/* Double array celtrac stores the sum of tracer composition */
-void compute_elemental_composition_ratio_method(E)
- struct All_variables *E;
+static void compute_elemental_composition_ratio_method(struct All_variables *E)
{
int kk;
@@ -300,123 +351,92 @@
double comp;
- static int been_here=0;
-
- if (been_here==0)
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- if ((E->trace.ieltrac[j]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 5u83a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((E->trace.celtrac[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 58hy8\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((E->trace.comp_el[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 8989y\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- been_here++;
-
- }
-
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* first zero arrays */
+ /* first zero arrays */
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- E->trace.ieltrac[j][kk]=0;
- E->trace.celtrac[j][kk]=0.0;
- }
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->composition.ieltrac[j][kk]=0;
+ E->composition.celtrac[j][kk]=0.0;
+ }
- numtracers=E->trace.itrac[j][0];
+ numtracers=E->trace.ntracers[j];
- /* Fill ieltrac and celtrac */
+ /* Fill ieltrac and celtrac */
- for (kk=1;kk<=numtracers;kk++)
- {
+ for (kk=1;kk<=numtracers;kk++)
+ {
- nelem=E->trace.itrac[j][kk];
- E->trace.ieltrac[j][nelem]++;
+ nelem=E->trace.ielement[j][kk];
+ E->composition.ieltrac[j][nelem]++;
- comp=E->trace.etrac[j][0][kk];
+ comp=E->trace.extraq[j][0][kk];
- if (comp>1.0000001)
- {
- fprintf(E->trace.fpt,"ERROR(compute elemental)-not ready for comp>1 yet (%f)(tr. %d) \n",comp,kk);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (comp>1.0000001)
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-not ready for comp>1 yet (%f)(tr. %d) \n",comp,kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- E->trace.celtrac[j][nelem]=E->trace.celtrac[j][nelem]+comp;
+ E->composition.celtrac[j][nelem]=E->composition.celtrac[j][nelem]+comp;
- }
+ }
- /* Check for empty entries and compute ratio. */
- /* If no tracers are in an element, use previous composition */
+ /* Check for empty entries and compute ratio. */
+ /* If no tracers are in an element, use previous composition */
- iempty=0;
+ iempty=0;
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
- if (E->trace.ieltrac[j][kk]==0)
- {
- iempty++;
- E->trace.comp_el[j][kk]=E->trace.oldel[j][kk];
- }
- else if (E->trace.ieltrac[j][kk]>0)
- {
- E->trace.comp_el[j][kk]=E->trace.celtrac[j][kk]/(1.0*E->trace.ieltrac[j][kk]);
- }
+ if (E->composition.ieltrac[j][kk]==0)
+ {
+ iempty++;
+ E->composition.comp_el[j][kk]=E->composition.oldel[j][kk];
+ }
+ else if (E->composition.ieltrac[j][kk]>0)
+ {
+ E->composition.comp_el[j][kk]=E->composition.celtrac[j][kk]/(1.0*E->composition.ieltrac[j][kk]);
+ }
- if (E->trace.comp_el[j][kk]>(1.000001) || E->trace.comp_el[j][kk]<(-0.000001))
- {
- fprintf(E->trace.fpt,"ERROR(compute elemental)-noway (3u5hd)\n");
- fprintf(E->trace.fpt,"COMPEL: %f (%d)(%d)\n",E->trace.comp_el[j][kk],kk,E->trace.ieltrac[j][kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- if (iempty>0)
- {
+ if (E->composition.comp_el[j][kk]>(1.000001) || E->composition.comp_el[j][kk]<(-0.000001))
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-noway (3u5hd)\n");
+ fprintf(E->trace.fpt,"COMPEL: %f (%d)(%d)\n",E->composition.comp_el[j][kk],kk,E->composition.ieltrac[j][kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ if (iempty>0)
+ {
- /*
- fprintf(E->trace.fpt,"%d empty elements filled with old values (%f percent)\n",iempty, (100.0*iempty/E->lmesh.nel));
- fflush(E->trace.fpt);
- */
+ /*
+ fprintf(E->trace.fpt,"%d empty elements filled with old values (%f percent)\n",iempty, (100.0*iempty/E->lmesh.nel));
+ fflush(E->trace.fpt);
+ */
- if ((1.0*iempty/E->lmesh.nel)>0.80)
- {
- fprintf(E->trace.fpt,"WARNING(compute_elemental...)-number of tracers is REALLY LOW\n");
- fflush(E->trace.fpt);
- if (E->trace.itracer_warnings==1) exit(10);
- }
- }
+ if ((1.0*iempty/E->lmesh.nel)>0.80)
+ {
+ fprintf(E->trace.fpt,"WARNING(compute_elemental...)-number of tracers is REALLY LOW\n");
+ fflush(E->trace.fpt);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
+ }
- /* Fill oldel */
+ /* Fill oldel */
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- E->trace.oldel[j][kk]=E->trace.comp_el[j][kk];
- }
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->composition.oldel[j][kk]=E->composition.comp_el[j][kk];
+ }
- } /* end j */
+ } /* end j */
E->trace.istat_iempty=E->trace.istat_iempty+iempty;
@@ -427,8 +447,7 @@
/* */
-void map_composition_to_nodes(E)
- struct All_variables *E;
+static void map_composition_to_nodes(struct All_variables *E)
{
int kk;
@@ -436,85 +455,52 @@
int j;
- static int been_here=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
- if (been_here==0)
- {
+ /* first, initialize node array */
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ E->composition.comp_node[j][kk]=0.0;
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ /* Loop through all elements */
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++) {
- if ((E->trace.comp_node[j]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(map_compostion_to_nodes)-no memory 983rk\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ /* for each element, loop through element nodes */
- been_here++;
- }
+ /* weight composition */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ for (nodenum=1;nodenum<=8;nodenum++) {
- /* first, initialize node array */
+ 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];
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
- E->trace.comp_node[j][kk]=0.0;
- }
+ }
- /* Loop through all elements */
+ } /* end nelem */
+ } /* end j */
- for (nelem=1;nelem<=E->lmesh.nel;nelem++)
- {
- /* for each element, loop through element nodes */
+ (E->exchange_node_d)(E,E->composition.comp_node,E->mesh.levmax);
- /* weight composition */
- for (nodenum=1;nodenum<=8;nodenum++)
- {
-
- E->trace.comp_node[j][E->ien[j][nelem].node[nodenum]] +=
- E->trace.comp_el[j][nelem]*
- E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
-
- }
-
- } /* end nelem */
- } /* end j */
-
- /* akm modified exchange node routine for doubles */
-
- (E->exchange_node_d)(E,E->trace.comp_node,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 (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
- E->trace.comp_node[j][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]);
+ }
- /* testing */
- /*
- for (kk=1;kk<=E->lmesh.nel;kk++)
- {
- fprintf(E->trace.fpt,"%d %f\n",kk,E->trace.comp_el[j][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 (kk=1;kk<=E->lmesh.nno;kk++)
- {
- fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->trace.comp_node[j][kk]);
- }
- */
+ } /* end j */
-
- } /* end j */
-
return;
}
@@ -555,46 +541,46 @@
/* if restart = 2 tracers may or may not be reset */
/* (read initial composition from file) */
-
+ //TODO: figure out how to remove been_here
if (been_here==0)
- {
- if (E->composition.ireset_initial_composition==1)
- {
- E->composition.initial_bulk_composition=volume;
- }
- else
- {
+ {
+ if (E->composition.ireset_initial_composition==1)
+ {
+ E->composition.initial_bulk_composition=volume;
+ }
+ else
+ {
- if (E->trace.ic_method!=2)
- {
- fprintf(E->trace.fpt,"ERROR(bulk composition)-wrong reset,restart combo\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (E->trace.ic_method!=2)
+ {
+ fprintf(E->trace.fpt,"ERROR(bulk composition)-wrong reset,restart combo\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,
- E->parallel.me,E->monitor.solution_cycles);
+ sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,
+ E->parallel.me,E->monitor.solution_cycles);
- fp=fopen(output_file,"r");
- fgets(input_s,200,fp);
- sscanf(input_s,"%d %d %lf %lf %lf",
- &istep,&idum1,&rdum1,&rdum2,&rdum3);
+ fp=fopen(output_file,"r");
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%d %d %lf %lf %lf",
+ &istep,&idum1,&rdum1,&rdum2,&rdum3);
- E->composition.initial_bulk_composition=rdum2;
- fclose(fp);
+ E->composition.initial_bulk_composition=rdum2;
+ fclose(fp);
- if (istep!=E->monitor.solution_cycles)
- {
- fprintf(E->trace.fpt,"ERROR(get_bulk_composition) %d %d\n",
- istep,E->monitor.solution_cycles);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ if (istep!=E->monitor.solution_cycles)
+ {
+ fprintf(E->trace.fpt,"ERROR(get_bulk_composition) %d %d\n",
+ istep,E->monitor.solution_cycles);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
E->composition.error_fraction=((volume-E->composition.initial_bulk_composition)/
- E->composition.initial_bulk_composition);
+ E->composition.initial_bulk_composition);
parallel_process_sync(E);
@@ -625,9 +611,9 @@
if ((fp=fopen(output_file,"r"))==NULL)
{
- fprintf(stderr,"ERROR(read_temp) - %s not found\n",output_file);
- fflush(stderr);
- exit(10);
+ fprintf(stderr,"ERROR(read_temp) - %s not found\n",output_file);
+ fflush(stderr);
+ exit(10);
}
fgets(input_s,1000,fp);
@@ -635,22 +621,22 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
{
- E->trace.comp_node[m]=(double *)malloc((E->lmesh.nno+1)*sizeof(double));
+ E->composition.comp_node[m]=(double *)malloc((E->lmesh.nno+1)*sizeof(double));
- fgets(input_s,1000,fp);
- sscanf(input_s,"%d %d",&ll,&mm);
- for(i=1;i<=E->lmesh.nno;i++)
- {
- if (fgets(input_s,1000,fp)==NULL)
- {
- fprintf(stderr,"ERROR(read_comp) -data for node %d not found\n",i);
- fflush(stderr);
- exit(10);
- }
- sscanf(input_s,"%lf",&g);
- E->trace.comp_node[m][i] = g;
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d",&ll,&mm);
+ for(i=1;i<=E->lmesh.nno;i++)
+ {
+ if (fgets(input_s,1000,fp)==NULL)
+ {
+ fprintf(stderr,"ERROR(read_comp) -data for node %d not found\n",i);
+ fflush(stderr);
+ exit(10);
+ }
+ sscanf(input_s,"%lf",&g);
+ E->composition.comp_node[m][i] = g;
- }
+ }
}
fclose (fp);
Modified: mc/3D/CitcomS/trunk/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -838,7 +838,7 @@
reltrac[j]=(double *) malloc((E->lmesh.nel+1)*sizeof(double));
for (kk=1;kk<=E->lmesh.nel;kk++)
{
- reltrac[j][kk]=(1.0*E->trace.ieltrac[j][kk]);
+ reltrac[j][kk]=(1.0*E->composition.ieltrac[j][kk]);
}
}
@@ -865,7 +865,7 @@
if (E->composition.chemical_buoyancy==1)
{
- return_horiz_ave(E,E->trace.comp_node,E->trace.Have_C);
+ return_horiz_ave(E,E->composition.comp_node,E->trace.Have_C);
if (E->parallel.me<E->parallel.nprocz)
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -33,6 +33,8 @@
#include "parallel_related.h"
#include "composition_related.h"
+void accumulate_tracers_in_element(struct All_variables *E);
+
static void get_neighboring_caps(struct All_variables *E);
static void pdebug(struct All_variables *E, int i);
@@ -44,6 +46,7 @@
static void predict_tracers(struct All_variables *E);
static void correct_tracers(struct All_variables *E);
static int isum_tracers(struct All_variables *E);
+static void make_tracer_array(struct All_variables *E);
@@ -66,18 +69,24 @@
input_int("tracer_ic_method",&(E->trace.ic_method),"0,0,nomax",m);
if (E->trace.ic_method==0)
- input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax",m);
+ input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax",m);
else if (E->trace.ic_method==1)
- input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
+ input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
else if (E->trace.ic_method==2) {
}
else {
- fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
- fflush(stderr);
- parallel_process_termination();
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
}
}
+
+ /* How many types of tracers, must be >= 1 */
+ input_int("tracer_types",&(E->trace.ntypes),"1,1,nomax",m);
+
+
+
/* Advection Scheme */
/* itracer_advection_scheme=1 (simple predictor corrector -uses only V(to)) */
@@ -85,20 +94,20 @@
E->trace.itracer_advection_scheme=2;
input_int("tracer_advection_scheme",&(E->trace.itracer_advection_scheme),
- "2,0,nomax",m);
+ "2,0,nomax",m);
if (E->trace.itracer_advection_scheme==1)
- {
- }
+ {
+ }
else if (E->trace.itracer_advection_scheme==2)
- {
- }
+ {
+ }
else
- {
- fprintf(stderr,"Sorry, only advection scheme 1 and 2 available (%d)\n",E->trace.itracer_advection_scheme);
- fflush(stderr);
- parallel_process_termination();
- }
+ {
+ fprintf(stderr,"Sorry, only advection scheme 1 and 2 available (%d)\n",E->trace.itracer_advection_scheme);
+ fflush(stderr);
+ parallel_process_termination();
+ }
/* Interpolation Scheme */
@@ -107,13 +116,13 @@
E->trace.itracer_interpolation_scheme=1;
input_int("tracer_interpolation_scheme",&(E->trace.itracer_interpolation_scheme),
- "1,0,nomax",m);
+ "1,0,nomax",m);
if (E->trace.itracer_interpolation_scheme<1 || E->trace.itracer_interpolation_scheme>2)
- {
- fprintf(stderr,"Sorry, only interpolation scheme 1 and 2 available\n");
- fflush(stderr);
- parallel_process_termination();
- }
+ {
+ fprintf(stderr,"Sorry, only interpolation scheme 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
/* Regular grid parameters */
/* (first fill uniform del[0] value) */
@@ -129,7 +138,7 @@
E->trace.ianalytical_tracer_test=0;
input_int("analytical_tracer_test",&(E->trace.ianalytical_tracer_test),
- "0,0,nomax",m);
+ "0,0,nomax",m);
composition_input(E);
@@ -146,7 +155,6 @@
int m;
void write_trace_instructions();
void viscosity_checks();
- void make_tracer_array();
void initialize_old_composition();
void find_tracers();
void make_regular_grid();
@@ -159,8 +167,15 @@
void tracer_post_processing();
void restart_tracers();
+ /* Some error control */
+
+ if (E->sphere.caps_per_proc>1) {
+ fprintf(stderr,"This code does not work for multiple caps per processor!\n");
+ parallel_process_termination();
+ }
+
+
m=E->parallel.me;
- E->trace.itracing=1;
/* some obscure initial parameters */
/* This parameter specifies how close a tracer can get to the boundary */
@@ -173,24 +188,25 @@
/* advection_quantites - those needed for advection */
// TODO: generalize it
- if (E->trace.itracer_advection_scheme==1) E->trace.number_of_advection_quantities=12;
- if (E->trace.itracer_advection_scheme==2) E->trace.number_of_advection_quantities=12;
+ if (E->trace.itracer_advection_scheme==1) E->trace.number_of_basic_quantities=12;
+ if (E->trace.itracer_advection_scheme==2) E->trace.number_of_basic_quantities=12;
/* extra_quantities - used for composition, etc. */
/* (can be increased for additional science i.e. tracing chemistry */
- // TODO: how to move this part to Composition_related.c?
- if (E->composition.ichemical_buoyancy==1) E->trace.number_of_extra_quantities=1;
- else E->trace.number_of_extra_quantities=0;
+ E->trace.number_of_extra_quantities = 0;
+ if (E->composition.ichemical_buoyancy==1)
+ E->trace.number_of_extra_quantities += 1;
- E->trace.number_of_tracer_quantities=E->trace.number_of_advection_quantities +
- E->trace.number_of_extra_quantities;
+ E->trace.number_of_tracer_quantities =
+ E->trace.number_of_basic_quantities +
+ E->trace.number_of_extra_quantities;
/* Fixed positions in tracer array */
- /* Comp is always in etrac position 0 */
- /* Current coordinates are always kept in rtrac positions 0-5 */
+ /* Comp is always in extraq position 0 */
+ /* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on advection scheme and/or science being done */
/* open tracing output file */
@@ -208,36 +224,24 @@
/* Some error control regarding size of pointer arrays */
- if (E->trace.number_of_advection_quantities>99)
- {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rtrac in global.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_extra_quantities>99)
- {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of etrac in global.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_tracer_quantities>99)
- {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in global.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
+ if (E->trace.number_of_basic_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ if (E->trace.number_of_extra_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ if (E->trace.number_of_tracer_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
write_trace_instructions(E);
- /* Some error control */
-
- if (E->sphere.caps_per_proc>1)
- {
- fprintf(E->trace.fpt,"This code does not work for multiple caps per processor!\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
get_neighboring_caps(E);
if (E->trace.ic_method==0) {
@@ -254,22 +258,21 @@
}
else if (E->trace.ic_method==2) restart_tracers(E);
else
- {
- fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
/* flush and wait for not real reason but it can't hurt */
fflush(E->trace.fpt);
parallel_process_sync(E);
- if (E->trace.itracer_interpolation_scheme==1)
- {
- define_uv_space(E);
- determine_shape_coefficients(E);
- }
+ if (E->trace.itracer_interpolation_scheme==1) {
+ define_uv_space(E);
+ determine_shape_coefficients(E);
+ }
/* flush and wait for not real reason but it can't hurt */
@@ -292,15 +295,13 @@
E->trace.ilast_tracer_count = isum_tracers(E);
fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
- if (E->trace.ianalytical_tracer_test==1)
- {
- //TODO: walk into this code...
- analytical_test(E);
- parallel_process_termination();
- }
+ if (E->trace.ianalytical_tracer_test==1) {
+ //TODO: walk into this code...
+ analytical_test(E);
+ parallel_process_termination();
+ }
composition_setup(E);
-
tracer_post_processing(E);
return;
@@ -320,19 +321,11 @@
void check_sum();
void tracer_post_processing();
- void write_tracers();
fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
fflush(E->trace.fpt);
- /* presave last timesteps tracers as preback */
- //TODO: use system("mv oldfile oldfile.bak\n") instead
- if ( (E->monitor.solution_cycles % E->control.record_every)==0) {
- //TODO: migrate to output_tracer
- //write_tracers(E,3);
- }
-
/* advect tracers */
predict_tracers(E);
@@ -342,7 +335,7 @@
//TODO: move
if (E->composition.ichemical_buoyancy==1) {
- fill_composition(E);
+ fill_composition(E);
}
tracer_post_processing(E);
@@ -365,38 +358,24 @@
double trace_fraction,total_time;
void get_bulk_composition();
- void write_tracers();
static int been_here=0;
//TODO: fix this function
//if (E->composition.ichemical_buoyancy==1) get_bulk_composition(E);
- if ( ((E->monitor.solution_cycles % E->control.record_every)==0) || (been_here==0))
- {
- //TODO: migrate to output_tracer
- //write_tracers(E,1);
- }
- //TODO: move to Output.c
- if ( ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
- {
- //TODO: migrate to output_tracer
- //write_tracers(E,2);
- }
-
-
fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
if (been_here!=0)
- {
- fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
- fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
- if (E->composition.ichemical_buoyancy==1)
- {
- fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
- E->trace.istat_iempty,(100.0*E->trace.istat_iempty/E->lmesh.nel));
- }
- }
+ {
+ fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
+ fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
+ if (E->composition.ichemical_buoyancy==1)
+ {
+ fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
+ E->trace.istat_iempty,(100.0*E->trace.istat_iempty/E->lmesh.nel));
+ }
+ }
/* fprintf(E->trace.fpt,"Error fraction: %f (composition: %f)\n",E->trace.error_fraction,E->trace.bulk_composition); */
@@ -412,30 +391,30 @@
/* compositional and error fraction data files */
//TODO: move
if (E->parallel.me==0)
- {
- if (been_here==0)
- {
- if (E->composition.ichemical_buoyancy==1)
- {
- sprintf(output_file,"%s.error_fraction.data",E->control.data_file);
- E->trace.fp_error_fraction=fopen(output_file,"w");
+ {
+ if (been_here==0)
+ {
+ if (E->composition.ichemical_buoyancy==1)
+ {
+ sprintf(output_file,"%s.error_fraction.data",E->control.data_file);
+ E->trace.fp_error_fraction=fopen(output_file,"w");
- sprintf(output_file,"%s.composition.data",E->control.data_file);
- E->trace.fp_composition=fopen(output_file,"w");
- }
- }
+ sprintf(output_file,"%s.composition.data",E->control.data_file);
+ E->trace.fp_composition=fopen(output_file,"w");
+ }
+ }
- if (E->composition.ichemical_buoyancy==1)
- {
+ if (E->composition.ichemical_buoyancy==1)
+ {
//TODO: to be init'd
- fprintf(E->trace.fp_error_fraction,"%e %e\n",E->monitor.elapsed_time,E->trace.error_fraction);
- fprintf(E->trace.fp_composition,"%e %e\n",E->monitor.elapsed_time,E->trace.bulk_composition);
+ fprintf(E->trace.fp_error_fraction,"%e %e\n",E->monitor.elapsed_time,E->trace.error_fraction);
+ fprintf(E->trace.fp_composition,"%e %e\n",E->monitor.elapsed_time,E->trace.bulk_composition);
- fflush(E->trace.fp_error_fraction);
- fflush(E->trace.fp_composition);
- }
+ fflush(E->trace.fp_error_fraction);
+ fflush(E->trace.fp_composition);
+ }
- }
+ }
@@ -446,90 +425,7 @@
return;
}
-/*********** WRITE TRACERS **********************************************/
-/* */
-/* if iflag=1, rewrite over last save array */
-/* if iflag=2, write periodic save file and gzip */
-
-
-void write_tracers(E,iflag)
- struct All_variables *E;
- int iflag;
-
-{
-
- char output_file[200];
- /* char doit_string[200]; */
- FILE *fpcomp;
-
- int kk,j;
- double comp;
- double theta,phi,rad;
-
- if (iflag==1)
- {
- sprintf(output_file,"%s.tracers.%d",E->control.data_file,
- E->parallel.me);
- }
- if (iflag==2)
- {
- sprintf(output_file,"%s.tracers.%d.%d",E->control.data_file,
- E->parallel.me, E->monitor.solution_cycles);
- }
- if (iflag==3)
- {
- sprintf(output_file,"%s.tracers.%d.preback",E->control.data_file,
- E->parallel.me);
- }
- fpcomp=fopen(output_file,"w");
-
- /* This may be slightly incompatible with previous version */
-
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- //TODO: move
-/* fprintf(fpcomp,"%6d %6d %.5e %.5e %.5e %d\n",E->trace.itrac[j][0], */
-/* E->monitor.solution_cycles,E->monitor.elapsed_time, */
-/* E->trace.bulk_composition, */
-/* E->trace.initial_bulk_composition,j); */
-
- comp=0.0;
- for (kk=1;kk<=E->trace.itrac[j][0];kk++)
- {
- theta=E->trace.rtrac[j][0][kk];
- phi=E->trace.rtrac[j][1][kk];
- rad=E->trace.rtrac[j][2][kk];
-
- if (E->composition.ichemical_buoyancy==1)
- {
- comp=E->trace.etrac[j][0][kk];
- fprintf(fpcomp,"%.12e %.12e %.12e %.12e \n",theta,phi,rad,comp);
- fflush(fpcomp);
- }
- else if (E->composition.ichemical_buoyancy==0)
- {
- fprintf(fpcomp,"%.12e %.12e %.12e \n",theta,phi,rad);
- fflush(fpcomp);
- }
- }
- fclose(fpcomp);
-
- /* maybe some problem with zipping on the fly? */
- /*
- if (iflag==2)
- {
- sprintf(doit_string,"gzip %s",output_file);
- system(doit_string);
- }
- */
-
- }
-
- return;
-}
-
-
/*********** PREDICT TRACERS **********************************************/
/* */
/* This function predicts tracers performing an euler step */
@@ -572,66 +468,66 @@
/* (already did after last stokes calculation, unless this is first step) */
if ((been_here==0) && (E->trace.itracer_advection_scheme==2))
- {
- get_cartesian_velocity_field(E);
- been_here++;
- }
+ {
+ get_cartesian_velocity_field(E);
+ been_here++;
+ }
if (E->trace.itracer_advection_scheme==1) get_cartesian_velocity_field(E);
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- numtracers=E->trace.itrac[j][0];
+ numtracers=E->trace.ntracers[j];
- for (kk=1;kk<=numtracers;kk++)
- {
+ for (kk=1;kk<=numtracers;kk++)
+ {
- theta0=E->trace.rtrac[j][0][kk];
- phi0=E->trace.rtrac[j][1][kk];
- rad0=E->trace.rtrac[j][2][kk];
- x0=E->trace.rtrac[j][3][kk];
- y0=E->trace.rtrac[j][4][kk];
- z0=E->trace.rtrac[j][5][kk];
+ theta0=E->trace.basicq[j][0][kk];
+ phi0=E->trace.basicq[j][1][kk];
+ rad0=E->trace.basicq[j][2][kk];
+ x0=E->trace.basicq[j][3][kk];
+ y0=E->trace.basicq[j][4][kk];
+ z0=E->trace.basicq[j][5][kk];
- nelem=E->trace.itrac[j][kk];
- get_velocity(E,j,nelem,theta0,phi0,rad0,velocity_vector);
+ nelem=E->trace.ielement[j][kk];
+ get_velocity(E,j,nelem,theta0,phi0,rad0,velocity_vector);
- x_pred=x0+velocity_vector[1]*dt;
- y_pred=y0+velocity_vector[2]*dt;
- z_pred=z0+velocity_vector[3]*dt;
+ x_pred=x0+velocity_vector[1]*dt;
+ y_pred=y0+velocity_vector[2]*dt;
+ z_pred=z0+velocity_vector[3]*dt;
- /* keep in box */
+ /* keep in box */
- cart_to_sphere(E,x_pred,y_pred,z_pred,&theta_pred,&phi_pred,&rad_pred);
- keep_in_sphere(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
+ cart_to_sphere(E,x_pred,y_pred,z_pred,&theta_pred,&phi_pred,&rad_pred);
+ keep_in_sphere(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
- /* Current Coordinates are always kept in positions 0-5. */
+ /* Current Coordinates are always kept in positions 0-5. */
- E->trace.rtrac[j][0][kk]=theta_pred;
- E->trace.rtrac[j][1][kk]=phi_pred;
- E->trace.rtrac[j][2][kk]=rad_pred;
- E->trace.rtrac[j][3][kk]=x_pred;
- E->trace.rtrac[j][4][kk]=y_pred;
- E->trace.rtrac[j][5][kk]=z_pred;
+ E->trace.basicq[j][0][kk]=theta_pred;
+ E->trace.basicq[j][1][kk]=phi_pred;
+ E->trace.basicq[j][2][kk]=rad_pred;
+ E->trace.basicq[j][3][kk]=x_pred;
+ E->trace.basicq[j][4][kk]=y_pred;
+ E->trace.basicq[j][5][kk]=z_pred;
- /* Fill in original coords (positions 6-8) */
+ /* Fill in original coords (positions 6-8) */
- E->trace.rtrac[j][6][kk]=x0;
- E->trace.rtrac[j][7][kk]=y0;
- E->trace.rtrac[j][8][kk]=z0;
+ E->trace.basicq[j][6][kk]=x0;
+ E->trace.basicq[j][7][kk]=y0;
+ E->trace.basicq[j][8][kk]=z0;
- /* Fill in original velocities (positions 9-11) */
+ /* Fill in original velocities (positions 9-11) */
- E->trace.rtrac[j][9][kk]=velocity_vector[1]; /* Vx */
- E->trace.rtrac[j][10][kk]=velocity_vector[2]; /* Vy */
- E->trace.rtrac[j][11][kk]=velocity_vector[3]; /* Vz */
+ E->trace.basicq[j][9][kk]=velocity_vector[1]; /* Vx */
+ E->trace.basicq[j][10][kk]=velocity_vector[2]; /* Vy */
+ E->trace.basicq[j][11][kk]=velocity_vector[3]; /* Vz */
- } /* end kk, predicting tracers */
- } /* end caps */
+ } /* end kk, predicting tracers */
+ } /* end caps */
/* find new tracer elements and caps */
@@ -683,67 +579,67 @@
dt=E->advection.timestep;
if ((E->parallel.me==0) && (E->trace.itracer_advection_scheme==2) )
- {
- fprintf(stderr,"AA:Correcting Tracers\n");
- fflush(stderr);
- }
+ {
+ fprintf(stderr,"AA:Correcting Tracers\n");
+ fflush(stderr);
+ }
/* If scheme==1, use same velocity (t=0) */
/* Else if scheme==2, use new velocity (t=dt) */
if (E->trace.itracer_advection_scheme==2)
- {
- get_cartesian_velocity_field(E);
- }
+ {
+ get_cartesian_velocity_field(E);
+ }
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->trace.itrac[j][0];kk++)
- {
+ {
+ for (kk=1;kk<=E->trace.ntracers[j];kk++)
+ {
- theta_pred=E->trace.rtrac[j][0][kk];
- phi_pred=E->trace.rtrac[j][1][kk];
- rad_pred=E->trace.rtrac[j][2][kk];
- x_pred=E->trace.rtrac[j][3][kk];
- y_pred=E->trace.rtrac[j][4][kk];
- z_pred=E->trace.rtrac[j][5][kk];
+ theta_pred=E->trace.basicq[j][0][kk];
+ phi_pred=E->trace.basicq[j][1][kk];
+ rad_pred=E->trace.basicq[j][2][kk];
+ x_pred=E->trace.basicq[j][3][kk];
+ y_pred=E->trace.basicq[j][4][kk];
+ z_pred=E->trace.basicq[j][5][kk];
- x0=E->trace.rtrac[j][6][kk];
- y0=E->trace.rtrac[j][7][kk];
- z0=E->trace.rtrac[j][8][kk];
+ x0=E->trace.basicq[j][6][kk];
+ y0=E->trace.basicq[j][7][kk];
+ z0=E->trace.basicq[j][8][kk];
- Vx0=E->trace.rtrac[j][9][kk];
- Vy0=E->trace.rtrac[j][10][kk];
- Vz0=E->trace.rtrac[j][11][kk];
+ Vx0=E->trace.basicq[j][9][kk];
+ Vy0=E->trace.basicq[j][10][kk];
+ Vz0=E->trace.basicq[j][11][kk];
- nelem=E->trace.itrac[j][kk];
+ nelem=E->trace.ielement[j][kk];
- get_velocity(E,j,nelem,theta_pred,phi_pred,rad_pred,velocity_vector);
+ get_velocity(E,j,nelem,theta_pred,phi_pred,rad_pred,velocity_vector);
- Vx_pred=velocity_vector[1];
- Vy_pred=velocity_vector[2];
- Vz_pred=velocity_vector[3];
+ Vx_pred=velocity_vector[1];
+ Vy_pred=velocity_vector[2];
+ Vz_pred=velocity_vector[3];
- x_cor=x0 + dt * 0.5*(Vx0+Vx_pred);
- y_cor=y0 + dt * 0.5*(Vy0+Vy_pred);
- z_cor=z0 + dt * 0.5*(Vz0+Vz_pred);
+ x_cor=x0 + dt * 0.5*(Vx0+Vx_pred);
+ y_cor=y0 + dt * 0.5*(Vy0+Vy_pred);
+ z_cor=z0 + dt * 0.5*(Vz0+Vz_pred);
- cart_to_sphere(E,x_cor,y_cor,z_cor,&theta_cor,&phi_cor,&rad_cor);
- keep_in_sphere(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
+ cart_to_sphere(E,x_cor,y_cor,z_cor,&theta_cor,&phi_cor,&rad_cor);
+ keep_in_sphere(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
- /* Fill in Current Positions (other positions are no longer important) */
+ /* Fill in Current Positions (other positions are no longer important) */
- E->trace.rtrac[j][0][kk]=theta_cor;
- E->trace.rtrac[j][1][kk]=phi_cor;
- E->trace.rtrac[j][2][kk]=rad_cor;
- E->trace.rtrac[j][3][kk]=x_cor;
- E->trace.rtrac[j][4][kk]=y_cor;
- E->trace.rtrac[j][5][kk]=z_cor;
+ E->trace.basicq[j][0][kk]=theta_cor;
+ E->trace.basicq[j][1][kk]=phi_cor;
+ E->trace.basicq[j][2][kk]=rad_cor;
+ E->trace.basicq[j][3][kk]=x_cor;
+ E->trace.basicq[j][4][kk]=y_cor;
+ E->trace.basicq[j][5][kk]=z_cor;
- } /* end kk, correcting tracers */
- } /* end caps */
+ } /* end kk, correcting tracers */
+ } /* end caps */
/* find new tracer elements and caps */
find_tracers(E);
@@ -765,21 +661,21 @@
/* gnomonic projection */
if (E->trace.itracer_interpolation_scheme==1)
- {
- gnomonic_interpolation(E,j,nelem,theta,phi,rad,velocity_vector);
- }
+ {
+ gnomonic_interpolation(E,j,nelem,theta,phi,rad,velocity_vector);
+ }
else if (E->trace.itracer_interpolation_scheme==2)
- {
- fprintf(E->trace.fpt,"Error(get velocity)-not ready for simple average interpolation scheme\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ fprintf(E->trace.fpt,"Error(get velocity)-not ready for simple average interpolation scheme\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
else
- {
- fprintf(E->trace.fpt,"Error(get velocity)-not ready for other interpolation schemes\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ fprintf(E->trace.fpt,"Error(get velocity)-not ready for other interpolation schemes\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
return;
}
@@ -874,47 +770,47 @@
/* if any shape functions are negative, goto next wedge */
if (shape2d[1]<eps||shape2d[2]<eps||shape2d[3]<eps)
- {
- inum=inum+1;
- /* AKMA clean this up */
- if (inum>3)
- {
- fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>3!\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if (inum>1 && itry==1)
- {
- fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>1\n");
- fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[1],shape2d[2],shape2d[3]);
- fprintf(E->trace.fpt,"u %f v %f element: %d \n",u,v, nelem);
- fprintf(E->trace.fpt,"Element uv boundaries: \n");
- for(kk=1;kk<=4;kk++)
- fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->trace.UV[j][1][E->ien[j][nelem].node[kk]],E->trace.UV[j][2][E->ien[j][nelem].node[kk]]);
- fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",theta,phi,rad);
- fprintf(E->trace.fpt,"Element theta-phi boundaries: \n");
- for(kk=1;kk<=4;kk++)
- fprintf(E->trace.fpt,"%d: Theta: %f Phi:%f\n",kk,E->sx[j][1][E->ien[j][nelem].node[kk]],E->sx[j][2][E->ien[j][nelem].node[kk]]);
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
- ival=icheck_element(E,j,nelem,x,y,z,rad);
- fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
- ival=iget_element(E,j,-99,x,y,z,theta,phi,rad);
- fprintf(E->trace.fpt,"New Element?: %d\n",ival);
- ival=icheck_column_neighbors(E,j,nelem,x,y,z,rad);
- fprintf(E->trace.fpt,"New Element (neighs)?: %d\n",ival);
- nelem=ival;
- ival=icheck_element(E,j,nelem,x,y,z,rad);
- fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
- itry++;
- if (ival>0) goto try_again;
- fprintf(E->trace.fpt,"NO LUCK\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ inum=inum+1;
+ /* AKMA clean this up */
+ if (inum>3)
+ {
+ fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>3!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if (inum>1 && itry==1)
+ {
+ fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>1\n");
+ fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[1],shape2d[2],shape2d[3]);
+ fprintf(E->trace.fpt,"u %f v %f element: %d \n",u,v, nelem);
+ fprintf(E->trace.fpt,"Element uv boundaries: \n");
+ for(kk=1;kk<=4;kk++)
+ fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->trace.UV[j][1][E->ien[j][nelem].node[kk]],E->trace.UV[j][2][E->ien[j][nelem].node[kk]]);
+ fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",theta,phi,rad);
+ fprintf(E->trace.fpt,"Element theta-phi boundaries: \n");
+ for(kk=1;kk<=4;kk++)
+ fprintf(E->trace.fpt,"%d: Theta: %f Phi:%f\n",kk,E->sx[j][1][E->ien[j][nelem].node[kk]],E->sx[j][2][E->ien[j][nelem].node[kk]]);
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+ ival=icheck_element(E,j,nelem,x,y,z,rad);
+ fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
+ ival=iget_element(E,j,-99,x,y,z,theta,phi,rad);
+ fprintf(E->trace.fpt,"New Element?: %d\n",ival);
+ ival=icheck_column_neighbors(E,j,nelem,x,y,z,rad);
+ fprintf(E->trace.fpt,"New Element (neighs)?: %d\n",ival);
+ nelem=ival;
+ ival=icheck_element(E,j,nelem,x,y,z,rad);
+ fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
+ itry++;
+ if (ival>0) goto try_again;
+ fprintf(E->trace.fpt,"NO LUCK\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- iwedge=2;
- goto next_wedge;
- }
+ iwedge=2;
+ goto next_wedge;
+ }
/* Determine radial shape functions */
/* There are 2 shape functions radially */
@@ -937,55 +833,55 @@
/* depending on wedge, set up velocity points */
if (iwedge==1)
- {
- vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
- vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
- vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[2]];
- vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
- vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
- vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[6]];
- vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
- vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
- vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[2]];
- vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
- vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
- vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[6]];
- vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
- vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
- vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[2]];
- vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
- vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
- vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[6]];
- vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
- }
+ {
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[2]];
+ vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
+ vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
+ vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[6]];
+ vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
+ vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
+ vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[2]];
+ vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
+ vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
+ vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[6]];
+ vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
+ vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
+ vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[2]];
+ vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
+ vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
+ vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[6]];
+ vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
+ }
if (iwedge==2)
- {
- vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
- vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
- vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[4]];
- vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
- vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
- vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[8]];
- vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
- vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
- vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[4]];
- vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
- vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
- vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[8]];
- vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
- vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
- vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[4]];
- vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
- vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
- vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[8]];
- }
+ {
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
+ vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[4]];
+ vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
+ vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
+ vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[8]];
+ vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
+ vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
+ vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[4]];
+ vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
+ vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
+ vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[8]];
+ vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
+ vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
+ vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[4]];
+ vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
+ vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
+ vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[8]];
+ }
velocity_vector[1]=vx[1]*shape[1]+vx[2]*shape[2]+shape[3]*vx[3]+
- vx[4]*shape[4]+vx[5]*shape[5]+shape[6]*vx[6];
+ vx[4]*shape[4]+vx[5]*shape[5]+shape[6]*vx[6];
velocity_vector[2]=vy[1]*shape[1]+vy[2]*shape[2]+shape[3]*vy[3]+
- vy[4]*shape[4]+vy[5]*shape[5]+shape[6]*vy[6];
+ vy[4]*shape[4]+vy[5]*shape[5]+shape[6]*vy[6];
velocity_vector[3]=vz[1]*shape[1]+vz[2]*shape[2]+shape[3]*vz[3]+
- vz[4]*shape[4]+vz[5]*shape[5]+shape[6]*vz[6];
+ vz[4]*shape[4]+vz[5]*shape[5]+shape[6]*vz[6];
@@ -1082,14 +978,14 @@
/* Some error control */
if (shaperad[1]>(top_bound)||shaperad[1]<(bottom_bound)||
- shaperad[2]>(top_bound)||shaperad[2]<(bottom_bound))
- {
- fprintf(E->trace.fpt,"ERROR(get_radial_shape)\n");
- fprintf(E->trace.fpt,"shaperad[1]: %f \n",shaperad[1]);
- fprintf(E->trace.fpt,"shaperad[2]: %f \n",shaperad[2]);
- fflush(E->trace.fpt);
- exit(10);
- }
+ shaperad[2]>(top_bound)||shaperad[2]<(bottom_bound))
+ {
+ fprintf(E->trace.fpt,"ERROR(get_radial_shape)\n");
+ fprintf(E->trace.fpt,"shaperad[1]: %f \n",shaperad[1]);
+ fprintf(E->trace.fpt,"shaperad[2]: %f \n",shaperad[2]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
return;
}
@@ -1157,43 +1053,38 @@
int kk;
- /* itracsize is physical size of tracer array */
+ /* max_ntracers is physical size of tracer array */
/* (initially make it 25% larger than required */
- E->trace.itracsize[j]=number_of_tracers+number_of_tracers/4;
+ E->trace.max_ntracers[j]=number_of_tracers+number_of_tracers/4;
+ E->trace.ntracers[j]=0;
/* make tracer arrays */
- if ((E->trace.itrac[j]=(int *) malloc(E->trace.itracsize[j]*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make tracer array)-no memory 1a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- E->trace.itrac[j][0]=0;
+ if ((E->trace.ielement[j]=(int *) malloc(E->trace.max_ntracers[j]*sizeof(int)))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(make tracer array)-no memory 1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
- {
- if ((E->trace.rtrac[j][kk]=(double *)malloc(E->trace.itracsize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1b.%d\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (kk=0;kk<E->trace.number_of_basic_quantities;kk++) {
+ if ((E->trace.basicq[j][kk]=(double *)malloc(E->trace.max_ntracers[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1b.%d\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- {
- if ((E->trace.etrac[j][kk]=(double *)malloc(E->trace.itracsize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1c.%d\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (kk=0;kk<E->trace.number_of_extra_quantities;kk++) {
+ if ((E->trace.extraq[j][kk]=(double *)malloc(E->trace.max_ntracers[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1c.%d\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- fprintf(E->trace.fpt,"Physical size of tracer arrays (itracsize): %d\n",
- E->trace.itracsize[j]);
+ fprintf(E->trace.fpt,"Physical size of tracer arrays (max_ntracers): %d\n",
+ E->trace.max_ntracers[j]);
fflush(E->trace.fpt);
return;
@@ -1205,7 +1096,7 @@
/* */
/* This function finds tracer elements and moves tracers to */
/* other processor domains if necessary. */
-/* Array itrac is filled with elemental values. */
+/* Array ielement is filled with elemental values. */
void find_tracers(E)
struct All_variables *E;
@@ -1237,77 +1128,74 @@
if (been_here==0)
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->trace.itrac[j][0];kk++)
- {
- E->trace.itrac[j][kk]=-99;
- }
- }
- been_here++;
- }
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->trace.ntracers[j];kk++)
+ {
+ E->trace.ielement[j][kk]=-99;
+ }
+ }
+ been_here++;
+ }
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* initialize arrays and statistical counters */
+ /* initialize arrays and statistical counters */
- E->trace.ilater[j]=0;
+ E->trace.ilater[j]=0;
- E->trace.istat1=0;
- for (kk=0;kk<=4;kk++)
- {
- E->trace.istat_ichoice[j][kk]=0;
- }
+ E->trace.istat1=0;
+ for (kk=0;kk<=4;kk++)
+ {
+ E->trace.istat_ichoice[j][kk]=0;
+ }
- //TODO: use while-loop instead of for-loop
- /* important to index by it, not kk */
+ //TODO: use while-loop instead of for-loop
+ /* important to index by it, not kk */
- it=0;
- num_tracers=E->trace.itrac[j][0];
+ it=0;
+ num_tracers=E->trace.ntracers[j];
- for (kk=1;kk<=num_tracers;kk++)
- {
+ for (kk=1;kk<=num_tracers;kk++)
+ {
- it++;
+ it++;
- theta=E->trace.rtrac[j][0][it];
- phi=E->trace.rtrac[j][1][it];
- rad=E->trace.rtrac[j][2][it];
- x=E->trace.rtrac[j][3][it];
- y=E->trace.rtrac[j][4][it];
- z=E->trace.rtrac[j][5][it];
+ theta=E->trace.basicq[j][0][it];
+ phi=E->trace.basicq[j][1][it];
+ rad=E->trace.basicq[j][2][it];
+ x=E->trace.basicq[j][3][it];
+ y=E->trace.basicq[j][4][it];
+ z=E->trace.basicq[j][5][it];
- iprevious_element=E->trace.itrac[j][it];
+ iprevious_element=E->trace.ielement[j][it];
- /* AKMA REMOVE */
- /*
- fprintf(E->trace.fpt,"BB. kk %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,x,y,z,theta,phi,rad);
- fflush(E->trace.fpt);
- */
- iel=iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad);
+ /* AKMA REMOVE */
+ /*
+ fprintf(E->trace.fpt,"BB. kk %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,x,y,z,theta,phi,rad);
+ fflush(E->trace.fpt);
+ */
+ iel=iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad);
- E->trace.itrac[j][it]=iel;
+ E->trace.ielement[j][it]=iel;
- if (iel<0)
- {
- put_away_later(E,j,it);
- eject_tracer(E,j,it);
- it--;
- }
+ if (iel<0)
+ {
+ put_away_later(E,j,it);
+ eject_tracer(E,j,it);
+ it--;
+ }
- } /* end tracers */
+ } /* end tracers */
- } /* end j */
+ } /* end j */
- /* fprintf(E->trace.fpt,"tracer# old:%d new:%d\n", */
- /* num_tracers, E->trace.itrac[1][0]); */
-
/* Now take care of tracers that exited cap */
/* REMOVE */
@@ -1320,15 +1208,15 @@
/* Free later arrays */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- if (E->trace.ilater[j]>0)
- {
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
- {
- free(E->trace.rlater[j][kk]);
- }
- }
- } /* end j */
+ {
+ if (E->trace.ilater[j]>0)
+ {
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
+ {
+ free(E->trace.rlater[j][kk]);
+ }
+ }
+ } /* end j */
/* Adjust Array Sizes */
@@ -1419,48 +1307,48 @@
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- E->trace.istat_isend=E->trace.ilater[j];
- }
+ {
+ E->trace.istat_isend=E->trace.ilater[j];
+ }
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- for (kk=1; kk<=E->trace.istat_isend; kk++) {
- fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
- E->trace.rlater[j][0][kk],
- E->trace.rlater[j][1][kk],
- E->trace.rlater[j][2][kk]);
- }
- fflush(E->trace.fpt);
+ for (kk=1; kk<=E->trace.istat_isend; kk++) {
+ fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
+ E->trace.rlater[j][0][kk],
+ E->trace.rlater[j][1][kk],
+ E->trace.rlater[j][2][kk]);
+ }
+ fflush(E->trace.fpt);
}
/* initialize isend and ireceive */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
/* # of neighbors in the horizontal plane */
num_ngb = E->parallel.TNUM_PASS[lev][j];
- isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
- for (kk=1;kk<=num_ngb;kk++) isend[j][kk]=0;
- for (kk=1;kk<=num_ngb;kk++) ireceive[j][kk]=0;
- }
+ isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
+ for (kk=1;kk<=num_ngb;kk++) isend[j][kk]=0;
+ for (kk=1;kk<=num_ngb;kk++) ireceive[j][kk]=0;
+ }
/* Allocate Maximum Memory to Send Arrays */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- ithiscap=0;
+ ithiscap=0;
- itemp_size=max(isize[j],1);
+ itemp_size=max(isize[j],1);
- if ((send[j][ithiscap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((send[j][ithiscap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
num_ngb = E->parallel.TNUM_PASS[lev][j];
- for (kk=1;kk<=num_ngb;kk++)
+ for (kk=1;kk<=num_ngb;kk++)
{
if ((send[j][kk]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
{
@@ -1469,7 +1357,7 @@
exit(10);
}
}
- }
+ }
/* For testing, remove this */
@@ -1479,9 +1367,9 @@
ithiscap=E->sphere.capid[j];
for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
{
- ithatcap=E->parallel.PROCESSOR[lev][j].pass[kk]+1;
- fprintf(E->trace.fpt,"cap: %d proc %d TNUM: %d ithatcap: %d\n",
- ithiscap,E->parallel.me,kk,ithatcap);
+ ithatcap=E->parallel.PROCESSOR[lev][j].pass[kk]+1;
+ fprintf(E->trace.fpt,"cap: %d proc %d TNUM: %d ithatcap: %d\n",
+ ithiscap,E->parallel.me,kk,ithatcap);
}
fflush(E->trace.fpt);
@@ -1492,68 +1380,68 @@
/* Pre communication */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* transfer tracers from rlater to send */
+ /* transfer tracers from rlater to send */
- numtracers=E->trace.ilater[j];
+ numtracers=E->trace.ilater[j];
- for (kk=1;kk<=numtracers;kk++)
- {
- rad=E->trace.rlater[j][2][kk];
- x=E->trace.rlater[j][3][kk];
- y=E->trace.rlater[j][4][kk];
- z=E->trace.rlater[j][5][kk];
+ for (kk=1;kk<=numtracers;kk++)
+ {
+ rad=E->trace.rlater[j][2][kk];
+ x=E->trace.rlater[j][3][kk];
+ y=E->trace.rlater[j][4][kk];
+ z=E->trace.rlater[j][5][kk];
- /* first check same cap if nprocz>1 */
+ /* first check same cap if nprocz>1 */
- if (E->parallel.nprocz>1)
- {
- ithatcap=0;
- icheck=icheck_cap(E,ithatcap,x,y,z,rad);
- if (icheck==1) goto foundit;
+ if (E->parallel.nprocz>1)
+ {
+ ithatcap=0;
+ icheck=icheck_cap(E,ithatcap,x,y,z,rad);
+ if (icheck==1) goto foundit;
- }
+ }
- /* check neighboring caps */
+ /* check neighboring caps */
- for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++)
- {
- ithatcap=pp;
- icheck=icheck_cap(E,ithatcap,x,y,z,rad);
- if (icheck==1) goto foundit;
- }
+ for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++)
+ {
+ ithatcap=pp;
+ icheck=icheck_cap(E,ithatcap,x,y,z,rad);
+ if (icheck==1) goto foundit;
+ }
- /* should not be here */
- if (icheck!=1)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
- fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
- icheck=icheck_cap(E,0,x,y,z,rad);
- if (icheck==1) fprintf(E->trace.fpt," icheck here!\n");
- else fprintf(E->trace.fpt,"icheck not here!\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ /* should not be here */
+ if (icheck!=1)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
+ fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
+ icheck=icheck_cap(E,0,x,y,z,rad);
+ if (icheck==1) fprintf(E->trace.fpt," icheck here!\n");
+ else fprintf(E->trace.fpt,"icheck not here!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- foundit:
+ foundit:
- isend[j][ithatcap]++;
+ isend[j][ithatcap]++;
- /* assign tracer to send */
+ /* assign tracer to send */
- isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
+ isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
- for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++)
- {
- ipos=isend_position+pp;
- send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
- }
+ for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++)
+ {
+ ipos=isend_position+pp;
+ send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
+ }
- } /* end kk, assigning tracers */
+ } /* end kk, assigning tracers */
- } /* end j */
+ } /* end j */
/* Send info to other processors regarding number of send tracers */
@@ -1563,58 +1451,58 @@
idb=0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- ithiscap=0;
+ ithiscap=0;
- /* if tracer is in same cap (nprocz>1) */
+ /* if tracer is in same cap (nprocz>1) */
- if (E->parallel.nprocz>1)
- {
- ireceive[j][ithiscap]=isend[j][ithiscap];
- }
+ if (E->parallel.nprocz>1)
+ {
+ ireceive[j][ithiscap]=isend[j][ithiscap];
+ }
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
- /* if neighbor cap is in another processor, send information via MPI */
+ /* if neighbor cap is in another processor, send information via MPI */
- idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
idb++;
MPI_Isend(&isend[j][ithatcap],1,MPI_INT,idestination_proc,
11,E->parallel.world,
&request[idb-1]);
- } /* end kk, number of neighbors */
+ } /* end kk, number of neighbors */
- } /* end j, caps per proc */
+ } /* end j, caps per proc */
/* Receive tracer count info */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
+ {
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
- /* if neighbor cap is in another processor, receive information via MPI */
+ /* if neighbor cap is in another processor, receive information via MPI */
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- if (idestination_proc!=E->parallel.me)
- {
+ if (idestination_proc!=E->parallel.me)
+ {
- idb++;
- MPI_Irecv(&ireceive[j][ithatcap],1,MPI_INT,isource_proc,
- 11,E->parallel.world,
- &request[idb-1]);
+ idb++;
+ MPI_Irecv(&ireceive[j][ithatcap],1,MPI_INT,isource_proc,
+ 11,E->parallel.world,
+ &request[idb-1]);
- } /* end if */
+ } /* end if */
- } /* end kk, number of neighbors */
- } /* end j, caps per proc */
+ } /* end kk, number of neighbors */
+ } /* end j, caps per proc */
/* Wait for non-blocking calls to complete */
@@ -1636,53 +1524,53 @@
/* Allocate memory in receive arrays */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
num_ngb = E->parallel.TNUM_PASS[lev][j];
- for (ithatcap=1;ithatcap<=num_ngb;ithatcap++)
- {
- isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ for (ithatcap=1;ithatcap<=num_ngb;ithatcap++)
+ {
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
- itemp_size=max(1,isize[j]);
+ itemp_size=max(1,isize[j]);
- if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
+ if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
/* Now, send the tracers to proper caps */
idb=0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- ithiscap=0;
+ {
+ ithiscap=0;
- /* same cap */
+ /* same cap */
- if (E->parallel.nprocz>1)
- {
+ if (E->parallel.nprocz>1)
+ {
- ithatcap=ithiscap;
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(isize[j]-1);mm++)
- {
- receive[j][ithatcap][mm]=send[j][ithatcap][mm];
- }
+ ithatcap=ithiscap;
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ for (mm=0;mm<=(isize[j]-1);mm++)
+ {
+ receive[j][ithatcap][mm]=send[j][ithatcap][mm];
+ }
- }
+ }
- /* neighbor caps */
+ /* neighbor caps */
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
- idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
idb++;
@@ -1690,22 +1578,22 @@
11,E->parallel.world,
&request[idb-1]);
- } /* end kk, number of neighbors */
+ } /* end kk, number of neighbors */
- } /* end j, caps per proc */
+ } /* end j, caps per proc */
/* Receive tracers */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- ithiscap=0;
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
+ ithiscap=0;
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
idb++;
@@ -1715,9 +1603,9 @@
11,E->parallel.world,
&request[idb-1]);
- } /* end kk, number of neighbors */
+ } /* end kk, number of neighbors */
- } /* end j, caps per proc */
+ } /* end j, caps per proc */
/* Wait for non-blocking calls to complete */
@@ -1730,96 +1618,96 @@
/* Sum up size of receive arrays (all tracers sent to this processor) */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isum[j]=0;
+ {
+ isum[j]=0;
- ithiscap=0;
+ ithiscap=0;
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
- isum[j]=isum[j]+ireceive[j][ithatcap];
- }
- if (E->parallel.nprocz>1) isum[j]=isum[j]+ireceive[j][ithiscap];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
+ isum[j]=isum[j]+ireceive[j][ithatcap];
+ }
+ if (E->parallel.nprocz>1) isum[j]=isum[j]+ireceive[j][ithiscap];
- itracers_subject_to_vertical_transport[j]=isum[j];
- }
+ itracers_subject_to_vertical_transport[j]=isum[j];
+ }
/* Allocate Memory for REC array */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
- if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- REC[j][0]=0.0;
- }
+ {
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
+ if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ REC[j][0]=0.0;
+ }
/* Put Received tracers in REC */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- irec[j]=0;
+ irec[j]=0;
- irec_position=0;
+ irec_position=0;
- ithiscap=0;
+ ithiscap=0;
- /* horizontal neighbors */
+ /* horizontal neighbors */
- for (ihorizontal_neighbor=1;ihorizontal_neighbor<=E->parallel.TNUM_PASS[lev][j];ihorizontal_neighbor++)
- {
+ for (ihorizontal_neighbor=1;ihorizontal_neighbor<=E->parallel.TNUM_PASS[lev][j];ihorizontal_neighbor++)
+ {
- ithatcap=ihorizontal_neighbor;
+ ithatcap=ihorizontal_neighbor;
- for (pp=1;pp<=ireceive[j][ithatcap];pp++)
- {
- irec[j]++;
- ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++)
+ {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos2=ipos+mm;
- REC[j][irec_position]=receive[j][ithatcap][ipos2];
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ ipos2=ipos+mm;
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
- irec_position++;
+ irec_position++;
- } /* end mm (cycling tracer quantities) */
- } /* end pp (cycling tracers) */
- } /* end ihorizontal_neighbors (cycling neighbors) */
+ } /* end mm (cycling tracer quantities) */
+ } /* end pp (cycling tracers) */
+ } /* end ihorizontal_neighbors (cycling neighbors) */
- /* for tracers in the same cap (nprocz>1) */
+ /* for tracers in the same cap (nprocz>1) */
- if (E->parallel.nprocz>1)
- {
- ithatcap=ithiscap;
- for (pp=1;pp<=ireceive[j][ithatcap];pp++)
- {
- irec[j]++;
- ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+ if (E->parallel.nprocz>1)
+ {
+ ithatcap=ithiscap;
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++)
+ {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos2=ipos+mm;
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ ipos2=ipos+mm;
- REC[j][irec_position]=receive[j][ithatcap][ipos2];
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
- irec_position++;
+ irec_position++;
- } /* end mm (cycling tracer quantities) */
+ } /* end mm (cycling tracer quantities) */
- } /* end pp (cycling tracers) */
+ } /* end pp (cycling tracers) */
- } /* endif nproc>1 */
+ } /* endif nproc>1 */
- } /* end j (cycling caps) */
+ } /* end j (cycling caps) */
/* Done filling REC */
@@ -1834,228 +1722,228 @@
/* CitcomS code. */
if (E->parallel.nprocz>1)
- {
+ {
- /* Allocate memory for send_z */
- /* Make send_z the size of receive array (max size) */
- /* (No dynamic reallocation of send_z necessary) */
+ /* Allocate memory for send_z */
+ /* Make send_z the size of receive array (max size) */
+ /* (No dynamic reallocation of send_z necessary) */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
- {
- isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
+ {
+ isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
- if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
+ if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
- ithat_processor=E->parallel.PROCESSORz[lev].pass[ivertical_neighbor];
+ ithat_processor=E->parallel.PROCESSORz[lev].pass[ivertical_neighbor];
- /* initialize isend_z and ireceive_z array */
+ /* initialize isend_z and ireceive_z array */
- isend_z[j][ivertical_neighbor]=0;
- ireceive_z[j][ivertical_neighbor]=0;
+ isend_z[j][ivertical_neighbor]=0;
+ ireceive_z[j][ivertical_neighbor]=0;
- /* sort through receive array and check radius */
+ /* sort through receive array and check radius */
- it=0;
- num_tracers=irec[j];
- for (kk=1;kk<=num_tracers;kk++)
- {
+ it=0;
+ num_tracers=irec[j];
+ for (kk=1;kk<=num_tracers;kk++)
+ {
- it++;
+ it++;
- ireceive_position=((it-1)*E->trace.number_of_tracer_quantities);
+ ireceive_position=((it-1)*E->trace.number_of_tracer_quantities);
- irad=ireceive_position+2;
+ irad=ireceive_position+2;
- rad=REC[j][irad];
+ rad=REC[j][irad];
- ival=icheck_that_processor_shell(E,j,ithat_processor,rad);
+ ival=icheck_that_processor_shell(E,j,ithat_processor,rad);
- /* if tracer is in other shell, take out of receive array and give to send_z*/
+ /* if tracer is in other shell, take out of receive array and give to send_z*/
- if (ival==1)
- {
+ if (ival==1)
+ {
- isend_z[j][ivertical_neighbor]++;
+ isend_z[j][ivertical_neighbor]++;
- isend_position=(isend_z[j][ivertical_neighbor]-1)*E->trace.number_of_tracer_quantities;
+ isend_position=(isend_z[j][ivertical_neighbor]-1)*E->trace.number_of_tracer_quantities;
- ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos=ireceive_position+mm;
- ipos2=isend_position+mm;
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ ipos=ireceive_position+mm;
+ ipos2=isend_position+mm;
- send_z[j][ivertical_neighbor][ipos2]=REC[j][ipos];
+ send_z[j][ivertical_neighbor][ipos2]=REC[j][ipos];
- /* eject tracer info from REC array, and replace with last tracer in array */
+ /* eject tracer info from REC array, and replace with last tracer in array */
- ipos3=ilast_receiver_position+mm;
- REC[j][ipos]=REC[j][ipos3];
+ ipos3=ilast_receiver_position+mm;
+ REC[j][ipos]=REC[j][ipos3];
- }
+ }
- it--;
- irec[j]--;
+ it--;
+ irec[j]--;
- } /* end if ival===1 */
+ } /* end if ival===1 */
- /* Otherwise, leave tracer */
+ /* Otherwise, leave tracer */
- } /* end kk (cycling through tracers) */
+ } /* end kk (cycling through tracers) */
- } /* end ivertical_neighbor */
+ } /* end ivertical_neighbor */
- } /* end j */
+ } /* end j */
- /* Send arrays are now filled. */
- /* Now send send information to vertical processor neighbor */
+ /* Send arrays are now filled. */
+ /* Now send send information to vertical processor neighbor */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
- MPI_Sendrecv(&isend_z[j][ivertical_neighbor],1,MPI_INT,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag,
- &ireceive_z[j][ivertical_neighbor],1,MPI_INT,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- itag,E->parallel.world,&status1);
+ MPI_Sendrecv(&isend_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag,
+ &ireceive_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag,E->parallel.world,&status1);
- /* for testing - remove */
- /*
- fprintf(E->trace.fpt,"PROC: %d IVN: %d (P: %d) SEND: %d REC: %d\n",
- E->parallel.me,ivertical_neighbor,E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- isend_z[j][ivertical_neighbor],ireceive_z[j][ivertical_neighbor]);
- fflush(E->trace.fpt);
- */
+ /* for testing - remove */
+ /*
+ fprintf(E->trace.fpt,"PROC: %d IVN: %d (P: %d) SEND: %d REC: %d\n",
+ E->parallel.me,ivertical_neighbor,E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ isend_z[j][ivertical_neighbor],ireceive_z[j][ivertical_neighbor]);
+ fflush(E->trace.fpt);
+ */
- } /* end ivertical_neighbor */
+ } /* end ivertical_neighbor */
- } /* end j */
+ } /* end j */
- /* Allocate memory to receive_z arrays */
+ /* Allocate memory to receive_z arrays */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
- {
- isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
+ {
+ isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
- if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
+ if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
- /* Send Tracers */
+ /* Send Tracers */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- isize_send=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
- isize_receive=ireceive_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+ isize_send=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ isize_receive=ireceive_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
- MPI_Sendrecv(send_z[j][ivertical_neighbor],isize_send,
- MPI_DOUBLE,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag+1,
- receive_z[j][ivertical_neighbor],isize_receive,
- MPI_DOUBLE,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- itag+1,E->parallel.world,&status2);
+ MPI_Sendrecv(send_z[j][ivertical_neighbor],isize_send,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag+1,
+ receive_z[j][ivertical_neighbor],isize_receive,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag+1,E->parallel.world,&status2);
- }
- }
+ }
+ }
- /* Put tracers into REC array */
+ /* Put tracers into REC array */
- /* First, reallocate memory to REC */
+ /* First, reallocate memory to REC */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isum[j]=0;
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- isum[j]=isum[j]+ireceive_z[j][ivertical_neighbor];
- }
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ isum[j]=0;
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+ isum[j]=isum[j]+ireceive_z[j][ivertical_neighbor];
+ }
- isum[j]=isum[j]+irec[j];
+ isum[j]=isum[j]+irec[j];
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
- if (isize[j]>0)
- {
- if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (i981)\n");
- fprintf(E->trace.fpt,"isize: %d\n",isize[j]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ if (isize[j]>0)
+ {
+ if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (i981)\n");
+ fprintf(E->trace.fpt,"isize: %d\n",isize[j]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
- for (kk=1;kk<=ireceive_z[j][ivertical_neighbor];kk++)
- {
- irec[j]++;
+ for (kk=1;kk<=ireceive_z[j][ivertical_neighbor];kk++)
+ {
+ irec[j]++;
- irec_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
- ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+ irec_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
- }
- }
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
+ }
+ }
- }
- }
+ }
+ }
- /* Free Vertical Arrays */
+ /* Free Vertical Arrays */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- free(send_z[j][ivertical_neighbor]);
- free(receive_z[j][ivertical_neighbor]);
- }
- }
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+ free(send_z[j][ivertical_neighbor]);
+ free(receive_z[j][ivertical_neighbor]);
+ }
+ }
- } /* endif nprocz>1 */
+ } /* endif nprocz>1 */
/* END OF VERTICAL TRANSPORT */
@@ -2063,50 +1951,50 @@
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=irec[j];kk++)
- {
- E->trace.itrac[j][0]++;
+ {
+ for (kk=1;kk<=irec[j];kk++)
+ {
+ E->trace.ntracers[j]++;
- if (E->trace.itrac[j][0]>(E->trace.itracsize[j]-5)) expand_tracer_arrays(E,j);
+ if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
- ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(E->trace.number_of_advection_quantities-1);mm++)
- {
- ipos=ireceive_position+mm;
+ for (mm=0;mm<=(E->trace.number_of_basic_quantities-1);mm++)
+ {
+ ipos=ireceive_position+mm;
- E->trace.rtrac[j][mm][E->trace.itrac[j][0]]=REC[j][ipos];
- }
- for (mm=0;mm<=(E->trace.number_of_extra_quantities-1);mm++)
- {
- ipos=ireceive_position+E->trace.number_of_advection_quantities+mm;
+ E->trace.basicq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
+ }
+ for (mm=0;mm<=(E->trace.number_of_extra_quantities-1);mm++)
+ {
+ ipos=ireceive_position+E->trace.number_of_basic_quantities+mm;
- E->trace.etrac[j][mm][E->trace.itrac[j][0]]=REC[j][ipos];
- }
+ E->trace.extraq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
+ }
- theta=E->trace.rtrac[j][0][E->trace.itrac[j][0]];
- phi=E->trace.rtrac[j][1][E->trace.itrac[j][0]];
- rad=E->trace.rtrac[j][2][E->trace.itrac[j][0]];
- x=E->trace.rtrac[j][3][E->trace.itrac[j][0]];
- y=E->trace.rtrac[j][4][E->trace.itrac[j][0]];
- z=E->trace.rtrac[j][5][E->trace.itrac[j][0]];
+ theta=E->trace.basicq[j][0][E->trace.ntracers[j]];
+ phi=E->trace.basicq[j][1][E->trace.ntracers[j]];
+ rad=E->trace.basicq[j][2][E->trace.ntracers[j]];
+ x=E->trace.basicq[j][3][E->trace.ntracers[j]];
+ y=E->trace.basicq[j][4][E->trace.ntracers[j]];
+ z=E->trace.basicq[j][5][E->trace.ntracers[j]];
- iel=iget_element(E,j,-99,x,y,z,theta,phi,rad);
+ iel=iget_element(E,j,-99,x,y,z,theta,phi,rad);
- if (iel<1)
- {
- fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
- fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",x,y,z,theta,phi,rad);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (iel<1)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
+ fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",x,y,z,theta,phi,rad);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- E->trace.itrac[j][E->trace.itrac[j][0]]=iel;
+ E->trace.ielement[j][E->trace.ntracers[j]]=iel;
- }
- }
+ }
+ }
fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
fflush(E->trace.fpt);
@@ -2115,24 +2003,24 @@
/* Free Arrays */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- ithiscap=0;
+ ithiscap=0;
- free(REC[j]);
+ free(REC[j]);
- free(send[j][ithiscap]);
+ free(send[j][ithiscap]);
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=kk;
- free(send[j][ithatcap]);
- free(receive[j][ithatcap]);
+ free(send[j][ithatcap]);
+ free(receive[j][ithatcap]);
- }
+ }
- }
+ }
fprintf(E->trace.fpt,"Leaving lost_souls()\n");
fflush(E->trace.fpt);
@@ -2155,71 +2043,71 @@
int icushion=100;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* if physical size is double tracer size, reduce it */
+ /* if physical size is double tracer size, reduce it */
- iempty_space=(E->trace.itracsize[j]-E->trace.itrac[j][0]);
+ iempty_space=(E->trace.max_ntracers[j]-E->trace.ntracers[j]);
- if (iempty_space>(E->trace.itrac[j][0]+icushion))
- {
+ if (iempty_space>(E->trace.ntracers[j]+icushion))
+ {
- inewsize=E->trace.itrac[j][0]+E->trace.itrac[j][0]/4+icushion;
+ inewsize=E->trace.ntracers[j]+E->trace.ntracers[j]/4+icushion;
- if (inewsize<1)
- {
- fprintf(E->trace.fpt,"Error(reduce tracer arrays)-something up (hdf3)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (inewsize<1)
+ {
+ fprintf(E->trace.fpt,"Error(reduce tracer arrays)-something up (hdf3)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- if ((E->trace.itrac[j]=(int *)realloc(E->trace.itrac[j],inewsize*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(reduce tracer arrays )-no memory (itracer)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((E->trace.ielement[j]=(int *)realloc(E->trace.ielement[j],inewsize*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(reduce tracer arrays )-no memory (ielement)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
- {
- if ((E->trace.rtrac[j][kk]=(double *)realloc(E->trace.rtrac[j][kk],inewsize*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
+ {
+ if ((E->trace.basicq[j][kk]=(double *)realloc(E->trace.basicq[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- {
- if ((E->trace.etrac[j][kk]=(double *)realloc(E->trace.etrac[j][kk],inewsize*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory 783 (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ if ((E->trace.extraq[j][kk]=(double *)realloc(E->trace.extraq[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory 783 (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- fprintf(E->trace.fpt,"Reducing physical memory of itrac, rtrac, and etrac to %d from %d\n",
- E->trace.itracsize[j],inewsize);
+ fprintf(E->trace.fpt,"Reducing physical memory of ielement, basicq, and extraq to %d from %d\n",
+ E->trace.max_ntracers[j],inewsize);
- E->trace.itracsize[j]=inewsize;
+ E->trace.max_ntracers[j]=inewsize;
- } /* end if */
+ } /* end if */
- } /* end j */
+ } /* end j */
return;
}
/********** PUT AWAY LATER ************************************/
/* */
-/* rlater has a similar structure to rtrac */
+/* rlater has a similar structure to basicq */
/* ilatersize is the physical memory and */
/* ilater is the number of tracers */
@@ -2241,20 +2129,20 @@
/* Memory is freed after parallel communications */
if (E->trace.ilater[j]==0)
- {
+ {
- E->trace.ilatersize[j]=E->trace.itracsize[j]/5;
+ E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
- {
- if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(put_away_later)-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end first particle initiating memory allocation */
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
+ {
+ if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(put_away_later)-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end first particle initiating memory allocation */
/* Put tracer in later array */
@@ -2263,16 +2151,16 @@
if (E->trace.ilater[j] >= (E->trace.ilatersize[j]-5)) expand_later_array(E,j);
- /* stack advection and extra quantities together (advection first) */
+ /* stack basic and extra quantities together (basic first) */
- for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
- {
- E->trace.rlater[j][kk][E->trace.ilater[j]]=E->trace.rtrac[j][kk][it];
- }
+ for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
+ {
+ E->trace.rlater[j][kk][E->trace.ilater[j]]=E->trace.basicq[j][kk][it];
+ }
for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- {
- E->trace.rlater[j][E->trace.number_of_advection_quantities+kk][E->trace.ilater[j]]=E->trace.etrac[j][kk][it];
- }
+ {
+ E->trace.rlater[j][E->trace.number_of_basic_quantities+kk][E->trace.ilater[j]]=E->trace.extraq[j][kk][it];
+ }
return;
}
@@ -2295,18 +2183,18 @@
inewsize=E->trace.ilatersize[j]+E->trace.ilatersize[j]/5+icushion;
for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
- {
- if ((E->trace.rlater[j][kk]=(double *)realloc(E->trace.rlater[j][kk],inewsize*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"AKM(expand later array )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ {
+ if ((E->trace.rlater[j][kk]=(double *)realloc(E->trace.rlater[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(expand later array )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
fprintf(E->trace.fpt,"Expanding physical memory of rlater to %d from %d\n",
- inewsize,E->trace.ilatersize[j]);
+ inewsize,E->trace.ilatersize[j]);
E->trace.ilatersize[j]=inewsize;
@@ -2327,23 +2215,23 @@
int kk;
- ilast_tracer=E->trace.itrac[j][0];
+ ilast_tracer=E->trace.ntracers[j];
/* put last tracer in ejected tracer position */
- E->trace.itrac[j][it]=E->trace.itrac[j][ilast_tracer];
+ E->trace.ielement[j][it]=E->trace.ielement[j][ilast_tracer];
- for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
- {
- E->trace.rtrac[j][kk][it]=E->trace.rtrac[j][kk][ilast_tracer];
- }
+ for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
+ {
+ E->trace.basicq[j][kk][it]=E->trace.basicq[j][kk][ilast_tracer];
+ }
for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- {
- E->trace.etrac[j][kk][it]=E->trace.etrac[j][kk][ilast_tracer];
- }
+ {
+ E->trace.extraq[j][kk][it]=E->trace.extraq[j][kk][ilast_tracer];
+ }
- E->trace.itrac[j][0]--;
+ E->trace.ntracers[j]--;
return;
}
@@ -2436,325 +2324,325 @@
numregnodes=0;
for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- thetamax=0.0;
- thetamin=M_PI;
+ thetamax=0.0;
+ thetamin=M_PI;
- phimax=two_pi;
- phimin=0.0;
+ phimax=two_pi;
+ phimin=0.0;
- for (kk=1;kk<=E->lmesh.nno;kk=kk+E->lmesh.noz)
- {
+ for (kk=1;kk<=E->lmesh.nno;kk=kk+E->lmesh.noz)
+ {
- theta=E->sx[j][1][kk];
- phi=E->sx[j][2][kk];
+ theta=E->sx[j][1][kk];
+ phi=E->sx[j][2][kk];
- thetamax=max(thetamax,theta);
- thetamin=min(thetamin,theta);
+ thetamax=max(thetamax,theta);
+ thetamin=min(thetamin,theta);
- }
+ }
- /* expand range slightly (should take care of poles) */
+ /* expand range slightly (should take care of poles) */
- thetamax=thetamax+expansion;
- thetamax=min(thetamax,M_PI);
+ thetamax=thetamax+expansion;
+ thetamax=min(thetamax,M_PI);
- thetamin=thetamin-expansion;
- thetamin=max(thetamin,0.0);
+ thetamin=thetamin-expansion;
+ thetamin=max(thetamin,0.0);
- /* Convert input data from degrees to radians */
+ /* Convert input data from degrees to radians */
- deltheta=E->trace.deltheta[0]*M_PI/180.0;
- delphi=E->trace.delphi[0]*M_PI/180.0;
+ deltheta=E->trace.deltheta[0]*M_PI/180.0;
+ delphi=E->trace.delphi[0]*M_PI/180.0;
- /* Adjust deltheta and delphi to fit a uniform number of regular elements */
+ /* Adjust deltheta and delphi to fit a uniform number of regular elements */
- numtheta=fabs(thetamax-thetamin)/deltheta;
- numphi=fabs(phimax-phimin)/delphi;
- nodestheta=numtheta+1;
- nodesphi=numphi+1;
- numregel=numtheta*numphi;
- numregnodes=nodestheta*nodesphi;
+ numtheta=fabs(thetamax-thetamin)/deltheta;
+ numphi=fabs(phimax-phimin)/delphi;
+ nodestheta=numtheta+1;
+ nodesphi=numphi+1;
+ numregel=numtheta*numphi;
+ numregnodes=nodestheta*nodesphi;
- if ((numtheta==0)||(numphi==0))
- {
- fprintf(E->trace.fpt,"Error(make_regular_grid): numtheta: %d numphi: %d\n",numtheta,numphi);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((numtheta==0)||(numphi==0))
+ {
+ fprintf(E->trace.fpt,"Error(make_regular_grid): numtheta: %d numphi: %d\n",numtheta,numphi);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- deltheta=fabs(thetamax-thetamin)/(1.0*numtheta);
- delphi=fabs(phimax-phimin)/(1.0*numphi);
+ deltheta=fabs(thetamax-thetamin)/(1.0*numtheta);
+ delphi=fabs(phimax-phimin)/(1.0*numphi);
- /* fill global variables */
+ /* fill global variables */
- E->trace.deltheta[j]=deltheta;
- E->trace.delphi[j]=delphi;
- E->trace.numtheta[j]=numtheta;
- E->trace.numphi[j]=numphi;
- E->trace.thetamax[j]=thetamax;
- E->trace.thetamin[j]=thetamin;
- E->trace.phimax[j]=phimax;
- E->trace.phimin[j]=phimin;
- E->trace.numregel[j]=numregel;
- E->trace.numregnodes[j]=numregnodes;
+ E->trace.deltheta[j]=deltheta;
+ E->trace.delphi[j]=delphi;
+ E->trace.numtheta[j]=numtheta;
+ E->trace.numphi[j]=numphi;
+ E->trace.thetamax[j]=thetamax;
+ E->trace.thetamin[j]=thetamin;
+ E->trace.phimax[j]=phimax;
+ E->trace.phimin[j]=phimin;
+ E->trace.numregel[j]=numregel;
+ E->trace.numregnodes[j]=numregnodes;
- if ( ((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)) < 0.5 )
- {
- fprintf(E->trace.fpt,"\n ! WARNING: regular/real ratio low: %f ! \n",
- ((1.0*numregel)/(1.0*E->lmesh.nel)) );
- fprintf(E->trace.fpt," Should reduce size of regular mesh\n");
- fprintf(stderr,"! WARNING: regular/real ratio low: %f ! \n",
- ((1.0*numregel)/(1.0*E->lmesh.nel)) );
- fprintf(stderr," Should reduce size of regular mesh\n");
- fflush(E->trace.fpt);
- fflush(stderr);
- if (E->trace.itracer_warnings==1) exit(10);
- }
+ if ( ((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)) < 0.5 )
+ {
+ fprintf(E->trace.fpt,"\n ! WARNING: regular/real ratio low: %f ! \n",
+ ((1.0*numregel)/(1.0*E->lmesh.nel)) );
+ fprintf(E->trace.fpt," Should reduce size of regular mesh\n");
+ fprintf(stderr,"! WARNING: regular/real ratio low: %f ! \n",
+ ((1.0*numregel)/(1.0*E->lmesh.nel)) );
+ fprintf(stderr," Should reduce size of regular mesh\n");
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
- /* print some output */
+ /* print some output */
- fprintf(E->trace.fpt,"\nRegular grid:\n");
- fprintf(E->trace.fpt,"Theta min: %f max: %f \n",thetamin,thetamax);
- fprintf(E->trace.fpt,"Phi min: %f max: %f \n",phimin,phimax);
- fprintf(E->trace.fpt,"Adjusted deltheta: %f delphi: %f\n",deltheta,delphi);
- fprintf(E->trace.fpt,"(numtheta: %d numphi: %d)\n",numtheta,numphi);
- fprintf(E->trace.fpt,"Number of regular elements: %d (nodes: %d)\n",numregel,numregnodes);
+ fprintf(E->trace.fpt,"\nRegular grid:\n");
+ fprintf(E->trace.fpt,"Theta min: %f max: %f \n",thetamin,thetamax);
+ fprintf(E->trace.fpt,"Phi min: %f max: %f \n",phimin,phimax);
+ fprintf(E->trace.fpt,"Adjusted deltheta: %f delphi: %f\n",deltheta,delphi);
+ fprintf(E->trace.fpt,"(numtheta: %d numphi: %d)\n",numtheta,numphi);
+ fprintf(E->trace.fpt,"Number of regular elements: %d (nodes: %d)\n",numregel,numregnodes);
- fprintf(E->trace.fpt,"regular/real ratio: %f\n",((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)));
- fflush(E->trace.fpt);
+ fprintf(E->trace.fpt,"regular/real ratio: %f\n",((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)));
+ fflush(E->trace.fpt);
- /* Allocate memory for regnodetoel */
- /* Regtoel is an integer array which represents nodes on */
- /* the regular mesh. Each node on the regular mesh contains */
- /* the real element value if one exists (-99 otherwise) */
+ /* Allocate memory for regnodetoel */
+ /* Regtoel is an integer array which represents nodes on */
+ /* the regular mesh. Each node on the regular mesh contains */
+ /* the real element value if one exists (-99 otherwise) */
- if ((E->trace.regnodetoel[j]=(int *)malloc((numregnodes+1)*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular) -no memory - uh3ud\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((E->trace.regnodetoel[j]=(int *)malloc((numregnodes+1)*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) -no memory - uh3ud\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- /* Initialize regnodetoel - reg elements not used =-99 */
+ /* Initialize regnodetoel - reg elements not used =-99 */
- for (kk=1;kk<=numregnodes;kk++)
- {
- E->trace.regnodetoel[j][kk]=-99;
- }
+ for (kk=1;kk<=numregnodes;kk++)
+ {
+ E->trace.regnodetoel[j][kk]=-99;
+ }
- /* Begin Mapping (only need to use surface elements) */
+ /* Begin Mapping (only need to use surface elements) */
- parallel_process_sync(E);
- if (E->parallel.me==0) fprintf(stderr,"Beginning Mapping\n");
- fflush(stderr);
+ parallel_process_sync(E);
+ if (E->parallel.me==0) fprintf(stderr,"Beginning Mapping\n");
+ fflush(stderr);
- /* Generate temporary arrays of max and min values for each surface element */
+ /* Generate temporary arrays of max and min values for each surface element */
- if ((tmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((tmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((fmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((fmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((tmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((tmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((fmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((fmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
- {
+ for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
+ {
- kk=mm/elz;
+ kk=mm/elz;
- theta_min=M_PI;
- theta_max=0.0;
- phi_min=two_pi;
- phi_max=0.0;
- for (pp=1;pp<=4;pp++)
- {
- node=E->ien[j][mm].node[pp];
- theta=E->sx[j][1][node];
- phi=E->sx[j][2][node];
+ theta_min=M_PI;
+ theta_max=0.0;
+ phi_min=two_pi;
+ phi_max=0.0;
+ for (pp=1;pp<=4;pp++)
+ {
+ node=E->ien[j][mm].node[pp];
+ theta=E->sx[j][1][node];
+ phi=E->sx[j][2][node];
- theta_min=min(theta_min,theta);
- theta_max=max(theta_max,theta);
- phi_min=min(phi_min,phi);
- phi_max=max(phi_max,phi);
- }
+ theta_min=min(theta_min,theta);
+ theta_max=max(theta_max,theta);
+ phi_min=min(phi_min,phi);
+ phi_max=max(phi_max,phi);
+ }
- /* add half difference to phi and expansion to theta to be safe */
+ /* add half difference to phi and expansion to theta to be safe */
- theta_max=theta_max+expansion;
- theta_min=theta_min-expansion;
+ theta_max=theta_max+expansion;
+ theta_min=theta_min-expansion;
- theta_max=min(M_PI,theta_max);
- theta_min=max(0.0,theta_min);
+ theta_max=min(M_PI,theta_max);
+ theta_min=max(0.0,theta_min);
- half_diff=0.5*(phi_max-phi_min);
- phi_max=phi_max+half_diff;
- phi_min=phi_min-half_diff;
+ half_diff=0.5*(phi_max-phi_min);
+ phi_max=phi_max+half_diff;
+ phi_min=phi_min-half_diff;
- fix_phi(&phi_max);
- fix_phi(&phi_min);
+ fix_phi(&phi_max);
+ fix_phi(&phi_min);
- if (phi_min>phi_max)
- {
- phi_min=0.0;
- phi_max=two_pi;
- }
+ if (phi_min>phi_max)
+ {
+ phi_min=0.0;
+ phi_max=two_pi;
+ }
- tmin[kk]=theta_min;
- tmax[kk]=theta_max;
- fmin[kk]=phi_min;
- fmax[kk]=phi_max;
- }
+ tmin[kk]=theta_min;
+ tmax[kk]=theta_max;
+ fmin[kk]=phi_min;
+ fmax[kk]=phi_max;
+ }
- /* end looking through elements */
+ /* end looking through elements */
- ifound_one=0;
+ ifound_one=0;
- rad=E->sphere.ro;
+ rad=E->sphere.ro;
- imap=0;
+ imap=0;
- for (kk=1;kk<=numregnodes;kk++)
- {
+ for (kk=1;kk<=numregnodes;kk++)
+ {
- E->trace.regnodetoel[j][kk]=-99;
+ E->trace.regnodetoel[j][kk]=-99;
- /* find theta and phi for a given regular node */
+ /* find theta and phi for a given regular node */
- idum1=(kk-1)/(numtheta+1);
- idum2=kk-1-idum1*(numtheta+1);
+ idum1=(kk-1)/(numtheta+1);
+ idum2=kk-1-idum1*(numtheta+1);
- theta=thetamin+(1.0*idum2*deltheta);
- phi=phimin+(1.0*idum1*delphi);
+ theta=thetamin+(1.0*idum2*deltheta);
+ phi=phimin+(1.0*idum1*delphi);
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
- ilast_el=1;
+ ilast_el=1;
- /* if previous element not found yet, check all surface elements */
+ /* if previous element not found yet, check all surface elements */
- /*
- if (ifound_one==0)
- {
- for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
- {
- pp=mm/elz;
- if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
- {
- ival=icheck_element_column(E,j,mm,x,y,z,rad);
- if (ival>0)
- {
- ilast_el=mm;
- ifound_one++;
- E->trace.regnodetoel[j][kk]=mm;
- goto foundit;
- }
- }
- }
- goto foundit;
- }
- */
+ /*
+ if (ifound_one==0)
+ {
+ for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
+ {
+ pp=mm/elz;
+ if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
+ {
+ ival=icheck_element_column(E,j,mm,x,y,z,rad);
+ if (ival>0)
+ {
+ ilast_el=mm;
+ ifound_one++;
+ E->trace.regnodetoel[j][kk]=mm;
+ goto foundit;
+ }
+ }
+ }
+ goto foundit;
+ }
+ */
- /* first check previous element */
+ /* first check previous element */
- ival=icheck_element_column(E,j,ilast_el,x,y,z,rad);
- if (ival>0)
- {
- E->trace.regnodetoel[j][kk]=ilast_el;
- goto foundit;
- }
+ ival=icheck_element_column(E,j,ilast_el,x,y,z,rad);
+ if (ival>0)
+ {
+ E->trace.regnodetoel[j][kk]=ilast_el;
+ goto foundit;
+ }
- /* check neighbors */
+ /* check neighbors */
- ival=icheck_column_neighbors(E,j,ilast_el,x,y,z,rad);
- if (ival>0)
- {
- E->trace.regnodetoel[j][kk]=ival;
- ilast_el=ival;
- goto foundit;
- }
+ ival=icheck_column_neighbors(E,j,ilast_el,x,y,z,rad);
+ if (ival>0)
+ {
+ E->trace.regnodetoel[j][kk]=ival;
+ ilast_el=ival;
+ goto foundit;
+ }
- /* check all */
+ /* check all */
- for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
- {
- pp=mm/elz;
- if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
- {
- ival=icheck_element_column(E,j,mm,x,y,z,rad);
- if (ival>0)
- {
- ilast_el=mm;
- E->trace.regnodetoel[j][kk]=mm;
- goto foundit;
- }
- }
- }
+ for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
+ {
+ pp=mm/elz;
+ if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
+ {
+ ival=icheck_element_column(E,j,mm,x,y,z,rad);
+ if (ival>0)
+ {
+ ilast_el=mm;
+ E->trace.regnodetoel[j][kk]=mm;
+ goto foundit;
+ }
+ }
+ }
- foundit:
+ foundit:
- if (E->trace.regnodetoel[j][kk]>0) imap++;
+ if (E->trace.regnodetoel[j][kk]>0) imap++;
- } /* end all regular nodes (kk) */
+ } /* end all regular nodes (kk) */
- fprintf(E->trace.fpt,"percentage mapped: %f\n", (1.0*imap)/(1.0*numregnodes)*100.0);
- fflush(E->trace.fpt);
+ fprintf(E->trace.fpt,"percentage mapped: %f\n", (1.0*imap)/(1.0*numregnodes)*100.0);
+ fflush(E->trace.fpt);
- /* free temporary arrays */
+ /* free temporary arrays */
- free(tmin);
- free(tmax);
- free(fmin);
- free(fmax);
+ free(tmin);
+ free(tmax);
+ free(fmin);
+ free(fmax);
- } /* end j */
+ } /* end j */
/* some error control */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=numregnodes;kk++)
- {
+ {
+ for (kk=1;kk<=numregnodes;kk++)
+ {
- if (E->trace.regnodetoel[j][kk]!=-99)
- {
- if ( (E->trace.regnodetoel[j][kk]<1)||(E->trace.regnodetoel[j][kk]>E->lmesh.nel) )
- {
- fprintf(stderr,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
- fprintf(E->trace.fpt,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
- fflush(E->trace.fpt);
- fflush(stderr);
- exit(10);
- }
- }
- }
- }
+ if (E->trace.regnodetoel[j][kk]!=-99)
+ {
+ if ( (E->trace.regnodetoel[j][kk]<1)||(E->trace.regnodetoel[j][kk]>E->lmesh.nel) )
+ {
+ fprintf(stderr,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
+ fprintf(E->trace.fpt,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ exit(10);
+ }
+ }
+ }
+ }
/* Now put regnodetoel information into regtoel */
@@ -2768,202 +2656,202 @@
for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* initialize statistical counter */
+ /* initialize statistical counter */
- for (pp=0;pp<=4;pp++) istat_ichoice[j][pp]=0;
+ for (pp=0;pp<=4;pp++) istat_ichoice[j][pp]=0;
- /* Allocate memory for regtoel */
- /* Regtoel consists of 4 positions for each regular element */
- /* Position[0] lists the number of element choices (later */
- /* referred to as ichoice), followed */
- /* by the possible element choices. */
- /* ex.) A regular element has 4 nodes. Each node resides in */
- /* a real element. The number of real elements a regular */
- /* element touches (one of its nodes are in) is ichoice. */
- /* Special ichoice notes: */
- /* ichoice=-1 all regular element nodes = -99 (no elements) */
- /* ichoice=0 all 4 corners within one element */
- /* ichoice=1 one element choice (diff from ichoice 0 in */
- /* that perhaps one reg node is in an element */
- /* and the rest are not (-99). */
- /* ichoice>1 Multiple elements to check */
+ /* Allocate memory for regtoel */
+ /* Regtoel consists of 4 positions for each regular element */
+ /* Position[0] lists the number of element choices (later */
+ /* referred to as ichoice), followed */
+ /* by the possible element choices. */
+ /* ex.) A regular element has 4 nodes. Each node resides in */
+ /* a real element. The number of real elements a regular */
+ /* element touches (one of its nodes are in) is ichoice. */
+ /* Special ichoice notes: */
+ /* ichoice=-1 all regular element nodes = -99 (no elements) */
+ /* ichoice=0 all 4 corners within one element */
+ /* ichoice=1 one element choice (diff from ichoice 0 in */
+ /* that perhaps one reg node is in an element */
+ /* and the rest are not (-99). */
+ /* ichoice>1 Multiple elements to check */
- numregel= E->trace.numregel[j];
+ numregel= E->trace.numregel[j];
- for (pp=0;pp<=4;pp++)
- {
- if ((E->trace.regtoel[j][pp]=(int *)malloc((numregel+1)*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(make regular)-no memory 98d (%d %d %d)\n",pp,numregel,j);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (pp=0;pp<=4;pp++)
+ {
+ if ((E->trace.regtoel[j][pp]=(int *)malloc((numregel+1)*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular)-no memory 98d (%d %d %d)\n",pp,numregel,j);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- numtheta=E->trace.numtheta[j];
- numphi=E->trace.numphi[j];
+ numtheta=E->trace.numtheta[j];
+ numphi=E->trace.numphi[j];
- for (nphi=1;nphi<=numphi;nphi++)
- {
- for (ntheta=1;ntheta<=numtheta;ntheta++)
- {
+ for (nphi=1;nphi<=numphi;nphi++)
+ {
+ for (ntheta=1;ntheta<=numtheta;ntheta++)
+ {
- iregel=ntheta+(nphi-1)*numtheta;
+ iregel=ntheta+(nphi-1)*numtheta;
- /* initialize regtoel (not necessary really) */
+ /* initialize regtoel (not necessary really) */
- for (pp=0;pp<=4;pp++) E->trace.regtoel[j][pp][iregel]=-33;
+ for (pp=0;pp<=4;pp++) E->trace.regtoel[j][pp][iregel]=-33;
- if ( (iregel>numregel)||(iregel<1) )
- {
- fprintf(E->trace.fpt,"ERROR(make_regular_grid)-weird iregel: %d (max: %d)\n",iregel,numregel);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ( (iregel>numregel)||(iregel<1) )
+ {
+ fprintf(E->trace.fpt,"ERROR(make_regular_grid)-weird iregel: %d (max: %d)\n",iregel,numregel);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- iregnode[1]=iregel+(nphi-1);
- iregnode[2]=iregel+nphi;
- iregnode[3]=iregel+nphi+E->trace.numtheta[j]+1;
- iregnode[4]=iregel+nphi+E->trace.numtheta[j];
+ iregnode[1]=iregel+(nphi-1);
+ iregnode[2]=iregel+nphi;
+ iregnode[3]=iregel+nphi+E->trace.numtheta[j]+1;
+ iregnode[4]=iregel+nphi+E->trace.numtheta[j];
- for (kk=1;kk<=4;kk++)
- {
- if ((iregnode[kk]<1)||(iregnode[kk]>numregnodes))
- {
- fprintf(E->trace.fpt,"ERROR(make regular)-bad regnode %d\n",iregnode[kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- if (E->trace.regnodetoel[j][iregnode[kk]]>E->lmesh.nel)
- {
- fprintf(E->trace.fpt,"AABB HERE %d %d %d %d\n",iregel,iregnode[kk],kk,E->trace.regnodetoel[j][iregnode[kk]]);
- fflush(E->trace.fpt);
- }
- }
+ for (kk=1;kk<=4;kk++)
+ {
+ if ((iregnode[kk]<1)||(iregnode[kk]>numregnodes))
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular)-bad regnode %d\n",iregnode[kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if (E->trace.regnodetoel[j][iregnode[kk]]>E->lmesh.nel)
+ {
+ fprintf(E->trace.fpt,"AABB HERE %d %d %d %d\n",iregel,iregnode[kk],kk,E->trace.regnodetoel[j][iregnode[kk]]);
+ fflush(E->trace.fpt);
+ }
+ }
- /* find number of choices */
+ /* find number of choices */
- ichoice=0;
- icount=0;
+ ichoice=0;
+ icount=0;
- for (kk=1;kk<=4;kk++)
- {
+ for (kk=1;kk<=4;kk++)
+ {
- if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
+ if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
- icount++;
- for (pp=1;pp<=(kk-1);pp++)
- {
- if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
- }
- ichoice++;
- itemp[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
+ icount++;
+ for (pp=1;pp<=(kk-1);pp++)
+ {
+ if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
+ }
+ ichoice++;
+ itemp[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
- if ((ichoice<0) || (ichoice>4) )
- {
- fprintf(E->trace.fpt,"ERROR(make regular) - weird ichoice %d \n",ichoice);
- fflush(E->trace.fpt);
- exit(10);
- }
- if ((itemp[ichoice]<0) || (itemp[ichoice]>E->lmesh.nel) )
- {
- fprintf(E->trace.fpt,"ERROR(make regular) - weird element choice %d %d\n",itemp[ichoice],ichoice);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((ichoice<0) || (ichoice>4) )
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) - weird ichoice %d \n",ichoice);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((itemp[ichoice]<0) || (itemp[ichoice]>E->lmesh.nel) )
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) - weird element choice %d %d\n",itemp[ichoice],ichoice);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- next_corner:
- ;
- } /* end kk */
+ next_corner:
+ ;
+ } /* end kk */
- istat_ichoice[j][ichoice]++;
+ istat_ichoice[j][ichoice]++;
- if ((ichoice<0) || (ichoice>4))
- {
- fprintf(E->trace.fpt,"ERROR(make_regular)-wierd ichoice %d\n",ichoice);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((ichoice<0) || (ichoice>4))
+ {
+ fprintf(E->trace.fpt,"ERROR(make_regular)-wierd ichoice %d\n",ichoice);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- if (ichoice==0)
- {
- E->trace.regtoel[j][0][iregel]=-1;
- /*
- fprintf(E->trace.fpt,"HH1: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
- */
- }
- else if ( (ichoice==1) && (icount==4) )
- {
- E->trace.regtoel[j][0][iregel]=0;
- E->trace.regtoel[j][1][iregel]=itemp[1];
+ if (ichoice==0)
+ {
+ E->trace.regtoel[j][0][iregel]=-1;
+ /*
+ fprintf(E->trace.fpt,"HH1: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
+ */
+ }
+ else if ( (ichoice==1) && (icount==4) )
+ {
+ E->trace.regtoel[j][0][iregel]=0;
+ E->trace.regtoel[j][1][iregel]=itemp[1];
- /*
- fprintf(E->trace.fpt,"HH2: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
- */
+ /*
+ fprintf(E->trace.fpt,"HH2: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
+ */
- if (itemp[1]<1 || itemp[1]>E->lmesh.nel)
- {
- fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- else if ( (ichoice>0) && (ichoice<5) )
- {
- E->trace.regtoel[j][0][iregel]=ichoice;
- for (pp=1;pp<=ichoice;pp++)
- {
- E->trace.regtoel[j][pp][iregel]=itemp[pp];
+ if (itemp[1]<1 || itemp[1]>E->lmesh.nel)
+ {
+ fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ else if ( (ichoice>0) && (ichoice<5) )
+ {
+ E->trace.regtoel[j][0][iregel]=ichoice;
+ for (pp=1;pp<=ichoice;pp++)
+ {
+ E->trace.regtoel[j][pp][iregel]=itemp[pp];
- /*
- fprintf(E->trace.fpt,"HH:(%p) iregel: %d ichoice: %d pp: %d value: %d %d\n",&E->trace.regtoel[j][pp][iregel],iregel,ichoice,pp,itemp[pp],E->trace.regtoel[j][pp][iregel]);
- */
- if (itemp[pp]<1 || itemp[pp]>E->lmesh.nel)
- {
- fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp 2 \n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
- else
- {
- fprintf(E->trace.fpt,"ERROR(make_regular)- should not be here! %d\n",ichoice);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ /*
+ fprintf(E->trace.fpt,"HH:(%p) iregel: %d ichoice: %d pp: %d value: %d %d\n",&E->trace.regtoel[j][pp][iregel],iregel,ichoice,pp,itemp[pp],E->trace.regtoel[j][pp][iregel]);
+ */
+ if (itemp[pp]<1 || itemp[pp]>E->lmesh.nel)
+ {
+ fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp 2 \n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"ERROR(make_regular)- should not be here! %d\n",ichoice);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- /* can now free regnodetoel */
+ /* can now free regnodetoel */
- free (E->trace.regnodetoel[j]);
+ free (E->trace.regnodetoel[j]);
- /* testing */
- for (kk=1;kk<=E->trace.numregel[j];kk++)
- {
- if ((E->trace.regtoel[j][0][kk]<-1)||(E->trace.regtoel[j][0][kk]>4))
- {
- fprintf(E->trace.fpt,"ERROR(make regular) regtoel ichoice0? %d %d \n",kk,E->trace.regtoel[j][pp][kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- for (pp=1;pp<=4;pp++)
- {
- if (((E->trace.regtoel[j][pp][kk]<1)&&(E->trace.regtoel[j][pp][kk]!=-33))||(E->trace.regtoel[j][pp][kk]>E->lmesh.nel))
- {
- fprintf(E->trace.fpt,"ERROR(make regular) (%p) regtoel? %d %d(%d) %d\n",&E->trace.regtoel[j][pp][kk],kk,pp,E->trace.regtoel[j][0][kk],E->trace.regtoel[j][pp][kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ /* testing */
+ for (kk=1;kk<=E->trace.numregel[j];kk++)
+ {
+ if ((E->trace.regtoel[j][0][kk]<-1)||(E->trace.regtoel[j][0][kk]>4))
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) regtoel ichoice0? %d %d \n",kk,E->trace.regtoel[j][pp][kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ for (pp=1;pp<=4;pp++)
+ {
+ if (((E->trace.regtoel[j][pp][kk]<1)&&(E->trace.regtoel[j][pp][kk]!=-33))||(E->trace.regtoel[j][pp][kk]>E->lmesh.nel))
+ {
+ fprintf(E->trace.fpt,"ERROR(make regular) (%p) regtoel? %d %d(%d) %d\n",&E->trace.regtoel[j][pp][kk],kk,pp,E->trace.regtoel[j][0][kk],E->trace.regtoel[j][pp][kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- } /* end j */
+ } /* end j */
fprintf(E->trace.fpt,"Mapping completed (%f seconds)\n",CPU_time0()-start_time);
@@ -2977,22 +2865,22 @@
/* Print out information regarding regular/real element coverage */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- isum=0;
- for (kk=0;kk<=4;kk++) isum=isum+istat_ichoice[j][kk];
- fprintf(E->trace.fpt,"\n\nInformation regarding number of real elements per regular elements\n");
- fprintf(E->trace.fpt," (stats done on regular elements that were used)\n");
- fprintf(E->trace.fpt,"Ichoice is number of real elements touched by a regular element\n");
- fprintf(E->trace.fpt," (ichoice=1 is optimal)\n");
- fprintf(E->trace.fpt," (if ichoice=0, no elements are touched by regular element)\n");
- fprintf(E->trace.fpt,"Ichoice=0: %f percent\n",(100.0*istat_ichoice[j][0])/(1.0*isum));
- fprintf(E->trace.fpt,"Ichoice=1: %f percent\n",(100.0*istat_ichoice[j][1])/(1.0*isum));
- fprintf(E->trace.fpt,"Ichoice=2: %f percent\n",(100.0*istat_ichoice[j][2])/(1.0*isum));
- fprintf(E->trace.fpt,"Ichoice=3: %f percent\n",(100.0*istat_ichoice[j][3])/(1.0*isum));
- fprintf(E->trace.fpt,"Ichoice=4: %f percent\n",(100.0*istat_ichoice[j][4])/(1.0*isum));
+ isum=0;
+ for (kk=0;kk<=4;kk++) isum=isum+istat_ichoice[j][kk];
+ fprintf(E->trace.fpt,"\n\nInformation regarding number of real elements per regular elements\n");
+ fprintf(E->trace.fpt," (stats done on regular elements that were used)\n");
+ fprintf(E->trace.fpt,"Ichoice is number of real elements touched by a regular element\n");
+ fprintf(E->trace.fpt," (ichoice=1 is optimal)\n");
+ fprintf(E->trace.fpt," (if ichoice=0, no elements are touched by regular element)\n");
+ fprintf(E->trace.fpt,"Ichoice=0: %f percent\n",(100.0*istat_ichoice[j][0])/(1.0*isum));
+ fprintf(E->trace.fpt,"Ichoice=1: %f percent\n",(100.0*istat_ichoice[j][1])/(1.0*isum));
+ fprintf(E->trace.fpt,"Ichoice=2: %f percent\n",(100.0*istat_ichoice[j][2])/(1.0*isum));
+ fprintf(E->trace.fpt,"Ichoice=3: %f percent\n",(100.0*istat_ichoice[j][3])/(1.0*isum));
+ fprintf(E->trace.fpt,"Ichoice=4: %f percent\n",(100.0*istat_ichoice[j][4])/(1.0*isum));
- } /* end j */
+ } /* end j */
return;
@@ -3065,17 +2953,17 @@
for (kk=1;kk<=number_of_neighbors;kk++)
- {
+ {
- if ((neighbor[kk]>=1)&&(neighbor[kk]<=E->lmesh.nel))
- {
- ival=icheck_element_column(E,j,neighbor[kk],x,y,z,rad);
- if (ival>0)
- {
- return neighbor[kk];
- }
- }
- }
+ if ((neighbor[kk]>=1)&&(neighbor[kk]<=E->lmesh.nel))
+ {
+ ival=icheck_element_column(E,j,neighbor[kk],x,y,z,rad);
+ if (ival>0)
+ {
+ return neighbor[kk];
+ }
+ }
+ }
return -99;
}
@@ -3100,13 +2988,13 @@
int numel=E->lmesh.nel;
for (nel=elz;nel<=numel;nel=nel+elz)
- {
- icheck=icheck_element_column(E,j,nel,x,y,z,rad);
- if (icheck==1)
- {
- return nel;
- }
- }
+ {
+ icheck=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (icheck==1)
+ {
+ return nel;
+ }
+ }
return -99;
@@ -3122,70 +3010,76 @@
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
if (E->trace.ic_method==0)
- {
- fprintf(E->trace.fpt,"Generating New Tracer Array\n");
- fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
- /* TODO: move
- if (E->composition.ichemical_buoyancy==1)
- {
- fprintf(E->trace.fpt,"Interface Height: %f\n",E->composition.z_interface);
- }
- */
- }
+ {
+ fprintf(E->trace.fpt,"Generating New Tracer Array\n");
+ fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
+ /* TODO: move
+ if (E->composition.ichemical_buoyancy==1)
+ {
+ fprintf(E->trace.fpt,"Interface Height: %f\n",E->composition.z_interface);
+ }
+ */
+ }
if (E->trace.ic_method==1)
- {
- fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
- }
+ {
+ fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
+ }
if (E->trace.ic_method==2)
- {
- fprintf(E->trace.fpt,"Restarting Tracers\n");
- }
+ {
+ fprintf(E->trace.fpt,"Restarting Tracers\n");
+ }
+ fprintf(E->trace.fpt,"Number of tracer types: %d", E->trace.ntypes);
+ if (E->trace.ntypes < 1) {
+ fprintf(E->trace.fpt, "Tracer types shouldn't be less than 1\n");
+ parallel_process_termination();
+ }
+
if (E->trace.itracer_advection_scheme==1)
- {
- fprintf(E->trace.fpt,"\nSimple predictor-corrector method used\n");
- fprintf(E->trace.fpt,"(Uses only velocity at to) \n");
- fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0) + v(xp,yp,zp))\n\n");
- }
+ {
+ fprintf(E->trace.fpt,"\nSimple predictor-corrector method used\n");
+ fprintf(E->trace.fpt,"(Uses only velocity at to) \n");
+ fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0) + v(xp,yp,zp))\n\n");
+ }
else if (E->trace.itracer_advection_scheme==2)
- {
- fprintf(E->trace.fpt,"\nPredictor-corrector method used\n");
- fprintf(E->trace.fpt,"(Uses only velocity at to and to+dt) \n");
- fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0,to) + v(xp,yp,zp,to+dt))\n\n");
- }
+ {
+ fprintf(E->trace.fpt,"\nPredictor-corrector method used\n");
+ fprintf(E->trace.fpt,"(Uses only velocity at to and to+dt) \n");
+ fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0,to) + v(xp,yp,zp,to+dt))\n\n");
+ }
else
- {
- fprintf(E->trace.fpt,"Sorry-Other Advection Schemes are Unavailable %d\n",E->trace.itracer_advection_scheme);
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
+ {
+ fprintf(E->trace.fpt,"Sorry-Other Advection Schemes are Unavailable %d\n",E->trace.itracer_advection_scheme);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
if (E->trace.itracer_interpolation_scheme==1)
- {
- fprintf(E->trace.fpt,"\nGreat Circle Projection Interpolation Scheme \n");
- }
+ {
+ fprintf(E->trace.fpt,"\nGreat Circle Projection Interpolation Scheme \n");
+ }
else if (E->trace.itracer_interpolation_scheme==2)
- {
- fprintf(E->trace.fpt,"\nSimple Averaging Interpolation Scheme \n");
- fprintf(E->trace.fpt,"\n(Not that great of a scheme!) \n");
+ {
+ fprintf(E->trace.fpt,"\nSimple Averaging Interpolation Scheme \n");
+ fprintf(E->trace.fpt,"\n(Not that great of a scheme!) \n");
- fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
+ fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
- }
+ }
else
- {
- fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
+ {
+ fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
/* regular grid stuff */
fprintf(E->trace.fpt,"Regular Grid-> deltheta: %f delphi: %f\n",
- E->trace.deltheta[0],E->trace.delphi[0]);
+ E->trace.deltheta[0],E->trace.delphi[0]);
@@ -3193,31 +3087,31 @@
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
- fprintf(E->trace.fpt,"Number of Advection Quantities: %d\n",
- E->trace.number_of_advection_quantities);
+ fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
+ E->trace.number_of_basic_quantities);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
- E->trace.number_of_extra_quantities);
+ E->trace.number_of_extra_quantities);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
- E->trace.number_of_tracer_quantities);
+ E->trace.number_of_tracer_quantities);
/* analytical test */
if (E->trace.ianalytical_tracer_test==1)
- {
- fprintf(E->trace.fpt,"\n\n ! Analytical Test Being Performed ! \n");
- fprintf(E->trace.fpt,"(some of the above parameters may not be used or applied\n");
- fprintf(E->trace.fpt,"Velocity functions given in main code\n");
- fflush(E->trace.fpt);
- }
+ {
+ fprintf(E->trace.fpt,"\n\n ! Analytical Test Being Performed ! \n");
+ fprintf(E->trace.fpt,"(some of the above parameters may not be used or applied\n");
+ fprintf(E->trace.fpt,"Velocity functions given in main code\n");
+ fflush(E->trace.fpt);
+ }
if (E->trace.itracer_warnings==0)
- {
- fprintf(E->trace.fpt,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
- fprintf(stderr,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
- fflush(E->trace.fpt);
- fflush(stderr);
- }
+ {
+ fprintf(E->trace.fpt,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fprintf(stderr,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ }
write_composition_instructions(E);
return;
@@ -3236,13 +3130,13 @@
char output_file[200];
char input_s[1000];
- int j,kk;
+ int i,j,kk;
int idum1,ncolumns;
int numtracers;
double rdum1;
double theta,phi,rad;
- double comp;
+ double extra[100];
double x,y,z;
void initialize_tracer_arrays();
@@ -3251,86 +3145,81 @@
FILE *fp1;
+ if (E->trace.number_of_extra_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(restart_tracers)-increase size of extra[]\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
- if ( (fp1=fopen(output_file,"r"))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ( (fp1=fopen(output_file,"r"))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
fprintf(stderr,"Restarting Tracers from %s\n",output_file);
fflush(stderr);
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- {
- fgets(input_s,200,fp1);
- sscanf(input_s,"%d %d %d %lf",
- &idum1, &numtracers, &ncolumns, &rdum1);
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fgets(input_s,200,fp1);
+ sscanf(input_s,"%d %d %d %lf",
+ &idum1, &numtracers, &ncolumns, &rdum1);
- /* some error control */
- if (E->composition.ichemical_buoyancy==0 && ncolumns!=3) {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
- fflush(E->trace.fpt);
- exit(10);
+ /* some error control */
+ if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ /* allocate memory for tracer arrays */
+
+ initialize_tracer_arrays(E,j,numtracers);
+ E->trace.ntracers[j]=numtracers;
+
+ for (kk=1;kk<=numtracers;kk++) {
+ fgets(input_s,200,fp1);
+ if (E->trace.number_of_extra_quantities==0) {
+ sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
}
- if (E->composition.ichemical_buoyancy==1 && ncolumns!=4) {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+ else if (E->trace.number_of_extra_quantities==1) {
+ sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&extra[0]);
+ }
+ /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
+ /* this part has to be changed... */
+ else {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
fflush(E->trace.fpt);
exit(10);
}
- /* allocate memory for tracer arrays */
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
- initialize_tracer_arrays(E,j,numtracers);
- E->trace.itrac[j][0]=numtracers;
+ /* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
- for (kk=1;kk<=numtracers;kk++)
- {
- comp=0.0;
- fgets(input_s,200,fp1);
- if (E->composition.ichemical_buoyancy==0)
- {
- sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
- }
- //TODO: move
- else if (E->composition.ichemical_buoyancy==1)
- {
- sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&comp);
- }
- else
- {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+ E->trace.basicq[j][0][kk]=theta;
+ E->trace.basicq[j][1][kk]=phi;
+ E->trace.basicq[j][2][kk]=rad;
+ E->trace.basicq[j][3][kk]=x;
+ E->trace.basicq[j][4][kk]=y;
+ E->trace.basicq[j][5][kk]=z;
- /* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
+ for (i=0; i<E->trace.number_of_extra_quantities; i++)
+ E->trace.extraq[j][i][kk]=extra[i];
- keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+ }
- E->trace.rtrac[j][0][kk]=theta;
- E->trace.rtrac[j][1][kk]=phi;
- E->trace.rtrac[j][2][kk]=rad;
- E->trace.rtrac[j][3][kk]=x;
- E->trace.rtrac[j][4][kk]=y;
- E->trace.rtrac[j][5][kk]=z;
- //TODO: move
- if (E->composition.ichemical_buoyancy==1) E->trace.etrac[j][0][kk]=comp;
+ fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
+ fflush(E->trace.fpt);
- }
+ }
- fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
- fflush(E->trace.fpt);
- }
-
-
return;
}
@@ -3338,8 +3227,7 @@
/* Here, each cap will generate tracers somewhere */
/* in the sphere - check if its in this cap - then check radial */
-void make_tracer_array(E)
- struct All_variables *E;
+static void make_tracer_array(struct All_variables *E)
{
int kk;
@@ -3365,86 +3253,81 @@
fflush(stderr);
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
- processor_fraction=( ( pow(E->sx[j][3][E->lmesh.noz],3.0)-pow(E->sx[j][3][1],3.0))/
- (pow(E->sphere.ro,3.0)-pow(E->sphere.ri,3.0)));
- tracers_cap=E->mesh.nel*E->trace.itperel*processor_fraction;
- /*
- fprintf(stderr,"AA: proc frac: %f (%d) %d %d %f %f\n",processor_fraction,tracers_cap,E->lmesh.nel,E->parallel.nprocz, E->sx[j][3][E->lmesh.noz],E->sx[j][3][1]);
- */
+ processor_fraction=( ( pow(E->sx[j][3][E->lmesh.noz],3.0)-pow(E->sx[j][3][1],3.0))/
+ (pow(E->sphere.ro,3.0)-pow(E->sphere.ri,3.0)));
+ tracers_cap=E->mesh.nel*E->trace.itperel*processor_fraction;
+ /*
+ fprintf(stderr,"AA: proc frac: %f (%d) %d %d %f %f\n",processor_fraction,tracers_cap,E->lmesh.nel,E->parallel.nprocz, E->sx[j][3][E->lmesh.noz],E->sx[j][3][1]);
+ */
- fprintf(E->trace.fpt,"\nGenerating %d Tracers\n",tracers_cap);
+ fprintf(E->trace.fpt,"\nGenerating %d Tracers\n",tracers_cap);
- initialize_tracer_arrays(E,j,tracers_cap);
+ initialize_tracer_arrays(E,j,tracers_cap);
- /* Tracers are placed randomly in cap */
- /* (intentionally using rand() instead of srand() )*/
+ /* Tracers are placed randomly in cap */
+ /* (intentionally using rand() instead of srand() )*/
- dmin=-1.0*E->sphere.ro;
- dmax=E->sphere.ro;
+ dmin=-1.0*E->sphere.ro;
+ dmax=E->sphere.ro;
- while (E->trace.itrac[j][0]<tracers_cap)
- {
+ while (E->trace.ntracers[j]<tracers_cap) {
- number_of_tries++;
- max_tries=500*tracers_cap;
+ number_of_tries++;
+ max_tries=500*tracers_cap;
- if (number_of_tries>max_tries)
- {
- fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
- fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (number_of_tries>max_tries) {
+ fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
+ fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- random1=(1.0*rand())/(1.0*RAND_MAX);
- random2=(1.0*rand())/(1.0*RAND_MAX);
- random3=(1.0*rand())/(1.0*RAND_MAX);
+ random1=(1.0*rand())/(1.0*RAND_MAX);
+ random2=(1.0*rand())/(1.0*RAND_MAX);
+ random3=(1.0*rand())/(1.0*RAND_MAX);
- x=dmin+random1*(dmax-dmin);
- y=dmin+random2*(dmax-dmin);
- z=dmin+random3*(dmax-dmin);
+ x=dmin+random1*(dmax-dmin);
+ y=dmin+random2*(dmax-dmin);
+ z=dmin+random3*(dmax-dmin);
- /* first check if within shell */
+ /* first check if within shell */
- rad=sqrt(x*x+y*y+z*z);
+ rad=sqrt(x*x+y*y+z*z);
- if (rad>=E->sx[j][3][E->lmesh.noz]) goto next_try;
- if (rad<E->sx[j][3][1]) goto next_try;
+ if (rad>=E->sx[j][3][E->lmesh.noz]) continue;
+ if (rad<E->sx[j][3][1]) continue;
- /* check if in current cap */
+ /* check if in current cap */
- ival=icheck_cap(E,0,x,y,z,rad);
+ ival=icheck_cap(E,0,x,y,z,rad);
- if (ival!=1) goto next_try;
+ if (ival!=1) continue;
- /* Made it, so record tracer information */
+ /* Made it, so record tracer information */
- cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
+ cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
- keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
- E->trace.itrac[j][0]++;
- kk=E->trace.itrac[j][0];
+ E->trace.ntracers[j]++;
+ kk=E->trace.ntracers[j];
- E->trace.rtrac[j][0][kk]=theta;
- E->trace.rtrac[j][1][kk]=phi;
- E->trace.rtrac[j][2][kk]=rad;
- E->trace.rtrac[j][3][kk]=x;
- E->trace.rtrac[j][4][kk]=y;
- E->trace.rtrac[j][5][kk]=z;
+ E->trace.basicq[j][0][kk]=theta;
+ E->trace.basicq[j][1][kk]=phi;
+ E->trace.basicq[j][2][kk]=rad;
+ E->trace.basicq[j][3][kk]=x;
+ E->trace.basicq[j][4][kk]=y;
+ E->trace.basicq[j][5][kk]=z;
- next_try:
- ;
- } /* end while */
+ } /* end while */
- }/* end j */
+ }/* end j */
fprintf(stderr,"DONE Making Tracer Array (%d)\n",E->parallel.me);
fflush(stderr);
@@ -3499,70 +3382,65 @@
iestimate=number_of_tracers/E->parallel.nproc + icushion;
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
- initialize_tracer_arrays(E,j,iestimate);
+ initialize_tracer_arrays(E,j,iestimate);
- for (kk=1;kk<=number_of_tracers;kk++)
- {
- fgets(input_s,200,fptracer);
- sscanf(input_s,"%lf %lf %lf",&rdum1,&rdum2,&rdum3);
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ fgets(input_s,200,fptracer);
+ sscanf(input_s,"%lf %lf %lf",&rdum1,&rdum2,&rdum3);
- theta=rdum1;
- phi=rdum2;
- rad=rdum3;
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+ theta=rdum1;
+ phi=rdum2;
+ rad=rdum3;
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
- /* make sure theta, phi is in range, and radius is within bounds */
+ /* make sure theta, phi is in range, and radius is within bounds */
- keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
- /* check whether tracer is within processor domain */
+ /* check whether tracer is within processor domain */
- icheck=1;
- if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,rad);
- if (icheck!=1) goto next_tracer;
+ icheck=1;
+ if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,rad);
+ if (icheck!=1) continue;
- icheck=icheck_cap(E,0,x,y,z,rad);
- if (icheck==0) goto next_tracer;
+ icheck=icheck_cap(E,0,x,y,z,rad);
+ if (icheck==0) continue;
- /* if still here, tracer is in processor domain */
+ /* if still here, tracer is in processor domain */
- E->trace.itrac[j][0]++;
+ E->trace.ntracers[j]++;
- if (E->trace.itrac[j][0]>=(E->trace.itracsize[j]-5)) expand_tracer_arrays(E,j);
+ if (E->trace.ntracers[j]>=(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
- E->trace.rtrac[j][0][E->trace.itrac[j][0]]=theta;
- E->trace.rtrac[j][1][E->trace.itrac[j][0]]=phi;
- E->trace.rtrac[j][2][E->trace.itrac[j][0]]=rad;
- E->trace.rtrac[j][3][E->trace.itrac[j][0]]=x;
- E->trace.rtrac[j][4][E->trace.itrac[j][0]]=y;
- E->trace.rtrac[j][5][E->trace.itrac[j][0]]=z;
+ E->trace.basicq[j][0][E->trace.ntracers[j]]=theta;
+ E->trace.basicq[j][1][E->trace.ntracers[j]]=phi;
+ E->trace.basicq[j][2][E->trace.ntracers[j]]=rad;
+ E->trace.basicq[j][3][E->trace.ntracers[j]]=x;
+ E->trace.basicq[j][4][E->trace.ntracers[j]]=y;
+ E->trace.basicq[j][5][E->trace.ntracers[j]]=z;
- next_tracer:
- ;
- } /* end kk, number of tracers */
+ } /* end kk, number of tracers */
- fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
- E->trace.itrac[j][0]);
+ fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
+ E->trace.ntracers[j]);
- } /* end j */
+ } /* end j */
fclose(fptracer);
icheck=isum_tracers(E);
- if (icheck!=number_of_tracers)
- {
- fprintf(E->trace.fpt,"ERROR(read_tracer_file) - tracers != number in file\n");
- fprintf(E->trace.fpt,"Tracers in system: %d\n", icheck);
- fprintf(E->trace.fpt,"Tracers in file: %d\n", number_of_tracers);
- fflush(E->trace.fpt);
- exit(10);
- }
+ if (icheck!=number_of_tracers) {
+ fprintf(E->trace.fpt,"ERROR(read_tracer_file) - tracers != number in file\n");
+ fprintf(E->trace.fpt,"Tracers in system: %d\n", icheck);
+ fprintf(E->trace.fpt,"Tracers in file: %d\n", number_of_tracers);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
return;
}
@@ -3716,15 +3594,15 @@
icheck=icheck_shell(E,nel,rad);
if (icheck==0)
- {
- return 0;
- }
+ {
+ return 0;
+ }
icheck=icheck_element_column(E,j,nel,x,y,z,rad);
if (icheck==0)
- {
- return 0;
- }
+ {
+ return 0;
+ }
return 1;
@@ -3788,23 +3666,23 @@
/* surface coords of element nodes */
for (kk=1;kk<=4;kk++)
- {
+ {
- node=E->ien[j][nel].node[kk+4];
+ node=E->ien[j][nel].node[kk+4];
- rnode[kk][1]=E->x[j][1][node];
- rnode[kk][2]=E->x[j][2][node];
- rnode[kk][3]=E->x[j][3][node];
+ rnode[kk][1]=E->x[j][1][node];
+ rnode[kk][2]=E->x[j][2][node];
+ rnode[kk][3]=E->x[j][3][node];
- rnode[kk][4]=E->sx[j][1][node];
- rnode[kk][5]=E->sx[j][2][node];
+ rnode[kk][4]=E->sx[j][1][node];
+ rnode[kk][5]=E->sx[j][2][node];
- rnode[kk][6]=E->SinCos[lev][j][2][node]; /* cos(theta) */
- rnode[kk][7]=E->SinCos[lev][j][0][node]; /* sin(theta) */
- rnode[kk][8]=E->SinCos[lev][j][3][node]; /* cos(phi) */
- rnode[kk][9]=E->SinCos[lev][j][1][node]; /* sin(phi) */
+ rnode[kk][6]=E->SinCos[lev][j][2][node]; /* cos(theta) */
+ rnode[kk][7]=E->SinCos[lev][j][0][node]; /* sin(theta) */
+ rnode[kk][8]=E->SinCos[lev][j][3][node]; /* cos(phi) */
+ rnode[kk][9]=E->SinCos[lev][j][1][node]; /* sin(phi) */
- }
+ }
/* test_point - project to outer radius */
@@ -3841,18 +3719,18 @@
for (kk=1;kk<=4;kk++)
- {
+ {
- rnode[kk][1]=E->trace.xcap[icap][kk];
- rnode[kk][2]=E->trace.ycap[icap][kk];
- rnode[kk][3]=E->trace.zcap[icap][kk];
- rnode[kk][4]=E->trace.theta_cap[icap][kk];
- rnode[kk][5]=E->trace.phi_cap[icap][kk];
- rnode[kk][6]=E->trace.cos_theta[icap][kk];
- rnode[kk][7]=E->trace.sin_theta[icap][kk];
- rnode[kk][8]=E->trace.cos_phi[icap][kk];
- rnode[kk][9]=E->trace.sin_phi[icap][kk];
- }
+ rnode[kk][1]=E->trace.xcap[icap][kk];
+ rnode[kk][2]=E->trace.ycap[icap][kk];
+ rnode[kk][3]=E->trace.zcap[icap][kk];
+ rnode[kk][4]=E->trace.theta_cap[icap][kk];
+ rnode[kk][5]=E->trace.phi_cap[icap][kk];
+ rnode[kk][6]=E->trace.cos_theta[icap][kk];
+ rnode[kk][7]=E->trace.sin_theta[icap][kk];
+ rnode[kk][8]=E->trace.cos_phi[icap][kk];
+ rnode[kk][9]=E->trace.sin_phi[icap][kk];
+ }
/* test_point - project to outer radius */
@@ -3961,46 +3839,46 @@
eps=1e-6;
if (number_of_tries>3)
- {
- fprintf(E->trace.fpt,"Error(icheck_bounds)-too many tries\n");
- fprintf(E->trace.fpt,"Rads: %f %f %f %f\n",rad1,rad2,rad3,rad4);
- fprintf(E->trace.fpt,"Test Point: %f %f %f \n",test_point[1],test_point[2],test_point[3]);
- fprintf(E->trace.fpt,"Nodal points: 1: %f %f %f\n",rnode1[1],rnode1[2],rnode1[3]);
- fprintf(E->trace.fpt,"Nodal points: 2: %f %f %f\n",rnode2[1],rnode2[2],rnode2[3]);
- fprintf(E->trace.fpt,"Nodal points: 3: %f %f %f\n",rnode3[1],rnode3[2],rnode3[3]);
- fprintf(E->trace.fpt,"Nodal points: 4: %f %f %f\n",rnode4[1],rnode4[2],rnode4[3]);
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ fprintf(E->trace.fpt,"Error(icheck_bounds)-too many tries\n");
+ fprintf(E->trace.fpt,"Rads: %f %f %f %f\n",rad1,rad2,rad3,rad4);
+ fprintf(E->trace.fpt,"Test Point: %f %f %f \n",test_point[1],test_point[2],test_point[3]);
+ fprintf(E->trace.fpt,"Nodal points: 1: %f %f %f\n",rnode1[1],rnode1[2],rnode1[3]);
+ fprintf(E->trace.fpt,"Nodal points: 2: %f %f %f\n",rnode2[1],rnode2[2],rnode2[3]);
+ fprintf(E->trace.fpt,"Nodal points: 3: %f %f %f\n",rnode3[1],rnode3[2],rnode3[3]);
+ fprintf(E->trace.fpt,"Nodal points: 4: %f %f %f\n",rnode4[1],rnode4[2],rnode4[3]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
if (fabs(rad1)<=tiny||fabs(rad2)<=tiny||fabs(rad3)<=tiny||fabs(rad4)<=tiny)
- {
- x=test_point[1];
- y=test_point[2];
- z=test_point[3];
- theta=myatan(sqrt(x*x+y*y),z);
- phi=myatan(y,x);
+ {
+ x=test_point[1];
+ y=test_point[2];
+ z=test_point[3];
+ theta=myatan(sqrt(x*x+y*y),z);
+ phi=myatan(y,x);
- if (theta<=M_PI/2.0)
- {
- theta=theta+eps;
- }
- else
- {
- theta=theta-eps;
- }
- phi=phi+eps;
- x=sin(theta)*cos(phi);
- y=sin(theta)*sin(phi);
- z=cos(theta);
- test_point[1]=x;
- test_point[2]=y;
- test_point[3]=z;
+ if (theta<=M_PI/2.0)
+ {
+ theta=theta+eps;
+ }
+ else
+ {
+ theta=theta-eps;
+ }
+ phi=phi+eps;
+ x=sin(theta)*cos(phi);
+ y=sin(theta)*sin(phi);
+ z=cos(theta);
+ test_point[1]=x;
+ test_point[2]=y;
+ test_point[3]=z;
- number_of_tries++;
- goto try_again;
+ number_of_tries++;
+ goto try_again;
- }
+ }
icheck=0;
if (rad1>0.0&&rad2>0.0&&rad3>0.0&&rad4>0.0) icheck=1;
@@ -4198,10 +4076,10 @@
/* check the radial range */
if (E->parallel.nprocz>1)
- {
- ival=icheck_processor_shell(E,j,rad);
- if (ival!=1) return -99;
- }
+ {
+ ival=icheck_processor_shell(E,j,rad);
+ if (ival!=1) return -99;
+ }
/* First check previous element */
@@ -4215,75 +4093,75 @@
iregel=iget_regel(E,j,theta,phi,&ntheta,&nphi);
if (iregel<=0)
- {
- return -99;
- }
+ {
+ return -99;
+ }
/* AKMA put safety here or in make grid */
if (E->trace.regtoel[j][0][iregel]==0)
- {
- iel=E->trace.regtoel[j][1][iregel];
- goto foundit;
- }
+ {
+ iel=E->trace.regtoel[j][1][iregel];
+ goto foundit;
+ }
/* first check previous element */
if (iprevious_element>0)
- {
- ival=icheck_element_column(E,j,iprevious_element,x,y,z,rad);
- if (ival==1)
- {
- iel=iprevious_element;
- goto foundit;
- }
- }
+ {
+ ival=icheck_element_column(E,j,iprevious_element,x,y,z,rad);
+ if (ival==1)
+ {
+ iel=iprevious_element;
+ goto foundit;
+ }
+ }
/* Check all regular mapping choices */
ichoice=0;
if (E->trace.regtoel[j][0][iregel]>0)
- {
+ {
- ichoice=E->trace.regtoel[j][0][iregel];
- for (kk=1;kk<=ichoice;kk++)
- {
- nelem=E->trace.regtoel[j][kk][iregel];
+ ichoice=E->trace.regtoel[j][0][iregel];
+ for (kk=1;kk<=ichoice;kk++)
+ {
+ nelem=E->trace.regtoel[j][kk][iregel];
- if (nelem!=iprevious_element)
- {
- ival=icheck_element_column(E,j,nelem,x,y,z,rad);
- if (ival==1)
- {
- iel=nelem;
- goto foundit;
- }
+ if (nelem!=iprevious_element)
+ {
+ ival=icheck_element_column(E,j,nelem,x,y,z,rad);
+ if (ival==1)
+ {
+ iel=nelem;
+ goto foundit;
+ }
- }
- }
- }
+ }
+ }
+ }
/* If here, it means that tracer could not be found quickly with regular element map */
/* First check previous element neighbors */
if (iprevious_element>0)
- {
- iel=icheck_column_neighbors(E,j,iprevious_element,x,y,z,rad);
- if (iel>0)
- {
- goto foundit;
- }
- }
+ {
+ iel=icheck_column_neighbors(E,j,iprevious_element,x,y,z,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
/* check if still in cap */
ival=icheck_cap(E,0,x,y,z,rad);
if (ival==0)
- {
- return -99;
- }
+ {
+ return -99;
+ }
/* if here, still in cap (hopefully, without a doubt) */
@@ -4295,45 +4173,45 @@
icorner[3]=elxz*(ely-1)+elz;
icorner[4]=elxz*ely;
for (kk=1;kk<=4;kk++)
- {
- ival=icheck_element_column(E,j,icorner[kk],x,y,z,rad);
- if (ival>0)
- {
- iel=icorner[kk];
- goto foundit;
- }
- }
+ {
+ ival=icheck_element_column(E,j,icorner[kk],x,y,z,rad);
+ if (ival>0)
+ {
+ iel=icorner[kk];
+ goto foundit;
+ }
+ }
/* if previous element is not known, check neighbors of those tried in iquick... */
if (iprevious_element<0)
- {
- if (ichoice>0)
- {
- for (kk=1;kk<=ichoice;kk++)
- {
- ineighbor=E->trace.regtoel[j][kk][iregel];
- iel=icheck_column_neighbors(E,j,ineighbor,x,y,z,rad);
- if (iel>0)
- {
- goto foundit;
- }
- }
- }
+ {
+ if (ichoice>0)
+ {
+ for (kk=1;kk<=ichoice;kk++)
+ {
+ ineighbor=E->trace.regtoel[j][kk][iregel];
+ iel=icheck_column_neighbors(E,j,ineighbor,x,y,z,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
+ }
- /* Decided to remove this part - not really needed and complicated */
- /*
- else
- {
- iel=icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad);
- if (iel>0)
- {
- goto foundit;
- }
- }
- */
- }
+ /* Decided to remove this part - not really needed and complicated */
+ /*
+ else
+ {
+ iel=icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
+ */
+ }
/* As a last resort, check all element columns */
@@ -4353,23 +4231,23 @@
*/
if (E->trace.istat1%100==0)
- {
- fprintf(E->trace.fpt,"Checked all elements %d times already this turn\n",E->trace.istat1);
- fprintf(stderr,"Checked all elements %d times already this turn\n",E->trace.istat1);
- fflush(E->trace.fpt);
- fflush(stderr);
- }
+ {
+ fprintf(E->trace.fpt,"Checked all elements %d times already this turn\n",E->trace.istat1);
+ fprintf(stderr,"Checked all elements %d times already this turn\n",E->trace.istat1);
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ }
if (iel>0)
- {
- goto foundit;
- }
+ {
+ goto foundit;
+ }
/* if still here, there is a problem */
fprintf(E->trace.fpt,"Error(iget_element) - element not found\n");
fprintf(E->trace.fpt,"x,y,z,theta,phi,iregel %f %f %f %f %f %d\n",
- x,y,z,theta,phi,iregel);
+ x,y,z,theta,phi,iregel);
fflush(E->trace.fpt);
exit(10);
@@ -4411,16 +4289,16 @@
iradial_element=ibottom_element;
for (kk=1;kk<=elz;kk++)
- {
+ {
- node=E->ien[j][iradial_element].node[8];
- top_rad=E->sx[j][3][node];
+ node=E->ien[j][iradial_element].node[8];
+ top_rad=E->sx[j][3][node];
- if (rad<top_rad) goto found_it;
+ if (rad<top_rad) goto found_it;
- iradial_element++;
+ iradial_element++;
- } /* end kk */
+ } /* end kk */
/* should not be here */
@@ -4464,22 +4342,22 @@
irange=2;
for (kk=-irange;kk<=irange;kk++)
- {
- for (pp=-irange;pp<=irange;pp++)
- {
- new_ntheta=ntheta+kk;
- new_nphi=nphi+pp;
- if ( (new_ntheta>0)&&(new_ntheta<=E->trace.numtheta[j])&&(new_nphi>0)&&(new_nphi<=E->trace.numphi[j]) )
- {
- iregel=new_ntheta+(new_nphi-1)*E->trace.numtheta[j];
- if ((iregel>0) && (iregel<=E->trace.numregel[j]))
- {
- ival=iquick_element_column_search(E,j,iregel,new_ntheta,new_nphi,x,y,z,theta,phi,rad,imap,&ichoice);
- if (ival>0) return ival;
- }
- }
- }
- }
+ {
+ for (pp=-irange;pp<=irange;pp++)
+ {
+ new_ntheta=ntheta+kk;
+ new_nphi=nphi+pp;
+ if ( (new_ntheta>0)&&(new_ntheta<=E->trace.numtheta[j])&&(new_nphi>0)&&(new_nphi<=E->trace.numphi[j]) )
+ {
+ iregel=new_ntheta+(new_nphi-1)*E->trace.numtheta[j];
+ if ((iregel>0) && (iregel<=E->trace.numregel[j]))
+ {
+ ival=iquick_element_column_search(E,j,iregel,new_ntheta,new_nphi,x,y,z,theta,phi,rad,imap,&ichoice);
+ if (ival>0) return ival;
+ }
+ }
+ }
+ }
return -99;
@@ -4547,35 +4425,35 @@
iregnode[4]=itemp2;
for (kk=1;kk<=4;kk++)
- {
- if ((iregnode[kk]<1) || (iregnode[kk]>E->trace.numregnodes[j]) )
- {
- fprintf(E->trace.fpt,"ERROR(iquick)-weird regnode %d\n",iregnode[kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ {
+ if ((iregnode[kk]<1) || (iregnode[kk]>E->trace.numregnodes[j]) )
+ {
+ fprintf(E->trace.fpt,"ERROR(iquick)-weird regnode %d\n",iregnode[kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
/* find number of choices */
ichoice=0;
icount=0;
for (kk=1;kk<=4;kk++)
- {
- if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
+ {
+ if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
- icount++;
- for (pp=1;pp<=(kk-1);pp++)
- {
- if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
- }
- ichoice++;
- imap[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
+ icount++;
+ for (pp=1;pp<=(kk-1);pp++)
+ {
+ if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
+ }
+ ichoice++;
+ imap[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
- next_corner:
- ;
- } /* end kk */
+ next_corner:
+ ;
+ } /* end kk */
*ich=ichoice;
@@ -4601,11 +4479,11 @@
/* check others */
for (kk=1;kk<=ichoice;kk++)
- {
- nel=imap[kk];
- ival=icheck_element_column(E,j,nel,x,y,z,rad);
- if (ival>0) return nel;
- }
+ {
+ nel=imap[kk];
+ ival=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (ival>0) return nel;
+ }
/* if still here, no element was found */
@@ -4670,44 +4548,44 @@
int kk;
int icushion;
- /* expand rtrac and itrac by 20% */
+ /* expand basicq and ielement by 20% */
icushion=100;
- inewsize=E->trace.itracsize[j]+E->trace.itracsize[j]/5+icushion;
+ inewsize=E->trace.max_ntracers[j]+E->trace.max_ntracers[j]/5+icushion;
- if ((E->trace.itrac[j]=(int *)realloc(E->trace.itrac[j],inewsize*sizeof(int)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (itracer)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ if ((E->trace.ielement[j]=(int *)realloc(E->trace.ielement[j],inewsize*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (ielement)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
- {
- if ((E->trace.rtrac[j][kk]=(double *)realloc(E->trace.rtrac[j][kk],inewsize*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
+ {
+ if ((E->trace.basicq[j][kk]=(double *)realloc(E->trace.basicq[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- {
- if ((E->trace.etrac[j][kk]=(double *)realloc(E->trace.etrac[j][kk],inewsize*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory 78 (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ {
+ if ((E->trace.extraq[j][kk]=(double *)realloc(E->trace.extraq[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory 78 (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- fprintf(E->trace.fpt,"Expanding physical memory of itrac, rtrac, and etrac to %d from %d\n",
- inewsize,E->trace.itracsize[j]);
+ fprintf(E->trace.fpt,"Expanding physical memory of ielement, basicq, and extraq to %d from %d\n",
+ inewsize,E->trace.max_ntracers[j]);
- E->trace.itracsize[j]=inewsize;
+ E->trace.max_ntracers[j]=inewsize;
return;
}
@@ -4741,54 +4619,54 @@
/* open memory for uv space coords */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- for (kk=1;kk<=2;kk++)
- {
+ for (kk=1;kk<=2;kk++)
+ {
//TODO: allocate for surface nodes only to save memory
- if ((E->trace.UV[j][kk]=(double *)malloc((numnodes+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(define uv)-not enough memory(a)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ if ((E->trace.UV[j][kk]=(double *)malloc((numnodes+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(define uv)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
- /* uv space requires a reference point */
- /* UV[j][1][0]=fixed theta */
- /* UV[j][2][0]=fixed phi */
+ /* uv space requires a reference point */
+ /* UV[j][1][0]=fixed theta */
+ /* UV[j][2][0]=fixed phi */
- midnode=numnodes/2;
+ midnode=numnodes/2;
- E->trace.UV[j][1][0]=E->sx[j][1][midnode];
- E->trace.UV[j][2][0]=E->sx[j][2][midnode];
+ E->trace.UV[j][1][0]=E->sx[j][1][midnode];
+ E->trace.UV[j][2][0]=E->sx[j][2][midnode];
- theta_f=E->sx[j][1][midnode];
- phi_f=E->sx[j][2][midnode];
+ theta_f=E->sx[j][1][midnode];
+ phi_f=E->sx[j][2][midnode];
- E->trace.cos_theta_f=cos(theta_f);
- E->trace.sin_theta_f=sin(theta_f);
+ E->trace.cos_theta_f=cos(theta_f);
+ E->trace.sin_theta_f=sin(theta_f);
- /* convert each nodal point to u and v */
+ /* convert each nodal point to u and v */
- for (node=1;node<=numnodes;node++)
- {
- theta=E->sx[j][1][node];
- phi=E->sx[j][2][node];
+ for (node=1;node<=numnodes;node++)
+ {
+ theta=E->sx[j][1][node];
+ phi=E->sx[j][2][node];
- cosc=cos(theta_f)*cos(theta)+sin(theta_f)*sin(theta)*
- cos(phi-phi_f);
- u=sin(theta)*sin(phi-phi_f)/cosc;
- v=(sin(theta_f)*cos(theta)-cos(theta_f)*sin(theta)*cos(phi-phi_f))/cosc;
+ cosc=cos(theta_f)*cos(theta)+sin(theta_f)*sin(theta)*
+ cos(phi-phi_f);
+ u=sin(theta)*sin(phi-phi_f)/cosc;
+ v=(sin(theta_f)*cos(theta)-cos(theta_f)*sin(theta)*cos(phi-phi_f))/cosc;
- E->trace.UV[j][1][node]=u;
- E->trace.UV[j][2][node]=v;
+ E->trace.UV[j][1][node]=u;
+ E->trace.UV[j][2][node]=v;
- }
+ }
- }/*end cap */
+ }/*end cap */
return;
}
@@ -4827,97 +4705,97 @@
fflush(stderr);
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ {
- /* first, allocate memory */
+ /* first, allocate memory */
- for(iwedge=1;iwedge<=2;iwedge++)
- {
- for (kk=1;kk<=9;kk++)
- {
+ for(iwedge=1;iwedge<=2;iwedge++)
+ {
+ for (kk=1;kk<=9;kk++)
+ {
//TODO: allocate for surface elements only to save memory
- if ((E->trace.shape_coefs[j][iwedge][kk]=
- (double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ if ((E->trace.shape_coefs[j][iwedge][kk]=
+ (double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- for (nelem=1;nelem<=E->lmesh.nel;nelem++)
- {
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++)
+ {
- /* find u,v of local nodes at one radius */
+ /* find u,v of local nodes at one radius */
- for(kk=1;kk<=4;kk++)
- {
- node=E->IEN[E->mesh.levmax][j][nelem].node[kk];
- u[kk]=E->trace.UV[j][1][node];
- v[kk]=E->trace.UV[j][2][node];
- }
+ for(kk=1;kk<=4;kk++)
+ {
+ node=E->IEN[E->mesh.levmax][j][nelem].node[kk];
+ u[kk]=E->trace.UV[j][1][node];
+ v[kk]=E->trace.UV[j][2][node];
+ }
- for(iwedge=1;iwedge<=2;iwedge++)
- {
+ for(iwedge=1;iwedge<=2;iwedge++)
+ {
- if (iwedge==1)
- {
- x1=u[1];
- x2=u[2];
- x3=u[3];
- y1=v[1];
- y2=v[2];
- y3=v[3];
- }
- if (iwedge==2)
- {
- x1=u[1];
- x2=u[3];
- x3=u[4];
- y1=v[1];
- y2=v[3];
- y3=v[4];
- }
+ if (iwedge==1)
+ {
+ x1=u[1];
+ x2=u[2];
+ x3=u[3];
+ y1=v[1];
+ y2=v[2];
+ y3=v[3];
+ }
+ if (iwedge==2)
+ {
+ x1=u[1];
+ x2=u[3];
+ x3=u[4];
+ y1=v[1];
+ y2=v[3];
+ y3=v[4];
+ }
- /* shape function 1 */
+ /* shape function 1 */
- delta=(x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
- a0=(x2*y3-x3*y2)/delta;
- a1=(y2-y3)/delta;
- a2=(x3-x2)/delta;
+ delta=(x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
+ a0=(x2*y3-x3*y2)/delta;
+ a1=(y2-y3)/delta;
+ a2=(x3-x2)/delta;
- E->trace.shape_coefs[j][iwedge][1][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][2][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][3][nelem]=a2;
+ E->trace.shape_coefs[j][iwedge][1][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][2][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][3][nelem]=a2;
- /* shape function 2 */
+ /* shape function 2 */
- delta=(x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
- a0=(x1*y3-x3*y1)/delta;
- a1=(y1-y3)/delta;
- a2=(x3-x1)/delta;
+ delta=(x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
+ a0=(x1*y3-x3*y1)/delta;
+ a1=(y1-y3)/delta;
+ a2=(x3-x1)/delta;
- E->trace.shape_coefs[j][iwedge][4][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][5][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][6][nelem]=a2;
+ E->trace.shape_coefs[j][iwedge][4][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][5][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][6][nelem]=a2;
- /* shape function 3 */
+ /* shape function 3 */
- delta=(x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
- a0=(x2*y1-x1*y2)/delta;
- a1=(y2-y1)/delta;
- a2=(x1-x2)/delta;
+ delta=(x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
+ a0=(x2*y1-x1*y2)/delta;
+ a1=(y2-y1)/delta;
+ a2=(x1-x2)/delta;
- E->trace.shape_coefs[j][iwedge][7][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][8][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][9][nelem]=a2;
+ E->trace.shape_coefs[j][iwedge][7][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][8][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][9][nelem]=a2;
- } /* end wedge */
- } /* end elem */
- } /* end cap */
+ } /* end wedge */
+ } /* end elem */
+ } /* end cap */
return;
@@ -4943,49 +4821,49 @@
static int been_here=0;
if (been_here==0)
- {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=3;kk++)
- {
- if ((E->trace.V0_cart[j][kk]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(get_cartesian_velocity)-no memory 82hd7\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=3;kk++)
+ {
+ if ((E->trace.V0_cart[j][kk]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(get_cartesian_velocity)-no memory 82hd7\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- been_here++;
- }
+ been_here++;
+ }
for (m=1;m<=E->sphere.caps_per_proc;m++)
- {
+ {
- for (i=1;i<=E->lmesh.nno;i++)
- {
- sint=E->SinCos[lev][m][0][i];
- sinf=E->SinCos[lev][m][1][i];
- cost=E->SinCos[lev][m][2][i];
- cosf=E->SinCos[lev][m][3][i];
+ for (i=1;i<=E->lmesh.nno;i++)
+ {
+ sint=E->SinCos[lev][m][0][i];
+ sinf=E->SinCos[lev][m][1][i];
+ cost=E->SinCos[lev][m][2][i];
+ cosf=E->SinCos[lev][m][3][i];
- v_theta=E->sphere.cap[m].V[1][i];
- v_phi=E->sphere.cap[m].V[2][i];
- v_rad=E->sphere.cap[m].V[3][i];
+ v_theta=E->sphere.cap[m].V[1][i];
+ v_phi=E->sphere.cap[m].V[2][i];
+ v_rad=E->sphere.cap[m].V[3][i];
- vx=v_theta*cost*cosf-v_phi*sinf+v_rad*sint*cosf;
- vy=v_theta*cost*sinf+v_phi*cosf+v_rad*sint*sinf;
- vz=-v_theta*sint+v_rad*cost;
+ vx=v_theta*cost*cosf-v_phi*sinf+v_rad*sint*cosf;
+ vy=v_theta*cost*sinf+v_phi*cosf+v_rad*sint*sinf;
+ vz=-v_theta*sint+v_rad*cost;
- E->trace.V0_cart[m][1][i]=vx;
- E->trace.V0_cart[m][2][i]=vy;
- E->trace.V0_cart[m][3][i]=vz;
- }
- }
+ E->trace.V0_cart[m][1][i]=vx;
+ E->trace.V0_cart[m][2][i]=vy;
+ E->trace.V0_cart[m][3][i]=vz;
+ }
+ }
return;
}
@@ -5025,12 +4903,12 @@
iold_number=E->trace.ilast_tracer_count;
if (number!=iold_number)
- {
- fprintf(E->trace.fpt,"ERROR(check_sum)-break in conservation %d %d\n",
- number,iold_number);
- fflush(E->trace.fpt);
- exit(10);
- }
+ {
+ fprintf(E->trace.fpt,"ERROR(check_sum)-break in conservation %d %d\n",
+ number,iold_number);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
E->trace.ilast_tracer_count=number;
@@ -5054,9 +4932,9 @@
imycount=0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- imycount=imycount+E->trace.itrac[j][0];
- }
+ {
+ imycount=imycount+E->trace.ntracers[j];
+ }
MPI_Allreduce(&imycount,&iallcount,1,MPI_INT,MPI_SUM,E->parallel.world);
@@ -5133,21 +5011,21 @@
/* Assign test velocity to Citcom nodes */
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->lmesh.nno;kk++)
- {
+ {
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
- theta=E->sx[j][1][kk];
- phi=E->sx[j][2][kk];
- rad=E->sx[j][3][kk];
+ theta=E->sx[j][1][kk];
+ phi=E->sx[j][2][kk];
+ rad=E->sx[j][3][kk];
- analytical_test_function(E,theta,phi,rad,vel_s,vel_c);
+ analytical_test_function(E,theta,phi,rad,vel_s,vel_c);
- E->sphere.cap[j].V[1][kk]=vel_s[1];
- E->sphere.cap[j].V[2][kk]=vel_s[2];
- E->sphere.cap[j].V[3][kk]=vel_s[3];
- }
- }
+ E->sphere.cap[j].V[1][kk]=vel_s[1];
+ E->sphere.cap[j].V[2][kk]=vel_s[2];
+ E->sphere.cap[j].V[3][kk]=vel_s[3];
+ }
+ }
time=0.0;
@@ -5159,72 +5037,72 @@
my_radf=0.0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- if (E->trace.itrac[j][0]>10)
- {
- fprintf(E->trace.fpt,"Warning(analytical)-too many tracers to print!\n");
- fflush(E->trace.fpt);
- if (E->trace.itracer_warnings==1) exit(10);
- }
- }
+ {
+ if (E->trace.ntracers[j]>10)
+ {
+ fprintf(E->trace.fpt,"Warning(analytical)-too many tracers to print!\n");
+ fflush(E->trace.fpt);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
+ }
/* print initial positions */
E->monitor.solution_cycles=0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (pp=1;pp<=E->trace.itrac[j][0];pp++)
- {
- theta=E->trace.rtrac[j][0][pp];
- phi=E->trace.rtrac[j][1][pp];
- rad=E->trace.rtrac[j][2][pp];
+ {
+ for (pp=1;pp<=E->trace.ntracers[j];pp++)
+ {
+ theta=E->trace.basicq[j][0][pp];
+ phi=E->trace.basicq[j][1][pp];
+ rad=E->trace.basicq[j][2][pp];
- fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+ fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
- if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+ if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
- if (pp==1)
- {
- my_theta0=theta;
- my_phi0=phi;
- my_rad0=rad;
- }
- }
- }
+ if (pp==1)
+ {
+ my_theta0=theta;
+ my_phi0=phi;
+ my_rad0=rad;
+ }
+ }
+ }
/* advect tracers */
for (kk=1;kk<=nsteps;kk++)
- {
- E->monitor.solution_cycles=kk;
+ {
+ E->monitor.solution_cycles=kk;
- time=time+dt;
+ time=time+dt;
- predict_tracers(E);
- correct_tracers(E);
+ predict_tracers(E);
+ correct_tracers(E);
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (pp=1;pp<=E->trace.itrac[j][0];pp++)
- {
- theta=E->trace.rtrac[j][0][pp];
- phi=E->trace.rtrac[j][1][pp];
- rad=E->trace.rtrac[j][2][pp];
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (pp=1;pp<=E->trace.ntracers[j];pp++)
+ {
+ theta=E->trace.basicq[j][0][pp];
+ phi=E->trace.basicq[j][1][pp];
+ rad=E->trace.basicq[j][2][pp];
- fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+ fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
- if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+ if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
- if ((kk==nsteps) && (pp==1))
- {
- my_thetaf=theta;
- my_phif=phi;
- my_radf=rad;
- }
- }
- }
+ if ((kk==nsteps) && (pp==1))
+ {
+ my_thetaf=theta;
+ my_phif=phi;
+ my_radf=rad;
+ }
+ }
+ }
- }
+ }
/* Get ready for comparison to Runge-Kutte (only works for one tracer) */
@@ -5236,9 +5114,9 @@
if (E->parallel.me==0) fprintf(stderr,"Comparison to Runge-Kutte\n");
for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- my_number=E->trace.itrac[j][0];
- }
+ {
+ my_number=E->trace.ntracers[j];
+ }
MPI_Allreduce(&my_number,&number,1,MPI_INT,MPI_SUM,E->parallel.world);
@@ -5248,13 +5126,13 @@
/* if more than 1 tracer, exit */
if (number!=1)
- {
- fprintf(E->trace.fpt,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
- if (E->parallel.me==0) fprintf(stderr,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
- fflush(E->trace.fpt);
- fflush(stderr);
- parallel_process_termination();
- }
+ {
+ fprintf(E->trace.fpt,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
+ if (E->parallel.me==0) fprintf(stderr,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ parallel_process_termination();
+ }
/* communicate starting and final positions */
@@ -5312,22 +5190,22 @@
fprintf(E->trace.fpt,"\n\n Difference between Citcom and RK: %e (diff per path length: %e)\n\n",difference,difperpath);
if (E->parallel.me==0)
- {
- fprintf(stderr,"Citcom calculation: steps: %d dt: %f\n",nsteps,dt);
- fprintf(stderr," (nodes per cap: %d x %d x %d)\n",E->lmesh.nox,E->lmesh.noy,(E->lmesh.noz-1)*E->parallel.nprocz+1);
- fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
- fprintf(stderr," final position: theta: %f phi: %f rad: %f\n", thetaf,phif,radf);
- fprintf(stderr," (final time: %f) \n",time );
+ {
+ fprintf(stderr,"Citcom calculation: steps: %d dt: %f\n",nsteps,dt);
+ fprintf(stderr," (nodes per cap: %d x %d x %d)\n",E->lmesh.nox,E->lmesh.noy,(E->lmesh.noz-1)*E->parallel.nprocz+1);
+ fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
+ fprintf(stderr," final position: theta: %f phi: %f rad: %f\n", thetaf,phif,radf);
+ fprintf(stderr," (final time: %f) \n",time );
- fprintf(stderr,"\n\nRunge-Kutte calculation: steps: %d dt: %f\n",nrunge_steps,runge_dt);
- fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
- fprintf(stderr," final position: theta: %f phi: %f rad: %f\n",xf_s[1],xf_s[2],xf_s[3]);
- fprintf(stderr," path length: %f \n",runge_path_length );
- fprintf(stderr," (final time: %f) \n",runge_time );
+ fprintf(stderr,"\n\nRunge-Kutte calculation: steps: %d dt: %f\n",nrunge_steps,runge_dt);
+ fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
+ fprintf(stderr," final position: theta: %f phi: %f rad: %f\n",xf_s[1],xf_s[2],xf_s[3]);
+ fprintf(stderr," path length: %f \n",runge_path_length );
+ fprintf(stderr," (final time: %f) \n",runge_time );
- fprintf(stderr,"\n\n Difference between Citcom and RK: %e (diff per path length: %e)\n\n",difference,difperpath);
+ fprintf(stderr,"\n\n Difference between Citcom and RK: %e (diff per path length: %e)\n\n",difference,difperpath);
- }
+ }
fflush(E->trace.fpt);
fflush(stderr);
@@ -5384,50 +5262,50 @@
path=0.0;
for (kk=1;kk<=nsteps;kk++)
- {
+ {
- /* get velocity at initial position */
+ /* get velocity at initial position */
- analytical_test_function(E,theta_0,phi_0,rad_0,vel0_s,vel0_c);
+ analytical_test_function(E,theta_0,phi_0,rad_0,vel0_s,vel0_c);
- /* Find predicted midpoint position */
+ /* Find predicted midpoint position */
- x_p=x_0+vel0_c[1]*dt*0.5;
- y_p=y_0+vel0_c[2]*dt*0.5;
- z_p=z_0+vel0_c[3]*dt*0.5;
+ x_p=x_0+vel0_c[1]*dt*0.5;
+ y_p=y_0+vel0_c[2]*dt*0.5;
+ z_p=z_0+vel0_c[3]*dt*0.5;
- /* convert to spherical */
+ /* convert to spherical */
- cart_to_sphere(E,x_p,y_p,z_p,&theta_p,&phi_p,&rad_p);
+ cart_to_sphere(E,x_p,y_p,z_p,&theta_p,&phi_p,&rad_p);
- /* get velocity at predicted midpoint position */
+ /* get velocity at predicted midpoint position */
- analytical_test_function(E,theta_p,phi_p,rad_p,velp_s,velp_c);
+ analytical_test_function(E,theta_p,phi_p,rad_p,velp_s,velp_c);
- /* Find corrected position using midpoint velocity */
+ /* Find corrected position using midpoint velocity */
- x_c=x_0+velp_c[1]*dt;
- y_c=y_0+velp_c[2]*dt;
- z_c=z_0+velp_c[3]*dt;
+ x_c=x_0+velp_c[1]*dt;
+ y_c=y_0+velp_c[2]*dt;
+ z_c=z_0+velp_c[3]*dt;
- /* convert to spherical */
+ /* convert to spherical */
- cart_to_sphere(E,x_c,y_c,z_c,&theta_c,&phi_c,&rad_c);
+ cart_to_sphere(E,x_c,y_c,z_c,&theta_c,&phi_c,&rad_c);
- /* compute path lenght */
+ /* compute path lenght */
- dtpath=sqrt((x_c-x_0)*(x_c-x_0)+(y_c-y_0)*(y_c-y_0)+(z_c-z_0)*(z_c-z_0));
- path=path+dtpath;
+ dtpath=sqrt((x_c-x_0)*(x_c-x_0)+(y_c-y_0)*(y_c-y_0)+(z_c-z_0)*(z_c-z_0));
+ path=path+dtpath;
- time=time+dt;
+ time=time+dt;
- x_0=x_c;
- y_0=y_c;
- z_0=z_c;
+ x_0=x_c;
+ y_0=y_c;
+ z_0=z_c;
- /* next time step */
+ /* next time step */
- }
+ }
/* fill final spherical and cartesian vectors to send back */
@@ -5445,6 +5323,46 @@
return;
}
+
+
+/********************************************************************/
+/* This function computes the number of tracers in each element. */
+/* Each tracer can be of different "type", which is the 0th index */
+/* of extraq. How to interprete "type" is left for the application. */
+
+void accumulate_tracers_in_element(struct All_variables *E)
+{
+ /* how many types of tracers? */
+ // TODO: fix to 1, generalized it later
+ const int itypes = 1;
+ int kk;
+ int numtracers;
+ int nelem;
+ int j;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ /* first zero arrays */
+ for (kk=1;kk<=E->lmesh.nel;kk++) {
+ }
+
+ numtracers=E->trace.ntracers[j];
+
+ /* Fill arrays */
+ for (kk=1;kk<=numtracers;kk++) {
+
+ nelem=E->trace.ielement[j][kk];
+
+ //E->composition.itypes[j][nelem]++;
+
+
+ }
+ }
+
+ return;
+}
+
+
/**************** ANALYTICAL TEST FUNCTION ******************/
/* */
/* vel_s[1] => velocity in theta direction */
@@ -5557,8 +5475,8 @@
}
/* Storing the current cap information */
- for (i=0; i<n; i++)
- rr[0][i] = xx[i];
+ for (i=0; i<n; i++)
+ rr[0][i] = xx[i];
/* Wait for non-blocking calls to complete */
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -393,7 +393,7 @@
void output_tracer(struct All_variables *E, int cycles)
{
- int j, n, ncolumns;
+ int i, j, n, ncolumns;
char output_file[255];
FILE *fp1;
@@ -414,32 +414,25 @@
}
else {
/* global model */
- if (E->composition.ichemical_buoyancy==1)
- ncolumns = 4;
- else if (E->composition.ichemical_buoyancy==0)
- ncolumns = 3;
+ ncolumns = 3 + E->trace.number_of_extra_quantities;
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.itrac[j][0],
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
ncolumns, E->monitor.elapsed_time);
- if (E->composition.ichemical_buoyancy==1) {
- for(n=1;n<=E->trace.itrac[j][0];n++) {
- fprintf(fp1,"%.12e %.12e %.12e %.12e\n",
- E->trace.rtrac[j][0][n],
- E->trace.rtrac[j][1][n],
- E->trace.rtrac[j][2][n],
- E->trace.etrac[j][0][n]);
+ for(n=1;n<=E->trace.ntracers[j];n++) {
+ /* write basic quantities (coordinate) */
+ fprintf(fp1,"%.12e %.12e %.12e",
+ E->trace.basicq[j][0][n],
+ E->trace.basicq[j][1][n],
+ E->trace.basicq[j][2][n]);
+
+ /* write extra quantities */
+ for (i=0; i<E->trace.number_of_extra_quantities; i++) {
+ fprintf(fp1," %.12e", E->trace.extraq[j][i][n]);
}
+ fprintf(fp1, "\n");
}
- else if (E->composition.ichemical_buoyancy==0) {
- for(n=1;n<=E->trace.itrac[j][0];n++) {
- fprintf(fp1,"%.12e %.12e %.12e\n",
- E->trace.rtrac[j][0][n],
- E->trace.rtrac[j][1][n],
- E->trace.rtrac[j][2][n]);
- }
- }
}
}
@@ -473,7 +466,7 @@
for(j=1;j<=E->sphere.caps_per_proc;j++) {
fprintf(fp1,"%3d %7d\n", j, E->lmesh.nno);
for(i=1;i<=E->lmesh.nno;i++) {
- fprintf(fp1,"%.6e\n",E->trace.comp_node[j][i]);
+ fprintf(fp1,"%.6e\n",E->composition.comp_node[j][i]);
}
}
@@ -508,7 +501,7 @@
for(j=1;j<=E->sphere.caps_per_proc;j++) {
fprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
for(i=1;i<=E->lmesh.nel;i++) {
- fprintf(fp1,"%.6e\n",E->trace.oldel[j][i]);
+ fprintf(fp1,"%.6e\n",E->composition.oldel[j][i]);
}
}
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -145,7 +145,7 @@
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->trace.comp_node[m][i];
+ buoy[m][i] -= temp2 * E->composition.comp_node[m][i];
}
phase_change_apply_410(E, buoy);
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-11 05:20:43 UTC (rev 6222)
@@ -674,6 +674,7 @@
int ireset_initial_composition;
int ibuoy_type;
int *ieltrac[13];
+ double *oldel[13];
double *celtrac[13];
double *comp_el[13];
double *comp_node[13];
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-11 05:20:43 UTC (rev 6222)
@@ -27,96 +27,89 @@
*/
struct Tracer {
-float *tracer_x;
-float *tracer_y;
-float *tracer_z;
-float *itcolor;
-float *x_space, *z_space, *y_space;
-int NUM_TRACERS;
-int LOCAL_NUM_TRACERS;
+ float *tracer_x;
+ float *tracer_y;
+ float *tracer_z;
+ float *itcolor;
+ float *x_space, *z_space, *y_space;
+ int NUM_TRACERS;
+ int LOCAL_NUM_TRACERS;
-int *LOCAL_ELEMENT;
+ int *LOCAL_ELEMENT;
-float *THETA_LOC_ELEM;
-float *FI_LOC_ELEM;
-float *R_LOC_ELEM;
+ float *THETA_LOC_ELEM;
+ float *FI_LOC_ELEM;
+ float *R_LOC_ELEM;
-float *THETA_LOC_ELEM_T;
-float *FI_LOC_ELEM_T;
-float *R_LOC_ELEM_T;
+ float *THETA_LOC_ELEM_T;
+ float *FI_LOC_ELEM_T;
+ float *R_LOC_ELEM_T;
};
struct TRACE{
-FILE *fpt;
+ FILE *fpt;
-FILE *fp_error_fraction;
-FILE *fp_composition;
+ FILE *fp_error_fraction;
+ FILE *fp_composition;
-char tracer_file[200];
+ char tracer_file[200];
-int number_of_advection_quantities;
-int number_of_extra_quantities;
-int number_of_tracer_quantities;
-int itracer_warnings;
+ int itracer_warnings;
+ int ianalytical_tracer_test;
+ int ic_method;
+ int itperel;
+ int itracer_advection_scheme;
+ int itracer_interpolation_scheme;
-int itracing;
-int ianalytical_tracer_test;
-int ic_method;
-int itperel;
-int itracer_advection_scheme;
-int itracer_interpolation_scheme;
+ double box_cushion;
-int icompositional_rheology;
+ /* tracer arrays */
+ int number_of_basic_quantities;
+ int number_of_extra_quantities;
+ int number_of_tracer_quantities;
-int itracsize[13];
-int *itrac[13];
+ double *basicq[13][100];
+ double *extraq[13][100];
-double box_cushion;
-double *rtrac[13][100];
-double *etrac[13][100];
-double *rlater[13][100];
-/*
-double **rtrac[13];
-double **etrac[13];
-double **rlater[13];
-*/
-double *oldel[13];
+ int ntracers[13];
+ int max_ntracers[13];
+ int *ielement[13];
-/* horizontally averaged */
-double *Have_C;
-double *Havel_tracers;
+ int ilatersize[13];
+ int ilater[13];
+ double *rlater[13][100];
+ /* tracer types */
+ int ntypes;
+ int *ntracer_in_element[13];
-int ilater[13];
-int ilatersize[13];
+ /* regular mesh parameters */
+ int numtheta[13];
+ int numphi[13];
+ unsigned int numregel[13];
+ unsigned int numregnodes[13];
+ double deltheta[13];
+ double delphi[13];
+ double thetamax[13];
+ double thetamin[13];
+ double phimax[13];
+ double phimin[13];
+ int *regnodetoel[13];
+ int *regtoel[13][5];
-/* regular mesh parameters */
-int numtheta[13];
-int numphi[13];
-unsigned int numregel[13];
-unsigned int numregnodes[13];
-double deltheta[13];
-double delphi[13];
-double thetamax[13];
-double thetamin[13];
-double phimax[13];
-double phimin[13];
-int *regnodetoel[13];
-int *regtoel[13][5];
+ /* statistical parameters */
+ int istat_ichoice[13][5];
+ int istat_isend;
+ int istat_iempty;
+ int istat1;
+ int istat_elements_checked;
+ int ilast_tracer_count;
-/* statistical parameters */
-int istat_ichoice[13][5];
-int istat_isend;
-int istat_iempty;
-int istat1;
-int istat_elements_checked;
-int ilast_tracer_count;
-
-/* Mesh information */
+ /* Mesh information */
double xcap[13][5];
double ycap[13][5];
double zcap[13][5];
@@ -130,22 +123,16 @@
double sin_phi[13][5];
-/* compositional information */
-int *ieltrac[13];
-double *celtrac[13];
-double *comp_el[13];
-double *comp_node[13];
+ /* gnomonic shape functions */
+ double *UV[13][3];
+ double cos_theta_f;
+ double sin_theta_f;
+ double *shape_coefs[13][3][10];
-/* gnomonic shape functions */
-double *UV[13][3];
-double cos_theta_f;
-double sin_theta_f;
-double *shape_coefs[13][3][10];
+ double *V0_cart[13][4];
-double *V0_cart[13][4];
+ double initial_bulk_composition;
+ double bulk_composition;
+ double error_fraction;
-double initial_bulk_composition;
-double bulk_composition;
-double error_fraction;
-
};
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-11 04:26:46 UTC (rev 6221)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-11 05:20:43 UTC (rev 6222)
@@ -600,9 +600,8 @@
}
}
+ getIntProperty(properties, "tracer_types", E->trace.ntypes, fp);
- getDoubleProperty(properties, "z_interface", E->composition.z_interface, fp);
-
getIntProperty(properties, "tracer_advection_scheme", E->trace.itracer_advection_scheme, fp);
getIntProperty(properties, "tracer_interpolation_scheme", E->trace.itracer_interpolation_scheme, fp);
@@ -612,14 +611,22 @@
E->trace.delphi[0] = tmp;
getIntProperty(properties, "analytical_tracer_test", E->trace.ianalytical_tracer_test, fp);
+
getIntProperty(properties, "chemical_buoyancy",
E->composition.ichemical_buoyancy, fp);
- getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
- getDoubleProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, fp);
- getIntProperty(properties, "reset_initial_composition", E->composition.ireset_initial_composition, fp);
- getIntProperty(properties, "compositional_rheology", E->trace.icompositional_rheology, fp);
- getDoubleProperty(properties, "compositional_prefactor", E->composition.compositional_rheology_prefactor, 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);
+ getIntProperty(properties, "reset_initial_composition", E->composition.ireset_initial_composition, fp);
+ getDoubleProperty(properties, "z_interface", E->composition.z_interface, fp);
+ }
+
+ getIntProperty(properties, "compositional_rheology", E->composition.icompositional_rheology, fp);
+
+ if (E->composition.icompositional_rheology==1) {
+ getDoubleProperty(properties, "compositional_prefactor", E->composition.compositional_rheology_prefactor, fp);
+ }
}
PUTS(("\n"));
More information about the cig-commits
mailing list