[cig-commits] r18782 - mc/3D/CitcomS/branches/eheien/lib
emheien at geodynamics.org
emheien at geodynamics.org
Tue Jul 19 11:34:06 PDT 2011
Author: emheien
Date: 2011-07-19 11:34:06 -0700 (Tue, 19 Jul 2011)
New Revision: 18782
Added:
mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp
Removed:
mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp
Modified:
mc/3D/CitcomS/branches/eheien/lib/Composition_related.c
mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/eheien/lib/Output.c
mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c
mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c
mc/3D/CitcomS/branches/eheien/lib/prototypes.h
mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h
Log:
Reverted development code to fix bugs
Modified: mc/3D/CitcomS/branches/eheien/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Composition_related.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Composition_related.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -35,7 +35,6 @@
static void allocate_composition_memory(struct All_variables *E);
static void compute_elemental_composition_ratio_method(struct All_variables *E);
-static void compute_elemental_composition_absolute_method(struct All_variables *E);
static void init_bulk_composition(struct All_variables *E);
static void check_initial_composition(struct All_variables *E);
static void fill_composition_from_neighbors(struct All_variables *E);
@@ -56,11 +55,11 @@
/* ibuoy_type=1 (ratio method) */
input_int("buoy_type",&(E->composition.ibuoy_type),"1,0,nomax",m);
- /* if (E->composition.ibuoy_type!=1) {
+ if (E->composition.ibuoy_type!=1) {
fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
fflush(stderr);
parallel_process_termination();
- } */
+ }
if (E->composition.ibuoy_type==0)
E->composition.ncomp = E->trace.nflavors;
@@ -162,6 +161,9 @@
void fill_composition(struct All_variables *E)
{
+ /* XXX: Currently, only the ratio method works here. */
+ /* Will have to come back here to include the absolute method. */
+
/* ratio method */
if (E->composition.ibuoy_type==1) {
@@ -169,15 +171,12 @@
}
/* absolute method */
- if (E->composition.ibuoy_type==0) {
- compute_elemental_composition_absolute_method(E);
- }
- /* if (E->composition.ibuoy_type!=1) {
+ if (E->composition.ibuoy_type!=1) {
fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
fflush(E->trace.fpt);
exit(10);
- } */
+ }
/* Map elemental composition to nodal points */
@@ -242,21 +241,20 @@
void init_composition(struct All_variables *E)
{
+ /* XXX: Currently, only the ratio method works here. */
+ /* Will have to come back here to include the absolute method. */
+
/* ratio method */
if (E->composition.ibuoy_type==1) {
compute_elemental_composition_ratio_method(E);
}
/* absolute method */
- if (E->composition.ibuoy_type==0) {
- compute_elemental_composition_absolute_method(E);
- }
-
- /* if (E->composition.ibuoy_type!=1) {
+ if (E->composition.ibuoy_type!=1) {
fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
fflush(E->trace.fpt);
exit(10);
- } */
+ }
/* for empty elements */
check_initial_composition(E);
@@ -303,7 +301,7 @@
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.num_tracer_flavors[j][flavor][e];
+ numtracers += E->trace.ntracer_flavor[j][flavor][e];
/* Check for empty entries and compute ratio. */
/* If no tracers are in an element, skip this element, */
@@ -317,7 +315,7 @@
for(i=0;i<E->composition.ncomp;i++) {
flavor = i + 1;
E->composition.comp_el[j][i][e] =
- E->trace.num_tracer_flavors[j][flavor][e] / (double)numtracers;
+ E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
}
}
@@ -338,70 +336,6 @@
return;
}
-
-
-/*********** COMPUTE ELEMENTAL ABSOLUTE METHOD *************************/
-/* */
-/* This function computes the composition per element. */
-/* The concentration of material i in an element is */
-/* defined as: */
-/* tracer i mass * (# of tracers of flavor i) / ( volume of element) */
-/* */
-/* where tracer i mass = tot tracer volume / tot num of tracers */
-/* */
-/* Added by DJB 01/14/10. Known caveats: (will address at some point) */
-/* 1) XXX: needs better error handling */
-/* 2) XXX: only tested with one flavor of tracer (flavor 0) */
-/* 3) XXX: only tested by reading in a file of tracer locations */
-/* 4) XXX: the initial volume of the tracer domain is HARD-CODED! */
-
-static void compute_elemental_composition_absolute_method(struct All_variables *E)
-{
- int i, j, e, flavor, numtracers;
- double domain_volume;
- float comp;
- float one = 1.0;
-
- /* This needs to be `manually' changed for the total volume */
- /* occupied by your tracers */
- domain_volume = 1e-2;
-
- for (j=1; j<=E->sphere.caps_per_proc; j++) {
- for (e=1; e<=E->lmesh.nel; e++) {
- numtracers = 0;
- for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.num_tracer_flavors[j][flavor][e];
-
- /* Check for empty entries */
- /* If no tracers are in an element, comp = 0.0 (i.e. is ambient) */
- if (numtracers == 0) {
- for(i=0;i<E->composition.ncomp;i++) {
- E->composition.comp_el[j][i][e] = 0.0;
- }
- continue;
- }
-
- /* Composition is proportional to (local) tracer density */
- for(i=0;i<E->composition.ncomp;i++) {
- flavor = i;
- comp =
- E->trace.num_tracer_flavors[j][flavor][e] / E->eco[j][e].area
- * domain_volume / E->trace.tracers[j].size();
-
- /* truncate composition at 1.0 */
- /* This violates mass conservation but prevents unphysical C */
- /* XXX: make truncation a switch for the user to specify */
- E->composition.comp_el[j][i][e] = fmin(comp,one);
-
- }
- }
- }
-
- return;
-}
-
-
-
/********** MAP COMPOSITION TO NODES ****************/
/* */
@@ -503,7 +437,7 @@
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.num_tracer_flavors[j][flavor][e];
+ numtracers += E->trace.ntracer_flavor[j][flavor][e];
if (numtracers == 0)
is_empty[e] = 1;
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -158,6 +158,24 @@
/* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
+ /* Determine number of tracer quantities */
+
+ /* advection_quantites - those needed for advection */
+ E->trace.number_of_basic_quantities=12;
+
+ /* extra_quantities - used for flavors, composition, etc. */
+ /* (can be increased for additional science i.e. tracing chemistry */
+
+ E->trace.number_of_extra_quantities = 0;
+ if (E->trace.nflavors > 0)
+ E->trace.number_of_extra_quantities += 1;
+
+
+ E->trace.number_of_tracer_quantities =
+ E->trace.number_of_basic_quantities +
+ E->trace.number_of_extra_quantities;
+
+
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
@@ -166,6 +184,22 @@
/* Some error control regarding size of pointer arrays */
+ 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);
@@ -254,6 +288,7 @@
double *receive_z[13][3];
double *REC[13];
+ void expand_tracer_arrays();
int icheck_that_processor_shell();
double CPU_time0();
@@ -277,7 +312,7 @@
fprintf(E->trace.fpt, "Entering lost_souls()\n");
- E->trace.istat_isend=E->trace.escaped_tracers[j].size();
+ E->trace.istat_isend=E->trace.ilater[j];
/** debug **
for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
@@ -292,7 +327,7 @@
/* initialize isend and ireceive */
/* # of neighbors in the horizontal plane */
- isize[j]=E->trace.escaped_tracers[j].size()*(12+1);
+ isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
for (kk=0;kk<=num_ngb;kk++) isend[j][kk]=0;
for (kk=0;kk<=num_ngb;kk++) ireceive[j][kk]=0;
@@ -372,7 +407,7 @@
/* Allocate memory in receive arrays */
for (ithatcap=0;ithatcap<=num_ngb;ithatcap++) {
- isize[j]=ireceive[j][ithatcap]*(12+1);
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
itemp_size=fmax(1,isize[j]);
@@ -393,7 +428,7 @@
if (E->parallel.nprocz>1) {
ithatcap=ithiscap;
- isize[j]=isend[j][ithatcap]*(12+1);
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
for (mm=0;mm<isize[j];mm++) {
receive[j][ithatcap][mm]=send[j][ithatcap][mm];
}
@@ -405,12 +440,12 @@
for (kk=1;kk<=num_ngb;kk++) {
idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- isize[j]=isend[j][kk]*(12+1);
+ isize[j]=isend[j][kk]*E->trace.number_of_tracer_quantities;
MPI_Isend(send[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
- isize[j]=ireceive[j][kk]*(12+1);
+ isize[j]=ireceive[j][kk]*E->trace.number_of_tracer_quantities;
MPI_Irecv(receive[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
@@ -440,7 +475,7 @@
/* Allocate Memory for REC array */
- isize[j]=isum[j]*(12+1);
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
isize[j]=fmax(isize[j],1);
if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
@@ -460,9 +495,9 @@
for (pp=0;pp<ireceive[j][ithatcap];pp++) {
irec[j]++;
- ipos=pp*(12+1);
+ ipos=pp*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<12+1;mm++) {
+ for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
ipos2=ipos+mm;
REC[j][irec_position]=receive[j][ithatcap][ipos2];
@@ -485,7 +520,7 @@
/* (No dynamic reallocation of send_z necessary) */
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=itracers_subject_to_vertical_transport[j]*(12+1);
+ isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
isize[j]=fmax(isize[j],1);
if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
@@ -511,7 +546,7 @@
num_tracers=irec[j];
for (kk=1;kk<=num_tracers;kk++) {
- ireceive_position=it*(12+1);
+ ireceive_position=it*E->trace.number_of_tracer_quantities;
it++;
irad=ireceive_position+2;
@@ -526,12 +561,12 @@
if (ival==1) {
- isend_position=isend_z[j][ivertical_neighbor]*(12+1);
+ isend_position=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
isend_z[j][ivertical_neighbor]++;
- ilast_receiver_position=(irec[j]-1)*(12+1);
+ ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<12+1;mm++) {
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
ipos=ireceive_position+mm;
ipos2=isend_position+mm;
@@ -591,7 +626,7 @@
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=ireceive_z[j][kk]*(12+1);
+ isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
isize[j]=fmax(isize[j],1);
if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
@@ -608,12 +643,12 @@
idestination_proc = E->parallel.PROCESSORz[lev].pass[kk];
- isize_send=isend_z[j][kk]*(12+1);
+ isize_send=isend_z[j][kk]*E->trace.number_of_tracer_quantities;
MPI_Isend(send_z[j][kk],isize_send,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
- isize_receive=ireceive_z[j][kk]*(12+1);
+ isize_receive=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
MPI_Irecv(receive_z[j][kk],isize_receive,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
@@ -635,7 +670,7 @@
isum[j]=isum[j]+irec[j];
- isize[j]=isum[j]*(12+1);
+ 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) {
@@ -651,11 +686,11 @@
for (kk=0;kk<ireceive_z[j][ivertical_neighbor];kk++) {
- irec_position=irec[j]*(12+1);
+ irec_position=irec[j]*E->trace.number_of_tracer_quantities;
irec[j]++;
- ireceive_position=kk*(12+1);
+ ireceive_position=kk*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<12+1;mm++) {
+ for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
}
}
@@ -676,22 +711,41 @@
for (kk=0;kk<irec[j];kk++) {
- Tracer new_tracer;
+ E->trace.ntracers[j]++;
- ireceive_position=kk*(12+1);
- new_tracer.readFromMem(&REC[j][ireceive_position]);
+ if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
- iel=(E->trace.iget_element)(E,j,-99,new_tracer.x(),new_tracer.y(),new_tracer.z(),new_tracer.theta(),new_tracer.phi(),new_tracer.rad());
+ ireceive_position=kk*E->trace.number_of_tracer_quantities;
+ for (mm=0;mm<E->trace.number_of_basic_quantities;mm++) {
+ ipos=ireceive_position+mm;
+
+ E->trace.basicq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
+ }
+ for (mm=0;mm<E->trace.number_of_extra_quantities;mm++) {
+ ipos=ireceive_position+E->trace.number_of_basic_quantities+mm;
+
+ E->trace.extraq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
+ }
+
+ 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=(E->trace.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",new_tracer.x(),new_tracer.y(),new_tracer.z(),new_tracer.theta(),new_tracer.phi(),new_tracer.rad());
+ 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);
}
- new_tracer.set_ielement(iel);
- E->trace.tracers[j].push_back(new_tracer);
+ E->trace.ielement[j][E->trace.ntracers[j]]=iel;
}
if(E->control.verbose){
@@ -724,21 +778,22 @@
{
const int j = 1;
int kk, pp;
- int ithatcap, icheck;
+ int numtracers, ithatcap, icheck;
int isend_position, ipos;
int lev = E->mesh.levmax;
double theta, phi, rad;
double x, y, z;
- TracerList::iterator tr;
/* transfer tracers from rlater to send */
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- rad=tr->rad();
- x=tr->x();
- y=tr->y();
- z=tr->z();
+ 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];
+
/* first check same cap if nprocz>1 */
if (E->parallel.nprocz>1) {
@@ -774,9 +829,12 @@
/* assign tracer to send */
- isend_position=(isend[j][ithatcap]-1)*(12+1);
+ isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
- tr->writeToMem(&send[j][ithatcap][isend_position]);
+ 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 */
@@ -833,6 +891,9 @@
const double eps=-1e-4;
+ void sphere_to_cart();
+
+
/* find u and v using spherical coordinates */
spherical_to_uv(E,j,theta,phi,&u,&v);
@@ -967,9 +1028,10 @@
/* 5 6 5 7 */
/* 6 7 6 8 */
-CartesianCoord full_get_velocity(struct All_variables *E,
+void full_get_velocity(struct All_variables *E,
int j, int nelem,
- double theta, double phi, double rad)
+ double theta, double phi, double rad,
+ double *velocity_vector)
{
int iwedge;
const int sphere_key = 0;
@@ -978,6 +1040,8 @@
double VV[4][9];
double vx[7],vy[7],vz[7];
+ void velo_from_element_d();
+
full_get_shape_functions(E, shape, nelem, theta, phi, rad);
iwedge=shape[0];
@@ -1029,11 +1093,16 @@
vz[6]=VV[3][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];
+ 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];
+ 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];
- return CartesianCoord(vx[1]*shape[1]+vx[2]*shape[2]+shape[3]*vx[3]+vx[4]*shape[4]+vx[5]*shape[5]+shape[6]*vx[6],
- vy[1]*shape[1]+vy[2]*shape[2]+shape[3]*vy[3]+vy[4]*shape[4]+vy[5]*shape[5]+shape[6]*vy[6],
- vz[1]*shape[1]+vz[2]*shape[2]+shape[3]*vz[3]+vz[4]*shape[4]+vz[5]*shape[5]+shape[6]*vz[6]);
+
+ return;
}
/***************************************************************/
@@ -1238,6 +1307,8 @@
double *fmin_array;
double *fmax_array;
+ void sphere_to_cart();
+
const double two_pi=2.0*M_PI;
elz=E->lmesh.elz;
@@ -1899,14 +1970,12 @@
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
- /*
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);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
E->trace.number_of_tracer_quantities);
- */
/* analytical test */
@@ -2452,11 +2521,34 @@
return;
}
+/******************************************************************/
+/* FIX THETA PHI */
+/* */
+/* This function constrains the value of theta to be */
+/* between 0 and PI, and */
+/* this function constrains the value of phi to be */
+/* between 0 and 2 PI */
+/* */
+static void fix_theta_phi(double *theta, double *phi)
+{
+ const double two_pi=2.0*M_PI;
+
+ fix_angle(theta);
+
+ if (*theta > M_PI) {
+ *theta = two_pi - *theta;
+ *phi += M_PI;
+ }
+
+ fix_angle(phi);
+
+ return;
+}
+
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
/* Returns -99 if not in this cap. */
-/* Returns -1 if in this cap but cannot find the element. */
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
@@ -2465,6 +2557,7 @@
double x, double y, double z,
double theta, double phi, double rad)
{
+ int icheck_processor_shell();
int iregel;
int iel;
int ntheta,nphi;
@@ -2645,7 +2738,7 @@
fprintf(E->trace.fpt,"x,y,z,theta,phi,iregel %.15e %.15e %.15e %.15e %.15e %d\n",
x,y,z,theta,phi,iregel);
fflush(E->trace.fpt);
- return -1;
+ exit(10);
foundit:
@@ -2952,10 +3045,14 @@
/* This function makes sure the particle is within the sphere, and */
/* phi and theta are within the proper degree range. */
-void full_keep_within_bounds(struct All_variables *E, CartesianCoord &cc, SphericalCoord &sc)
+void full_keep_within_bounds(struct All_variables *E,
+ double *x, double *y, double *z,
+ double *theta, double *phi, double *rad)
{
- sc.constrainThetaPhi();
- //fix_radius(E,rad,theta,phi,x,y,z);
+ fix_theta_phi(theta, phi);
+ fix_radius(E,rad,theta,phi,x,y,z);
+
+ return;
}
@@ -2967,6 +3064,7 @@
/* a test function (in "analytical_test_function"). */
void analytical_test(struct All_variables *E)
+
{
#if 0
int kk,pp;
@@ -2999,6 +3097,7 @@
void predict_tracers();
void correct_tracers();
void analytical_runge_kutte();
+ void sphere_to_cart();
fprintf(E->trace.fpt,"Starting Analytical Test\n");
@@ -3223,7 +3322,6 @@
/*************** ANALYTICAL RUNGE KUTTE ******************/
/* */
void analytical_runge_kutte(struct All_variables *E, int nsteps, double dt, double *x0_s, double *x0_c, double *xf_s, double *xf_c, double *vec)
-
{
int kk;
@@ -3241,6 +3339,8 @@
double time;
double path,dtpath;
+ void sphere_to_cart();
+ void cart_to_sphere();
void analytical_test_function();
theta_0=x0_s[1];
Modified: mc/3D/CitcomS/branches/eheien/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Output.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Output.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -58,14 +58,7 @@
extern void get_STD_topo(struct All_variables *, float**, float**,
float**, float**, int);
extern void get_CBF_topo(struct All_variables *, float**, float**);
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
-#include "anisotropic_viscosity.h"
-void output_avisc(struct All_variables *, int);
-#endif
-
-
-
/**********************************************************************/
void output_common_input(struct All_variables *E)
@@ -112,11 +105,7 @@
output_velo(E, cycles);
output_visc(E, cycles);
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
- output_avisc(E, cycles);
-#endif
-
output_surf_botm(E, cycles);
/* optional output below */
@@ -322,30 +311,7 @@
return;
}
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
-void output_avisc(struct All_variables *E, int cycles)
-{
- int i, j;
- char output_file[255];
- FILE *fp1;
- int lev = E->mesh.levmax;
- if(E->viscosity.allow_anisotropic_viscosity){
- sprintf(output_file,"%s.avisc.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
- fp1 = output_open(output_file, "w");
- 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,"%.4e %.4e %.4e %.4e\n",E->VI2[lev][j][i],E->VIn1[lev][j][i],E->VIn2[lev][j][i],E->VIn3[lev][j][i]);
- }
-
- fclose(fp1);
- }
- return;
-}
-#endif
-
void output_velo(struct All_variables *E, int cycles)
{
int i, j;
@@ -552,8 +518,7 @@
void output_seismic(struct All_variables *E, int cycles)
{
void get_prem(double, double*, double*, double*);
- void compute_seismic_model(struct All_variables*,
- double*, double*, double*);
+ void compute_seismic_model(const struct All_variables*, double*, double*, double*);
char output_file[255];
FILE* fp;
@@ -562,9 +527,9 @@
double *rho, *vp, *vs;
const int len = E->lmesh.nno;
- rho = (double *)malloc(len * sizeof(double));
- vp = (double *)malloc(len * sizeof(double));
- vs = (double *)malloc(len * sizeof(double));
+ rho = (double*)malloc(len * sizeof(double));
+ vp = (double*)malloc(len * sizeof(double));
+ vs = (double*)malloc(len * sizeof(double));
if(rho==NULL || vp==NULL || vs==NULL) {
fprintf(stderr, "Error while allocating memory\n");
abort();
@@ -659,28 +624,28 @@
int i, j, n, ncolumns;
char output_file[255];
FILE *fp1;
- TracerList::iterator tr;
-
sprintf(output_file,"%s.tracer.%d.%d", E->control.data_file,
E->parallel.me, cycles);
fp1 = output_open(output_file, "w");
- //ncolumns = 3 + E->trace.number_of_extra_quantities;
- ncolumns = 3 + 1;
+ 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.tracers[j].size(),
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
ncolumns, E->monitor.elapsed_time);
- for(tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ for(n=1;n<=E->trace.ntracers[j];n++) {
/* write basic quantities (coordinate) */
- fprintf(fp1,"%.12e %.12e %.12e %.12e",
- tr->theta(),
- tr->phi(),
- tr->rad(),
- tr->flavor());
+ 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");
}
Modified: mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -117,10 +117,6 @@
int open_file_zipped(char *, FILE **,struct All_variables *);
void gzip_file(char *);
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
-#include "anisotropic_viscosity.h"
-void gzdir_output_avisc(struct All_variables *, int);
-#endif
extern void temperatures_conform_bcs(struct All_variables *);
extern void myerror(struct All_variables *,char *);
@@ -167,9 +163,6 @@
else new VTK output won't
work */
gzdir_output_visc(E, out_cycles);
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
- gzdir_output_avisc(E, out_cycles);
-#endif
gzdir_output_surf_botm(E, out_cycles);
@@ -286,7 +279,6 @@
/* write nodal coordinate to file, big endian */
for(j=1;j <= E->sphere.caps_per_proc;j++) {
for(i=1;i <= E->lmesh.nno;i++) {
- /* cartesian coordinates */
x[0]=E->x[j][1][i];x[1]=E->x[j][2][i];x[2]=E->x[j][3][i];
if(be_write_float_to_file(x,3,fp1) != 3)
BE_WERROR;
@@ -761,70 +753,8 @@
return;
}
-#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
-/*
- anisotropic viscosity
-*/
-void gzdir_output_avisc(struct All_variables *E, int cycles)
-{
- int i, j;
- char output_file[255];
- gzFile *gz1;
- FILE *fp1;
- int lev = E->mesh.levmax;
- float ftmp;
- /* for dealing with several processors */
- MPI_Status mpi_stat;
- int mpi_rc;
- int mpi_inmsg, mpi_success_message = 1;
- if(E->viscosity.allow_anisotropic_viscosity){
-
- if(E->output.gzdir.vtk_io < 2){
- snprintf(output_file,255,
- "%s/%d/avisc.%d.%d.gz", E->control.data_dir,
- cycles,E->parallel.me, cycles);
- gz1 = gzdir_output_open(output_file,"w");
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(gz1,"%3d %7d\n",j,E->lmesh.nno);
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(gz1,"%.4e %.4e %.4e %.4e\n",E->VI2[lev][j][i],E->VIn1[lev][j][i],E->VIn2[lev][j][i],E->VIn3[lev][j][i]);
- }
-
- gzclose(gz1);
- }else{
- if(E->output.gzdir.vtk_io == 2)
- parallel_process_sync(E);
- /* new legacy VTK */
- get_vtk_filename(output_file,0,E,cycles);
- if((E->parallel.me == 0) || (E->output.gzdir.vtk_io == 3)){
- fp1 = output_open(output_file,"a");
- myfprintf(fp1,"SCALARS vis2 float 1\n");
- myfprintf(fp1,"LOOKUP_TABLE default\n");
- }else{
- /* if not first, wait for previous */
- mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, E->parallel.world, &mpi_stat);
- /* open for append */
- fp1 = output_open(output_file,"a");
- }
- for(j=1; j<= E->sphere.caps_per_proc;j++)
- for(i=1;i<=E->lmesh.nno;i++){
- ftmp = E->VI2[lev][j][i];
- if(be_write_float_to_file(&ftmp,1,fp1)!=1)
- BE_WERROR;
- }
- fclose(fp1);fflush(fp1); /* close file and flush buffer */
- if(E->output.gzdir.vtk_io == 2)
- if(E->parallel.me < E->parallel.nproc-1){
- mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, E->parallel.world);
- }
- }
- }
- return;
-}
-#endif
-
void gzdir_output_surf_botm(struct All_variables *E, int cycles)
{
int i, j, s;
@@ -984,7 +914,7 @@
snprintf(output_file,255,"%s/%d/horiz_avg.%d.%d.gz", E->control.data_dir,
cycles,E->parallel.me, cycles);
fp1=gzdir_output_open(output_file,"w");
- for(j=1;j<=E->lmesh.noz;j++) { /* format: r <T> <vh> <vr> (<C>) */
+ for(j=1;j<=E->lmesh.noz;j++) {
gzprintf(fp1,"%.4e %.4e %.4e %.4e",E->sx[1][3][j],E->Have.T[j],E->Have.V[1][j],E->Have.V[2][j]);
if (E->composition.on) {
@@ -1078,28 +1008,29 @@
int i, j, n, ncolumns;
char output_file[255];
gzFile *fp1;
- TracerList::iterator tr;
snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz",
E->control.data_dir,cycles,
E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
- ncolumns = 3 + 1;
- //ncolumns = 3 + E->trace.number_of_extra_quantities;
+ ncolumns = 3 + E->trace.number_of_extra_quantities;
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.tracers[j].size(),
+ gzprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
ncolumns, E->monitor.elapsed_time);
- for(tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ for(n=1;n<=E->trace.ntracers[j];n++) {
/* write basic quantities (coordinate) */
- gzprintf(fp1,"%9.5e %9.5e %9.5e %9.5e",
- tr->theta(),
- tr->phi(),
- tr->rad(),
- tr->flavor());
+ gzprintf(fp1,"%9.5e %9.5e %9.5e",
+ 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++) {
+ gzprintf(fp1," %9.5e", E->trace.extraq[j][i][n]);
+ }
gzprintf(fp1, "\n");
}
@@ -1513,8 +1444,8 @@
/* this should not be called with (i,i,size i) */
void be_flipit(void *d, void *s, size_t len)
{
- unsigned char *dest = (unsigned char *)d;
- unsigned char *src = (unsigned char *)s;
+ unsigned char *dest = (unsigned char*)d;
+ unsigned char *src = (unsigned char*)s;
src += len - 1;
for (; len; len--)
*dest++ = *src--;
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -49,7 +49,7 @@
static void make_mesh_ijk(struct All_variables *E);
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
- Tracer send_tracer, int j);
+ int kk, int j);
static void put_found_tracers(struct All_variables *E,
int recv_size, double *recv,
int j);
@@ -92,12 +92,48 @@
/* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
+ /* Determine number of tracer quantities */
+
+ /* advection_quantites - those needed for advection */
+ E->trace.number_of_basic_quantities=12;
+
+ /* extra_quantities - used for flavors, composition, etc. */
+ /* (can be increased for additional science i.e. tracing chemistry */
+
+ E->trace.number_of_extra_quantities = 0;
+ if (E->trace.nflavors > 0)
+ E->trace.number_of_extra_quantities += 1;
+
+
+ E->trace.number_of_tracer_quantities =
+ E->trace.number_of_basic_quantities +
+ E->trace.number_of_extra_quantities;
+
+
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on science being done */
+ /* Some error control regarding size of pointer arrays */
+
+ 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);
/* The bounding box of neiboring processors */
@@ -173,14 +209,12 @@
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
- /*
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);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
E->trace.number_of_tracer_quantities);
- */
@@ -486,14 +520,16 @@
/******** GET VELOCITY ***************************************/
-CartesianCoord regional_get_velocity(struct All_variables *E,
+void regional_get_velocity(struct All_variables *E,
int m, int nelem,
- double theta, double phi, double rad)
+ double theta, double phi, double rad,
+ double *velocity_vector)
{
+ void velo_from_element_d();
+
double shp[9], VV[4][9], tmp;
int n, d, node;
const int sphere_key = 0;
- CartesianCoord res_vector;
/* get shape functions at (theta, phi, rad) */
regional_get_shape_functions(E, shp, nelem, theta, phi, rad);
@@ -507,10 +543,13 @@
Interpolate the velocity on the tracer position
***/
- for(n=1; n<=8; n++) {
- res_vector._x += VV[1][n] * shp[n];
- res_vector._y += VV[2][n] * shp[n];
- res_vector._z += VV[3][n] * shp[n];
+ for(d=1; d<=3; d++)
+ velocity_vector[d] = 0;
+
+
+ for(d=1; d<=3; d++) {
+ for(n=1; n<=8; n++)
+ velocity_vector[d] += VV[d][n] * shp[n];
}
@@ -531,46 +570,52 @@
fflush(E->trace.fpt);
/**/
- return res_vector;
+ return;
}
-void regional_keep_within_bounds(struct All_variables *E, CartesianCoord &cc, SphericalCoord &sc)
+void regional_keep_within_bounds(struct All_variables *E,
+ double *x, double *y, double *z,
+ double *theta, double *phi, double *rad)
{
+ void sphere_to_cart();
int changed = 0;
- if (sc._theta > E->control.theta_max - E->trace.box_cushion) {
- sc._theta = E->control.theta_max - E->trace.box_cushion;
+ if (*theta > E->control.theta_max - E->trace.box_cushion) {
+ *theta = E->control.theta_max - E->trace.box_cushion;
changed = 1;
}
- if (sc._theta < E->control.theta_min + E->trace.box_cushion) {
- sc._theta = E->control.theta_min + E->trace.box_cushion;
+ if (*theta < E->control.theta_min + E->trace.box_cushion) {
+ *theta = E->control.theta_min + E->trace.box_cushion;
changed = 1;
}
- if (sc._phi > E->control.fi_max - E->trace.box_cushion) {
- sc._phi = E->control.fi_max - E->trace.box_cushion;
+ if (*phi > E->control.fi_max - E->trace.box_cushion) {
+ *phi = E->control.fi_max - E->trace.box_cushion;
changed = 1;
}
- if (sc._phi < E->control.fi_min + E->trace.box_cushion) {
- sc._phi = E->control.fi_min + E->trace.box_cushion;
+ if (*phi < E->control.fi_min + E->trace.box_cushion) {
+ *phi = E->control.fi_min + E->trace.box_cushion;
changed = 1;
}
- if (sc._rad > E->sphere.ro - E->trace.box_cushion) {
- sc._rad = E->sphere.ro - E->trace.box_cushion;
+ if (*rad > E->sphere.ro - E->trace.box_cushion) {
+ *rad = E->sphere.ro - E->trace.box_cushion;
changed = 1;
}
- if (sc._rad < E->sphere.ri + E->trace.box_cushion) {
- sc._rad = E->sphere.ri + E->trace.box_cushion;
+ if (*rad < E->sphere.ri + E->trace.box_cushion) {
+ *rad = E->sphere.ri + E->trace.box_cushion;
changed = 1;
}
- if (changed) cc = sc.toCartesian();
+ if (changed)
+ sphere_to_cart(E, *theta, *phi, *rad, x, y, z);
+
+ return;
}
@@ -589,6 +634,7 @@
double *send[2];
double *recv[2];
+ void expand_tracer_arrays();
int icheck_that_processor_shell();
int ipass;
@@ -599,7 +645,7 @@
double CPU_time0();
double begin_time = CPU_time0();
- E->trace.istat_isend = E->trace.escaped_tracers[j].size();
+ E->trace.istat_isend = E->trace.ilater[j];
/* the bounding box */
for (d=0; d<E->mesh.nsd; d++) {
@@ -643,8 +689,8 @@
/* Allocate Maximum Memory to Send Arrays */
- max_send_size = fmax(2*E->trace.escaped_tracers[j].size(), E->trace.tracers[j].size()/100);
- itemp_size = max_send_size * (12+1);
+ max_send_size = fmax(2*E->trace.ilater[j], E->trace.ntracers[j]/100);
+ itemp_size = max_send_size * E->trace.number_of_tracer_quantities;
if ((send[0] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
@@ -660,74 +706,70 @@
}
- TracerList::iterator tr;
-
for (d=0; d<E->mesh.nsd; d++) {
- int original_size = E->trace.escaped_tracers[j].size();
+ int original_size = E->trace.ilater[j];
int idb;
+ int kk = 1;
int isend[2], irecv[2];
isend[0] = isend[1] = 0;
+
/* move out-of-bound tracers to send array */
- for (tr=E->trace.escaped_tracers[j].begin();tr!=E->trace.escaped_tracers[j].end();++tr) {
+ while (kk<=E->trace.ilater[j]) {
double coord;
/* Is the tracer within the bounds in the d-th dimension */
- switch(d) {
- case 0:
- coord=tr->theta();
- break;
- case 1:
- coord=tr->phi();
- break;
- case 2:
- coord=tr->rad();
- break;
- }
+ coord = E->trace.rlater[j][d][kk];
if (coord < bounds[d][0]) {
- put_lost_tracers(E, &(isend[0]), send[0], *tr, j);
- tr = E->trace.escaped_tracers[j].erase(tr);
+ put_lost_tracers(E, &(isend[0]), send[0], kk, j);
}
else if (coord >= bounds[d][1]) {
- put_lost_tracers(E, &(isend[1]), send[1], *tr, j);
- tr = E->trace.escaped_tracers[j].erase(tr);
+ put_lost_tracers(E, &(isend[1]), send[1], kk, j);
}
else {
/* check next tracer */
+ kk++;
}
/* reallocate send if size too small */
- if ((isend[0] > max_send_size - 5) || (isend[1] > max_send_size - 5)) {
+ if ((isend[0] > max_send_size - 5) ||
+ (isend[1] > max_send_size - 5)) {
isize = max_send_size + max_send_size/4 + 10;
- itemp_size = isize * (12+1);
+ itemp_size = isize * E->trace.number_of_tracer_quantities;
- if ((send[0] = (double *)realloc(send[0], itemp_size*sizeof(double))) == NULL) {
+ if ((send[0] = (double *)realloc(send[0],
+ itemp_size*sizeof(double)))
+ == NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s4)\n");
fflush(E->trace.fpt);
exit(10);
}
- if ((send[1] = (double *)realloc(send[1], itemp_size*sizeof(double))) == NULL) {
+ if ((send[1] = (double *)realloc(send[1],
+ itemp_size*sizeof(double)))
+ == NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s5)\n");
fflush(E->trace.fpt);
exit(10);
}
- fprintf(E->trace.fpt,"Expanding physical memory of send to %d from %d\n",
+ fprintf(E->trace.fpt,"Expanding physical memory of send to "
+ "%d from %d\n",
isize, max_send_size);
max_send_size = isize;
}
- }
+ } /* end of while kk */
+
/* check the total # of tracers is conserved */
- if ((isend[0] + isend[1] + E->trace.escaped_tracers[j].size()) != original_size) {
+ if ((isend[0] + isend[1] + E->trace.ilater[j]) != original_size) {
fprintf(E->trace.fpt, "original_size: %d, rlater_size: %d, "
"send_size: %d\n",
- original_size, (int)E->trace.escaped_tracers[j].size(), kk);
+ original_size, E->trace.ilater[j], kk);
}
@@ -785,7 +827,7 @@
/* Allocate memory in receive arrays */
for (i=0; i<2; i++) {
- isize = irecv[i] * (12+1);
+ isize = irecv[i] * E->trace.number_of_tracer_quantities;
itemp_size = fmax(1, isize);
if ((recv[i] = (double *)malloc(itemp_size*sizeof(double)))
@@ -804,11 +846,11 @@
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
- isize = isend[i] * (12+1);
+ isize = isend[i] * E->trace.number_of_tracer_quantities;
MPI_Isend(send[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
- isize = irecv[i] * (12+1);
+ isize = irecv[i] * E->trace.number_of_tracer_quantities;
MPI_Irecv(recv[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
@@ -847,16 +889,14 @@
/* rlater should be empty by now */
- if (E->trace.escaped_tracers[j].size() > 0) {
+ if (E->trace.ilater[j] > 0) {
fprintf(E->trace.fpt, "Error(regional_lost_souls) lost tracers\n");
- /*
for (kk=1; kk<=E->trace.ilater[j]; kk++) {
fprintf(E->trace.fpt, "lost #%d xx=(%e, %e, %e)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
}
- */
fflush(E->trace.fpt);
exit(10);
}
@@ -874,18 +914,27 @@
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
- Tracer send_tracer, int j)
+ int kk, int j)
{
int ilast_tracer, isend_position, ipos;
int pp;
/* move the tracer from rlater to send */
- isend_position = (*send_size) * (12+1);
+ isend_position = (*send_size) * E->trace.number_of_tracer_quantities;
- send_tracer.writeToMem(&send[isend_position]);
-
+ for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
+ ipos = isend_position + pp;
+ send[ipos] = E->trace.rlater[j][pp][kk];
+ }
(*send_size)++;
+ /* eject the tracer from rlater */
+ ilast_tracer = E->trace.ilater[j];
+ for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
+ E->trace.rlater[j][pp][kk] = E->trace.rlater[j][pp][ilast_tracer];
+ }
+ E->trace.ilater[j]--;
+
return;
}
@@ -898,24 +947,25 @@
int recv_size, double *recv,
int j)
{
+ void expand_tracer_arrays();
+ void expand_later_array();
int icheck_processor_shell();
int kk, pp;
- int ipos, inside, iel;
+ int ipos, ilast, inside, iel;
+ double theta, phi, rad;
for (kk=0; kk<recv_size; kk++) {
- Tracer new_tracer;
+ ipos = kk * E->trace.number_of_tracer_quantities;
+ theta = recv[ipos];
+ phi = recv[ipos + 1];
+ rad = recv[ipos + 2];
- ipos = kk * new_tracer.size();
-
- /* found the element */
- new_tracer.readFromMem(&recv[ipos]);
-
/* check whether this tracer is inside this proc */
/* check radius first, since it is cheaper */
- inside = icheck_processor_shell(E, j, new_tracer.rad());
+ inside = icheck_processor_shell(E, j, rad);
if (inside == 1)
- inside = regional_icheck_cap(E, 0, new_tracer.theta(), new_tracer.phi(), new_tracer.rad(), new_tracer.rad());
+ inside = regional_icheck_cap(E, 0, theta, phi, rad, rad);
else
inside = 0;
@@ -927,21 +977,58 @@
/**/
if (inside) {
- iel = regional_iget_element(E, j, -99, 0, 0, 0, new_tracer.theta(), new_tracer.phi(), new_tracer.rad());
+ E->trace.ntracers[j]++;
+ ilast = E->trace.ntracers[j];
+
+ if (E->trace.ntracers[j] > (E->trace.max_ntracers[j]-5))
+ expand_tracer_arrays(E, j);
+
+ for (pp=0; pp<E->trace.number_of_basic_quantities; pp++)
+ E->trace.basicq[j][pp][ilast] = recv[ipos+pp];
+
+ ipos += E->trace.number_of_basic_quantities;
+ for (pp=0; pp<E->trace.number_of_extra_quantities; pp++)
+ E->trace.extraq[j][pp][ilast] = recv[ipos+pp];
+
+
+ /* found the element */
+ iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
+
if (iel<1) {
- fprintf(E->trace.fpt, "Error(regional lost souls) - element not here?\n");
- fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n", new_tracer.theta(), new_tracer.phi(), new_tracer.rad());
+ fprintf(E->trace.fpt, "Error(regional lost souls) - "
+ "element not here?\n");
+ fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n",
+ theta, phi, rad);
fflush(E->trace.fpt);
exit(10);
}
- new_tracer.set_ielement(iel);
- E->trace.tracers[j].push_back(new_tracer);
+ E->trace.ielement[j][ilast] = iel;
}
else {
- E->trace.escaped_tracers[j].push_back(new_tracer);
+ if (E->trace.ilatersize[j]==0) {
+
+ E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
+
+ for (kk=0;kk<E->trace.number_of_tracer_quantities;kk++) {
+ if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(put_found_tracers)-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end first particle initiating memory allocation */
+
+ E->trace.ilater[j]++;
+ ilast = E->trace.ilater[j];
+
+ if (E->trace.ilater[j] > (E->trace.ilatersize[j]-5))
+ expand_later_array(E, j);
+
+ for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++)
+ E->trace.rlater[j][pp][ilast] = recv[ipos+pp];
} /* end of if-else */
/** debug **
Deleted: mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp 2011-07-19 18:34:06 UTC (rev 18782)
@@ -1,1507 +0,0 @@
-/*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
- *<LicenseText>
- *
- * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
- * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
- * Copyright (C) 1994-2005, California Institute of Technology.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- *
- *</LicenseText>
- *
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- */
-/*
-
- Tracer_setup.c
-
- A program which initiates the distribution of tracers
- and advects those tracers in a time evolving velocity field.
- Called and used from the CitCOM finite element code.
- Written 2/96 M. Gurnis for Citcom in cartesian geometry
- Modified by Lijie in 1998 and by Vlad and Eh in 2005 for the
- regional version of CitcomS. In 2003, Allen McNamara wrote the
- tracer module for the global version of CitcomS. In 2007, Eh Tan
- merged the two versions of tracer codes together.
-*/
-
-#include <math.h>
-#include "global_defs.h"
-#include "parsing.h"
-#include "parallel_related.h"
-#include "composition_related.h"
-
-#ifdef USE_GGRD
-#include "ggrd_handling.h"
-#endif
-
-#ifdef USE_GZDIR
-int open_file_zipped(char *, FILE **,struct All_variables *);
-void gzip_file(char *);
-#endif
-
-int icheck_that_processor_shell(struct All_variables *E,
- int j, int nprocessor, double rad);
-void tracer_post_processing(struct All_variables *E);
-void count_tracers_of_flavors(struct All_variables *E);
-
-int full_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad);
-int regional_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad);
-
-static void find_tracers(struct All_variables *E);
-static void predict_tracers(struct All_variables *E);
-static void correct_tracers(struct All_variables *E);
-static void make_tracer_array(struct All_variables *E);
-static void generate_random_tracers(struct All_variables *E,
- int tracers_cap, int j);
-static void read_tracer_file(struct All_variables *E);
-static void read_old_tracer_file(struct All_variables *E);
-static void check_sum(struct All_variables *E);
-static int isum_tracers(struct All_variables *E);
-static void init_tracer_flavors(struct All_variables *E);
-static void put_away_later(struct All_variables *E, int j, int it);
-int read_double_vector(FILE *, int , double *);
-void cart_to_sphere(struct All_variables *,
- double , double , double ,
- double *, double *, double *);
-void sphere_to_cart(struct All_variables *,
- double , double , double ,
- double *, double *, double *);
-int icheck_processor_shell(struct All_variables *,
- int , double );
-
-
-void SphericalCoord::writeToMem(double *mem) const {
- mem[0] = _theta;
- mem[1] = _phi;
- mem[2] = _rad;
-}
-
-void SphericalCoord::readFromMem(const double *mem) {
- _theta = mem[0];
- _phi = mem[1];
- _rad = mem[2];
-}
-
-void CartesianCoord::writeToMem(double *mem) const {
- mem[0] = _x;
- mem[1] = _y;
- mem[2] = _z;
-}
-
-void CartesianCoord::readFromMem(const double *mem) {
- _x = mem[0];
- _y = mem[1];
- _z = mem[2];
-}
-
-
-size_t Tracer::size(void) {
- return _sc.size() // spherical coords
- + _cc.size() // Cartesian coords
- + _cc0.size() // Original coords
- + _Vc.size() // Velocity
- + 1; // flavor
-}
-
-void Tracer::writeToMem(double *mem) const {
- _sc.writeToMem(&mem[0]);
- _cc.writeToMem(&mem[3]);
- _cc0.writeToMem(&mem[6]);
- _Vc.writeToMem(&mem[9]);
- mem[12] = _flavor;
-}
-
-void Tracer::readFromMem(const double *mem) {
- _sc.readFromMem(&mem[0]);
- _cc.readFromMem(&mem[3]);
- _cc0.readFromMem(&mem[6]);
- _Vc.readFromMem(&mem[9]);
- _flavor = mem[12];
-}
-
-
-SphericalCoord CartesianCoord::toSpherical(void) const
-{
- double temp;
- double theta, phi, rad;
-
- temp=_x*_x+_y*_y;
-
- theta=atan2(sqrt(temp),_z);
- phi=myatan(_y,_x);
- rad=sqrt(temp+_z*_z);
-
- return SphericalCoord(theta, phi, rad);
-}
-
-CartesianCoord SphericalCoord::toCartesian(void) const
-{
- double sint,cost,sinf,cosf;
- double x, y, z;
-
- sint=sin(_theta);
- cost=cos(_theta);
- sinf=sin(_phi);
- cosf=cos(_phi);
-
- x=_rad*sint*cosf;
- y=_rad*sint*sinf;
- z=_rad*cost;
-
- return CartesianCoord(x,y,z);
-}
-
-
-
-const CartesianCoord CartesianCoord::operator+(const CartesianCoord &other) const {
- return CartesianCoord( this->_x+other._x,
- this->_y+other._y,
- this->_z+other._z);
-}
-
-// Multiply each element by a constant factor
-const CartesianCoord CartesianCoord::operator*(const double &val) const {
- return CartesianCoord( this->_x*val,
- this->_y*val,
- this->_z*val);
-}
-
-
-// Returns a constrained angle between 0 and 2 PI
-double SphericalCoord::constrainAngle(const double angle) const {
- return 2.0 * M_PI * floor(angle/(2.0*M_PI));
-}
-
-// Constrains theta to be between 0 and PI while simultaneously constraining phi between 0 and 2 PI
-void SphericalCoord::constrainThetaPhi(void) {
- _theta = constrainAngle(_theta);
- if (_theta > M_PI) {
- _theta = 2*M_PI - _theta;
- _phi += M_PI;
- }
- _phi = constrainAngle(_phi);
-}
-
-
-
-
-void tracer_input(struct All_variables *E)
-{
- void full_tracer_input();
- void myerror();
- void report();
- char message[100];
- int m=E->parallel.me;
- int i;
-
- input_boolean("tracer",&(E->control.tracer),"off",m);
- input_boolean("tracer_enriched",
- &(E->control.tracer_enriched),"off",m);
- if(E->control.tracer_enriched){
- if(!E->control.tracer) /* check here so that we can get away
- with only one if statement in
- Advection_diffusion */
- myerror(E,"need to switch on tracers for tracer_enriched");
-
- input_float("Q0_enriched",&(E->control.Q0ER),"0.0",m);
- snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g (only one composition!)",
- E->control.Q0,E->control.Q0ER);
- report(E,message);
- //
- // this check doesn't work at this point in the code, and we didn't want to put it into every call to
- // Advection diffusion
- //
- //if(E->composition.ncomp != 1)
- //myerror(E,"enriched tracers cannot deal with more than one composition yet");
-
- }
- if(E->control.tracer) {
-
- /* tracer_ic_method=0 (random generated array) */
- /* tracer_ic_method=1 (all proc read the same file) */
- /* tracer_ic_method=2 (each proc reads its restart file) */
- 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);
- }
- else if (E->trace.ic_method==1)
- input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
- else if (E->trace.ic_method==2) {
- /* Use 'datadir_old', 'datafile_old', and 'solution_cycles_init' */
- /* to form the filename */
- }
- else {
- fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
- parallel_process_termination();
- }
-
-
- /* How many flavors of tracers */
- /* If tracer_flavors > 0, each element will report the number of
- * tracers of each flavor inside it. This information can be used
- * later for many purposes. One of it is to compute composition,
- * either using absolute method or ratio method. */
- input_int("tracer_flavors",&(E->trace.nflavors),"0,0,nomax",m);
-
- /* 0: default from layers
- 1: from netcdf grds
-
-
- 99: from grds, overriding checkpoints during restart
- (1 and 99 require ggrd)
- */
-
- input_int("ic_method_for_flavors",
- &(E->trace.ic_method_for_flavors),"0,0,nomax",m);
-
-
- if (E->trace.nflavors > 1) {
- switch(E->trace.ic_method_for_flavors){
- /* default method */
- case 0:
- /* flavors initialized from layers */
- E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
- *sizeof(double));
- for(i=0; i<E->trace.nflavors-1; i++)
- E->trace.z_interface[i] = 0.7;
-
- input_double_vector("z_interface", E->trace.nflavors-1,
- E->trace.z_interface, m);
- break;
- /*
- two grd init method, second will override restart
- */
-#ifdef USE_GGRD
- case 1:
- case 99: /* will override restart */
- /* from grid in top n materials, this will override
- the checkpoint input */
- input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
- input_int("ictracer_grd_layers",&(E->trace.ggrd_layers),"2",m); /*
-
- >0 : which top layers to use, layer <= ictracer_grd_layers
- <0 : only use one layer layer == -ictracer_grd_layers
-
- */
- break;
-
-#endif
- default:
- fprintf(stderr,"ic_method_for_flavors %i undefined (1 and 99 only for ggrd mode)\n",E->trace.ic_method_for_flavors);
- parallel_process_termination();
- break;
- }
- }
-
- /* Warning level */
- input_boolean("itracer_warnings",&(E->trace.itracer_warnings),"on",m);
-
-
- if(E->parallel.nprocxy == 12)
- full_tracer_input(E);
-
-
- composition_input(E);
-
- }
-
- return;
-}
-
-
-void tracer_initial_settings(struct All_variables *E)
-{
- E->trace.advection_time = 0;
- E->trace.find_tracers_time = 0;
- E->trace.lost_souls_time = 0;
-
- if(E->parallel.nprocxy == 1) {
- E->problem_tracer_setup = regional_tracer_setup;
-
- E->trace.keep_within_bounds = regional_keep_within_bounds;
- E->trace.get_velocity = regional_get_velocity;
- E->trace.iget_element = regional_iget_element;
- }
- else {
- E->problem_tracer_setup = full_tracer_setup;
-
- E->trace.keep_within_bounds = full_keep_within_bounds;
- E->trace.get_velocity = full_get_velocity;
- E->trace.iget_element = full_iget_element;
- }
-}
-
-
-
-/*****************************************************************************/
-/* This function is the primary tracing routine called from Citcom.c */
-/* In this code, unlike the original 3D cartesian code, force is filled */
-/* during Stokes solution. No need to call thermal_buoyancy() after tracing. */
-
-
-void tracer_advection(struct All_variables *E)
-{
- double CPU_time0();
- double begin_time = CPU_time0();
-
- /* advect tracers */
- predict_tracers(E);
- correct_tracers(E);
-
- /* check that the number of tracers is conserved */
- check_sum(E);
-
- /* count # of tracers of each flavor */
- if (E->trace.nflavors > 0)
- count_tracers_of_flavors(E);
-
- /* update the composition field */
- if (E->composition.on) {
- fill_composition(E);
- }
-
- E->trace.advection_time += CPU_time0() - begin_time;
-
- tracer_post_processing(E);
-
- return;
-}
-
-
-
-/********* TRACER POST PROCESSING ****************************************/
-
-void tracer_post_processing(struct All_variables *E)
-{
- int i;
-
- /* reset statistical counters */
-
- E->trace.istat_isend=0;
- E->trace.istat_elements_checked=0;
- E->trace.istat1=0;
-
- /* write timing information every 20 steps */
- if ((E->monitor.solution_cycles % 20) == 0) {
- fprintf(E->trace.fpt, "STEP %d\n", E->monitor.solution_cycles);
-
- fprintf(E->trace.fpt, "Advecting tracers takes %f seconds.\n",
- E->trace.advection_time - E->trace.find_tracers_time);
- fprintf(E->trace.fpt, "Finding element takes %f seconds.\n",
- E->trace.find_tracers_time - E->trace.lost_souls_time);
- fprintf(E->trace.fpt, "Exchanging lost tracers takes %f seconds.\n",
- E->trace.lost_souls_time);
- }
-
- if(E->control.verbose){
- fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
-
- fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
-
- fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
-
- /* compositional and error fraction data files */
- //TODO: move
- if (E->composition.on) {
- 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);
- E->trace.istat_iempty=0;
-
-
- get_bulk_composition(E);
-
- if (E->parallel.me==0) {
-
- fprintf(E->fp,"composition: %e",E->monitor.elapsed_time);
- for (i=0; i<E->composition.ncomp; i++)
- fprintf(E->fp," %e", E->composition.bulk_composition[i]);
- fprintf(E->fp,"\n");
-
- fprintf(E->fp,"composition_error_fraction: %e",E->monitor.elapsed_time);
- for (i=0; i<E->composition.ncomp; i++)
- fprintf(E->fp," %e", E->composition.error_fraction[i]);
- fprintf(E->fp,"\n");
-
- }
- }
- fflush(E->trace.fpt);
- }
-
- return;
-}
-
-
-/*********** PREDICT TRACERS **********************************************/
-/* */
-/* This function predicts tracers performing an euler step */
-/* */
-/* */
-/* Note positions used in tracer array */
-/* [positions 0-5 are always fixed with current coordinates */
-/* Positions 6-8 contain original Cartesian coordinates. */
-/* Positions 9-11 contain original Cartesian velocities. */
-/* */
-
-
-static void predict_tracers(struct All_variables *E)
-{
-
- int j;
- int nelem;
-
- double dt;
- //double x0,y0,z0;
- //double theta_pred,phi_pred,rad_pred;
- //double x_pred,y_pred,z_pred;
- //double velocity_vector[4];
- TracerList::iterator tr;
-
- dt=E->advection.timestep;
-
- // For each cap
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- // For each tracer listed for the given cap
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- CartesianCoord pred, vel_vec, start_pos;
- SphericalCoord pred_sph;
-
- // Get the starting position of the tracer
- start_pos = tr->getCartesianPos();
-
- // Calculate the velocity vector where the tracer is (interpolating if necessary)
- nelem=tr->ielement();
- vel_vec = (E->trace.get_velocity)(E,j,nelem,tr->theta(),tr->phi(),tr->rad());
-
- // Use linear extrapolation to predict tracers new position
- pred = start_pos + vel_vec * dt;
- pred_sph = pred.toSpherical();
-
- /* keep in box */
-
- (E->trace.keep_within_bounds)(E, pred, pred_sph);
-
- // Update the tracer with the bounded positions
- tr->setCoords(pred_sph, pred);
-
- /* Fill in original coords (positions 6-8) */
- /* Fill in original velocities (positions 9-11) */
-
- tr->setOrigVals(start_pos, vel_vec);
-
- }
- } /* end caps */
-
- /* find new tracer elements and caps */
-
- find_tracers(E);
-
- return;
-
-}
-
-
-/*********** CORRECT TRACERS **********************************************/
-/* */
-/* This function corrects tracers using both initial and */
-/* predicted velocities */
-/* */
-/* */
-/* Note positions used in tracer array */
-/* [positions 0-5 are always fixed with current coordinates */
-/* Positions 6-8 contain original Cartesian coordinates. */
-/* Positions 9-11 contain original Cartesian velocities. */
-/* */
-
-
-static void correct_tracers(struct All_variables *E)
-{
-
- int j;
- int kk;
- int nelem;
-
-
- double dt;
- //double x0,y0,z0;
- //double theta_cor,phi_cor,rad_cor;
- //double x_cor,y_cor,z_cor;
- //double velocity_vector[4];
- //double Vx0,Vy0,Vz0;
- //double Vx_pred,Vy_pred,Vz_pred;
-
- TracerList::iterator tr;
-
- dt=E->advection.timestep;
-
-
- // For each cap
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- // For each tracer listed for the given cap
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- CartesianCoord orig_pos, orig_vel, vel_vec, new_coord;
- SphericalCoord new_coord_sph;
-
- orig_pos = tr->getOrigCartesianPos();
- orig_vel = tr->getCartesianVel();
-
- nelem=tr->ielement();
-
- vel_vec = (E->trace.get_velocity)(E,j,nelem,tr->theta(),tr->phi(),tr->rad());
-
- new_coord = orig_pos + (orig_vel + vel_vec) * (dt * 0.5);
- new_coord_sph = new_coord.toSpherical();
-
- (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
-
- /* Fill in Current Positions (other positions are no longer important) */
-
- tr->setCoords(new_coord_sph, new_coord);
-
- } /* end kk, correcting tracers */
- } /* end caps */
-
- /* find new tracer elements and caps */
-
- find_tracers(E);
-
- return;
-}
-
-
-/************ FIND TRACERS *************************************/
-/* */
-/* This function finds tracer elements and moves tracers to */
-/* other processor domains if necessary. */
-/* Array ielement is filled with elemental values. */
-
-static void find_tracers(struct All_variables *E)
-{
-
- int iel;
- int kk;
- int j;
- int iprevious_element;
-
- double theta,phi,rad;
- double x,y,z;
- double time_stat1;
- double time_stat2;
-
- void put_away_later();
- void full_lost_souls();
- void regional_lost_souls();
-
- double begin_time = CPU_time0();
-
- TracerList::iterator tr;
-
- // For each cap
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- /* initialize arrays and statistical counters */
-
- E->trace.istat1=0;
- for (kk=0;kk<=4;kk++) {
- E->trace.istat_ichoice[j][kk]=0;
- }
-
- // For each tracer listed for the given cap
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
-
- theta=tr->theta();
- phi=tr->phi();
- rad=tr->rad();
- x=tr->x();
- y=tr->y();
- z=tr->z();
-
- iprevious_element=tr->ielement();
-
- iel=(E->trace.iget_element)(E,j,iprevious_element,x,y,z,theta,phi,rad);
- /* debug *
- fprintf(E->trace.fpt,"BB. kk %d %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,iel,x,y,z,theta,phi,rad);
- fflush(E->trace.fpt);
- /**/
-
- tr->set_ielement(iel);
-
- if (iel == -99) {
- /* tracer is inside other processors */
- E->trace.escaped_tracers[j].push_back(*tr);
- tr = E->trace.tracers[j].erase(tr);
- } else if (iel == -1) {
- /* tracer is inside this processor,
- * but cannot find its element.
- * Throw away the tracer. */
-
- if (E->trace.itracer_warnings) exit(10);
-
- tr = E->trace.tracers[j].erase(tr);
- }
-
- } /* end tracers */
-
- } /* end j */
-
-
- /* Now take care of tracers that exited cap */
-
- /* REMOVE */
- /*
- parallel_process_termination();
- */
-
- if (E->parallel.nprocxy == 12)
- full_lost_souls(E);
- else
- regional_lost_souls(E);
-
- E->trace.find_tracers_time += CPU_time0() - begin_time;
-
- return;
-}
-
-
-/***********************************************************************/
-/* This function computes the number of tracers in each element. */
-/* Each tracer can be of different "flavors", which is the 0th index */
-/* of extraq. How to interprete "flavor" is left for the application. */
-
-void count_tracers_of_flavors(struct All_variables *E)
-{
- TracerList::iterator tr;
-
- int j, flavor, e;
-
- // For each cap
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- /* first zero arrays */
- for (flavor=0; flavor<E->trace.nflavors; flavor++)
- for (e=1; e<=E->lmesh.nel; e++)
- E->trace.num_tracer_flavors[j][flavor][e] = 0;
-
- /* Fill arrays */
- // For each tracer listed for the given cap
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- e = tr->ielement();
- flavor = tr->flavor();
- E->trace.num_tracer_flavors[j][flavor][e]++;
- }
- }
-
- /* debug */
- /**
- for (j=1; j<=E->sphere.caps_per_proc; j++) {
- for (e=1; e<=E->lmesh.nel; e++) {
- fprintf(E->trace.fpt, "element=%d ntracer_flaver =", e);
- for (flavor=0; flavor<E->trace.nflavors; flavor++) {
- fprintf(E->trace.fpt, " %d",
- E->trace.ntracer_flavor[j][flavor][e]);
- }
- fprintf(E->trace.fpt, "\n");
- }
- }
- fflush(E->trace.fpt);
- /**/
-
- return;
-}
-
-
-
-void initialize_tracers(struct All_variables *E)
-{
-
- if (E->trace.ic_method==0)
- make_tracer_array(E);
- else if (E->trace.ic_method==1)
- read_tracer_file(E);
- else if (E->trace.ic_method==2)
- read_old_tracer_file(E);
- else {
- fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
- /* total number of tracers */
-
- E->trace.ilast_tracer_count = isum_tracers(E);
- fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
- if(E->parallel.me==0)
- fprintf(stderr, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
-
-
- /* find elements */
-
- find_tracers(E);
-
-
- /* count # of tracers of each flavor */
-
- if (E->trace.nflavors > 0)
- count_tracers_of_flavors(E);
-
- return;
-}
-
-
-/************** MAKE TRACER ARRAY ********************************/
-/* Here, each processor will generate tracers somewhere */
-/* in the sphere - check if its in this cap - then check radial */
-
-static void make_tracer_array(struct All_variables *E)
-{
-
- int tracers_cap;
- int j;
- double processor_fraction;
-
- void init_tracer_flavors();
-
- if (E->parallel.me==0) fprintf(stderr,"Making Tracer Array\n");
-
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- processor_fraction=E->lmesh.volume/E->mesh.volume;
- 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);
-
- generate_random_tracers(E, tracers_cap, j);
-
-
-
- }/* end j */
-
-
- /* Initialize tracer flavors */
- if (E->trace.nflavors) init_tracer_flavors(E);
-
- return;
-}
-
-
-
-static void generate_random_tracers(struct All_variables *E,
- int tracers_cap, int j)
-{
- int kk;
- int ival;
- int number_of_tries=0;
- int max_tries;
-
- double x,y,z;
- //double theta,phi,rad;
- double xmin,xmax,ymin,ymax,zmin,zmax;
- double random1,random2,random3;
-
-
- /* Finding the min/max of the cartesian coordinates. */
- /* One must loop over E->X to find the min/max, since the 8 corner */
- /* nodes may not be the min/max. */
- xmin = ymin = zmin = E->sphere.ro;
- xmax = ymax = zmax = -E->sphere.ro;
- for (kk=1; kk<=E->lmesh.nno; kk++) {
- x = E->x[j][1][kk];
- y = E->x[j][2][kk];
- z = E->x[j][3][kk];
-
- xmin = ((xmin < x) ? xmin : x);
- xmax = ((xmax > x) ? xmax : x);
- ymin = ((ymin < y) ? ymin : y);
- ymax = ((ymax > y) ? ymax : y);
- zmin = ((zmin < z) ? zmin : z);
- zmax = ((zmax > z) ? zmax : z);
- }
-
- /* Tracers are placed randomly in cap */
- /* (intentionally using rand() instead of srand() )*/
- while (E->trace.tracers[j].size()<tracers_cap) {
-
- number_of_tries++;
- max_tries=100*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 1
- random1=drand48();
- random2=drand48();
- random3=drand48();
-#else /* never called */
- random1=(1.0*rand())/(1.0*RAND_MAX);
- random2=(1.0*rand())/(1.0*RAND_MAX);
- random3=(1.0*rand())/(1.0*RAND_MAX);
-#endif
-
- x=xmin+random1*(xmax-xmin);
- y=ymin+random2*(ymax-ymin);
- z=zmin+random3*(zmax-zmin);
-
- CartesianCoord new_coord(x,y,z);
- SphericalCoord new_coord_sph;
-
- /* first check if within shell */
- new_coord_sph = new_coord.toSpherical();
-
- if (new_coord_sph._rad>=E->sx[j][3][E->lmesh.noz]) continue;
- if (new_coord_sph._rad<E->sx[j][3][1]) continue;
-
-
- /* check if in current cap */
- if (E->parallel.nprocxy==1)
- ival=regional_icheck_cap(E,0,new_coord_sph._theta,new_coord_sph._phi,new_coord_sph._rad,new_coord_sph._rad);
- else
- ival=full_icheck_cap(E,0,new_coord._x,new_coord._y,new_coord._z,new_coord_sph._rad);
-
- if (ival!=1) continue;
-
- /* Made it, so record tracer information */
-
- (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
-
- E->trace.tracers[j].push_back(Tracer(new_coord_sph, new_coord));
-
- } /* end while */
-
- return;
-}
-
-
-/******** READ TRACER ARRAY *********************************************/
-/* */
-/* This function reads tracers from input file. */
-/* All processors read the same input file, then sort out which ones */
-/* belong. */
-
-static void read_tracer_file(struct All_variables *E)
-{
-
- char input_s[1000];
-
- int number_of_tracers, ncolumns;
- int kk;
- int icheck;
- int iestimate;
- int icushion;
- int i, j;
-
-
- int icheck_processor_shell();
-
- //double x,y,z;
- //double theta,phi,rad;
- double buffer[100];
-
- FILE *fptracer;
-
- fptracer=fopen(E->trace.tracer_file,"r");
-
- fgets(input_s,200,fptracer);
- if(sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns) != 2) {
- fprintf(stderr,"Error while reading file '%s'\n", E->trace.tracer_file);
- exit(8);
- }
- fprintf(E->trace.fpt,"%d Tracers, %d columns in file \n",
- number_of_tracers, ncolumns);
-
- /* some error control */
- //if (E->trace.number_of_extra_quantities+3 != ncolumns) {
- if (1+3 != ncolumns) {
- fprintf(E->trace.fpt,"ERROR(read tracer file)-wrong # of columns\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
-
- /* initially size tracer arrays to number of tracers divided by processors */
-
- icushion=100;
-
- /* for absolute tracer method */
- iestimate=number_of_tracers/E->parallel.nproc + icushion;
-
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- for (kk=1;kk<=number_of_tracers;kk++) {
- SphericalCoord new_coord_sph;
- CartesianCoord new_coord;
- int len, ncol;
- ncol = 3 + 1;
- //ncol = 3 + E->trace.number_of_extra_quantities;
-
- len = read_double_vector(fptracer, ncol, buffer);
- if (len != ncol) {
- fprintf(E->trace.fpt,"ERROR(read tracer file) - wrong input file format: %s\n", E->trace.tracer_file);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- new_coord_sph.readFromMem(buffer);
- new_coord = new_coord_sph.toCartesian();
-
- /* make sure theta, phi is in range, and radius is within bounds */
-
- (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
-
- /* check whether tracer is within processor domain */
-
- icheck=1;
- if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,new_coord_sph._rad);
- if (icheck!=1) continue;
-
- if (E->parallel.nprocxy==1)
- icheck=regional_icheck_cap(E,0,new_coord_sph._theta,new_coord_sph._phi,new_coord_sph._rad,new_coord_sph._rad);
- else
- icheck=full_icheck_cap(E,0,new_coord._x,new_coord._y,new_coord._z,new_coord_sph._rad);
-
- if (icheck==0) continue;
-
- /* if still here, tracer is in processor domain */
-
- E->trace.tracers[j].push_back(Tracer(new_coord_sph, new_coord));
-
- } /* end kk, number of tracers */
-
- fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
- (int)E->trace.tracers[j].size());
-
- /** debug **
- for (kk=1; kk<=E->trace.ntracers[j]; kk++) {
- fprintf(E->trace.fpt, "tracer#=%d sph_coord=(%g,%g,%g)", kk,
- E->trace.basicq[j][0][kk],
- E->trace.basicq[j][1][kk],
- E->trace.basicq[j][2][kk]);
- fprintf(E->trace.fpt, " extraq=");
- for (i=0; i<E->trace.number_of_extra_quantities; i++)
- fprintf(E->trace.fpt, " %g", E->trace.extraq[j][i][kk]);
- fprintf(E->trace.fpt, "\n");
- }
- fflush(E->trace.fpt);
- /**/
-
- } /* 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);
- }
-
- return;
-}
-
-
-/************** READ OLD TRACER FILE *************************************/
-/* */
-/* This function read tracers written from previous calculation */
-/* and the tracers are read as seperate files for each processor domain. */
-
-static void read_old_tracer_file(struct All_variables *E)
-{
-
- char output_file[200];
- char input_s[1000];
-
- int i,j,kk,rezip;
- int idum1,ncolumns;
- int numtracers;
-
- double rdum1;
- //double theta,phi,rad;
- //double x,y,z;
- double buffer[100];
-
- FILE *fp1;
-
- /*if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-increase size of extra[]\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }*/
-
-
- /* deal with different output formats */
-#ifdef USE_GZDIR
- if(strcmp(E->output.format, "ascii-gz") == 0){
- sprintf(output_file,"%s/%d/tracer.%d.%d",
- E->control.data_dir_old,E->monitor.solution_cycles_init,E->parallel.me,E->monitor.solution_cycles_init);
- rezip = open_file_zipped(output_file,&fp1,E);
- }else{
- 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(read_old_tracer_file)-gziped file not found %s\n",output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-#else
- 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(read_old_tracer_file)-file not found %s\n",output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
-#endif
-
- fprintf(stderr,"Read old tracers from %s\n",output_file);
-
-
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fgets(input_s,200,fp1);
- if(sscanf(input_s,"%d %d %d %lf",
- &idum1, &numtracers, &ncolumns, &rdum1) != 4) {
- fprintf(stderr,"Error while reading file '%s'\n", output_file);
- exit(8);
- }
-
-
- /* some error control */
- //if (E->trace.number_of_extra_quantities+3 != ncolumns) {
- if (1+3 != ncolumns) {
- fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-wrong # of columns\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- for (kk=1;kk<=numtracers;kk++) {
- SphericalCoord new_coord_sph;
- CartesianCoord new_coord;
- int len, ncol;
- ncol = 3 + 1;
- //ncol = 3 + E->trace.number_of_extra_quantities;
-
- len = read_double_vector(fp1, ncol, buffer);
- if (len != ncol) {
- fprintf(E->trace.fpt,"ERROR(read_old_tracer_file) - wrong input file format: %s\n", output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- new_coord_sph.readFromMem(buffer);
- new_coord = new_coord_sph.toCartesian();
-
- /* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
-
- (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
-
- E->trace.tracers[j].push_back(Tracer(new_coord_sph, new_coord));
-
- }
-
- /** debug **
- for (kk=1; kk<=E->trace.ntracers[j]; kk++) {
- fprintf(E->trace.fpt, "tracer#=%d sph_coord=(%g,%g,%g)", kk,
- E->trace.basicq[j][0][kk],
- E->trace.basicq[j][1][kk],
- E->trace.basicq[j][2][kk]);
- fprintf(E->trace.fpt, " extraq=");
- for (i=0; i<E->trace.number_of_extra_quantities; i++)
- fprintf(E->trace.fpt, " %g", E->trace.extraq[j][i][kk]);
- fprintf(E->trace.fpt, "\n");
- }
- fflush(E->trace.fpt);
- /**/
-
- fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
- fflush(E->trace.fpt);
-
- }
- fclose(fp1);
-#ifdef USE_GZDIR
- if(strcmp(E->output.format, "ascii-gz") == 0)
- if(rezip) /* rezip */
- gzip_file(output_file);
-#endif
-
- return;
-}
-
-
-
-
-
-/*********** CHECK SUM **************************************************/
-/* */
-/* This functions checks to make sure number of tracers is preserved */
-
-static void check_sum(struct All_variables *E)
-{
-
- int number, iold_number;
-
- number = isum_tracers(E);
-
- 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);
- if (E->trace.itracer_warnings)
- parallel_process_termination();
- }
-
- E->trace.ilast_tracer_count = number;
-
- return;
-}
-
-
-/************* ISUM TRACERS **********************************************/
-/* */
-/* This function uses MPI to sum all tracers and returns number of them. */
-
-static int isum_tracers(struct All_variables *E)
-{
- int imycount;
- int iallcount;
- int j;
-
- iallcount = 0;
-
- imycount = 0;
- for (j=1; j<=E->sphere.caps_per_proc; j++)
- imycount = imycount + E->trace.tracers[j].size();
-
- MPI_Allreduce(&imycount,&iallcount,1,MPI_INT,MPI_SUM,E->parallel.world);
-
- return iallcount;
-}
-
-
-
-/********** CART TO SPHERE ***********************/
-void cart_to_sphere(struct All_variables *E,
- double x, double y, double z,
- double *theta, double *phi, double *rad)
-{
-
- double temp;
-
- temp=x*x+y*y;
-
- *rad=sqrt(temp+z*z);
- *theta=atan2(sqrt(temp),z);
- *phi=myatan(y,x);
-
-
- return;
-}
-
-/********** SPHERE TO CART ***********************/
-void sphere_to_cart(struct All_variables *E,
- double theta, double phi, double rad,
- double *x, double *y, double *z)
-{
-
- double sint,cost,sinf,cosf;
- double temp;
-
- sint=sin(theta);
- cost=cos(theta);
- sinf=sin(phi);
- cosf=cos(phi);
-
- temp=rad*sint;
-
- *x=temp*cosf;
- *y=temp*sinf;
- *z=rad*cost;
-
- return;
-}
-
-
-
-static void init_tracer_flavors(struct All_variables *E)
-{
- int j, kk;
- int i;
- double flavor;
- double rad;
- TracerList::iterator tr;
-
- switch(E->trace.ic_method_for_flavors){
- case 0:
- /* ic_method_for_flavors == 0 (layered structure) */
- /* any tracer above z_interface[i] is of flavor i */
- /* any tracer below z_interface is of flavor (nflavors-1) */
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- rad = tr->rad();
-
- flavor = E->trace.nflavors - 1;
- for (i=0; i<E->trace.nflavors-1; i++) {
- if (rad > E->trace.z_interface[i]) {
- flavor = i;
- break;
- }
- }
- tr->set_flavor(flavor);
- }
- }
- break;
-
- case 1: /* from grd in top n layers */
- case 99: /* (will override restart) */
-#ifndef USE_GGRD
- fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
- E->trace.ic_method_for_flavors);
- parallel_process_termination();
-#else
- ggrd_init_tracer_flavors(E);
-#endif
- break;
-
-
- default:
-
- fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
- parallel_process_termination();
- break;
- }
-
- return;
-}
-
-
-/******************* get_neighboring_caps ************************************/
-/* */
-/* Communicate with neighboring processors to get their cap boundaries, */
-/* which is later used by (E->trace.icheck_cap)() */
-/* */
-
-void get_neighboring_caps(struct All_variables *E)
-{
- void sphere_to_cart();
-
- const int ncorners = 4; /* # of top corner nodes */
- int i, j, n, d, kk, lev, idb;
- int num_ngb, neighbor_proc, tag;
- MPI_Status status[200];
- MPI_Request request[200];
-
- int node[ncorners];
- double xx[ncorners*2], rr[12][ncorners*2];
- int nox,noy,noz;
- double x,y,z;
- double theta,phi,rad;
-
- nox=E->lmesh.nox;
- noy=E->lmesh.noy;
- noz=E->lmesh.noz;
-
- node[0]=nox*noz*(noy-1)+noz;
- node[1]=noz;
- node[2]=noz*nox;
- node[3]=noz*nox*noy;
-
- lev = E->mesh.levmax;
- tag = 45;
-
- for (j=1; j<=E->sphere.caps_per_proc; j++) {
-
- /* loop over top corners to get their coordinates */
- n = 0;
- for (i=0; i<ncorners; i++) {
- for (d=0; d<2; d++) {
- xx[n] = E->sx[j][d+1][node[i]];
- n++;
- }
- }
-
- idb = 0;
- num_ngb = E->parallel.TNUM_PASS[lev][j];
- for (kk=1; kk<=num_ngb; kk++) {
- neighbor_proc = E->parallel.PROCESSOR[lev][j].pass[kk];
-
- MPI_Isend(xx, n, MPI_DOUBLE, neighbor_proc,
- tag, E->parallel.world, &request[idb]);
- idb++;
-
- MPI_Irecv(rr[kk], n, MPI_DOUBLE, neighbor_proc,
- tag, E->parallel.world, &request[idb]);
- idb++;
- }
-
- /* Storing the current cap information */
- for (i=0; i<n; i++)
- rr[0][i] = xx[i];
-
- /* Wait for non-blocking calls to complete */
-
- MPI_Waitall(idb, request, status);
-
- /* Storing the received cap information
- * XXX: this part assumes:
- * 1) E->sphere.caps_per_proc==1
- * 2) E->mesh.nsd==3
- */
- for (kk=0; kk<=num_ngb; kk++) {
- n = 0;
- for (i=1; i<=ncorners; i++) {
- theta = rr[kk][n++];
- phi = rr[kk][n++];
- rad = E->sphere.ro;
-
- sphere_to_cart(E, theta, phi, rad, &x, &y, &z);
-
- E->trace.xcap[kk][i] = x;
- E->trace.ycap[kk][i] = y;
- E->trace.zcap[kk][i] = z;
- E->trace.theta_cap[kk][i] = theta;
- E->trace.phi_cap[kk][i] = phi;
- E->trace.rad_cap[kk][i] = rad;
- E->trace.cos_theta[kk][i] = cos(theta);
- E->trace.sin_theta[kk][i] = sin(theta);
- E->trace.cos_phi[kk][i] = cos(phi);
- E->trace.sin_phi[kk][i] = sin(phi);
- }
- } /* end kk, number of neighbors */
-
- /* debugging output *
- for (kk=0; kk<=num_ngb; kk++) {
- if (kk==0)
- neighbor_proc = E->parallel.me;
- else
- neighbor_proc = E->parallel.PROCESSOR[lev][1].pass[kk];
-
- for (i=1; i<=ncorners; i++) {
- fprintf(E->trace.fpt, "pass=%d rank=%d corner=%d "
- "sx=(%g, %g, %g)\n",
- kk, neighbor_proc, i,
- E->trace.theta_cap[kk][i],
- E->trace.phi_cap[kk][i],
- E->trace.rad_cap[kk][i]);
- }
- }
- fflush(E->trace.fpt);
- /**/
- }
-
- return;
-}
-
-
-
-
-/********** ICHECK PROCESSOR SHELL *************/
-/* returns -99 if rad is below current shell */
-/* returns 0 if rad is above current shell */
-/* returns 1 if rad is within current shell */
-/* */
-/* Shell, here, refers to processor shell */
-/* */
-/* shell is defined as bottom boundary up to */
-/* and not including the top boundary unless */
-/* the shell in question is the top shell */
-
-int icheck_processor_shell(struct All_variables *E,
- int j, double rad)
-{
-
- const int noz = E->lmesh.noz;
- const int nprocz = E->parallel.nprocz;
- double top_r, bottom_r;
-
- if (nprocz==1) return 1;
-
- top_r = E->sx[j][3][noz];
- bottom_r = E->sx[j][3][1];
-
- /* First check bottom */
-
- if (rad<bottom_r) return -99;
-
-
- /* Check top */
-
- if (rad<top_r) return 1;
-
- /* top processor */
-
- if ( (rad<=top_r) && (E->parallel.me_loc[3]==nprocz-1) ) return 1;
-
- /* If here, means point is above processor */
- return 0;
-}
-
-
-/********* ICHECK THAT PROCESSOR SHELL ********/
-/* */
-/* Checks whether a given radius is within */
-/* a given processors radial domain. */
-/* Returns 0 if not, 1 if so. */
-/* The domain is defined as including the bottom */
-/* radius, but excluding the top radius unless */
-/* we the processor domain is the one that */
-/* is at the surface (then both boundaries are */
-/* included). */
-
-int icheck_that_processor_shell(struct All_variables *E,
- int j, int nprocessor, double rad)
-{
- int icheck_processor_shell();
- int me = E->parallel.me;
-
- /* nprocessor is right on top of me */
- if (nprocessor == me+1) {
- if (icheck_processor_shell(E, j, rad) == 0) return 1;
- else return 0;
- }
-
- /* nprocessor is right on bottom of me */
- if (nprocessor == me-1) {
- if (icheck_processor_shell(E, j, rad) == -99) return 1;
- else return 0;
- }
-
- /* Shouldn't be here */
- fprintf(E->trace.fpt, "Should not be here\n");
- fprintf(E->trace.fpt, "Error(check_shell) nprocessor: %d, radius: %f\n",
- nprocessor, rad);
- fflush(E->trace.fpt);
- exit(10);
-
- return 0;
-}
-
-
Copied: mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp (from rev 18754, mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.c)
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp (rev 0)
+++ mc/3D/CitcomS/branches/eheien/lib/Tracer_setup.cpp 2011-07-19 18:34:06 UTC (rev 18782)
@@ -0,0 +1,1797 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+/*
+
+ Tracer_setup.c
+
+ A program which initiates the distribution of tracers
+ and advects those tracers in a time evolving velocity field.
+ Called and used from the CitCOM finite element code.
+ Written 2/96 M. Gurnis for Citcom in cartesian geometry
+ Modified by Lijie in 1998 and by Vlad and Eh in 2005 for the
+ regional version of CitcomS. In 2003, Allen McNamara wrote the
+ tracer module for the global version of CitcomS. In 2007, Eh Tan
+ merged the two versions of tracer codes together.
+*/
+
+#include <math.h>
+#include "global_defs.h"
+#include "parsing.h"
+#include "parallel_related.h"
+#include "composition_related.h"
+
+#ifdef USE_GGRD
+#include "ggrd_handling.h"
+#endif
+
+#ifdef USE_GZDIR
+int open_file_zipped(char *, FILE **,struct All_variables *);
+void gzip_file(char *);
+#endif
+
+int icheck_that_processor_shell(struct All_variables *E,
+ int j, int nprocessor, double rad);
+void expand_later_array(struct All_variables *E, int j);
+void expand_tracer_arrays(struct All_variables *E, int j);
+void tracer_post_processing(struct All_variables *E);
+void allocate_tracer_arrays(struct All_variables *E,
+ int j, int number_of_tracers);
+void count_tracers_of_flavors(struct All_variables *E);
+
+int full_icheck_cap(struct All_variables *E, int icap,
+ double x, double y, double z, double rad);
+int regional_icheck_cap(struct All_variables *E, int icap,
+ double x, double y, double z, double rad);
+
+static void find_tracers(struct All_variables *E);
+static void predict_tracers(struct All_variables *E);
+static void correct_tracers(struct All_variables *E);
+static void make_tracer_array(struct All_variables *E);
+static void generate_random_tracers(struct All_variables *E,
+ int tracers_cap, int j);
+static void read_tracer_file(struct All_variables *E);
+static void read_old_tracer_file(struct All_variables *E);
+static void check_sum(struct All_variables *E);
+static int isum_tracers(struct All_variables *E);
+static void init_tracer_flavors(struct All_variables *E);
+static void reduce_tracer_arrays(struct All_variables *E);
+static void put_away_later(struct All_variables *E, int j, int it);
+static void eject_tracer(struct All_variables *E, int j, int it);
+int read_double_vector(FILE *, int , double *);
+void cart_to_sphere(struct All_variables *,
+ double , double , double ,
+ double *, double *, double *);
+void sphere_to_cart(struct All_variables *,
+ double , double , double ,
+ double *, double *, double *);
+int icheck_processor_shell(struct All_variables *,
+ int , double );
+
+
+
+void tracer_input(struct All_variables *E)
+{
+ void full_tracer_input();
+ void myerror();
+ void report();
+ char message[100];
+ int m=E->parallel.me;
+ int i;
+
+ input_boolean("tracer",&(E->control.tracer),"off",m);
+ input_boolean("tracer_enriched",
+ &(E->control.tracer_enriched),"off",m);
+ if(E->control.tracer_enriched){
+ if(!E->control.tracer) /* check here so that we can get away
+ with only one if statement in
+ Advection_diffusion */
+ myerror(E,"need to switch on tracers for tracer_enriched");
+
+ input_float("Q0_enriched",&(E->control.Q0ER),"0.0",m);
+ snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g (only one composition!)",
+ E->control.Q0,E->control.Q0ER);
+ report(E,message);
+ //
+ // this check doesn't work at this point in the code, and we didn't want to put it into every call to
+ // Advection diffusion
+ //
+ //if(E->composition.ncomp != 1)
+ //myerror(E,"enriched tracers cannot deal with more than one composition yet");
+
+ }
+ if(E->control.tracer) {
+
+ /* tracer_ic_method=0 (random generated array) */
+ /* tracer_ic_method=1 (all proc read the same file) */
+ /* tracer_ic_method=2 (each proc reads its restart file) */
+ 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);
+ }
+ else if (E->trace.ic_method==1)
+ input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
+ else if (E->trace.ic_method==2) {
+ /* Use 'datadir_old', 'datafile_old', and 'solution_cycles_init' */
+ /* to form the filename */
+ }
+ else {
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ parallel_process_termination();
+ }
+
+
+ /* How many flavors of tracers */
+ /* If tracer_flavors > 0, each element will report the number of
+ * tracers of each flavor inside it. This information can be used
+ * later for many purposes. One of it is to compute composition,
+ * either using absolute method or ratio method. */
+ input_int("tracer_flavors",&(E->trace.nflavors),"0,0,nomax",m);
+
+ /* 0: default from layers
+ 1: from netcdf grds
+
+
+ 99: from grds, overriding checkpoints during restart
+ (1 and 99 require ggrd)
+ */
+
+ input_int("ic_method_for_flavors",
+ &(E->trace.ic_method_for_flavors),"0,0,nomax",m);
+
+
+ if (E->trace.nflavors > 1) {
+ switch(E->trace.ic_method_for_flavors){
+ /* default method */
+ case 0:
+ /* flavors initialized from layers */
+ E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
+ *sizeof(double));
+ for(i=0; i<E->trace.nflavors-1; i++)
+ E->trace.z_interface[i] = 0.7;
+
+ input_double_vector("z_interface", E->trace.nflavors-1,
+ E->trace.z_interface, m);
+ break;
+ /*
+ two grd init method, second will override restart
+ */
+#ifdef USE_GGRD
+ case 1:
+ case 99: /* will override restart */
+ /* from grid in top n materials, this will override
+ the checkpoint input */
+ input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
+ input_int("ictracer_grd_layers",&(E->trace.ggrd_layers),"2",m); /*
+
+ >0 : which top layers to use, layer <= ictracer_grd_layers
+ <0 : only use one layer layer == -ictracer_grd_layers
+
+ */
+ break;
+
+#endif
+ default:
+ fprintf(stderr,"ic_method_for_flavors %i undefined (1 and 99 only for ggrd mode)\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
+ }
+ }
+
+ /* Warning level */
+ input_boolean("itracer_warnings",&(E->trace.itracer_warnings),"on",m);
+
+
+ if(E->parallel.nprocxy == 12)
+ full_tracer_input(E);
+
+
+ composition_input(E);
+
+ }
+
+ return;
+}
+
+
+void tracer_initial_settings(struct All_variables *E)
+{
+ E->trace.advection_time = 0;
+ E->trace.find_tracers_time = 0;
+ E->trace.lost_souls_time = 0;
+
+ if(E->parallel.nprocxy == 1) {
+ E->problem_tracer_setup = regional_tracer_setup;
+
+ E->trace.keep_within_bounds = regional_keep_within_bounds;
+ E->trace.get_velocity = regional_get_velocity;
+ E->trace.iget_element = regional_iget_element;
+ }
+ else {
+ E->problem_tracer_setup = full_tracer_setup;
+
+ E->trace.keep_within_bounds = full_keep_within_bounds;
+ E->trace.get_velocity = full_get_velocity;
+ E->trace.iget_element = full_iget_element;
+ }
+}
+
+
+
+/*****************************************************************************/
+/* This function is the primary tracing routine called from Citcom.c */
+/* In this code, unlike the original 3D cartesian code, force is filled */
+/* during Stokes solution. No need to call thermal_buoyancy() after tracing. */
+
+
+void tracer_advection(struct All_variables *E)
+{
+ double CPU_time0();
+ double begin_time = CPU_time0();
+
+ /* advect tracers */
+ predict_tracers(E);
+ correct_tracers(E);
+
+ /* check that the number of tracers is conserved */
+ check_sum(E);
+
+ /* count # of tracers of each flavor */
+ if (E->trace.nflavors > 0)
+ count_tracers_of_flavors(E);
+
+ /* update the composition field */
+ if (E->composition.on) {
+ fill_composition(E);
+ }
+
+ E->trace.advection_time += CPU_time0() - begin_time;
+
+ tracer_post_processing(E);
+
+ return;
+}
+
+
+
+/********* TRACER POST PROCESSING ****************************************/
+
+void tracer_post_processing(struct All_variables *E)
+{
+ int i;
+
+ /* reset statistical counters */
+
+ E->trace.istat_isend=0;
+ E->trace.istat_elements_checked=0;
+ E->trace.istat1=0;
+
+ /* write timing information every 20 steps */
+ if ((E->monitor.solution_cycles % 20) == 0) {
+ fprintf(E->trace.fpt, "STEP %d\n", E->monitor.solution_cycles);
+
+ fprintf(E->trace.fpt, "Advecting tracers takes %f seconds.\n",
+ E->trace.advection_time - E->trace.find_tracers_time);
+ fprintf(E->trace.fpt, "Finding element takes %f seconds.\n",
+ E->trace.find_tracers_time - E->trace.lost_souls_time);
+ fprintf(E->trace.fpt, "Exchanging lost tracers takes %f seconds.\n",
+ E->trace.lost_souls_time);
+ }
+
+ if(E->control.verbose){
+ fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
+
+ fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
+
+ fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
+
+ /* compositional and error fraction data files */
+ //TODO: move
+ if (E->composition.on) {
+ 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);
+ E->trace.istat_iempty=0;
+
+
+ get_bulk_composition(E);
+
+ if (E->parallel.me==0) {
+
+ fprintf(E->fp,"composition: %e",E->monitor.elapsed_time);
+ for (i=0; i<E->composition.ncomp; i++)
+ fprintf(E->fp," %e", E->composition.bulk_composition[i]);
+ fprintf(E->fp,"\n");
+
+ fprintf(E->fp,"composition_error_fraction: %e",E->monitor.elapsed_time);
+ for (i=0; i<E->composition.ncomp; i++)
+ fprintf(E->fp," %e", E->composition.error_fraction[i]);
+ fprintf(E->fp,"\n");
+
+ }
+ }
+ fflush(E->trace.fpt);
+ }
+
+ return;
+}
+
+
+/*********** PREDICT TRACERS **********************************************/
+/* */
+/* This function predicts tracers performing an euler step */
+/* */
+/* */
+/* Note positions used in tracer array */
+/* [positions 0-5 are always fixed with current coordinates */
+/* Positions 6-8 contain original Cartesian coordinates. */
+/* Positions 9-11 contain original Cartesian velocities. */
+/* */
+
+
+static void predict_tracers(struct All_variables *E)
+{
+
+ int numtracers;
+ int j;
+ int kk;
+ int nelem;
+
+ double dt;
+ double theta0,phi0,rad0;
+ double x0,y0,z0;
+ double theta_pred,phi_pred,rad_pred;
+ double x_pred,y_pred,z_pred;
+ double velocity_vector[4];
+
+ void cart_to_sphere();
+
+
+ dt=E->advection.timestep;
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ numtracers=E->trace.ntracers[j];
+
+ for (kk=1;kk<=numtracers;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.ielement[j][kk];
+ (E->trace.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;
+
+
+ /* keep in box */
+
+ cart_to_sphere(E,x_pred,y_pred,z_pred,&theta_pred,&phi_pred,&rad_pred);
+ (E->trace.keep_within_bounds)(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
+
+ /* Current Coordinates are always kept in positions 0-5. */
+
+ 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) */
+
+ 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) */
+
+ 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 */
+
+ /* find new tracer elements and caps */
+
+ find_tracers(E);
+
+ return;
+
+}
+
+
+/*********** CORRECT TRACERS **********************************************/
+/* */
+/* This function corrects tracers using both initial and */
+/* predicted velocities */
+/* */
+/* */
+/* Note positions used in tracer array */
+/* [positions 0-5 are always fixed with current coordinates */
+/* Positions 6-8 contain original Cartesian coordinates. */
+/* Positions 9-11 contain original Cartesian velocities. */
+/* */
+
+
+static void correct_tracers(struct All_variables *E)
+{
+
+ int j;
+ int kk;
+ int nelem;
+
+
+ double dt;
+ double x0,y0,z0;
+ double theta_pred,phi_pred,rad_pred;
+ double x_pred,y_pred,z_pred;
+ double theta_cor,phi_cor,rad_cor;
+ double x_cor,y_cor,z_cor;
+ double velocity_vector[4];
+ double Vx0,Vy0,Vz0;
+ double Vx_pred,Vy_pred,Vz_pred;
+
+ void cart_to_sphere();
+
+
+ dt=E->advection.timestep;
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=E->trace.ntracers[j];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.basicq[j][6][kk];
+ y0=E->trace.basicq[j][7][kk];
+ z0=E->trace.basicq[j][8][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.ielement[j][kk];
+
+ (E->trace.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];
+
+ 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);
+ (E->trace.keep_within_bounds)(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
+
+ /* Fill in Current Positions (other positions are no longer important) */
+
+ 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 */
+
+ /* find new tracer elements and caps */
+
+ find_tracers(E);
+
+ return;
+}
+
+
+/************ FIND TRACERS *************************************/
+/* */
+/* This function finds tracer elements and moves tracers to */
+/* other processor domains if necessary. */
+/* Array ielement is filled with elemental values. */
+
+static void find_tracers(struct All_variables *E)
+{
+
+ int iel;
+ int kk;
+ int j;
+ int it;
+ int iprevious_element;
+ int num_tracers;
+
+ double theta,phi,rad;
+ double x,y,z;
+ double time_stat1;
+ double time_stat2;
+
+ void put_away_later();
+ void eject_tracer();
+ void reduce_tracer_arrays();
+ void sphere_to_cart();
+ void full_lost_souls();
+ void regional_lost_souls();
+
+ double CPU_time0();
+ double begin_time = CPU_time0();
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+
+ /* initialize arrays and statistical counters */
+
+ E->trace.ilater[j]=E->trace.ilatersize[j]=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 */
+
+ it=0;
+ num_tracers=E->trace.ntracers[j];
+
+ for (kk=1;kk<=num_tracers;kk++) {
+
+ 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.ielement[j][it];
+
+ iel=(E->trace.iget_element)(E,j,iprevious_element,x,y,z,theta,phi,rad);
+ /* debug *
+ fprintf(E->trace.fpt,"BB. kk %d %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,iel,x,y,z,theta,phi,rad);
+ fflush(E->trace.fpt);
+ /**/
+
+ E->trace.ielement[j][it]=iel;
+
+ if (iel<0) {
+ put_away_later(E,j,it);
+ eject_tracer(E,j,it);
+ it--;
+ }
+
+ } /* end tracers */
+
+ } /* end j */
+
+
+ /* Now take care of tracers that exited cap */
+
+ /* REMOVE */
+ /*
+ parallel_process_termination();
+ */
+
+ if (E->parallel.nprocxy == 12)
+ full_lost_souls(E);
+ else
+ regional_lost_souls(E);
+
+ /* Free later arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ if (E->trace.ilatersize[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 */
+
+ reduce_tracer_arrays(E);
+
+ E->trace.find_tracers_time += CPU_time0() - begin_time;
+
+ return;
+}
+
+
+/***********************************************************************/
+/* This function computes the number of tracers in each element. */
+/* Each tracer can be of different "flavors", which is the 0th index */
+/* of extraq. How to interprete "flavor" is left for the application. */
+
+void count_tracers_of_flavors(struct All_variables *E)
+{
+
+ int j, flavor, e, kk;
+ int numtracers;
+
+ for (j=1; j<=E->sphere.caps_per_proc; j++) {
+
+ /* first zero arrays */
+ for (flavor=0; flavor<E->trace.nflavors; flavor++)
+ for (e=1; e<=E->lmesh.nel; e++)
+ E->trace.ntracer_flavor[j][flavor][e] = 0;
+
+ numtracers=E->trace.ntracers[j];
+
+ /* Fill arrays */
+ for (kk=1; kk<=numtracers; kk++) {
+ e = E->trace.ielement[j][kk];
+ flavor = E->trace.extraq[j][0][kk];
+ E->trace.ntracer_flavor[j][flavor][e]++;
+ }
+ }
+
+ /* debug */
+ /**
+ for (j=1; j<=E->sphere.caps_per_proc; j++) {
+ for (e=1; e<=E->lmesh.nel; e++) {
+ fprintf(E->trace.fpt, "element=%d ntracer_flaver =", e);
+ for (flavor=0; flavor<E->trace.nflavors; flavor++) {
+ fprintf(E->trace.fpt, " %d",
+ E->trace.ntracer_flavor[j][flavor][e]);
+ }
+ fprintf(E->trace.fpt, "\n");
+ }
+ }
+ fflush(E->trace.fpt);
+ /**/
+
+ return;
+}
+
+
+
+void initialize_tracers(struct All_variables *E)
+{
+
+ if (E->trace.ic_method==0)
+ make_tracer_array(E);
+ else if (E->trace.ic_method==1)
+ read_tracer_file(E);
+ else if (E->trace.ic_method==2)
+ read_old_tracer_file(E);
+ else {
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ /* total number of tracers */
+
+ E->trace.ilast_tracer_count = isum_tracers(E);
+ fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
+ if(E->parallel.me==0)
+ fprintf(stderr, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
+
+
+ /* find elements */
+
+ find_tracers(E);
+
+
+ /* count # of tracers of each flavor */
+
+ if (E->trace.nflavors > 0)
+ count_tracers_of_flavors(E);
+
+ return;
+}
+
+
+/************** MAKE TRACER ARRAY ********************************/
+/* Here, each processor will generate tracers somewhere */
+/* in the sphere - check if its in this cap - then check radial */
+
+static void make_tracer_array(struct All_variables *E)
+{
+
+ int tracers_cap;
+ int j;
+ double processor_fraction;
+
+ void generate_random_tracers();
+ void init_tracer_flavors();
+
+ if (E->parallel.me==0) fprintf(stderr,"Making Tracer Array\n");
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ processor_fraction=E->lmesh.volume/E->mesh.volume;
+ 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);
+
+ generate_random_tracers(E, tracers_cap, j);
+
+
+
+ }/* end j */
+
+
+ /* Initialize tracer flavors */
+ if (E->trace.nflavors) init_tracer_flavors(E);
+
+ return;
+}
+
+
+
+static void generate_random_tracers(struct All_variables *E,
+ int tracers_cap, int j)
+{
+ void cart_to_sphere();
+ int kk;
+ int ival;
+ int number_of_tries=0;
+ int max_tries;
+
+ double x,y,z;
+ double theta,phi,rad;
+ double xmin,xmax,ymin,ymax,zmin,zmax;
+ double random1,random2,random3;
+
+
+ allocate_tracer_arrays(E,j,tracers_cap);
+
+ /* Finding the min/max of the cartesian coordinates. */
+ /* One must loop over E->X to find the min/max, since the 8 corner */
+ /* nodes may not be the min/max. */
+ xmin = ymin = zmin = E->sphere.ro;
+ xmax = ymax = zmax = -E->sphere.ro;
+ for (kk=1; kk<=E->lmesh.nno; kk++) {
+ x = E->x[j][1][kk];
+ y = E->x[j][2][kk];
+ z = E->x[j][3][kk];
+
+ xmin = ((xmin < x) ? xmin : x);
+ xmax = ((xmax > x) ? xmax : x);
+ ymin = ((ymin < y) ? ymin : y);
+ ymax = ((ymax > y) ? ymax : y);
+ zmin = ((zmin < z) ? zmin : z);
+ zmax = ((zmax > z) ? zmax : z);
+ }
+
+ /* Tracers are placed randomly in cap */
+ /* (intentionally using rand() instead of srand() )*/
+
+ while (E->trace.ntracers[j]<tracers_cap) {
+
+ number_of_tries++;
+ max_tries=100*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 1
+ random1=drand48();
+ random2=drand48();
+ random3=drand48();
+#else
+ random1=(1.0*rand())/(1.0*RAND_MAX);
+ random2=(1.0*rand())/(1.0*RAND_MAX);
+ random3=(1.0*rand())/(1.0*RAND_MAX);
+#endif
+
+ x=xmin+random1*(xmax-xmin);
+ y=ymin+random2*(ymax-ymin);
+ z=zmin+random3*(zmax-zmin);
+
+ /* first check if within shell */
+
+ cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
+
+ if (rad>=E->sx[j][3][E->lmesh.noz]) continue;
+ if (rad<E->sx[j][3][1]) continue;
+
+
+ /* check if in current cap */
+ if (E->parallel.nprocxy==1)
+ ival=regional_icheck_cap(E,0,theta,phi,rad,rad);
+ else
+ ival=full_icheck_cap(E,0,x,y,z,rad);
+
+ if (ival!=1) continue;
+
+ /* Made it, so record tracer information */
+
+ (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+
+ E->trace.ntracers[j]++;
+ kk=E->trace.ntracers[j];
+
+ 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;
+
+ } /* end while */
+
+ return;
+}
+
+
+/******** READ TRACER ARRAY *********************************************/
+/* */
+/* This function reads tracers from input file. */
+/* All processors read the same input file, then sort out which ones */
+/* belong. */
+
+static void read_tracer_file(struct All_variables *E)
+{
+
+ char input_s[1000];
+
+ int number_of_tracers, ncolumns;
+ int kk;
+ int icheck;
+ int iestimate;
+ int icushion;
+ int i, j;
+
+
+ int icheck_processor_shell();
+ void sphere_to_cart();
+ void cart_to_sphere();
+ void expand_tracer_arrays();
+
+ double x,y,z;
+ double theta,phi,rad;
+ double buffer[100];
+
+ FILE *fptracer;
+
+ fptracer=fopen(E->trace.tracer_file,"r");
+
+ fgets(input_s,200,fptracer);
+ if(sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns) != 2) {
+ fprintf(stderr,"Error while reading file '%s'\n", E->trace.tracer_file);
+ exit(8);
+ }
+ fprintf(E->trace.fpt,"%d Tracers, %d columns in file \n",
+ number_of_tracers, ncolumns);
+
+ /* some error control */
+ if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ fprintf(E->trace.fpt,"ERROR(read tracer file)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ /* initially size tracer arrays to number of tracers divided by processors */
+
+ icushion=100;
+
+ iestimate=number_of_tracers/E->parallel.nproc + icushion;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ allocate_tracer_arrays(E,j,iestimate);
+
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ int len, ncol;
+ ncol = 3 + E->trace.number_of_extra_quantities;
+
+ len = read_double_vector(fptracer, ncol, buffer);
+ if (len != ncol) {
+ fprintf(E->trace.fpt,"ERROR(read tracer file) - wrong input file format: %s\n", E->trace.tracer_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ theta = buffer[0];
+ phi = buffer[1];
+ rad = buffer[2];
+
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+
+
+ /* make sure theta, phi is in range, and radius is within bounds */
+
+ (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+
+ /* check whether tracer is within processor domain */
+
+ icheck=1;
+ if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,rad);
+ if (icheck!=1) continue;
+
+ if (E->parallel.nprocxy==1)
+ icheck=regional_icheck_cap(E,0,theta,phi,rad,rad);
+ else
+ icheck=full_icheck_cap(E,0,x,y,z,rad);
+
+ if (icheck==0) continue;
+
+ /* if still here, tracer is in processor domain */
+
+
+ E->trace.ntracers[j]++;
+
+ if (E->trace.ntracers[j]>=(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
+
+ 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;
+
+ for (i=0; i<E->trace.number_of_extra_quantities; i++)
+ E->trace.extraq[j][i][E->trace.ntracers[j]]=buffer[i+3];
+
+ } /* end kk, number of tracers */
+
+ fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
+ E->trace.ntracers[j]);
+
+ /** debug **
+ for (kk=1; kk<=E->trace.ntracers[j]; kk++) {
+ fprintf(E->trace.fpt, "tracer#=%d sph_coord=(%g,%g,%g)", kk,
+ E->trace.basicq[j][0][kk],
+ E->trace.basicq[j][1][kk],
+ E->trace.basicq[j][2][kk]);
+ fprintf(E->trace.fpt, " extraq=");
+ for (i=0; i<E->trace.number_of_extra_quantities; i++)
+ fprintf(E->trace.fpt, " %g", E->trace.extraq[j][i][kk]);
+ fprintf(E->trace.fpt, "\n");
+ }
+ fflush(E->trace.fpt);
+ /**/
+
+ } /* 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);
+ }
+
+ return;
+}
+
+
+/************** READ OLD TRACER FILE *************************************/
+/* */
+/* This function read tracers written from previous calculation */
+/* and the tracers are read as seperate files for each processor domain. */
+
+static void read_old_tracer_file(struct All_variables *E)
+{
+
+ char output_file[200];
+ char input_s[1000];
+
+ int i,j,kk,rezip;
+ int idum1,ncolumns;
+ int numtracers;
+
+ double rdum1;
+ double theta,phi,rad;
+ double x,y,z;
+ double buffer[100];
+
+ void sphere_to_cart();
+
+ FILE *fp1;
+
+ if (E->trace.number_of_extra_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-increase size of extra[]\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+
+
+ /* deal with different output formats */
+#ifdef USE_GZDIR
+ if(strcmp(E->output.format, "ascii-gz") == 0){
+ sprintf(output_file,"%s/%d/tracer.%d.%d",
+ E->control.data_dir_old,E->monitor.solution_cycles_init,E->parallel.me,E->monitor.solution_cycles_init);
+ rezip = open_file_zipped(output_file,&fp1,E);
+ }else{
+ 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(read_old_tracer_file)-gziped file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+#else
+ 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(read_old_tracer_file)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+#endif
+
+ fprintf(stderr,"Read old tracers from %s\n",output_file);
+
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fgets(input_s,200,fp1);
+ if(sscanf(input_s,"%d %d %d %lf",
+ &idum1, &numtracers, &ncolumns, &rdum1) != 4) {
+ fprintf(stderr,"Error while reading file '%s'\n", output_file);
+ exit(8);
+ }
+
+
+ /* some error control */
+ if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ /* allocate memory for tracer arrays */
+
+ allocate_tracer_arrays(E,j,numtracers);
+ E->trace.ntracers[j]=numtracers;
+
+ for (kk=1;kk<=numtracers;kk++) {
+ int len, ncol;
+ ncol = 3 + E->trace.number_of_extra_quantities;
+
+ len = read_double_vector(fp1, ncol, buffer);
+ if (len != ncol) {
+ fprintf(E->trace.fpt,"ERROR(read_old_tracer_file) - wrong input file format: %s\n", output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ theta = buffer[0];
+ phi = buffer[1];
+ rad = buffer[2];
+
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+
+ /* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
+
+ (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+
+ 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;
+
+ for (i=0; i<E->trace.number_of_extra_quantities; i++)
+ E->trace.extraq[j][i][kk]=buffer[i+3];
+
+ }
+
+ /** debug **
+ for (kk=1; kk<=E->trace.ntracers[j]; kk++) {
+ fprintf(E->trace.fpt, "tracer#=%d sph_coord=(%g,%g,%g)", kk,
+ E->trace.basicq[j][0][kk],
+ E->trace.basicq[j][1][kk],
+ E->trace.basicq[j][2][kk]);
+ fprintf(E->trace.fpt, " extraq=");
+ for (i=0; i<E->trace.number_of_extra_quantities; i++)
+ fprintf(E->trace.fpt, " %g", E->trace.extraq[j][i][kk]);
+ fprintf(E->trace.fpt, "\n");
+ }
+ fflush(E->trace.fpt);
+ /**/
+
+ fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
+ fflush(E->trace.fpt);
+
+ }
+ fclose(fp1);
+#ifdef USE_GZDIR
+ if(strcmp(E->output.format, "ascii-gz") == 0)
+ if(rezip) /* rezip */
+ gzip_file(output_file);
+#endif
+
+ return;
+}
+
+
+
+
+
+/*********** CHECK SUM **************************************************/
+/* */
+/* This functions checks to make sure number of tracers is preserved */
+
+static void check_sum(struct All_variables *E)
+{
+
+ int number, iold_number;
+
+ number = isum_tracers(E);
+
+ 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);
+ if (E->trace.itracer_warnings)
+ parallel_process_termination();
+ }
+
+ E->trace.ilast_tracer_count = number;
+
+ return;
+}
+
+
+/************* ISUM TRACERS **********************************************/
+/* */
+/* This function uses MPI to sum all tracers and returns number of them. */
+
+static int isum_tracers(struct All_variables *E)
+{
+ int imycount;
+ int iallcount;
+ int j;
+
+ iallcount = 0;
+
+ imycount = 0;
+ for (j=1; j<=E->sphere.caps_per_proc; j++)
+ imycount = imycount + E->trace.ntracers[j];
+
+ MPI_Allreduce(&imycount,&iallcount,1,MPI_INT,MPI_SUM,E->parallel.world);
+
+ return iallcount;
+}
+
+
+
+/********** CART TO SPHERE ***********************/
+void cart_to_sphere(struct All_variables *E,
+ double x, double y, double z,
+ double *theta, double *phi, double *rad)
+{
+
+ double temp;
+
+ temp=x*x+y*y;
+
+ *rad=sqrt(temp+z*z);
+ *theta=atan2(sqrt(temp),z);
+ *phi=myatan(y,x);
+
+
+ return;
+}
+
+/********** SPHERE TO CART ***********************/
+void sphere_to_cart(struct All_variables *E,
+ double theta, double phi, double rad,
+ double *x, double *y, double *z)
+{
+
+ double sint,cost,sinf,cosf;
+ double temp;
+
+ sint=sin(theta);
+ cost=cos(theta);
+ sinf=sin(phi);
+ cosf=cos(phi);
+
+ temp=rad*sint;
+
+ *x=temp*cosf;
+ *y=temp*sinf;
+ *z=rad*cost;
+
+ return;
+}
+
+
+
+static void init_tracer_flavors(struct All_variables *E)
+{
+ int j, kk, number_of_tracers;
+ int i;
+ double flavor;
+ double rad;
+
+ switch(E->trace.ic_method_for_flavors){
+ case 0:
+ /* ic_method_for_flavors == 0 (layered structure) */
+ /* any tracer above z_interface[i] is of flavor i */
+ /* any tracer below z_interface is of flavor (nflavors-1) */
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ number_of_tracers = E->trace.ntracers[j];
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ rad = E->trace.basicq[j][2][kk];
+
+ flavor = E->trace.nflavors - 1;
+ for (i=0; i<E->trace.nflavors-1; i++) {
+ if (rad > E->trace.z_interface[i]) {
+ flavor = i;
+ break;
+ }
+ }
+ E->trace.extraq[j][0][kk] = flavor;
+ }
+ }
+ break;
+
+ case 1: /* from grd in top n layers */
+ case 99: /* (will override restart) */
+#ifndef USE_GGRD
+ fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
+ E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+#else
+ ggrd_init_tracer_flavors(E);
+#endif
+ break;
+
+
+ default:
+
+ fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
+ }
+
+ return;
+}
+
+
+/******************* get_neighboring_caps ************************************/
+/* */
+/* Communicate with neighboring processors to get their cap boundaries, */
+/* which is later used by (E->trace.icheck_cap)() */
+/* */
+
+void get_neighboring_caps(struct All_variables *E)
+{
+ void sphere_to_cart();
+
+ const int ncorners = 4; /* # of top corner nodes */
+ int i, j, n, d, kk, lev, idb;
+ int num_ngb, neighbor_proc, tag;
+ MPI_Status status[200];
+ MPI_Request request[200];
+
+ int node[ncorners];
+ double xx[ncorners*2], rr[12][ncorners*2];
+ int nox,noy,noz;
+ double x,y,z;
+ double theta,phi,rad;
+
+ nox=E->lmesh.nox;
+ noy=E->lmesh.noy;
+ noz=E->lmesh.noz;
+
+ node[0]=nox*noz*(noy-1)+noz;
+ node[1]=noz;
+ node[2]=noz*nox;
+ node[3]=noz*nox*noy;
+
+ lev = E->mesh.levmax;
+ tag = 45;
+
+ for (j=1; j<=E->sphere.caps_per_proc; j++) {
+
+ /* loop over top corners to get their coordinates */
+ n = 0;
+ for (i=0; i<ncorners; i++) {
+ for (d=0; d<2; d++) {
+ xx[n] = E->sx[j][d+1][node[i]];
+ n++;
+ }
+ }
+
+ idb = 0;
+ num_ngb = E->parallel.TNUM_PASS[lev][j];
+ for (kk=1; kk<=num_ngb; kk++) {
+ neighbor_proc = E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ MPI_Isend(xx, n, MPI_DOUBLE, neighbor_proc,
+ tag, E->parallel.world, &request[idb]);
+ idb++;
+
+ MPI_Irecv(rr[kk], n, MPI_DOUBLE, neighbor_proc,
+ tag, E->parallel.world, &request[idb]);
+ idb++;
+ }
+
+ /* Storing the current cap information */
+ for (i=0; i<n; i++)
+ rr[0][i] = xx[i];
+
+ /* Wait for non-blocking calls to complete */
+
+ MPI_Waitall(idb, request, status);
+
+ /* Storing the received cap information
+ * XXX: this part assumes:
+ * 1) E->sphere.caps_per_proc==1
+ * 2) E->mesh.nsd==3
+ */
+ for (kk=0; kk<=num_ngb; kk++) {
+ n = 0;
+ for (i=1; i<=ncorners; i++) {
+ theta = rr[kk][n++];
+ phi = rr[kk][n++];
+ rad = E->sphere.ro;
+
+ sphere_to_cart(E, theta, phi, rad, &x, &y, &z);
+
+ E->trace.xcap[kk][i] = x;
+ E->trace.ycap[kk][i] = y;
+ E->trace.zcap[kk][i] = z;
+ E->trace.theta_cap[kk][i] = theta;
+ E->trace.phi_cap[kk][i] = phi;
+ E->trace.rad_cap[kk][i] = rad;
+ E->trace.cos_theta[kk][i] = cos(theta);
+ E->trace.sin_theta[kk][i] = sin(theta);
+ E->trace.cos_phi[kk][i] = cos(phi);
+ E->trace.sin_phi[kk][i] = sin(phi);
+ }
+ } /* end kk, number of neighbors */
+
+ /* debugging output *
+ for (kk=0; kk<=num_ngb; kk++) {
+ if (kk==0)
+ neighbor_proc = E->parallel.me;
+ else
+ neighbor_proc = E->parallel.PROCESSOR[lev][1].pass[kk];
+
+ for (i=1; i<=ncorners; i++) {
+ fprintf(E->trace.fpt, "pass=%d rank=%d corner=%d "
+ "sx=(%g, %g, %g)\n",
+ kk, neighbor_proc, i,
+ E->trace.theta_cap[kk][i],
+ E->trace.phi_cap[kk][i],
+ E->trace.rad_cap[kk][i]);
+ }
+ }
+ fflush(E->trace.fpt);
+ /**/
+ }
+
+ return;
+}
+
+
+/**************** INITIALIZE TRACER ARRAYS ************************************/
+/* */
+/* This function allocates memories to tracer arrays. */
+
+void allocate_tracer_arrays(struct All_variables *E,
+ int j, int number_of_tracers)
+{
+
+ int kk;
+
+ /* max_ntracers is physical size of tracer array */
+ /* (initially make it 25% larger than required */
+
+ E->trace.max_ntracers[j]=number_of_tracers+number_of_tracers/4;
+ E->trace.ntracers[j]=0;
+
+ /* make tracer arrays */
+
+ 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=1;kk<E->trace.max_ntracers[j];kk++)
+ E->trace.ielement[j][kk]=-99;
+
+
+ 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;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);
+ }
+ }
+
+ if (E->trace.nflavors > 0) {
+ E->trace.ntracer_flavor[j]=(int **)malloc(E->trace.nflavors*sizeof(int*));
+ for (kk=0;kk<E->trace.nflavors;kk++) {
+ if ((E->trace.ntracer_flavor[j][kk]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==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 (max_ntracers): %d\n",
+ E->trace.max_ntracers[j]);
+ fflush(E->trace.fpt);
+
+ return;
+}
+
+
+
+/****** EXPAND TRACER ARRAYS *****************************************/
+
+void expand_tracer_arrays(struct All_variables *E, int j)
+{
+
+ int inewsize;
+ int kk;
+ int icushion;
+
+ /* expand basicq and ielement by 20% */
+
+ icushion=100;
+
+ inewsize=E->trace.max_ntracers[j]+E->trace.max_ntracers[j]/5+icushion;
+
+ 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_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.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 ielement, basicq, and extraq to %d from %d\n",
+ inewsize,E->trace.max_ntracers[j]);
+
+ E->trace.max_ntracers[j]=inewsize;
+
+ return;
+}
+
+
+
+
+/****** REDUCE TRACER ARRAYS *****************************************/
+
+static void reduce_tracer_arrays(struct All_variables *E)
+{
+
+ int inewsize;
+ int kk;
+ int iempty_space;
+ int j;
+
+ int icushion=100;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+
+ /* if physical size is double tracer size, reduce it */
+
+ iempty_space=(E->trace.max_ntracers[j]-E->trace.ntracers[j]);
+
+ if (iempty_space>(E->trace.ntracers[j]+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 ((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_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.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 ielement, basicq, and extraq to %d from %d\n",
+ E->trace.max_ntracers[j],inewsize);
+
+ E->trace.max_ntracers[j]=inewsize;
+
+ } /* end if */
+
+ } /* end j */
+
+ return;
+}
+
+
+/********** PUT AWAY LATER ************************************/
+/* */
+/* rlater has a similar structure to basicq */
+/* ilatersize is the physical memory and */
+/* ilater is the number of tracers */
+
+static void put_away_later(struct All_variables *E, int j, int it)
+{
+ int kk;
+ void expand_later_array();
+
+
+ /* The first tracer in initiates memory allocation. */
+ /* Memory is freed after parallel communications */
+
+ if (E->trace.ilatersize[j]==0) {
+
+ 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 */
+
+
+ /* Put tracer in later array */
+
+ E->trace.ilater[j]++;
+
+ if (E->trace.ilater[j] >= (E->trace.ilatersize[j]-5)) expand_later_array(E,j);
+
+ /* stack basic and extra quantities together (basic first) */
+
+ 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_basic_quantities+kk][E->trace.ilater[j]]=E->trace.extraq[j][kk][it];
+
+
+ return;
+}
+
+
+/****** EXPAND LATER ARRAY *****************************************/
+
+void expand_later_array(struct All_variables *E, int j)
+{
+
+ int inewsize;
+ int kk;
+ int icushion;
+
+ /* expand rlater by 20% */
+
+ icushion=100;
+
+ 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);
+ }
+ }
+
+
+ fprintf(E->trace.fpt,"Expanding physical memory of rlater to %d from %d\n",
+ inewsize,E->trace.ilatersize[j]);
+
+ E->trace.ilatersize[j]=inewsize;
+
+ return;
+}
+
+
+/***** EJECT TRACER ************************************************/
+
+static void eject_tracer(struct All_variables *E, int j, int it)
+{
+
+ int ilast_tracer;
+ int kk;
+
+
+ ilast_tracer=E->trace.ntracers[j];
+
+ /* put last tracer in ejected tracer position */
+
+ E->trace.ielement[j][it]=E->trace.ielement[j][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.extraq[j][kk][it]=E->trace.extraq[j][kk][ilast_tracer];
+
+
+
+ E->trace.ntracers[j]--;
+
+ return;
+}
+
+
+
+/********** ICHECK PROCESSOR SHELL *************/
+/* returns -99 if rad is below current shell */
+/* returns 0 if rad is above current shell */
+/* returns 1 if rad is within current shell */
+/* */
+/* Shell, here, refers to processor shell */
+/* */
+/* shell is defined as bottom boundary up to */
+/* and not including the top boundary unless */
+/* the shell in question is the top shell */
+
+int icheck_processor_shell(struct All_variables *E,
+ int j, double rad)
+{
+
+ const int noz = E->lmesh.noz;
+ const int nprocz = E->parallel.nprocz;
+ double top_r, bottom_r;
+
+ if (nprocz==1) return 1;
+
+ top_r = E->sx[j][3][noz];
+ bottom_r = E->sx[j][3][1];
+
+ /* First check bottom */
+
+ if (rad<bottom_r) return -99;
+
+
+ /* Check top */
+
+ if (rad<top_r) return 1;
+
+ /* top processor */
+
+ if ( (rad<=top_r) && (E->parallel.me_loc[3]==nprocz-1) ) return 1;
+
+ /* If here, means point is above processor */
+ return 0;
+}
+
+
+/********* ICHECK THAT PROCESSOR SHELL ********/
+/* */
+/* Checks whether a given radius is within */
+/* a given processors radial domain. */
+/* Returns 0 if not, 1 if so. */
+/* The domain is defined as including the bottom */
+/* radius, but excluding the top radius unless */
+/* we the processor domain is the one that */
+/* is at the surface (then both boundaries are */
+/* included). */
+
+int icheck_that_processor_shell(struct All_variables *E,
+ int j, int nprocessor, double rad)
+{
+ int icheck_processor_shell();
+ int me = E->parallel.me;
+
+ /* nprocessor is right on top of me */
+ if (nprocessor == me+1) {
+ if (icheck_processor_shell(E, j, rad) == 0) return 1;
+ else return 0;
+ }
+
+ /* nprocessor is right on bottom of me */
+ if (nprocessor == me-1) {
+ if (icheck_processor_shell(E, j, rad) == -99) return 1;
+ else return 0;
+ }
+
+ /* Shouldn't be here */
+ fprintf(E->trace.fpt, "Should not be here\n");
+ fprintf(E->trace.fpt, "Error(check_shell) nprocessor: %d, radius: %f\n",
+ nprocessor, rad);
+ fflush(E->trace.fpt);
+ exit(10);
+
+ return 0;
+}
+
+
Modified: mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c 2011-07-19 18:34:06 UTC (rev 18782)
@@ -937,7 +937,8 @@
return;
}
-void visc_from_P(struct All_variables *E, float **EEta) /* "plasticity" implementation
+void visc_from_P(struct All_variables *E, float **EEta)
+ /* "plasticity" implementation
psrw = FALSE
@@ -1379,7 +1380,7 @@
E->sx[m][3][E->ien[m][ee].node[8]]);
/* if ee has tracers in it and is within the channel */
- if((E->trace.num_tracer_flavors[m][flavor][ee] > 0) &&
+ if((E->trace.ntracer_flavor[m][flavor][ee] > 0) &&
(rad_mean <= rr + E->viscosity.lv_channel_thickness)) {
F[e] = E->viscosity.lv_reduction;
break;
@@ -1440,7 +1441,7 @@
ee = (k-1)*E->lmesh.elz + ii;
/* if ee has tracers in it */
- if(E->trace.num_tracer_flavors[m][flavor][ee] > 0) {
+ if(E->trace.ntracer_flavor[m][flavor][ee] > 0) {
F[e] = E->viscosity.lv_reduction;
break;
}
Modified: mc/3D/CitcomS/branches/eheien/lib/prototypes.h
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/prototypes.h 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/prototypes.h 2011-07-19 18:34:06 UTC (rev 18782)
@@ -156,10 +156,10 @@
void full_lost_souls(struct All_variables *);
void full_get_shape_functions(struct All_variables *, double [9], int, double, double, double);
double full_interpolate_data(struct All_variables *, double [9], double [9]);
-CartesianCoord full_get_velocity(struct All_variables *, int, int, double, double, double);
-int full_icheck_cap(struct All_variables *, int, Tracer);
+void full_get_velocity(struct All_variables *, int, int, double, double, double, double *);
+int full_icheck_cap(struct All_variables *, int, double, double, double, double);
int full_iget_element(struct All_variables *, int, int, double, double, double, double, double, double);
-void full_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
+void full_keep_within_bounds(struct All_variables *, double *, double *, double *, double *, double *, double *);
void analytical_test(struct All_variables *);
void analytical_runge_kutte(struct All_variables *, int, double, double *, double *, double *, double *, double *);
void analytical_test_function(struct All_variables *, double, double, double, double *, double *);
@@ -429,11 +429,11 @@
int regional_iget_element(struct All_variables *, int, int, double, double, double, double, double, double);
int isearch_all(double *, int, double);
int isearch_neighbors(double *, int, double, int);
-int regional_icheck_cap(struct All_variables *, int, Tracer);
+int regional_icheck_cap(struct All_variables *, int, double, double, double, double);
void regional_get_shape_functions(struct All_variables *, double [9], int, double, double, double);
double regional_interpolate_data(struct All_variables *, double [9], double [9]);
-CartesianCoord regional_get_velocity(struct All_variables *, int, int, double, double, double);
-void regional_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
+void regional_get_velocity(struct All_variables *, int, int, double, double, double, double *);
+void regional_keep_within_bounds(struct All_variables *, double *, double *, double *, double *, double *, double *);
void regional_lost_souls(struct All_variables *);
/* Regional_version_dependent.c */
void regional_node_locations(struct All_variables *);
Modified: mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h 2011-07-18 06:12:07 UTC (rev 18781)
+++ mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h 2011-07-19 18:34:06 UTC (rev 18782)
@@ -26,106 +26,12 @@
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
-#include <vector>
-#include <list>
+
/* forward declaration */
struct All_variables;
-#ifndef _TRACER_DEFS_H_
-#define _TRACER_DEFS_H_
-class CartesianCoord;
-
-// Position or vector in spherical coordinates
-class SphericalCoord {
-public:
- double _theta, _phi, _rad;
- SphericalCoord(void) : _theta(0), _phi(0), _rad(0) {};
- SphericalCoord(double theta, double phi, double rad) : _theta(theta), _phi(phi), _rad(rad) {};
-
- size_t size(void) const { return 3; };
- void writeToMem(double *mem) const;
- void readFromMem(const double *mem);
- CartesianCoord toCartesian(void) const;
-
- void constrainThetaPhi(void);
- double constrainAngle(const double angle) const;
-};
-
-// Position or vector in Cartesian coordinates
-class CartesianCoord {
-public:
- double _x, _y, _z;
- CartesianCoord(void) : _x(0), _y(0), _z(0) {};
- CartesianCoord(double x, double y, double z) : _x(x), _y(y), _z(z) {};
-
- size_t size(void) const { return 3; };
- void writeToMem(double *mem) const;
- void readFromMem(const double *mem);
- SphericalCoord toSpherical(void) const;
-
- const CartesianCoord operator+(const CartesianCoord &other) const;
- const CartesianCoord operator*(const double &val) const;
-};
-
-class Tracer {
-private:
- // Tracer position in spherical coordinates
- SphericalCoord _sc;
- // Tracer position in Cartesian coordinates
- CartesianCoord _cc;
- // Previous Cartesian position
- CartesianCoord _cc0;
- // Previous Cartesian velocity
- CartesianCoord _Vc;
-
- // Tracer flavor (meaning should be application dependent)
- double _flavor;
-
- int _ielement;
-
-public:
- Tracer(void) : _sc(), _cc(), _cc0(), _Vc() {};
- Tracer(SphericalCoord new_sc, CartesianCoord new_cc) : _sc(new_sc), _cc(new_cc), _cc0(), _Vc() {};
-
- CartesianCoord getCartesianPos(void) const { return _cc; };
- SphericalCoord getSphericalPos(void) const { return _sc; };
- CartesianCoord getOrigCartesianPos(void) const { return _cc0; };
- CartesianCoord getCartesianVel(void) const { return _Vc; };
-
- void setCoords(SphericalCoord new_sc, CartesianCoord new_cc) {
- _sc = new_sc;
- _cc = new_cc;
- }
- void setOrigVals(CartesianCoord new_cc0, CartesianCoord new_vc) {
- _cc0 = new_cc0;
- _Vc = new_vc;
- }
-
- double theta(void) { return _sc._theta; };
- double phi(void) { return _sc._phi; };
- double rad(void) { return _sc._rad; };
-
- double x(void) { return _cc._x; };
- double y(void) { return _cc._y; };
- double z(void) { return _cc._z; };
-
- int ielement(void) { return _ielement; };
- void set_ielement(int ielement) { _ielement = ielement; };
-
- double flavor(void) { return _flavor; };
- void set_flavor(double flavor) { _flavor = flavor; };
-
- size_t size(void);
- void writeToMem(double *mem) const;
- void readFromMem(const double *mem);
-};
-
-typedef std::list<Tracer> TracerList;
-
-#endif
-
struct TRACE{
FILE *fpt;
@@ -141,14 +47,24 @@
double box_cushion;
/* tracer arrays */
- TracerList tracers[13];
+ int number_of_basic_quantities;
+ int number_of_extra_quantities;
+ int number_of_tracer_quantities;
- TracerList escaped_tracers[13];
+ double *basicq[13][100];
+ double *extraq[13][100];
+ int ntracers[13];
+ int max_ntracers[13];
+ int *ielement[13];
+
+ int ilatersize[13];
+ int ilater[13];
+ double *rlater[13][100];
+
/* tracer flavors */
int nflavors;
- std::map<int, std::map<int, int> > num_tracer_flavors[13];
- //int **ntracer_flavor[13];
+ int **ntracer_flavor[13];
int ic_method_for_flavors;
double *z_interface;
@@ -230,8 +146,10 @@
int (* iget_element)(struct All_variables*, int, int,
double, double, double, double, double, double);
- CartesianCoord (* get_velocity)(struct All_variables*, int, int,
- double, double, double);
+ void (* get_velocity)(struct All_variables*, int, int,
+ double, double, double, double*);
- void (* keep_within_bounds)(struct All_variables*, CartesianCoord &, SphericalCoord &);
+ void (* keep_within_bounds)(struct All_variables*,
+ double*, double*, double*,
+ double*, double*, double*);
};
More information about the CIG-COMMITS
mailing list