[cig-commits] r6273 - in mc/3D/CitcomS/trunk: bin lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Mar 16 17:21:02 PDT 2007
Author: tan2
Date: 2007-03-16 17:21:01 -0700 (Fri, 16 Mar 2007)
New Revision: 6273
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/lib/Full_obsolete.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
mc/3D/CitcomS/trunk/module/misc.c
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Merged the regional and global tracer codes.
Basically, Vlad's code is mostly gone and is replaced by Allen's and mine code.
A few functions, eg. analytical_test(), low_viscosity_*_factor(), are disabled,
need to work on them later.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -59,6 +59,7 @@
void open_time();
void output_finalize();
void PG_timestep_init();
+ void tracer_advection();
float dot();
float cpu_time_on_vp_it;
@@ -116,7 +117,7 @@
if (E->control.stokes) {
if(E->control.tracer==1)
- (E->problem_tracer_advection)(E);
+ tracer_advection(E);
parallel_process_termination();
}
@@ -155,7 +156,7 @@
}
if(E->control.tracer==1)
- (E->problem_tracer_advection)(E);
+ tracer_advection(E);
general_stokes_solver(E);
Modified: mc/3D/CitcomS/trunk/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -887,6 +887,182 @@
}
+/****** ICHECK REGULAR NEIGHBORS *****************************/
+/* */
+/* This function searches the regular element neighborhood. */
+
+/* This function is no longer used! */
+
+int icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad)
+ struct All_variables *E;
+ int j,ntheta,nphi;
+ double x,y,z;
+ double theta,phi,rad;
+{
+
+ int new_ntheta,new_nphi;
+ int kk,pp;
+ int iregel;
+ int ival;
+ int imap[5];
+ int ichoice;
+ int irange;
+
+ int iquick_element_column_search();
+
+ fprintf(E->trace.fpt,"ERROR(icheck_regular_neighbors)-this subroutine is no longer used !\n");
+ fflush(E->trace.fpt);
+ exit(10);
+
+ irange=2;
+
+ for (kk=-irange;kk<=irange;kk++)
+ {
+ for (pp=-irange;pp<=irange;pp++)
+ {
+ new_ntheta=ntheta+kk;
+ new_nphi=nphi+pp;
+ if ( (new_ntheta>0)&&(new_ntheta<=E->trace.numtheta[j])&&(new_nphi>0)&&(new_nphi<=E->trace.numphi[j]) )
+ {
+ iregel=new_ntheta+(new_nphi-1)*E->trace.numtheta[j];
+ if ((iregel>0) && (iregel<=E->trace.numregel[j]))
+ {
+ ival=iquick_element_column_search(E,j,iregel,new_ntheta,new_nphi,x,y,z,theta,phi,rad,imap,&ichoice);
+ if (ival>0) return ival;
+ }
+ }
+ }
+ }
+
+
+ return -99;
+}
+
+
+/****** IQUICK ELEMENT SEARCH *****************************/
+/* */
+/* This function does a quick regular to real element */
+/* map check. Element number, if found, is returned. */
+/* Otherwise, -99 is returned. */
+/* Pointers to imap and ichoice are used because they may */
+/* prove to be convenient. */
+/* This routine is no longer used */
+
+int iquick_element_column_search(E,j,iregel,ntheta,nphi,x,y,z,theta,phi,rad,imap,ich)
+ struct All_variables *E;
+ int j,iregel;
+ int ntheta,nphi;
+ double x,y,z,theta,phi,rad;
+ int *imap;
+ int *ich;
+{
+
+ int iregnode[5];
+ int kk,pp;
+ int nel,ival;
+ int ichoice;
+ int icount;
+ int itemp1;
+ int itemp2;
+
+ int icheck_element_column();
+
+ fprintf(E->trace.fpt,"ERROR(iquick element)-this routine is no longer used!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+
+ /* REMOVE*/
+ /*
+ ichoice=*ich;
+
+ fprintf(E->trace.fpt,"AA: ichoice: %d\n",ichoice);
+ fflush(E->trace.fpt);
+ */
+
+ /* find regular nodes on regular element */
+
+ /*
+ iregnode[1]=iregel+(nphi-1);
+ iregnode[2]=iregel+nphi;
+ iregnode[3]=iregel+nphi+E->trace.numtheta[j]+1;
+ iregnode[4]=iregel+nphi+E->trace.numtheta[j];
+ */
+
+ itemp1=iregel+nphi;
+ itemp2=itemp1+E->trace.numtheta[j];
+
+ iregnode[1]=itemp1-1;
+ iregnode[2]=itemp1;
+ iregnode[3]=itemp2+1;
+ iregnode[4]=itemp2;
+
+ for (kk=1;kk<=4;kk++)
+ {
+ if ((iregnode[kk]<1) || (iregnode[kk]>E->trace.numregnodes[j]) )
+ {
+ fprintf(E->trace.fpt,"ERROR(iquick)-weird regnode %d\n",iregnode[kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ /* find number of choices */
+
+ ichoice=0;
+ icount=0;
+ for (kk=1;kk<=4;kk++)
+ {
+ if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
+
+ icount++;
+ for (pp=1;pp<=(kk-1);pp++)
+ {
+ if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
+ }
+ ichoice++;
+ imap[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
+
+
+ next_corner:
+ ;
+ } /* end kk */
+
+ *ich=ichoice;
+
+ /* statistical counter */
+
+ E->trace.istat_ichoice[j][ichoice]++;
+
+ if (ichoice==0) return -99;
+
+ /* Here, no check is performed if all 4 corners */
+ /* lie within a given element. */
+ /* It may be possible (not sure) but unlikely */
+ /* that the tracer is still not in that element */
+
+ /* Decided to comment this out. */
+ /* May not be valid for large regular grids. */
+ /*
+ */
+ /* AKMA */
+
+ if ((ichoice==1)&&(icount==4)) return imap[1];
+
+ /* check others */
+
+ for (kk=1;kk<=ichoice;kk++)
+ {
+ nel=imap[kk];
+ ival=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (ival>0) return nel;
+ }
+
+ /* if still here, no element was found */
+
+ return -99;
+}
+
+
/*************************************************************************/
/* from */
/*************************************************************************/
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -33,9 +33,7 @@
#include "parallel_related.h"
#include "composition_related.h"
-void count_tracers_of_flavors(struct All_variables *E);
-
-static void get_neighboring_caps(struct All_variables *E);
+void get_neighboring_caps(struct All_variables *E);
static void pdebug(struct All_variables *E, int i);
static void fix_radius(struct All_variables *E,
@@ -43,82 +41,18 @@
double *x, double *y, double *z);
static void fix_theta_phi(double *theta, double *phi);
static void fix_phi(double *phi);
-static void predict_tracers(struct All_variables *E);
-static void correct_tracers(struct All_variables *E);
-static int isum_tracers(struct All_variables *E);
static void make_tracer_array(struct All_variables *E);
-static void init_tracer_flavors(struct All_variables *E);
+void init_tracer_flavors(struct All_variables *E);
/******* FULL TRACER INPUT *********************/
-void full_tracer_input(E)
- struct All_variables *E;
-
+void full_tracer_input(struct All_variables *E)
{
int m = E->parallel.me;
- /* Initial condition, this option is ignored if E->control.restart is 1,
- * ie. restarted from a previous run */
- /* 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) */
- if(E->control.restart)
- E->trace.ic_method = 2;
- else {
- 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) {
- }
- else {
- fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
- fflush(stderr);
- 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);
-
-
- input_int("ic_method_for_flavors",&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
- if (E->trace.ic_method_for_flavors == 0)
- input_double("z_interface",&(E->trace.z_interface),"0.7",m);
-
-
- /* Advection Scheme */
-
- /* itracer_advection_scheme=1 (simple predictor corrector -uses only V(to)) */
- /* itracer_advection_scheme=2 (predictor-corrector - uses V(to) and V(to+dt)) */
-
- E->trace.itracer_advection_scheme=2;
- input_int("tracer_advection_scheme",&(E->trace.itracer_advection_scheme),
- "2,0,nomax",m);
-
- if (E->trace.itracer_advection_scheme==1)
- {
- }
- else if (E->trace.itracer_advection_scheme==2)
- {
- }
- else
- {
- fprintf(stderr,"Sorry, only advection scheme 1 and 2 available (%d)\n",E->trace.itracer_advection_scheme);
- fflush(stderr);
- parallel_process_termination();
- }
-
-
/* Interpolation Scheme */
/* itracer_interpolation_scheme=1 (gnometric projection) */
/* itracer_interpolation_scheme=2 (simple average) */
@@ -150,14 +84,12 @@
"0,0,nomax",m);
- composition_input(E);
return;
}
/***** FULL TRACER SETUP ************************/
-void full_tracer_setup(E)
- struct All_variables *E;
+void full_tracer_setup(struct All_variables *E)
{
char output_file[200];
@@ -169,11 +101,11 @@
void initialize_tracer_elements();
void define_uv_space();
void determine_shape_coefficients();
- void check_sum();
void read_tracer_file();
void analytical_test();
void tracer_post_processing();
void restart_tracers();
+ int isum_tracers();
/* Some error control */
@@ -270,12 +202,11 @@
read_tracer_file(E);
else if (E->trace.ic_method==2)
restart_tracers(E);
- else
- {
- fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
- fflush(E->trace.fpt);
- exit(10);
- }
+ else {
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
/* total number of tracers */
@@ -283,11 +214,6 @@
fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
- /* flush and wait for not real reason but it can't hurt */
- fflush(E->trace.fpt);
- parallel_process_sync(E);
-
-
/* find elements */
find_tracers(E);
@@ -312,290 +238,83 @@
}
-/******* TRACING *************************************************************/
-/* */
-/* 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 full_tracer_advection(E)
- struct All_variables *E;
+void full_put_lost_tracers(struct All_variables *E,
+ int isend[13][13], double *send[13][13])
{
+ int j, kk, pp;
+ int numtracers, ithatcap, icheck;
+ int isend_position, ipos;
+ int lev = E->mesh.levmax;
+ double theta, phi, rad;
+ double x, y, z;
- void check_sum();
- void tracer_post_processing();
-
-
- fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
- fflush(E->trace.fpt);
-
- /* advect tracers */
-
- predict_tracers(E);
- correct_tracers(E);
-
- check_sum(E);
-
- if (E->trace.nflavors > 0)
- count_tracers_of_flavors(E);
-
- if (E->composition.on) {
- fill_composition(E);
- }
-
- tracer_post_processing(E);
-
- return;
-}
-
-
-
-/********* TRACER POST PROCESSING ****************************************/
-
-void tracer_post_processing(E)
- struct All_variables *E;
-
-{
-
- char output_file[200];
-
- double convection_time,tracer_time;
- double trace_fraction,total_time;
-
-
- fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
-
- fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
-
- fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
-
- if (E->composition.ichemical_buoyancy==1) {
- fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
- E->trace.istat_iempty,(100.0*E->trace.istat_iempty)/E->lmesh.nel);
- }
-
-
- /* reset statistical counters */
-
- E->trace.istat_isend=0;
- E->trace.istat_iempty=0;
- E->trace.istat_elements_checked=0;
- E->trace.istat1=0;
-
- /* compositional and error fraction data files */
- //TODO: move
- if (E->parallel.me==0) {
- if (E->composition.ichemical_buoyancy==1) {
- fprintf(E->fp,"composition: %e %e\n",E->monitor.elapsed_time,E->composition.bulk_composition);
- fprintf(E->fp,"composition_error_fraction: %e %e\n",E->monitor.elapsed_time,E->composition.error_fraction);
-
- }
-
- }
-
- 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 */
-/* regardless of advection scheme]. */
-/* 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 get_velocity();
- void cart_to_sphere();
- void keep_in_sphere();
- void find_tracers();
-
-
- dt=E->advection.timestep;
-
-
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- numtracers=E->trace.ntracers[j];
+ /* transfer tracers from rlater to send */
+ 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];
- 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];
+ /* first check same cap if nprocz>1 */
- nelem=E->trace.ielement[j][kk];
- get_velocity(E,j,nelem,theta0,phi0,rad0,velocity_vector);
+ if (E->parallel.nprocz>1) {
+ ithatcap=0;
+ icheck=full_icheck_cap(E,ithatcap,x,y,z,rad);
+ if (icheck==1) goto foundit;
- x_pred=x0+velocity_vector[1]*dt;
- y_pred=y0+velocity_vector[2]*dt;
- z_pred=z0+velocity_vector[3]*dt;
+ }
+ /* check neighboring caps */
- /* keep in box */
+ for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++) {
+ ithatcap=pp;
+ icheck=full_icheck_cap(E,ithatcap,x,y,z,rad);
+ if (icheck==1) goto foundit;
+ }
- cart_to_sphere(E,x_pred,y_pred,z_pred,&theta_pred,&phi_pred,&rad_pred);
- keep_in_sphere(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
- /* Current Coordinates are always kept in positions 0-5. */
+ /* should not be here */
+ if (icheck!=1) {
+ fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
+ fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
+ icheck=full_icheck_cap(E,0,x,y,z,rad);
+ if (icheck==1) fprintf(E->trace.fpt," icheck here!\n");
+ else fprintf(E->trace.fpt,"icheck not here!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
- 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;
+ foundit:
- /* Fill in original coords (positions 6-8) */
+ isend[j][ithatcap]++;
- E->trace.basicq[j][6][kk]=x0;
- E->trace.basicq[j][7][kk]=y0;
- E->trace.basicq[j][8][kk]=z0;
+ /* assign tracer to send */
- /* Fill in original velocities (positions 9-11) */
+ isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
- 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 */
+ for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++) {
+ ipos=isend_position+pp;
+ send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
+ }
+ } /* end kk, assigning tracers */
- } /* end kk, predicting tracers */
- } /* end caps */
-
- /* find new tracer elements and caps */
-
- find_tracers(E);
-
+ } /* end j */
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 */
-/* regardless of advection scheme]. */
-/* 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 get_velocity();
- void cart_to_sphere();
- void keep_in_sphere();
- void find_tracers();
-
-
- 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];
-
- 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);
- keep_in_sphere(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
-
- /* Fill in Current Positions (other positions are no longer important) */
-
- 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;
-}
-
/******** GET VELOCITY ***************************************/
-void get_velocity(E,j,nelem,theta,phi,rad,velocity_vector)
- struct All_variables *E;
- int j,nelem;
- double theta,phi,rad;
- double *velocity_vector;
+void full_get_velocity(struct All_variables *E,
+ int j, int nelem,
+ double theta, double phi, double rad,
+ double *velocity_vector)
{
void gnomonic_interpolation();
@@ -687,7 +406,6 @@
void spherical_to_uv();
void get_2dshape();
void velo_from_element_d();
- int iget_element();
int icheck_element();
int icheck_column_neighbors();
@@ -739,7 +457,7 @@
sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
ival=icheck_element(E,j,nelem,x,y,z,rad);
fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
- ival=iget_element(E,j,-99,x,y,z,theta,phi,rad);
+ ival=(E->trace.iget_element)(E,j,-99,x,y,z,theta,phi,rad);
fprintf(E->trace.fpt,"New Element?: %d\n",ival);
ival=icheck_column_neighbors(E,j,nelem,x,y,z,rad);
fprintf(E->trace.fpt,"New Element (neighs)?: %d\n",ival);
@@ -989,1201 +707,6 @@
}
-/**************** INITIALIZE TRACER ARRAYS ************************************/
-/* */
-/* This function allocates memories to tracer arrays. */
-
-void initialize_tracer_arrays(E,j,number_of_tracers)
- struct All_variables *E;
- int j, 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*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;
-}
-
-
-
-/************ FIND TRACERS *************************************/
-/* */
-/* This function finds tracer elements and moves tracers to */
-/* other processor domains if necessary. */
-/* Array ielement is filled with elemental values. */
-
-void find_tracers(E)
- 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;
-
- int iget_element();
- void put_away_later();
- void eject_tracer();
- void reduce_tracer_arrays();
- void lost_souls();
- void sphere_to_cart();
-
- time_stat1=CPU_time0();
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
-
- /* initialize arrays and statistical counters */
-
- E->trace.ilater[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];
-
- /* AKMA REMOVE */
- /*
- fprintf(E->trace.fpt,"BB. kk %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,x,y,z,theta,phi,rad);
- fflush(E->trace.fpt);
- */
- iel=iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad);
-
- E->trace.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();
- */
-
- lost_souls(E);
-
- /* Free later arrays */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- if (E->trace.ilater[j]>0)
- {
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
- {
- free(E->trace.rlater[j][kk]);
- }
- }
- } /* end j */
-
-
- /* Adjust Array Sizes */
-
- reduce_tracer_arrays(E);
-
- time_stat2=CPU_time0();
-
- fprintf(E->trace.fpt,"AA: time for find tracers: %f\n", time_stat2-time_stat1);
-
- return;
-}
-
-/************** LOST SOULS ****************************************************/
-/* */
-/* This function is used to transport tracers to proper processor domains. */
-/* (MPI parallel) */
-/* All of the tracers that were sent to rlater arrays are destined to another */
-/* cap and sent there. Then they are raised up or down for multiple z procs. */
-/* isend[j][n]=number of tracers this processor cap is sending to cap n */
-/* ireceive[j][n]=number of tracers this processor cap is receiving from cap n */
-
-
-void lost_souls(E)
- struct All_variables *E;
-{
- int ithiscap;
- int ithatcap=1;
- int isend[13][13];
- int ireceive[13][13];
- int isize[13];
- int kk,pp,j;
- int mm;
- int numtracers;
- int icheck=0;
- int isend_position;
- int ipos,ipos2,ipos3;
- int idb;
- int idestination_proc=0;
- int isource_proc;
- int isend_z[13][3];
- int ireceive_z[13][3];
- int isum[13];
- int irad;
- int ival;
- int ithat_processor;
- int ireceive_position;
- int ihorizontal_neighbor;
- int ivertical_neighbor;
- int ilast_receiver_position;
- int it;
- int irec[13];
- int irec_position;
- int iel;
- int num_tracers;
- int isize_send;
- int isize_receive;
- int itemp_size;
- int itracers_subject_to_vertical_transport[13];
-
- double x,y,z;
- double theta,phi,rad;
- double *send[13][13];
- double *receive[13][13];
- double *send_z[13][3];
- double *receive_z[13][3];
- double *REC[13];
-
- int iget_element();
- int icheck_cap();
- void expand_tracer_arrays();
-
- int number_of_caps=12;
- int lev=E->mesh.levmax;
- int num_ngb;
-
- /* Note, if for some reason, the number of neighbors exceeds */
- /* 50, which is unlikely, the MPI arrays must be increased. */
- MPI_Status status[200];
- MPI_Request request[200];
- MPI_Status status1;
- MPI_Status status2;
- int itag=1;
-
-
- parallel_process_sync(E);
- fprintf(E->trace.fpt, "Entering lost_souls()\n");
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- E->trace.istat_isend=E->trace.ilater[j];
- }
-
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
- for (kk=1; kk<=E->trace.istat_isend; kk++) {
- fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
- E->trace.rlater[j][0][kk],
- E->trace.rlater[j][1][kk],
- E->trace.rlater[j][2][kk]);
- }
- fflush(E->trace.fpt);
- }
-
- /* initialize isend and ireceive */
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- /* # of neighbors in the horizontal plane */
- num_ngb = E->parallel.TNUM_PASS[lev][j];
- isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
- for (kk=1;kk<=num_ngb;kk++) isend[j][kk]=0;
- for (kk=1;kk<=num_ngb;kk++) ireceive[j][kk]=0;
- }
-
- /* Allocate Maximum Memory to Send Arrays */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- ithiscap=0;
-
- itemp_size=max(isize[j],1);
-
- if ((send[j][ithiscap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- num_ngb = E->parallel.TNUM_PASS[lev][j];
- for (kk=1;kk<=num_ngb;kk++)
- {
- if ((send[j][kk]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
-
-
- /* For testing, remove this */
- /*
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- ithiscap=E->sphere.capid[j];
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=E->parallel.PROCESSOR[lev][j].pass[kk]+1;
- fprintf(E->trace.fpt,"cap: %d proc %d TNUM: %d ithatcap: %d\n",
- ithiscap,E->parallel.me,kk,ithatcap);
-
- }
- fflush(E->trace.fpt);
- }
- */
-
-
- /* Pre communication */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- /* transfer tracers from rlater to send */
-
- 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)
- {
- ithatcap=0;
- icheck=icheck_cap(E,ithatcap,x,y,z,rad);
- if (icheck==1) goto foundit;
-
- }
-
- /* check neighboring caps */
-
- for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++)
- {
- ithatcap=pp;
- icheck=icheck_cap(E,ithatcap,x,y,z,rad);
- if (icheck==1) goto foundit;
- }
-
-
- /* should not be here */
- if (icheck!=1)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
- fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
- icheck=icheck_cap(E,0,x,y,z,rad);
- if (icheck==1) fprintf(E->trace.fpt," icheck here!\n");
- else fprintf(E->trace.fpt,"icheck not here!\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- foundit:
-
- isend[j][ithatcap]++;
-
- /* assign tracer to send */
-
- isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
-
- for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++)
- {
- ipos=isend_position+pp;
- send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
- }
-
- } /* end kk, assigning tracers */
-
- } /* end j */
-
-
- /* Send info to other processors regarding number of send tracers */
-
- /* idb is the request array index variable */
- /* Each send and receive has a request variable */
-
- idb=0;
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- ithiscap=0;
-
- /* if tracer is in same cap (nprocz>1) */
-
- if (E->parallel.nprocz>1)
- {
- ireceive[j][ithiscap]=isend[j][ithiscap];
- }
-
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
-
- /* if neighbor cap is in another processor, send information via MPI */
-
- idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
-
- idb++;
- MPI_Isend(&isend[j][ithatcap],1,MPI_INT,idestination_proc,
- 11,E->parallel.world,
- &request[idb-1]);
-
- } /* end kk, number of neighbors */
-
- } /* end j, caps per proc */
-
- /* Receive tracer count info */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
-
- /* if neighbor cap is in another processor, receive information via MPI */
-
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
-
- if (idestination_proc!=E->parallel.me)
- {
-
- idb++;
- MPI_Irecv(&ireceive[j][ithatcap],1,MPI_INT,isource_proc,
- 11,E->parallel.world,
- &request[idb-1]);
-
- } /* end if */
-
- } /* end kk, number of neighbors */
- } /* end j, caps per proc */
-
- /* Wait for non-blocking calls to complete */
-
- MPI_Waitall(idb,request,status);
-
- /* Testing, should remove */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- fprintf(E->trace.fpt,"j: %d send %d to cap %d\n",j,isend[j][kk],isource_proc);
- fprintf(E->trace.fpt,"j: %d rec %d from cap %d\n",j,ireceive[j][kk],isource_proc);
- }
- }
-
-
- /* Allocate memory in receive arrays */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- num_ngb = E->parallel.TNUM_PASS[lev][j];
- for (ithatcap=1;ithatcap<=num_ngb;ithatcap++)
- {
- isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
-
- itemp_size=max(1,isize[j]);
-
- if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
-
- /* Now, send the tracers to proper caps */
-
- idb=0;
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- ithiscap=0;
-
- /* same cap */
-
- if (E->parallel.nprocz>1)
- {
-
- ithatcap=ithiscap;
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
- for (mm=0;mm<=(isize[j]-1);mm++)
- {
- receive[j][ithatcap][mm]=send[j][ithatcap][mm];
- }
-
- }
-
- /* neighbor caps */
-
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
-
- idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
-
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
-
- idb++;
-
- MPI_Isend(send[j][ithatcap],isize[j],MPI_DOUBLE,idestination_proc,
- 11,E->parallel.world,
- &request[idb-1]);
-
- } /* end kk, number of neighbors */
-
- } /* end j, caps per proc */
-
-
- /* Receive tracers */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- ithiscap=0;
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
-
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
-
- idb++;
-
- isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
-
- MPI_Irecv(receive[j][ithatcap],isize[j],MPI_DOUBLE,isource_proc,
- 11,E->parallel.world,
- &request[idb-1]);
-
- } /* end kk, number of neighbors */
-
- } /* end j, caps per proc */
-
- /* Wait for non-blocking calls to complete */
-
- MPI_Waitall(idb,request,status);
-
-
- /* Put all received tracers in array REC[j] */
- /* This makes things more convenient. */
-
- /* Sum up size of receive arrays (all tracers sent to this processor) */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isum[j]=0;
-
- ithiscap=0;
-
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
- isum[j]=isum[j]+ireceive[j][ithatcap];
- }
- if (E->parallel.nprocz>1) isum[j]=isum[j]+ireceive[j][ithiscap];
-
- itracers_subject_to_vertical_transport[j]=isum[j];
- }
-
- /* Allocate Memory for REC array */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
- if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- REC[j][0]=0.0;
- }
-
- /* Put Received tracers in REC */
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- irec[j]=0;
-
- irec_position=0;
-
- ithiscap=0;
-
- /* horizontal neighbors */
-
- for (ihorizontal_neighbor=1;ihorizontal_neighbor<=E->parallel.TNUM_PASS[lev][j];ihorizontal_neighbor++)
- {
-
- ithatcap=ihorizontal_neighbor;
-
- for (pp=1;pp<=ireceive[j][ithatcap];pp++)
- {
- irec[j]++;
- ipos=(pp-1)*E->trace.number_of_tracer_quantities;
-
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos2=ipos+mm;
- REC[j][irec_position]=receive[j][ithatcap][ipos2];
-
- irec_position++;
-
- } /* end mm (cycling tracer quantities) */
- } /* end pp (cycling tracers) */
- } /* end ihorizontal_neighbors (cycling neighbors) */
-
- /* for tracers in the same cap (nprocz>1) */
-
- if (E->parallel.nprocz>1)
- {
- ithatcap=ithiscap;
- for (pp=1;pp<=ireceive[j][ithatcap];pp++)
- {
- irec[j]++;
- ipos=(pp-1)*E->trace.number_of_tracer_quantities;
-
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos2=ipos+mm;
-
- REC[j][irec_position]=receive[j][ithatcap][ipos2];
-
- irec_position++;
-
- } /* end mm (cycling tracer quantities) */
-
- } /* end pp (cycling tracers) */
-
- } /* endif nproc>1 */
-
- } /* end j (cycling caps) */
-
- /* Done filling REC */
-
-
-
- /* VERTICAL COMMUNICATION */
-
- /* Note: For generality, I assume that both multiple */
- /* caps per processor as well as multiple processors */
- /* in the radial direction. These are probably */
- /* inconsistent parameter choices for the regular */
- /* CitcomS code. */
-
- if (E->parallel.nprocz>1)
- {
-
- /* Allocate memory for send_z */
- /* Make send_z the size of receive array (max size) */
- /* (No dynamic reallocation of send_z necessary) */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
- {
- isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
-
- if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
-
- ithat_processor=E->parallel.PROCESSORz[lev].pass[ivertical_neighbor];
-
- /* initialize isend_z and ireceive_z array */
-
- isend_z[j][ivertical_neighbor]=0;
- ireceive_z[j][ivertical_neighbor]=0;
-
- /* sort through receive array and check radius */
-
- it=0;
- num_tracers=irec[j];
- for (kk=1;kk<=num_tracers;kk++)
- {
-
- it++;
-
- ireceive_position=((it-1)*E->trace.number_of_tracer_quantities);
-
- irad=ireceive_position+2;
-
- rad=REC[j][irad];
-
- ival=icheck_that_processor_shell(E,j,ithat_processor,rad);
-
-
- /* if tracer is in other shell, take out of receive array and give to send_z*/
-
- if (ival==1)
- {
-
- isend_z[j][ivertical_neighbor]++;
-
- isend_position=(isend_z[j][ivertical_neighbor]-1)*E->trace.number_of_tracer_quantities;
-
- ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
-
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- ipos=ireceive_position+mm;
- ipos2=isend_position+mm;
-
- send_z[j][ivertical_neighbor][ipos2]=REC[j][ipos];
-
-
- /* eject tracer info from REC array, and replace with last tracer in array */
-
- ipos3=ilast_receiver_position+mm;
- REC[j][ipos]=REC[j][ipos3];
-
- }
-
- it--;
- irec[j]--;
-
- } /* end if ival===1 */
-
- /* Otherwise, leave tracer */
-
- } /* end kk (cycling through tracers) */
-
- } /* end ivertical_neighbor */
-
- } /* end j */
-
-
- /* Send arrays are now filled. */
- /* Now send send information to vertical processor neighbor */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
-
- MPI_Sendrecv(&isend_z[j][ivertical_neighbor],1,MPI_INT,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag,
- &ireceive_z[j][ivertical_neighbor],1,MPI_INT,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- itag,E->parallel.world,&status1);
-
- /* for testing - remove */
- /*
- fprintf(E->trace.fpt,"PROC: %d IVN: %d (P: %d) SEND: %d REC: %d\n",
- E->parallel.me,ivertical_neighbor,E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- isend_z[j][ivertical_neighbor],ireceive_z[j][ivertical_neighbor]);
- fflush(E->trace.fpt);
- */
-
- } /* end ivertical_neighbor */
-
- } /* end j */
-
-
- /* Allocate memory to receive_z arrays */
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
- {
- isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
-
- if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end j */
-
- /* Send Tracers */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- isize_send=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
- isize_receive=ireceive_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
-
- MPI_Sendrecv(send_z[j][ivertical_neighbor],isize_send,
- MPI_DOUBLE,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag+1,
- receive_z[j][ivertical_neighbor],isize_receive,
- MPI_DOUBLE,
- E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
- itag+1,E->parallel.world,&status2);
-
- }
- }
-
- /* Put tracers into REC array */
-
- /* First, reallocate memory to REC */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- isum[j]=0;
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- isum[j]=isum[j]+ireceive_z[j][ivertical_neighbor];
- }
-
- isum[j]=isum[j]+irec[j];
-
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
-
- if (isize[j]>0)
- {
- if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(lost souls)-no memory (i981)\n");
- fprintf(E->trace.fpt,"isize: %d\n",isize[j]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
-
- for (kk=1;kk<=ireceive_z[j][ivertical_neighbor];kk++)
- {
- irec[j]++;
-
- irec_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
- ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
-
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
- {
- REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
- }
- }
-
- }
- }
-
- /* Free Vertical Arrays */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
- {
- free(send_z[j][ivertical_neighbor]);
- free(receive_z[j][ivertical_neighbor]);
- }
- }
-
- } /* endif nprocz>1 */
-
- /* END OF VERTICAL TRANSPORT */
-
- /* Put away tracers */
-
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
- for (kk=1;kk<=irec[j];kk++)
- {
- E->trace.ntracers[j]++;
-
- if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
-
- ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
-
- for (mm=0;mm<=(E->trace.number_of_basic_quantities-1);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-1);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=iget_element(E,j,-99,x,y,z,theta,phi,rad);
-
- if (iel<1)
- {
- fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
- fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",x,y,z,theta,phi,rad);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- E->trace.ielement[j][E->trace.ntracers[j]]=iel;
-
- }
- }
-
- fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
- fflush(E->trace.fpt);
- parallel_process_sync(E);
-
- /* Free Arrays */
-
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
-
- ithiscap=0;
-
- free(REC[j]);
-
- free(send[j][ithiscap]);
-
- for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
- {
- ithatcap=kk;
-
- free(send[j][ithatcap]);
- free(receive[j][ithatcap]);
-
- }
-
- }
- fprintf(E->trace.fpt,"Leaving lost_souls()\n");
- fflush(E->trace.fpt);
-
- return;
-}
-
-
-/****** REDUCE TRACER ARRAYS *****************************************/
-
-void reduce_tracer_arrays(E)
- 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 */
-
-
-void put_away_later(E,j,it)
- struct All_variables *E;
- int it,j;
-
-{
-
-
- int kk;
-
-
- void expand_later_array();
-
-
- /* The first tracer in initiates memory allocation. */
- /* Memory is freed after parallel communications */
-
- if (E->trace.ilater[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(E,j)
- 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 ************************************************/
-
-
-void eject_tracer(E,j,it)
- struct All_variables *E;
- int it,j;
-
-{
-
- 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;
-}
-
-
-
/*********** MAKE REGULAR GRID ********************************/
/* */
/* This function generates the finer regular grid which is */
@@ -2833,121 +1356,6 @@
}
-
-/********* ICHECK COLUMN NEIGHBORS ***************************/
-/* */
-/* This function check whether a point is in a neighboring */
-/* column. Neighbor surface element number is returned */
-
-int icheck_column_neighbors(E,j,nel,x,y,z,rad)
- struct All_variables *E;
- int j,nel;
- double x,y,z,rad;
-{
-
- int ival;
- int neighbor[25];
- int elx,ely,elz;
- int elxz;
- int kk;
-
- int icheck_element_column();
-
- /*
- const int number_of_neighbors=24;
- */
-
- /* maybe faster to only check inner ring */
-
- const int number_of_neighbors=8;
-
- elx=E->lmesh.elx;
- ely=E->lmesh.ely;
- elz=E->lmesh.elz;
-
- elxz=elx*elz;
-
- /* inner ring */
-
- neighbor[1]=nel-elxz-elz;
- neighbor[2]=nel-elxz;
- neighbor[3]=nel-elxz+elz;
- neighbor[4]=nel-elz;
- neighbor[5]=nel+elz;
- neighbor[6]=nel+elxz-elz;
- neighbor[7]=nel+elxz;
- neighbor[8]=nel+elxz+elz;
-
- /* outer ring */
-
- neighbor[9]=nel+2*elxz-elz;
- neighbor[10]=nel+2*elxz;
- neighbor[11]=nel+2*elxz+elz;
- neighbor[12]=nel+2*elxz+2*elz;
- neighbor[13]=nel+elxz+2*elz;
- neighbor[14]=nel+2*elz;
- neighbor[15]=nel-elxz+2*elz;
- neighbor[16]=nel-2*elxz+2*elz;
- neighbor[17]=nel-2*elxz+elz;
- neighbor[18]=nel-2*elxz;
- neighbor[19]=nel-2*elxz-elz;
- neighbor[20]=nel-2*elxz-2*elz;
- neighbor[21]=nel-elxz-2*elz;
- neighbor[22]=nel-2*elz;
- neighbor[23]=nel+elxz-2*elz;
- neighbor[24]=nel+2*elxz-2*elz;
-
-
- for (kk=1;kk<=number_of_neighbors;kk++)
- {
-
- if ((neighbor[kk]>=1)&&(neighbor[kk]<=E->lmesh.nel))
- {
- ival=icheck_element_column(E,j,neighbor[kk],x,y,z,rad);
- if (ival>0)
- {
- return neighbor[kk];
- }
- }
- }
-
- return -99;
-}
-
-/********** ICHECK ALL COLUMNS ********************************/
-/* */
-/* This function check all columns until the proper one for */
-/* a point (x,y,z) is found. The surface element is returned */
-/* else -99 if can't be found. */
-
-int icheck_all_columns(E,j,x,y,z,rad)
- struct All_variables *E;
- int j;
- double x,y,z,rad;
-{
-
- int icheck;
- int nel;
- int icheck_element_column();
-
- int elz=E->lmesh.elz;
- int numel=E->lmesh.nel;
-
- for (nel=elz;nel<=numel;nel=nel+elz)
- {
- icheck=icheck_element_column(E,j,nel,x,y,z,rad);
- if (icheck==1)
- {
- return nel;
- }
- }
-
-
- return -99;
-}
-
-
-
/**** WRITE TRACE INSTRUCTIONS ***************/
void write_trace_instructions(E)
struct All_variables *E;
@@ -3067,112 +1475,6 @@
}
-/************** RESTART TRACERS ******************************************/
-/* */
-/* This function restarts tracers written from previous calculation */
-/* and the tracers are read as seperate files for each processor domain. */
-
-void restart_tracers(E)
- struct All_variables *E;
-{
-
- char output_file[200];
- char input_s[1000];
-
- int i,j,kk;
- int idum1,ncolumns;
- int numtracers;
-
- double rdum1;
- double theta,phi,rad;
- double extra[100];
- double x,y,z;
-
- void initialize_tracer_arrays();
- void sphere_to_cart();
- void keep_in_sphere();
-
- FILE *fp1;
-
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(restart_tracers)-increase size of extra[]\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
- sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
-
- if ( (fp1=fopen(output_file,"r"))==NULL) {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- fprintf(stderr,"Restarting Tracers from %s\n",output_file);
- fflush(stderr);
-
-
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fgets(input_s,200,fp1);
- sscanf(input_s,"%d %d %d %lf",
- &idum1, &numtracers, &ncolumns, &rdum1);
-
- /* some error control */
- if (E->trace.number_of_extra_quantities+3 != ncolumns) {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- /* allocate memory for tracer arrays */
-
- initialize_tracer_arrays(E,j,numtracers);
- E->trace.ntracers[j]=numtracers;
-
- for (kk=1;kk<=numtracers;kk++) {
- fgets(input_s,200,fp1);
- if (E->trace.number_of_extra_quantities==0) {
- sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
- }
- else if (E->trace.number_of_extra_quantities==1) {
- sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&extra[0]);
- }
- /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
- /* this part has to be changed... */
- else {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- 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 */
-
- keep_in_sphere(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]=extra[i];
-
- }
-
- fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
- fflush(E->trace.fpt);
-
- }
- fclose(fp1);
-
-
- return;
-}
-
/************** MAKE TRACER ARRAY ********************************/
/* Here, each cap will generate tracers somewhere */
/* in the sphere - check if its in this cap - then check radial */
@@ -3180,25 +1482,12 @@
static void make_tracer_array(struct All_variables *E)
{
- int kk;
int tracers_cap;
int j;
- int ival;
- int number_of_tries=0;
- int max_tries;
-
- double x,y,z;
- double theta,phi,rad;
- double dmin,dmax;
- double random1,random2,random3;
double processor_fraction;
- void cart_to_sphere();
- void keep_in_sphere();
- void initialize_tracer_arrays();
+ void generate_random_tracers();
- int icheck_cap();
-
if (E->parallel.me==0) fprintf(stderr,"Making Tracer Array\n");
fflush(stderr);
@@ -3214,69 +1503,10 @@
fprintf(E->trace.fpt,"\nGenerating %d Tracers\n",tracers_cap);
- initialize_tracer_arrays(E,j,tracers_cap);
+ generate_random_tracers(E, tracers_cap, j);
- /* Tracers are placed randomly in cap */
- /* (intentionally using rand() instead of srand() )*/
- dmin=-1.0*E->sphere.ro;
- dmax=E->sphere.ro;
-
- while (E->trace.ntracers[j]<tracers_cap) {
-
- number_of_tries++;
- max_tries=500*tracers_cap;
-
- if (number_of_tries>max_tries) {
- fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
- fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
- fflush(E->trace.fpt);
- exit(10);
- }
-
-
- random1=(1.0*rand())/(1.0*RAND_MAX);
- random2=(1.0*rand())/(1.0*RAND_MAX);
- random3=(1.0*rand())/(1.0*RAND_MAX);
-
- x=dmin+random1*(dmax-dmin);
- y=dmin+random2*(dmax-dmin);
- z=dmin+random3*(dmax-dmin);
-
- /* first check if within shell */
-
- rad=sqrt(x*x+y*y+z*z);
-
- if (rad>=E->sx[j][3][E->lmesh.noz]) continue;
- if (rad<E->sx[j][3][1]) continue;
-
-
- /* check if in current cap */
-
- ival=icheck_cap(E,0,x,y,z,rad);
-
- if (ival!=1) continue;
-
- /* Made it, so record tracer information */
-
- cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
-
- keep_in_sphere(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 */
-
-
}/* end j */
@@ -3290,269 +1520,123 @@
return;
}
-/******** READ TRACER ARRAY *********************************************/
-/* */
-/* This function reads tracers from input file. */
-/* All processors read the same input file, then sort out which ones */
-/* belong. */
-void read_tracer_file(E)
- struct All_variables *E;
-{
- char input_s[1000];
- int number_of_tracers, ncolumns;
- int kk;
- int icheck;
- int iestimate;
- int icushion;
- int i, j;
+/********* ICHECK COLUMN NEIGHBORS ***************************/
+/* */
+/* This function check whether a point is in a neighboring */
+/* column. Neighbor surface element number is returned */
- int icheck_cap();
- int icheck_processor_shell();
- int isum_tracers();
- void initialize_tracer_arrays();
- void keep_in_sphere();
- void sphere_to_cart();
- void cart_to_sphere();
- void expand_tracer_arrays();
-
- double x,y,z;
- double theta,phi,rad;
- double extra[100];
-
- FILE *fptracer;
-
- fptracer=fopen(E->trace.tracer_file,"r");
- fprintf(E->trace.fpt,"Opening %s\n",E->trace.tracer_file);
-
- fgets(input_s,200,fptracer);
- sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns);
- 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++) {
-
- initialize_tracer_arrays(E,j,iestimate);
-
- for (kk=1;kk<=number_of_tracers;kk++) {
- fgets(input_s,200,fptracer);
- if (E->trace.number_of_extra_quantities==0) {
- sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
- }
- else if (E->trace.number_of_extra_quantities==1) {
- sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&extra[0]);
- }
- /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
- /* this part has to be changed... */
- else {
- fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
-
-
- /* make sure theta, phi is in range, and radius is within bounds */
-
- keep_in_sphere(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;
-
- icheck=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]]=extra[i];
-
- } /* end kk, number of tracers */
-
- fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
- E->trace.ntracers[j]);
-
- } /* end j */
-
- fclose(fptracer);
-
- icheck=isum_tracers(E);
-
- if (icheck!=number_of_tracers) {
- fprintf(E->trace.fpt,"ERROR(read_tracer_file) - tracers != number in file\n");
- fprintf(E->trace.fpt,"Tracers in system: %d\n", icheck);
- fprintf(E->trace.fpt,"Tracers in file: %d\n", number_of_tracers);
- fflush(E->trace.fpt);
- exit(10);
- }
-
- return;
-}
-
-
-/********** CART TO SPHERE ***********************/
-void cart_to_sphere(E,x,y,z,theta,phi,rad)
+int icheck_column_neighbors(E,j,nel,x,y,z,rad)
struct All_variables *E;
- double x,y,z;
- double *theta,*phi,*rad;
+ int j,nel;
+ double x,y,z,rad;
{
- double temp;
- double myatan();
+ int ival;
+ int neighbor[25];
+ int elx,ely,elz;
+ int elxz;
+ int kk;
- temp=x*x+y*y;
+ int icheck_element_column();
- *rad=sqrt(temp+z*z);
- *theta=atan2(sqrt(temp),z);
- *phi=myatan(y,x);
+ /*
+ const int number_of_neighbors=24;
+ */
+ /* maybe faster to only check inner ring */
- return;
-}
+ const int number_of_neighbors=8;
-/********** SPHERE TO CART ***********************/
-void sphere_to_cart(E,theta,phi,rad,x,y,z)
- struct All_variables *E;
- double theta,phi,rad;
- double *x,*y,*z;
-{
+ elx=E->lmesh.elx;
+ ely=E->lmesh.ely;
+ elz=E->lmesh.elz;
- double sint,cost,sinf,cosf;
- double temp;
+ elxz=elx*elz;
- sint=sin(theta);
- cost=cos(theta);
- sinf=sin(phi);
- cosf=cos(phi);
+ /* inner ring */
- temp=rad*sint;
+ neighbor[1]=nel-elxz-elz;
+ neighbor[2]=nel-elxz;
+ neighbor[3]=nel-elxz+elz;
+ neighbor[4]=nel-elz;
+ neighbor[5]=nel+elz;
+ neighbor[6]=nel+elxz-elz;
+ neighbor[7]=nel+elxz;
+ neighbor[8]=nel+elxz+elz;
- *x=temp*cosf;
- *y=temp*sinf;
- *z=rad*cost;
+ /* outer ring */
+ neighbor[9]=nel+2*elxz-elz;
+ neighbor[10]=nel+2*elxz;
+ neighbor[11]=nel+2*elxz+elz;
+ neighbor[12]=nel+2*elxz+2*elz;
+ neighbor[13]=nel+elxz+2*elz;
+ neighbor[14]=nel+2*elz;
+ neighbor[15]=nel-elxz+2*elz;
+ neighbor[16]=nel-2*elxz+2*elz;
+ neighbor[17]=nel-2*elxz+elz;
+ neighbor[18]=nel-2*elxz;
+ neighbor[19]=nel-2*elxz-elz;
+ neighbor[20]=nel-2*elxz-2*elz;
+ neighbor[21]=nel-elxz-2*elz;
+ neighbor[22]=nel-2*elz;
+ neighbor[23]=nel+elxz-2*elz;
+ neighbor[24]=nel+2*elxz-2*elz;
- return;
-}
+ for (kk=1;kk<=number_of_neighbors;kk++)
+ {
-/********* 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). */
+ if ((neighbor[kk]>=1)&&(neighbor[kk]<=E->lmesh.nel))
+ {
+ ival=icheck_element_column(E,j,neighbor[kk],x,y,z,rad);
+ if (ival>0)
+ {
+ return neighbor[kk];
+ }
+ }
+ }
-int icheck_that_processor_shell(E,j,nprocessor,rad)
- 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 -99;
}
-/********** 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 */
+/********** ICHECK ALL COLUMNS ********************************/
+/* */
+/* This function check all columns until the proper one for */
+/* a point (x,y,z) is found. The surface element is returned */
+/* else -99 if can't be found. */
-int icheck_processor_shell(E,j,rad)
+int icheck_all_columns(E,j,x,y,z,rad)
struct All_variables *E;
- double rad;
int j;
-
+ double x,y,z,rad;
{
- const int noz = E->lmesh.noz;
- const int nprocz = E->parallel.nprocz;
- double top_r, bottom_r;
+ int icheck;
+ int nel;
+ int icheck_element_column();
- if (nprocz==1) return 1;
+ int elz=E->lmesh.elz;
+ int numel=E->lmesh.nel;
- top_r = E->sx[j][3][noz];
- bottom_r = E->sx[j][3][1];
+ for (nel=elz;nel<=numel;nel=nel+elz)
+ {
+ icheck=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (icheck==1)
+ {
+ return nel;
+ }
+ }
- /* 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;
+ return -99;
}
+
/******* ICHECK ELEMENT *************************************/
/* */
/* This function serves to determine if a point lies within */
@@ -3678,10 +1762,8 @@
/* This function serves to determine if a point lies within */
/* a given cap */
/* */
-int icheck_cap(E,icap,x,y,z,rad)
- struct All_variables *E;
- int icap;
- double x,y,z,rad;
+int full_icheck_cap(struct All_variables *E, int icap,
+ double x, double y, double z, double rad)
{
double test_point[4];
@@ -4016,12 +2098,10 @@
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
-int iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad)
- struct All_variables *E;
- int j;
- int iprevious_element;
- double x,y,z;
- double theta,phi,rad;
+int full_iget_element(struct All_variables *E,
+ int j, int iprevious_element,
+ double x, double y, double z,
+ double theta, double phi, double rad)
{
int iregel;
@@ -4037,9 +2117,6 @@
int nelem;
int iget_regel();
- int iquick_element_column_search();
- int icheck_cap();
- int icheck_regular_neighbors();
int iget_radial_element();
elx=E->lmesh.elx;
@@ -4133,7 +2210,7 @@
/* check if still in cap */
- ival=icheck_cap(E,0,x,y,z,rad);
+ ival=full_icheck_cap(E,0,x,y,z,rad);
if (ival==0)
{
return -99;
@@ -4176,17 +2253,6 @@
}
}
- /* Decided to remove this part - not really needed and complicated */
- /*
- else
- {
- iel=icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad);
- if (iel>0)
- {
- goto foundit;
- }
- }
- */
}
/* As a last resort, check all element columns */
@@ -4196,7 +2262,7 @@
iel=icheck_all_columns(E,j,x,y,z,rad);
/*
- fprintf(E->trace.fpt,"WARNING(iget_element)-doing a full search!\n");
+ fprintf(E->trace.fpt,"WARNING(full_iget_element)-doing a full search!\n");
fprintf(E->trace.fpt," Most often means tracers have moved more than 1 element away\n");
fprintf(E->trace.fpt," or regular element resolution is way too low.\n");
fprintf(E->trace.fpt," COLUMN: %d \n",iel);
@@ -4221,7 +2287,7 @@
/* if still here, there is a problem */
- fprintf(E->trace.fpt,"Error(iget_element) - element not found\n");
+ fprintf(E->trace.fpt,"Error(full_iget_element) - element not found\n");
fprintf(E->trace.fpt,"x,y,z,theta,phi,iregel %f %f %f %f %f %d\n",
x,y,z,theta,phi,iregel);
fflush(E->trace.fpt);
@@ -4288,185 +2354,7 @@
return iradial_element;
}
-/****** ICHECK REGULAR NEIGHBORS *****************************/
-/* */
-/* This function searches the regular element neighborhood. */
-/* This function is no longer used! */
-
-int icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad)
- struct All_variables *E;
- int j,ntheta,nphi;
- double x,y,z;
- double theta,phi,rad;
-{
-
- int new_ntheta,new_nphi;
- int kk,pp;
- int iregel;
- int ival;
- int imap[5];
- int ichoice;
- int irange;
-
- int iquick_element_column_search();
-
- fprintf(E->trace.fpt,"ERROR(icheck_regular_neighbors)-this subroutine is no longer used !\n");
- fflush(E->trace.fpt);
- exit(10);
-
- irange=2;
-
- for (kk=-irange;kk<=irange;kk++)
- {
- for (pp=-irange;pp<=irange;pp++)
- {
- new_ntheta=ntheta+kk;
- new_nphi=nphi+pp;
- if ( (new_ntheta>0)&&(new_ntheta<=E->trace.numtheta[j])&&(new_nphi>0)&&(new_nphi<=E->trace.numphi[j]) )
- {
- iregel=new_ntheta+(new_nphi-1)*E->trace.numtheta[j];
- if ((iregel>0) && (iregel<=E->trace.numregel[j]))
- {
- ival=iquick_element_column_search(E,j,iregel,new_ntheta,new_nphi,x,y,z,theta,phi,rad,imap,&ichoice);
- if (ival>0) return ival;
- }
- }
- }
- }
-
-
- return -99;
-}
-
-
-
-
-
-/****** IQUICK ELEMENT SEARCH *****************************/
-/* */
-/* This function does a quick regular to real element */
-/* map check. Element number, if found, is returned. */
-/* Otherwise, -99 is returned. */
-/* Pointers to imap and ichoice are used because they may */
-/* prove to be convenient. */
-/* This routine is no longer used */
-
-int iquick_element_column_search(E,j,iregel,ntheta,nphi,x,y,z,theta,phi,rad,imap,ich)
- struct All_variables *E;
- int j,iregel;
- int ntheta,nphi;
- double x,y,z,theta,phi,rad;
- int *imap;
- int *ich;
-{
-
- int iregnode[5];
- int kk,pp;
- int nel,ival;
- int ichoice;
- int icount;
- int itemp1;
- int itemp2;
-
- int icheck_element_column();
-
- fprintf(E->trace.fpt,"ERROR(iquick element)-this routine is no longer used!\n");
- fflush(E->trace.fpt);
- exit(10);
-
- /* REMOVE*/
- /*
- ichoice=*ich;
-
- fprintf(E->trace.fpt,"AA: ichoice: %d\n",ichoice);
- fflush(E->trace.fpt);
- */
-
- /* find regular nodes on regular element */
-
- /*
- iregnode[1]=iregel+(nphi-1);
- iregnode[2]=iregel+nphi;
- iregnode[3]=iregel+nphi+E->trace.numtheta[j]+1;
- iregnode[4]=iregel+nphi+E->trace.numtheta[j];
- */
-
- itemp1=iregel+nphi;
- itemp2=itemp1+E->trace.numtheta[j];
-
- iregnode[1]=itemp1-1;
- iregnode[2]=itemp1;
- iregnode[3]=itemp2+1;
- iregnode[4]=itemp2;
-
- for (kk=1;kk<=4;kk++)
- {
- if ((iregnode[kk]<1) || (iregnode[kk]>E->trace.numregnodes[j]) )
- {
- fprintf(E->trace.fpt,"ERROR(iquick)-weird regnode %d\n",iregnode[kk]);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- /* find number of choices */
-
- ichoice=0;
- icount=0;
- for (kk=1;kk<=4;kk++)
- {
- if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
-
- icount++;
- for (pp=1;pp<=(kk-1);pp++)
- {
- if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
- }
- ichoice++;
- imap[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
-
-
- next_corner:
- ;
- } /* end kk */
-
- *ich=ichoice;
-
- /* statistical counter */
-
- E->trace.istat_ichoice[j][ichoice]++;
-
- if (ichoice==0) return -99;
-
- /* Here, no check is performed if all 4 corners */
- /* lie within a given element. */
- /* It may be possible (not sure) but unlikely */
- /* that the tracer is still not in that element */
-
- /* Decided to comment this out. */
- /* May not be valid for large regular grids. */
- /*
- */
- /* AKMA */
-
- if ((ichoice==1)&&(icount==4)) return imap[1];
-
- /* check others */
-
- for (kk=1;kk<=ichoice;kk++)
- {
- nel=imap[kk];
- ival=icheck_element_column(E,j,nel,x,y,z,rad);
- if (ival>0) return nel;
- }
-
- /* if still here, no element was found */
-
- return -99;
-}
-
-
/*********** IGET REGEL ******************************************/
/* */
/* This function returns the regular element in which a point */
@@ -4512,60 +2400,7 @@
}
-/****** EXPAND TRACER ARRAYS *****************************************/
-void expand_tracer_arrays(E,j)
- 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;
-}
-
/****************************************************************/
/* DEFINE UV SPACE */
/* */
@@ -4798,63 +2633,7 @@
return;
}
-/*********** CHECK SUM **************************************************/
-/* */
-/* This functions checks to make sure number of tracers is preserved */
-void check_sum(E)
- 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);
- exit(10);
- }
-
- 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 num;
- 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);
-
- num=iallcount;
-
- return num;
-}
-
-
-
/* &&&&&&&&&&&&&&&&&&&& ANALYTICAL TESTS &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&**************/
/**************** ANALYTICAL TEST *********************************************************/
@@ -4866,7 +2645,7 @@
struct All_variables *E;
{
-
+#if 0
int kk,pp;
int nsteps;
int j;
@@ -5119,7 +2898,7 @@
fflush(E->trace.fpt);
fflush(stderr);
-
+#endif
return;
}
@@ -5235,79 +3014,6 @@
-static void init_tracer_flavors(struct All_variables *E)
-{
- int j, kk, number_of_tracers;
- double rad;
-
-
- /* ic_method_for_flavors == 0 (layerd structure) */
- /* any tracer above z_interface is of flavor 0 */
- /* any tracer below z_interface is of flavor 1 */
- if (E->trace.ic_method_for_flavors == 0) {
- 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];
-
- if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
- if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
- }
- }
- }
-
- 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;
-}
-
-
/**************** ANALYTICAL TEST FUNCTION ******************/
/* */
/* vel_s[1] => velocity in theta direction */
@@ -5359,118 +3065,6 @@
}
-/******************* get_neighboring_caps ************************************/
-/* */
-/* Communicate with neighboring processors to get their cap boundaries, */
-/* which is later used by icheck_cap() */
-/* */
-
-static 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*3], rr[12][ncorners*3];
- int nox,noy,noz,dims;
- double x,y,z;
- double theta,phi,rad;
-
- dims=E->mesh.nsd;
- 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<dims; 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++) {
- for (i=1; i<=ncorners; i++) {
- theta = rr[kk][(i-1)*dims];
- phi = rr[kk][(i-1)*dims+1];
- rad = rr[kk][(i-1)*dims+2];
-
- 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++) {
- for (i=1; i<=ncorners; i++) {
- fprintf(E->trace.fpt, "pass=%d corner=%d sx=(%g, %g, %g)\n",
- kk, i,
- E->trace.theta_cap[kk][i],
- E->trace.phi_cap[kk][i],
- E->trace.rad_cap[kk][i]);
- }
- }
- fflush(E->trace.fpt);
-
- }
-
- return;
-}
-
-
/**** PDEBUG ***********************************************************/
static void pdebug(struct All_variables *E, int i)
{
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -401,40 +401,26 @@
E->parallel.me, cycles);
fp1 = output_open(output_file);
- //TODO: unify
+ ncolumns = 3 + E->trace.number_of_extra_quantities;
- if(E->parallel.nprocxy==1) {
- /* regional model */
- fprintf(fp1,"%d %d %.5e\n", cycles, E->Tracer.LOCAL_NUM_TRACERS,
- E->monitor.elapsed_time);
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
+ ncolumns, E->monitor.elapsed_time);
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
- fprintf(fp1,"%.4e %.4e %.4e %.4e\n", E->Tracer.itcolor[n], E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n]);
- }
- }
- else {
- /* global model */
- ncolumns = 3 + E->trace.number_of_extra_quantities;
+ for(n=1;n<=E->trace.ntracers[j];n++) {
+ /* write basic quantities (coordinate) */
+ fprintf(fp1,"%.12e %.12e %.12e",
+ E->trace.basicq[j][0][n],
+ E->trace.basicq[j][1][n],
+ E->trace.basicq[j][2][n]);
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
- ncolumns, E->monitor.elapsed_time);
-
- for(n=1;n<=E->trace.ntracers[j];n++) {
- /* write basic quantities (coordinate) */
- fprintf(fp1,"%.12e %.12e %.12e",
- E->trace.basicq[j][0][n],
- E->trace.basicq[j][1][n],
- E->trace.basicq[j][2][n]);
-
- /* write extra quantities */
- for (i=0; i<E->trace.number_of_extra_quantities; i++) {
- fprintf(fp1," %.12e", E->trace.extraq[j][i][n]);
- }
- fprintf(fp1, "\n");
+ /* 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");
}
+
}
fclose(fp1);
@@ -452,24 +438,18 @@
E->parallel.me, cycles);
fp1 = output_open(output_file);
- //TODO: unify
+ fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
+ cycles, E->lmesh.nno,
+ E->monitor.elapsed_time,
+ E->composition.initial_bulk_composition,
+ E->composition.bulk_composition);
- if(E->parallel.nprocxy==1) {
- }
- else {
- fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
- cycles, E->lmesh.nno,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n", j, E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++) {
+ fprintf(fp1,"%.6e\n",E->composition.comp_node[j][i]);
+ }
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%3d %7d\n", j, E->lmesh.nno);
- for(i=1;i<=E->lmesh.nno;i++) {
- fprintf(fp1,"%.6e\n",E->composition.comp_node[j][i]);
- }
- }
-
}
fclose(fp1);
@@ -487,26 +467,19 @@
E->parallel.me, cycles);
fp1 = output_open(output_file);
- //TODO: unify
+ fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
+ cycles, E->lmesh.nel,
+ E->monitor.elapsed_time,
+ E->composition.initial_bulk_composition,
+ E->composition.bulk_composition);
- if(E->parallel.nprocxy==1) {
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
+ for(i=1;i<=E->lmesh.nel;i++) {
+ fprintf(fp1,"%.6e\n",E->composition.comp_el[j][i]);
+ }
}
- else {
- fprintf(fp1,"%d %d %.5e %.5e %.5e\n",
- cycles, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
- for(i=1;i<=E->lmesh.nel;i++) {
- fprintf(fp1,"%.6e\n",E->composition.comp_el[j][i]);
- }
- }
-
- }
-
fclose(fp1);
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -1536,7 +1536,7 @@
*/
status = set_attribute_int(input, "tracer", E->control.tracer);
- status = set_attribute_string(input, "tracer_file", E->control.tracer_file);
+ status = set_attribute_string(input, "tracer_file", E->trace.tracer_file);
/*
* Visc.inventory
Modified: mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -26,18 +26,8 @@
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
-/*
-Regional_tracer_advection.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 CitcomS
-
-*/
-
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
@@ -52,771 +42,530 @@
#include "element_definitions.h"
#include "global_defs.h"
-static void mean_elem_coord(struct All_variables *E);
+static void make_mesh_ijk(struct All_variables *E);
+
void regional_tracer_setup(struct All_variables *E)
{
- int i,j,k,node;
- int m,ntr;
- int n_x,n_y,n_z;
- int node1,node2,node3,node4,node5,node6,node7,node8;
- int local_element;
- float THETA_LOC_ELEM,FI_LOC_ELEM,R_LOC_ELEM;
- float idummy,xdummy,ydummy,zdummy;
- FILE *fp;
- int nox,noy,noz;
char output_file[255];
- MPI_Comm world;
- MPI_Status status;
- n_x=0;
- n_y=0;
- n_z=0;
+ /* Some error control */
- nox=E->lmesh.nox;
- noy=E->lmesh.noy;
- noz=E->lmesh.noz;
-
- sprintf(output_file,"%s",E->control.tracer_file);
- fp=fopen(output_file,"r");
- if (fp == NULL) {
- fprintf(E->fp,"(Tracer_advection #1) Cannot open %s\n", output_file);
- exit(8);
+ if (E->sphere.caps_per_proc>1) {
+ fprintf(stderr,"This code does not work for multiple caps per processor!\n");
+ parallel_process_termination();
}
- fscanf(fp,"%d",&(E->Tracer.NUM_TRACERS));
- E->Tracer.tracer_x=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.tracer_y=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.tracer_z=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.itcolor=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.x_space=(float*) malloc((nox+1)*sizeof(float));
- E->Tracer.y_space=(float*) malloc((noy+1)*sizeof(float));
- E->Tracer.z_space=(float*) malloc((noz+1)*sizeof(float));
+ /* open tracing output file */
- E->Tracer.LOCAL_ELEMENT=(int*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(int));
+ sprintf(output_file,"%s.tracer_log.%d",E->control.data_file,E->parallel.me);
+ E->trace.fpt=fopen(output_file,"w");
- /* for rheology stuff */
- E->Tracer.THETA_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
- E->Tracer.FI_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
- E->Tracer.R_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
+ /* reset statistical counters */
- E->Tracer.THETA_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.FI_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.R_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->trace.istat_isend=0;
+ E->trace.istat_iempty=0;
+ E->trace.istat_elements_checked=0;
+ E->trace.istat1=0;
+ /* some obscure initial parameters */
+ /* This parameter specifies how close a tracer can get to the boundary */
+ E->trace.box_cushion=0.00001;
- /***comment by Vlad 03/15/2005
- each processor holds its own number of tracers
- ***/
+ /* AKMA turn this back on after debugging */
+ E->trace.itracer_warnings=1;
- ntr=1;
- for(i=1;i<=E->Tracer.NUM_TRACERS;i++) {
- fscanf(fp,"%f %f %f %f", &idummy, &xdummy, &ydummy, &zdummy);
- if(xdummy >= E->sx[1][1][1] && xdummy <= E->sx[1][1][nox*noy*noz]) {
- if(ydummy >= E->sx[1][2][1] && ydummy <= E->sx[1][2][nox*noy*noz]) {
- if(zdummy >= E->sx[1][3][1] && zdummy <= E->sx[1][3][nox*noy*noz]) {
- E->Tracer.itcolor[ntr]=idummy;
- E->Tracer.tracer_x[ntr]=xdummy;
- E->Tracer.tracer_y[ntr]=ydummy;
- E->Tracer.tracer_z[ntr]=zdummy;
- ntr++;
- }
- }
- }
- }
+ /* Determine number of tracer quantities */
- /***comment by Vlad 3/30/2005
- E->Tracer.LOCAL_NUM_TRACERS is the initial number
- of tracers in each processor
- ***/
+ /* advection_quantites - those needed for advection */
+ E->trace.number_of_basic_quantities=12;
- E->Tracer.LOCAL_NUM_TRACERS=ntr-1;
+ /* 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;
- /***comment by Vlad 1/26/2005
- reading the local mesh coordinate
- ***/
+ 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 advection scheme and/or science being done */
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=1;i<=nox;i++)
- {
- j=1;
- k=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.x_space[i]=E->sx[m][1][node];
- }
- for(j=1;j<=noy;j++)
- {
- i=1;
- k=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.y_space[j]=E->sx[m][2][node];
- }
+ /* Some error control regarding size of pointer arrays */
- for(k=1;k<=noz;k++)
- {
- i=1;
- j=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.z_space[k]=E->sx[m][3][node];
- }
+ 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();
+ }
- } /* end of m */
+ /* The bounding box of neiboring processors */
+ get_neighboring_caps(E);
- /***comment by Vlad 04/20/2006
- reading the local element eight nodes coordinate,
- then computing the mean element coordinates
- ***/
+ if (E->trace.ic_method==0) {
+ fprintf(E->trace.fpt,"Not ready for this inputs yet, tracer_ic_method=%d\n", E->trace.ic_method);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ else if (E->trace.ic_method==1)
+ read_tracer_file(E);
+ else if (E->trace.ic_method==2)
+ restart_tracers(E);
+ else {
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ /* total number of tracers */
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(j=1;j<=(noy-1);j++) {
- for(i=1;i<=(nox-1);i++) {
- for(k=1;k<=(noz-1);k++) {
- n_x=i;
- n_y=j;
- n_z=k;
+ E->trace.ilast_tracer_count = isum_tracers(E);
+ fprintf(E->trace.fpt, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
- node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
- node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
- node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
- node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
- node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
- node6 = n_z + n_x*noz + n_y*noz*nox;
- node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
- node8 = n_z+1 + n_x*noz + n_y*noz*nox;
- /* for rheology stuff */
- local_element=node1-(n_x-1)-(n_y-1)*(nox+noz-1);
+ make_mesh_ijk(E);
- E->Tracer.THETA_LOC_ELEM[local_element]=(E->sx[m][1][node1]+E->sx[m][1][node2])/2;
- E->Tracer.FI_LOC_ELEM[local_element]=(E->sx[m][2][node1]+E->sx[m][2][node5])/2;
- E->Tracer.R_LOC_ELEM[local_element]=(E->sx[m][3][node1]+E->sx[m][3][node3])/2;
- }
- }
- }
- } /* end of m */
+ /* find elements */
- mean_elem_coord(E);
+ find_tracers(E);
- return;
-}
+ /* count # of tracers of each flavor */
+ if (E->trace.nflavors > 0)
+ count_tracers_of_flavors(E);
+ composition_setup(E);
+ tracer_post_processing(E);
-void regional_tracer_advection(struct All_variables *E)
-{
- int i,j,k,l,m,n,o,p;
- int n_x,n_y,n_z;
- int node1,node2,node3,node4,node5,node6,node7,node8;
- int nno,nox,noy,noz;
- int iteration;
- float x_tmp, y_tmp, z_tmp;
- float volume, tr_dx, tr_dy, tr_dz, dx, dy, dz;
- float w1,w2,w3,w4,w5,w6,w7,w8;
- float tr_v[NCS][4];
- MPI_Comm world;
- MPI_Status status[4];
- MPI_Status status_count;
- MPI_Request request[4];
- float xmin,xmax,ymin,ymax,zmin,zmax;
+ return;
+}
- float x_global_min,x_global_max,y_global_min,y_global_max,z_global_min,z_global_max;
+static void make_mesh_ijk(struct All_variables *E)
+{
+ int m,i,j,k,node;
+ int nox,noy,noz;
- float *tr_color_1,*tr_x_1,*tr_y_1,*tr_z_1;
- float *tr_color_new[13],*tr_x_new[13],*tr_y_new[13],*tr_z_new[13];
- int *Left_proc_list,*Right_proc_list;
- int *jump_new,*count_new;
- int *jj;
-
- int proc;
- int Previous_num_tracers,Current_num_tracers;
- int locx,locy,locz;
- int left,right,up,down,back,front;
- int temp_tracers;
-
-
- world=E->parallel.world;
-
-
nox=E->lmesh.nox;
noy=E->lmesh.noy;
noz=E->lmesh.noz;
- nno=nox*noy*noz;
- Left_proc_list=(int*) malloc(6*sizeof(int));
- Right_proc_list=(int*) malloc(6*sizeof(int));
- jump_new=(int*) malloc(6*sizeof(int));
- count_new=(int*) malloc(6*sizeof(int));
- jj=(int*) malloc(6*sizeof(int));
+ E->trace.x_space=(double*) malloc(nox*sizeof(double));
+ E->trace.y_space=(double*) malloc(noy*sizeof(double));
+ E->trace.z_space=(double*) malloc(noz*sizeof(double));
- tr_x_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_y_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_z_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_color_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
-
- for(i=0;i<=11;i++){
- tr_color_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_x_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_y_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_z_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- }
-
-
- /*** comment by Vlad 3/25/2005
- This part of code gets the bounding box of the local mesh.
+ /***comment by Vlad 1/26/2005
+ reading the local mesh coordinate
***/
- xmin=E->Tracer.x_space[1];
- xmax=E->Tracer.x_space[nox];
- ymin=E->Tracer.y_space[1];
- ymax=E->Tracer.y_space[noy];
- zmin=E->Tracer.z_space[1];
- zmax=E->Tracer.z_space[noz];
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(i=0;i<nox;i++)
+ E->trace.x_space[i]=E->sx[m][1][i*noz+1];
- /*fprintf(stderr,"%d %d\n", E->parallel.nprocx, E->parallel.loc2proc_map[0][0][1][0]);*/
+ for(j=0;j<noy;j++)
+ E->trace.y_space[j]=E->sx[m][2][j*nox*noz+1];
- /*** comment by Tan2 1/25/2005
- Copy the velocity array.
- ***/
+ for(k=0;k<noz;k++)
+ E->trace.z_space[k]=E->sx[m][3][k+1];
+ } /* end of m */
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=1;i<=nno;i++)
- for(j=1;j<=3;j++) {
- E->GV[m][j][i]=E->sphere.cap[m].V[j][i];
- }
- }
- /*** comment by vlad 03/17/2005
- advecting tracers in each processor
- ***/
+ /* debug *
+ for(i=0;i<nox;i++)
+ fprintf(E->trace.fpt, "i=%d x=%e\n", i, E->trace.x_space[i]);
+ for(j=0;j<noy;j++)
+ fprintf(E->trace.fpt, "j=%d y=%e\n", j, E->trace.y_space[j]);
+ for(k=0;k<noz;k++)
+ fprintf(E->trace.fpt, "k=%d z=%e\n", k, E->trace.z_space[k]);
+ /**
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 0));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 1));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 2));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 3));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 4));
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 0));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 1));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 2));
- n_x=0;
- n_y=0;
- n_z=0;
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 2));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 3));
+ fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 4));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.5));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 1.1));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.55));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 1.0));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.551));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.99));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.7));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.75));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.775));
+ fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.7750001));
+ parallel_process_termination();
+ /**/
- /* mid point method uses 2 iterations */
+ return;
+}
- x_tmp=E->Tracer.tracer_x[n];
- y_tmp=E->Tracer.tracer_y[n];
- z_tmp=E->Tracer.tracer_z[n];
- for(iteration=1;iteration<=2;iteration++)
- {
+/********** IGET ELEMENT *****************************************/
+/* */
+/* This function returns the the real element for a given point. */
+/* Returns -99 in not in this cap. */
+/* iprevious_element, if known, is the last known element. If */
+/* it is not known, input a negative number. */
+int regional_iget_element(struct All_variables *E,
+ int m, int iprevious_element,
+ double x, double y, double z,
+ double theta, double phi, double rad)
+{
+ int e, i, j, k;
+ int ii, jj, kk;
+ int elx, ely, elz;
- /*** comment by Tan2 1/25/2005
- Find the element that contains the tracer.
+ elx = E->lmesh.elx;
+ ely = E->lmesh.ely;
+ elz = E->lmesh.elz;
- nodex n_x n_x+1
- | * |
- <----------->
- tr_dx
+ //TODO: take care of upper bound
- <-------------------------------->
- dx
- ***/
- for(i=1;i<nox;i++) {
- if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
- tr_dx=x_tmp-E->Tracer.x_space[i];
- dx=E->Tracer.x_space[i+1]-E->Tracer.x_space[i];
- n_x=i;
+ /* Search neighboring elements if the previous element is known */
+ if (iprevious_element > 0) {
+ e = iprevious_element - 1;
+ k = e % elz;
+ i = (e / elz) % elx;
+ j = e / (elz*elx);
- E->Tracer.THETA_LOC_ELEM_T[n]=(E->Tracer.x_space[i+1]+E->Tracer.x_space[i])/2;
- }
- }
+ ii = isearch_neighbors(E->trace.x_space, elx+1, theta, i);
+ jj = isearch_neighbors(E->trace.y_space, ely+1, phi, j);
+ kk = isearch_neighbors(E->trace.z_space, elz+1, rad, k);
- for(j=1;j<noy;j++) {
- if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
- tr_dy=y_tmp-E->Tracer.y_space[j];
- dy=E->Tracer.y_space[j+1]-E->Tracer.y_space[j];
- n_y=j;
+ if (ii>=0 && jj>=0 && kk>=0)
+ return jj*elx*elz + ii*elz + kk + 1;
+ }
- E->Tracer.FI_LOC_ELEM_T[n]=(E->Tracer.y_space[j+1]+E->Tracer.y_space[j])/2;
+ /* Search all elements if either the previous element is unknown */
+ /* or failed to find in the neighboring elements */
+ ii = isearch_all(E->trace.x_space, elx+1, theta);
+ jj = isearch_all(E->trace.y_space, ely+1, phi);
+ kk = isearch_all(E->trace.z_space, elz+1, rad);
- }
- }
+ if (ii<0 || jj<0 || kk<0)
+ return -99;
- for(k=1;k<noz;k++) {
- if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
- tr_dz=z_tmp-E->Tracer.z_space[k];
- dz=E->Tracer.z_space[k+1]-E->Tracer.z_space[k];
- n_z=k;
+ return jj*elx*elz + ii*elz + kk + 1;
+}
- E->Tracer.R_LOC_ELEM_T[n]=(E->Tracer.z_space[k+1]+E->Tracer.z_space[k])/2;
- }
- }
+/* array is an ordered array of length nsize */
+/* return an index i, such that array[i] <= a < array[i+1] */
+/* return -1 if not found. */
+/* Note that -1 is returned if a == array[nsize-1] */
+int isearch_all(double *array, int nsize, double a)
+{
+ int high, i, low;
- //fprintf(stderr,"tracer: %d %f %f %f\n",n,E->Tracer.THETA_LOC_ELEM_T[n],E->Tracer.FI_LOC_ELEM_T[n],E->Tracer.R_LOC_ELEM_T[n]);
+ /* check the min/max bound */
+ if ((a < array[0]) || (a >= array[nsize-1]))
+ return -1;
+ /* binary search */
+ for (low=0, high=nsize-1; high-low>1;) {
+ i = (high+low) / 2;
+ if ( a < array[i] ) high = i;
+ else low = i;
+ }
- /*** comment by Tan2 1/25/2005
- Calculate shape functions from tr_dx, tr_dy, tr_dz
- This assumes linear element
- ***/
+ return low;
+}
- /* compute volumetic weighting functions */
- w1=tr_dx*tr_dz*tr_dy;
- w2=(dx-tr_dx)*tr_dz*tr_dy;
- w3=tr_dx*(dz-tr_dz)*tr_dy;
- w4=(dx-tr_dx)*(dz-tr_dz)*tr_dy;
- w5=tr_dx*tr_dz*(dy-tr_dy);
- w6=(dx-tr_dx)*tr_dz*(dy-tr_dy);
- w7=tr_dx*(dz-tr_dz)*(dy-tr_dy);
- w8=(dx-tr_dx)*(dz-tr_dz)*(dy-tr_dy);
+/* Similar the isearch_all(), but with a hint */
+int isearch_neighbors(double *array, int nsize,
+ double a, int hint)
+{
+ /* search the nearest neighbors only */
+ const int number_of_neighbors = 3;
+ int neighbors[5];
+ int n, i;
+ neighbors[0] = hint;
+ neighbors[1] = hint-1;
+ neighbors[2] = hint+1;
+ neighbors[3] = hint-2;
+ neighbors[4] = hint+2;
- volume=dx*dz*dy;
+ /**/
+ for (n=0; n<number_of_neighbors; n++) {
+ i = neighbors[n];
+ if ((i >= 0) && (i < nsize-1) &&
+ (a >= array[i]) && (a < array[i+1]))
+ return i;
+ }
- /*** comment by Tan2 1/25/2005
- Calculate the 8 node numbers of current element
- ***/
+ return -1;
+}
- node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
- node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
- node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
- node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
- node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
- node6 = n_z + n_x*noz + n_y*noz*nox;
- node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
- node8 = n_z+1 + n_x*noz + n_y*noz*nox;
+/* */
+/* This function serves to determine if a point lies within */
+/* a given cap */
+/* */
+int regional_icheck_cap(struct All_variables *E, int icap,
+ double theta, double phi, double rad, double junk)
+{
+ double theta_min, theta_max;
+ double phi_min, phi_max;
- /* printf("%d %d %d %d %d %d %d %d %d\n",E->parallel.me, node1,node2,node3,node4,node5,node6,node7,node8);
- //printf("%d %f %f %f %f %f %f %f %f\n", E->parallel.me, E->GV[1][2][node1], E->GV[1][2][node2], E->GV[1][2][node3], E->GV[1][2][node4], E->GV[1][2][node5], E->GV[1][2][node6], E->GV[1][2][node7], E->GV[1][2][node8]);
- */
+ /* corner 2 is the lower-left corner */
+ /* corner 4 is the upper-right corner */
- /*** comment by Tan2 1/25/2005
- Interpolate the velocity on the tracer position
- ***/
+ theta_min = E->trace.theta_cap[icap][2];
+ theta_max = E->trace.theta_cap[icap][4];
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(j=1;j<=3;j++) {
- tr_v[m][j]=w8*E->GV[m][j][node1]
- +w7*E->GV[m][j][node2]
- +w6*E->GV[m][j][node3]
- +w5*E->GV[m][j][node4]
- +w4*E->GV[m][j][node5]
- +w3*E->GV[m][j][node6]
- +w2*E->GV[m][j][node7]
- +w1*E->GV[m][j][node8];
- tr_v[m][j]=tr_v[m][j]/volume;
+ phi_min = E->trace.phi_cap[icap][2];
+ phi_max = E->trace.phi_cap[icap][4];
- }
+ if ((theta >= theta_min) && (theta < theta_max) &&
+ (phi >= phi_min) && (phi < phi_max))
+ return 1;
+ //TODO: deal with upper right bounds
+ return 0;
+}
- E->Tracer.LOCAL_ELEMENT[n]=node1-(n_x-1)-(n_y-1)*(nox+noz-1);
+static void get_shape_functions(struct All_variables *E,
+ double w[9], int nelem,
+ double theta, double phi, double rad)
+{
+ int e, i, j, k;
+ int elx, ely, elz;
+ double tr_dx, tr_dy, tr_dz;
+ double dx, dy, dz;
+ double volume;
+ elx = E->lmesh.elx;
+ ely = E->lmesh.ely;
+ elz = E->lmesh.elz;
+ e = nelem - 1;
+ k = e % elz;
+ i = (e / elz) % elx;
+ j = e / (elz*elx);
- //fprintf(stderr,"%s %d %s %d %f %f %f %f\n", "The tracer no:", n,"is in element no:", E->Tracer.LOCAL_ELEMENT[n], E->Tracer.y_space[n_y], E->Tracer.tracer_y[n], E->Tracer.FI_LOC_ELEM[515], E->Tracer.y_space[n_y+1]);
+ /*** comment by Tan2 1/25/2005
+ Find the element that contains the tracer.
+ node(i) tracer node(i+1)
+ | * |
+ <----------->
+ tr_dx
+ <-------------------------------->
+ dx
+ ***/
+ tr_dx = theta - E->trace.x_space[i];
+ dx = E->trace.x_space[i+1] - E->trace.x_space[i];
+ tr_dy = phi - E->trace.y_space[j];
+ dy = E->trace.y_space[j+1] - E->trace.y_space[j];
- /*** comment by Tan2 1/25/2005
- advect tracer using mid-point method (2nd order accuracy)
- ***/
+ tr_dz = rad - E->trace.z_space[k];
+ dz = E->trace.z_space[k+1] - E->trace.z_space[k];
- /* mid point method */
- if(iteration == 1) {
- x_tmp = x_tmp + (E->advection.timestep/2.0)*tr_v[m][1]/E->Tracer.z_space[n_z];
- y_tmp = y_tmp + (E->advection.timestep/2.0)*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
- z_tmp = z_tmp + (E->advection.timestep/2.0)*tr_v[m][3];
- }
- if( iteration == 2) {
- E->Tracer.tracer_x[n] += E->advection.timestep*tr_v[m][1]/E->Tracer.z_space[n_z];
- E->Tracer.tracer_y[n] += E->advection.timestep*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
- E->Tracer.tracer_z[n] += E->advection.timestep*tr_v[m][3];
+ /*** comment by Tan2 1/25/2005
+ Calculate shape functions from tr_dx, tr_dy, tr_dz
+ This assumes linear element
+ ***/
- //fprintf(stderr,"%d %d %f %f %f %f %f %f\n", E->parallel.me, E->monitor.solution_cycles, E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n], tr_v[m][1],tr_v[m][2],tr_v[m][3]);
+ /* compute volumetic weighting functions */
+ volume = dx*dz*dy;
- }
+ w[1] = (dx-tr_dx) * (dy-tr_dy) * (dz-tr_dz) / volume;
+ w[2] = tr_dx * (dy-tr_dy) * (dz-tr_dz) / volume;
+ w[3] = tr_dx * tr_dy * (dz-tr_dz) / volume;
+ w[4] = (dx-tr_dx) * tr_dy * (dz-tr_dz) / volume;
+ w[5] = (dx-tr_dx) * (dy-tr_dy) * tr_dz / volume;
+ w[6] = tr_dx * (dy-tr_dy) * tr_dz / volume;
+ w[7] = tr_dx * tr_dy * tr_dz / volume;
+ w[8] = (dx-tr_dx) * tr_dy * tr_dz / volume;
+ /** debug **
+ fprintf(E->trace.fpt, "dr=(%e,%e,%e) tr_dr=(%e,%e,%e)\n",
+ dx, dy, dz, tr_dx, tr_dy, tr_dz);
+ fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e %e %e\n",
+ w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8]);
+ fprintf(E->trace.fpt, "sum(shp): %e\n",
+ w[1]+ w[2]+ w[3]+ w[4]+ w[5]+ w[6]+ w[7]+ w[8]);
+ fflush(E->trace.fpt);
+ /**/
+ return;
+}
- } /* end of m */
+/******** GET VELOCITY ***************************************/
- } /* end of iteration loop */
+void regional_get_velocity(struct All_variables *E,
+ int m, int nelem,
+ double theta, double phi, double rad,
+ double *velocity_vector)
+{
+ void velo_from_element_d();
+ double weight[9], VV[4][9], tmp;
+ int n, d, node;
+ const int sphere_key = 0;
- /*** Comment by Vlad 12/15/2005
- Put the tracers back in the box if they go out
- ***/
+ /* get shape functions at (theta, phi, rad) */
+ get_shape_functions(E, weight, nelem, theta, phi, rad);
- /*** Comment by Vlad 12/15/2005
- get the bounding box of the global mesh
- ***/
- x_global_min = E->control.theta_min;
- x_global_max = E->control.theta_max;
- y_global_min = E->control.fi_min;
- y_global_max = E->control.fi_max;
- z_global_min = E->sphere.ri;
- z_global_max = E->sphere.ro;
+ /* get cartesian velocity */
+ velo_from_element_d(E, VV, m, nelem, sphere_key);
- //printf("%f %f %f %f %f %f\n", E->sphere.cap[1].theta[1],E->sphere.cap[1].theta[3],E->sphere.cap[1].fi[1],E->sphere.cap[1].fi[3],E->sphere.ri,E->sphere.ro);
- if(E->Tracer.tracer_x[n] > x_global_max)
- E->Tracer.tracer_x[n] = x_global_max;
- if(E->Tracer.tracer_x[n] < x_global_min)
- E->Tracer.tracer_x[n] = x_global_min;
- if(E->Tracer.tracer_y[n] > y_global_max)
- E->Tracer.tracer_y[n] = y_global_max;
- if(E->Tracer.tracer_y[n] < y_global_min)
- E->Tracer.tracer_y[n] = y_global_min;
- if(E->Tracer.tracer_z[n] > z_global_max)
- E->Tracer.tracer_z[n] = z_global_max;
- if(E->Tracer.tracer_z[n] < z_global_min)
- E->Tracer.tracer_z[n] = z_global_min;
-
-
-
- }/* end of tracer loop */
-
- /*** Comment by Vlad 3/25/2005
- MPI for the tracer-advection code
+ /*** comment by Tan2 1/25/2005
+ Interpolate the velocity on the tracer position
***/
+ for(d=1; d<=3; d++)
+ velocity_vector[d] = 0;
- m = 0;
- locx = E->parallel.me_loc[1];
- locy = E->parallel.me_loc[2];
- locz = E->parallel.me_loc[3];
-
- /* Am I the left-most proc.? If not, who is on my left? */
- if (locy == 0)
- left = -1;
- else
- left = E->parallel.loc2proc_map[m][locx][locy-1][locz];
-
- /* Am I the right-most proc.? If not, who is on my right? */
- if (locy == E->parallel.nprocy-1)
- right = -1;
- else
- right = E->parallel.loc2proc_map[m][locx][locy+1][locz];
-
- /* Am I the lower-most proc.? If not, who is beneath me? */
- if (locz == 0)
- down = -1;
- else
- down = E->parallel.loc2proc_map[m][locx][locy][locz-1];
-
- /* Am I the upper-most proc.? If not, who is above me? */
- if (locz == E->parallel.nprocz-1)
- up = -1;
- else
- up = E->parallel.loc2proc_map[m][locx][locy][locz+1];
-
- /* Am I the back-most proc.? If not, who is behind me? */
- if (locx == 0)
- back = -1;
- else
- back = E->parallel.loc2proc_map[m][locx-1][locy][locz];
-
- /* Am I the front-most proc.? If not, who is in front of me? */
- if (locx == E->parallel.nprocx-1)
- front = -1;
- else
- front = E->parallel.loc2proc_map[m][locx+1][locy][locz];
-
-
- Left_proc_list[0]=left;
- Left_proc_list[1]=right;
- Left_proc_list[2]=down;
- Left_proc_list[3]=up;
- Left_proc_list[4]=back;
- Left_proc_list[5]=front;
-
- Right_proc_list[0]=right;
- Right_proc_list[1]=left;
- Right_proc_list[2]=up;
- Right_proc_list[3]=down;
- Right_proc_list[4]=front;
- Right_proc_list[5]=back;
-
- jump_new[0]=0;
- jump_new[1]=0;
- jump_new[2]=0;
- jump_new[3]=0;
- jump_new[4]=0;
- jump_new[5]=0;
-
- count_new[0]=0;
- count_new[1]=0;
- count_new[2]=0;
- count_new[3]=0;
- count_new[4]=0;
- count_new[5]=0;
-
- jj[0]=1;
- jj[1]=0;
- jj[2]=3;
- jj[3]=2;
- jj[4]=5;
- jj[5]=4;
-
- temp_tracers=0;
- Current_num_tracers=0;
-
- for(i=0;i<=11;i++){
- for(j=0;j<=E->Tracer.NUM_TRACERS;j++){
- tr_color_new[i][j]=999;
- tr_x_new[i][j]=999;
- tr_y_new[i][j]=999;
- tr_z_new[i][j]=999;
-
- tr_color_1[j]=999;
- tr_x_1[j]=999;
- tr_y_1[j]=999;
- tr_z_1[j]=999;
- }
+ for(d=1; d<=3; d++) {
+ for(n=1; n<=8; n++)
+ velocity_vector[d] += VV[d][n] * weight[n];
}
- i=0;
- j=0;
- k=0;
- l=0;
- m=0;
- o=0;
- p=0;
-
-
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++){
-
- if(E->Tracer.tracer_y[n]>ymax) {
- /* excluding Nan */
- if(E->Tracer.tracer_y[n]+100 != 100) {
- tr_color_new[0][i]=E->Tracer.itcolor[n];
- tr_x_new[0][i]=E->Tracer.tracer_x[n];
- tr_y_new[0][i]=E->Tracer.tracer_y[n];
- tr_z_new[0][i]=E->Tracer.tracer_z[n];
- i++;
- jump_new[0]=i;
- }
- }
- else if(E->Tracer.tracer_y[n]<ymin) {
- if(E->Tracer.tracer_y[n]+100 != 100) {
- tr_color_new[1][j]=E->Tracer.itcolor[n];
- tr_x_new[1][j]=E->Tracer.tracer_x[n];
- tr_y_new[1][j]=E->Tracer.tracer_y[n];
- tr_z_new[1][j]=E->Tracer.tracer_z[n];
- j++;
- jump_new[1]=j;
- }
- }
- else if(E->Tracer.tracer_z[n]>zmax) {
- if(E->Tracer.tracer_z[n]+100 != 100) {
- tr_color_new[2][k]=E->Tracer.itcolor[n];
- tr_x_new[2][k]=E->Tracer.tracer_x[n];
- tr_y_new[2][k]=E->Tracer.tracer_y[n];
- tr_z_new[2][k]=E->Tracer.tracer_z[n];
- k++;
- jump_new[2]=k;
- }
- }
- else if(E->Tracer.tracer_z[n]<zmin) {
- if(E->Tracer.tracer_z[n]+100 != 100) {
- tr_color_new[3][l]=E->Tracer.itcolor[n];
- tr_x_new[3][l]=E->Tracer.tracer_x[n];
- tr_y_new[3][l]=E->Tracer.tracer_y[n];
- tr_z_new[3][l]=E->Tracer.tracer_z[n];
- l++;
- jump_new[3]=l;
- }
- }
-
- else if(E->Tracer.tracer_x[n]>xmax) {
- if(E->Tracer.tracer_x[n]+100 != 100) {
- tr_color_new[4][m]=E->Tracer.itcolor[n];
- tr_x_new[4][m]=E->Tracer.tracer_x[n];
- tr_y_new[4][m]=E->Tracer.tracer_y[n];
- tr_z_new[4][m]=E->Tracer.tracer_z[n];
- m++;
- jump_new[4]=m;
- }
- }
- else if(E->Tracer.tracer_x[n]<xmin) {
- if(E->Tracer.tracer_x[n]+100 != 100) {
- tr_color_new[5][o]=E->Tracer.itcolor[n];
- tr_x_new[5][o]=E->Tracer.tracer_x[n];
- tr_y_new[5][o]=E->Tracer.tracer_y[n];
- tr_z_new[5][o]=E->Tracer.tracer_z[n];
- o++;
- jump_new[5]=o;
- }
- }
-
- else {
- tr_color_1[p]=E->Tracer.itcolor[n];
- tr_x_1[p]=E->Tracer.tracer_x[n];
- tr_y_1[p]=E->Tracer.tracer_y[n];
- tr_z_1[p]=E->Tracer.tracer_z[n];
- p++;
- }
+ /** debug **
+ for(d=1; d<=3; d++) {
+ fprintf(E->trace.fpt, "VV: %e %e %e %e %e %e %e %e: %e\n",
+ VV[d][1], VV[d][2], VV[d][3], VV[d][4],
+ VV[d][5], VV[d][6], VV[d][7], VV[d][8],
+ velocity_vector[d]);
}
- Previous_num_tracers=E->Tracer.LOCAL_NUM_TRACERS;
- Current_num_tracers=Previous_num_tracers-jump_new[0]-jump_new[1]-jump_new[2]-jump_new[3]-jump_new[4]-jump_new[5];
+ tmp = 0;
+ for(n=1; n<=8; n++)
+ tmp += E->sx[m][1][E->ien[m][nelem].node[n]] * weight[n];
- /* compact the remaining tracer */
- for(p=1;p<=Current_num_tracers;p++){
- E->Tracer.itcolor[p]=tr_color_1[p-1];
- E->Tracer.tracer_x[p]=tr_x_1[p-1];
- E->Tracer.tracer_y[p]=tr_y_1[p-1];
- E->Tracer.tracer_z[p]=tr_z_1[p-1];
- }
+ fprintf(E->trace.fpt, "THETA: %e -> %e\n", theta, tmp);
+ fflush(E->trace.fpt);
+ /**/
- for(i=0;i<=5;i++){
+ return;
+}
- j=jj[i];
- if (Left_proc_list[i] >= 0) {
- proc=Left_proc_list[i];
- MPI_Irecv(tr_color_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 11+i, world, &request[0]);
- MPI_Irecv(tr_x_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 12+i, world, &request[1]);
- MPI_Irecv(tr_y_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 13+i, world, &request[2]);
- MPI_Irecv(tr_z_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 14+i, world, &request[3]);
- }
+void regional_put_lost_tracers(struct All_variables *E,
+ int isend[13][13], double *send[13][13])
+{
+ int j, kk, pp;
+ int numtracers, ithatcap, icheck;
+ int isend_position, ipos;
+ int lev = E->mesh.levmax;
+ double theta, phi, rad;
- if (Right_proc_list[i] >= 0) {
- proc=Right_proc_list[i];
- MPI_Send(tr_color_new[i], jump_new[i], MPI_FLOAT, proc, 11+i, world);
- MPI_Send(tr_x_new[i], jump_new[i], MPI_FLOAT, proc, 12+i, world);
- MPI_Send(tr_y_new[i], jump_new[i], MPI_FLOAT, proc, 13+i, world);
- MPI_Send(tr_z_new[i], jump_new[i], MPI_FLOAT, proc, 14+i, world);
- }
- if (Left_proc_list[i] >= 0) {
- MPI_Waitall(4, request, status);
- status_count = status[0];
- MPI_Get_count(&status_count, MPI_FLOAT, &count_new[i]);
- }
+ for (j=1; j<=E->sphere.caps_per_proc; j++) {
+ /* transfer tracers from rlater to send */
- temp_tracers=temp_tracers+count_new[i]-jump_new[i];
- E->Tracer.LOCAL_NUM_TRACERS=Previous_num_tracers+temp_tracers;
+ numtracers = E->trace.ilater[j];
+ for (kk=1; kk<=numtracers; kk++) {
+ theta = E->trace.rlater[j][0][kk];
+ phi = E->trace.rlater[j][1][kk];
+ rad = E->trace.rlater[j][2][kk];
- /* append the tracers */
- if(i <= 1){
- for(n=Current_num_tracers+count_new[j];n<=Current_num_tracers+count_new[i]+count_new[j]-1;n++) {
- m=Current_num_tracers+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
+ /* first check same cap if nprocz>1 */
+ if (E->parallel.nprocz>1) {
+ ithatcap = 0;
+ icheck = regional_icheck_cap(E, ithatcap, theta, phi, rad, rad);
+ if (icheck == 1) goto foundit;
+
}
- }
+ /* check neighboring caps */
- else if (i <= 3) {
- for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];n<=Current_num_tracers+count_new[0]+count_new[1]+count_new[i]+count_new[j]-1;n++) {
- m=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
-
+ for (pp=1; pp<=E->parallel.TNUM_PASS[lev][j]; pp++) {
+ ithatcap = pp;
+ icheck = regional_icheck_cap(E, ithatcap, theta, phi, rad, rad);
+ if (icheck == 1) goto foundit;
}
- }
- else {
- for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];n<=E->Tracer.LOCAL_NUM_TRACERS-1;n++) {
- m=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
+ /* should not be here */
+ if (icheck!=1) {
+ fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
+ fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",
+ theta,phi,rad);
+ icheck = regional_icheck_cap(E, 0, theta, phi, rad,rad);
+ if (icheck == 1) fprintf(E->trace.fpt," icheck here!\n");
+ else fprintf(E->trace.fpt,"icheck not here!\n");
+ fflush(E->trace.fpt);
+ exit(10);
}
- }
+ foundit:
- }
+ isend[j][ithatcap]++;
+ /* assign tracer to send */
- free (tr_color_1);
- free (tr_x_1);
- free (tr_y_1);
- free (tr_z_1);
- for(i=0;i<=11;i++) {
- free (tr_color_new[i]);
- free (tr_x_new[i]);
- free (tr_y_new[i]);
- free (tr_z_new[i]);
- }
+ isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
-
- return;
-}
-
-
-
-/*avoid 1st time step problem with E->Tracer.THETA_LOC_ELEM_T[n],E->Tracer.FI_LOC_ELEM_T[n],E->Tracer.R_LOC_ELEM_T[n] */
-
-static void mean_elem_coord(struct All_variables *E)
-{
- int n, i, j, k;
- float x_tmp, y_tmp, z_tmp;
-
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
- x_tmp=E->Tracer.tracer_x[n];
- y_tmp=E->Tracer.tracer_y[n];
- z_tmp=E->Tracer.tracer_z[n];
-
- for(i=1;i<E->lmesh.nox;i++) {
- if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
- E->Tracer.THETA_LOC_ELEM_T[n]=(E->Tracer.x_space[i+1]+E->Tracer.x_space[i])/2;
+ for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++) {
+ ipos=isend_position+pp;
+ send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
}
- }
- for(j=1;j<E->lmesh.noy;j++) {
- if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
- E->Tracer.FI_LOC_ELEM_T[n]=(E->Tracer.y_space[j+1]+E->Tracer.y_space[j])/2;
+ } /* end kk, assigning tracers */
- }
- }
+ } /* end j */
- for(k=1;k<E->lmesh.noz;k++) {
- if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
- E->Tracer.R_LOC_ELEM_T[n]=(E->Tracer.z_space[k+1]+E->Tracer.z_space[k])/2;
- }
- }
- }
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -33,47 +33,2153 @@
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 CitcomS
-
+ 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.
*/
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
+#include <math.h>
#include "global_defs.h"
#include "parsing.h"
+void tracer_post_processing(struct All_variables *E);
+
+static void predict_tracers(struct All_variables *E);
+static void correct_tracers(struct All_variables *E);
+void count_tracers_of_flavors(struct All_variables *E);
+static void check_sum(struct All_variables *E);
+int isum_tracers(struct All_variables *E);
+static void put_lost_tracers(struct All_variables *E,
+ int isend[13][13], double *send[13][13]);
+int icheck_that_processor_shell(struct All_variables *E,
+ int j, int nprocessor, double rad);
+
+
void tracer_input(struct All_variables *E)
{
- void full_tracer_input();
- int m=E->parallel.me;
+ void full_tracer_input();
+ int m=E->parallel.me;
- input_int("tracer",&(E->control.tracer),"0",m);
- if(E->control.tracer) {
- if(E->parallel.nprocxy == 1)
- input_string("tracer_file",E->control.tracer_file,"tracer.dat",m);
- else
- full_tracer_input(E);
- }
+ input_int("tracer",&(E->control.tracer),"0",m);
+ if(E->control.tracer) {
+
+ /* Initial condition, this option is ignored if E->control.restart is 1,
+ * ie. restarted from a previous run */
+ /* 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) */
+ if(E->control.restart)
+ E->trace.ic_method = 2;
+ else {
+ 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) {
+ }
+ else {
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ fflush(stderr);
+ 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);
+
+
+ input_int("ic_method_for_flavors",&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
+ if (E->trace.ic_method_for_flavors == 0)
+ input_double("z_interface",&(E->trace.z_interface),"0.7",m);
+
+
+ /* Advection Scheme */
+
+ /* itracer_advection_scheme=1 (simple predictor corrector -uses only V(to)) */
+ /* itracer_advection_scheme=2 (predictor-corrector - uses V(to) and V(to+dt)) */
+
+ E->trace.itracer_advection_scheme=2;
+ input_int("tracer_advection_scheme",&(E->trace.itracer_advection_scheme),
+ "2,0,nomax",m);
+
+ if (E->trace.itracer_advection_scheme==1)
+ {}
+ else if (E->trace.itracer_advection_scheme==2)
+ {}
+ else {
+ fprintf(stderr,"Sorry, only advection scheme 1 and 2 available (%d)\n",E->trace.itracer_advection_scheme);
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+
+
+
+ if(E->parallel.nprocxy == 12)
+ full_tracer_input(E);
+
+
+ composition_input(E);
+ }
+
+ return;
}
-void tracer_initial_settings(E)
- struct All_variables *E;
+void tracer_initial_settings(struct All_variables *E)
{
void full_tracer_setup();
- void full_tracer_advection();
+ void full_get_velocity();
+ int full_iget_element();
void regional_tracer_setup();
- void regional_tracer_advection();
+ void regional_get_velocity();
+ int regional_iget_element();
if(E->parallel.nprocxy == 1) {
E->problem_tracer_setup = regional_tracer_setup;
- E->problem_tracer_advection = regional_tracer_advection;
+
+ E->trace.get_velocity = regional_get_velocity;
+ E->trace.iget_element = regional_iget_element;
}
else {
E->problem_tracer_setup = full_tracer_setup;
- E->problem_tracer_advection = full_tracer_advection;
+
+ 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)
+{
+
+ fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
+
+ /* 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);
+ }
+
+ tracer_post_processing(E);
+
+ return;
+}
+
+
+
+/********* TRACER POST PROCESSING ****************************************/
+
+void tracer_post_processing(struct All_variables *E)
+{
+ void get_bulk_composition();
+ char output_file[200];
+
+ double convection_time,tracer_time;
+ double trace_fraction,total_time;
+
+
+ fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
+
+ fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
+
+ fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
+
+ if (E->composition.ichemical_buoyancy==1) {
+ fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
+ E->trace.istat_iempty,(100.0*E->trace.istat_iempty)/E->lmesh.nel);
+ }
+
+
+ /* reset statistical counters */
+
+ E->trace.istat_isend=0;
+ E->trace.istat_iempty=0;
+ E->trace.istat_elements_checked=0;
+ E->trace.istat1=0;
+
+ /* compositional and error fraction data files */
+ //TODO: move
+ if (E->composition.ichemical_buoyancy==1) {
+ get_bulk_composition(E);
+ if (E->parallel.me==0) {
+ fprintf(E->fp,"composition: %e %e\n",E->monitor.elapsed_time,E->composition.bulk_composition);
+ fprintf(E->fp,"composition_error_fraction: %e %e\n",E->monitor.elapsed_time,E->composition.error_fraction);
+
+ }
+
+ }
+
+ 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 */
+/* regardless of advection scheme]. */
+/* 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();
+ void keep_in_sphere();
+ void find_tracers();
+
+
+ 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);
+ keep_in_sphere(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
+
+ /* Current Coordinates are always kept in positions 0-5. */
+
+ 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 */
+/* regardless of advection scheme]. */
+/* 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();
+ void keep_in_sphere();
+ void find_tracers();
+
+
+ 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);
+ keep_in_sphere(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
+
+ /* Fill in Current Positions (other positions are no longer important) */
+
+ 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. */
+
+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 lost_souls();
+ void sphere_to_cart();
+
+ time_stat1=CPU_time0();
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+
+ /* initialize arrays and statistical counters */
+
+ E->trace.ilater[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();
+ */
+
+ lost_souls(E);
+
+ /* Free later arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ if (E->trace.ilater[j]>0) {
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++) {
+ free(E->trace.rlater[j][kk]);
+ }
+ }
+ } /* end j */
+
+
+ /* Adjust Array Sizes */
+
+ reduce_tracer_arrays(E);
+
+ time_stat2=CPU_time0();
+
+ fprintf(E->trace.fpt,"AA: time for find tracers: %f\n", time_stat2-time_stat1);
+
+ return;
+}
+
+/************** LOST SOULS ****************************************************/
+/* */
+/* This function is used to transport tracers to proper processor domains. */
+/* (MPI parallel) */
+/* All of the tracers that were sent to rlater arrays are destined to another */
+/* cap and sent there. Then they are raised up or down for multiple z procs. */
+/* isend[j][n]=number of tracers this processor cap is sending to cap n */
+/* ireceive[j][n]=number of tracers this processor cap is receiving from cap n */
+
+
+void lost_souls(struct All_variables *E)
+{
+ int ithiscap;
+ int ithatcap=1;
+ int isend[13][13];
+ int ireceive[13][13];
+ int isize[13];
+ int kk,pp,j;
+ int mm;
+ int numtracers;
+ int icheck=0;
+ int isend_position;
+ int ipos,ipos2,ipos3;
+ int idb;
+ int idestination_proc=0;
+ int isource_proc;
+ int isend_z[13][3];
+ int ireceive_z[13][3];
+ int isum[13];
+ int irad;
+ int ival;
+ int ithat_processor;
+ int ireceive_position;
+ int ihorizontal_neighbor;
+ int ivertical_neighbor;
+ int ilast_receiver_position;
+ int it;
+ int irec[13];
+ int irec_position;
+ int iel;
+ int num_tracers;
+ int isize_send;
+ int isize_receive;
+ int itemp_size;
+ int itracers_subject_to_vertical_transport[13];
+
+ double x,y,z;
+ double theta,phi,rad;
+ double *send[13][13];
+ double *receive[13][13];
+ double *send_z[13][3];
+ double *receive_z[13][3];
+ double *REC[13];
+
+ void expand_tracer_arrays();
+
+ int number_of_caps=12;
+ int lev=E->mesh.levmax;
+ int num_ngb;
+
+ /* Note, if for some reason, the number of neighbors exceeds */
+ /* 50, which is unlikely, the MPI arrays must be increased. */
+ MPI_Status status[200];
+ MPI_Request request[200];
+ MPI_Status status1;
+ MPI_Status status2;
+ int itag=1;
+
+
+ parallel_process_sync(E);
+ fprintf(E->trace.fpt, "Entering lost_souls()\n");
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ E->trace.istat_isend=E->trace.ilater[j];
+ }
+
+
+ /* debug *
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1; kk<=E->trace.istat_isend; kk++) {
+ fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
+ E->trace.rlater[j][0][kk],
+ E->trace.rlater[j][1][kk],
+ E->trace.rlater[j][2][kk]);
+ }
+ fflush(E->trace.fpt);
+ }
+ /**/
+
+
+ /* initialize isend and ireceive */
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ /* # of neighbors in the horizontal plane */
+ num_ngb = E->parallel.TNUM_PASS[lev][j];
+ isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
+ for (kk=0;kk<=num_ngb;kk++) isend[j][kk]=0;
+ for (kk=0;kk<=num_ngb;kk++) ireceive[j][kk]=0;
+ }
+
+ /* Allocate Maximum Memory to Send Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ itemp_size=max(isize[j],1);
+
+ num_ngb = E->parallel.TNUM_PASS[lev][j];
+ for (kk=0;kk<=num_ngb;kk++) {
+ if ((send[j][kk]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+
+ /* debug *
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ ithiscap=E->sphere.capid[j];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=E->parallel.PROCESSOR[lev][j].pass[kk]+1;
+ fprintf(E->trace.fpt,"cap: %d proc %d TNUM: %d ithatcap: %d\n",
+ ithiscap,E->parallel.me,kk,ithatcap);
+
+ }
+ fflush(E->trace.fpt);
+ }
+ /**/
+
+
+ /* Pre communication */
+ if (E->parallel.nprocxy == 12)
+ full_put_lost_tracers(E, isend, send);
+ else
+ regional_put_lost_tracers(E, isend, send);
+
+ /* Send info to other processors regarding number of send tracers */
+
+ /* idb is the request array index variable */
+ /* Each send and receive has a request variable */
+
+ idb=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ ithiscap=0;
+
+ /* if tracer is in same cap (nprocz>1) */
+
+ if (E->parallel.nprocz>1) {
+ ireceive[j][ithiscap]=isend[j][ithiscap];
+ }
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+
+ /* if neighbor cap is in another processor, send information via MPI */
+
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ idb++;
+ MPI_Isend(&isend[j][ithatcap],1,MPI_INT,idestination_proc,
+ 11,E->parallel.world,
+ &request[idb-1]);
+
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+ /* Receive tracer count info */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+
+ /* if neighbor cap is in another processor, receive information via MPI */
+
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ if (idestination_proc!=E->parallel.me) {
+
+ idb++;
+ MPI_Irecv(&ireceive[j][ithatcap],1,MPI_INT,isource_proc,
+ 11,E->parallel.world,
+ &request[idb-1]);
+
+ } /* end if */
+
+ } /* end kk, number of neighbors */
+ } /* end j, caps per proc */
+
+ /* Wait for non-blocking calls to complete */
+
+ MPI_Waitall(idb,request,status);
+
+ /* Testing, should remove */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+ fprintf(E->trace.fpt,"j: %d send %d to cap %d\n",j,isend[j][kk],isource_proc);
+ fprintf(E->trace.fpt,"j: %d rec %d from cap %d\n",j,ireceive[j][kk],isource_proc);
+ }
+ }
+
+
+ /* Allocate memory in receive arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ num_ngb = E->parallel.TNUM_PASS[lev][j];
+ for (ithatcap=1;ithatcap<=num_ngb;ithatcap++) {
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ itemp_size=max(1,isize[j]);
+
+ if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
+
+ /* Now, send the tracers to proper caps */
+
+ idb=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ ithiscap=0;
+
+ /* same cap */
+
+ if (E->parallel.nprocz>1) {
+
+ ithatcap=ithiscap;
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ for (mm=0;mm<=(isize[j]-1);mm++) {
+ receive[j][ithatcap][mm]=send[j][ithatcap][mm];
+ }
+
+ }
+
+ /* neighbor caps */
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ idb++;
+
+ MPI_Isend(send[j][ithatcap],isize[j],MPI_DOUBLE,idestination_proc,
+ 11,E->parallel.world,
+ &request[idb-1]);
+
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+
+ /* Receive tracers */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ ithiscap=0;
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ idb++;
+
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ MPI_Irecv(receive[j][ithatcap],isize[j],MPI_DOUBLE,isource_proc,
+ 11,E->parallel.world,
+ &request[idb-1]);
+
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+ /* Wait for non-blocking calls to complete */
+
+ MPI_Waitall(idb,request,status);
+
+
+ /* Put all received tracers in array REC[j] */
+ /* This makes things more convenient. */
+
+ /* Sum up size of receive arrays (all tracers sent to this processor) */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ isum[j]=0;
+
+ ithiscap=0;
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+ isum[j]=isum[j]+ireceive[j][ithatcap];
+ }
+ if (E->parallel.nprocz>1) isum[j]=isum[j]+ireceive[j][ithiscap];
+
+ itracers_subject_to_vertical_transport[j]=isum[j];
+ }
+
+ /* Allocate Memory for REC array */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
+ if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ REC[j][0]=0.0;
+ }
+
+ /* Put Received tracers in REC */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ irec[j]=0;
+
+ irec_position=0;
+
+ ithiscap=0;
+
+ /* horizontal neighbors */
+
+ for (ihorizontal_neighbor=1;ihorizontal_neighbor<=E->parallel.TNUM_PASS[lev][j];ihorizontal_neighbor++) {
+
+ ithatcap=ihorizontal_neighbor;
+
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++) {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ ipos2=ipos+mm;
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
+
+ irec_position++;
+
+ } /* end mm (cycling tracer quantities) */
+ } /* end pp (cycling tracers) */
+ } /* end ihorizontal_neighbors (cycling neighbors) */
+
+ /* for tracers in the same cap (nprocz>1) */
+
+ if (E->parallel.nprocz>1) {
+ ithatcap=ithiscap;
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++) {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ ipos2=ipos+mm;
+
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
+
+ irec_position++;
+
+ } /* end mm (cycling tracer quantities) */
+
+ } /* end pp (cycling tracers) */
+
+ } /* endif nproc>1 */
+
+ } /* end j (cycling caps) */
+
+ /* Done filling REC */
+
+
+
+ /* VERTICAL COMMUNICATION */
+
+ /* Note: For generality, I assume that both multiple */
+ /* caps per processor as well as multiple processors */
+ /* in the radial direction. These are probably */
+ /* inconsistent parameter choices for the regular */
+ /* CitcomS code. */
+
+ if (E->parallel.nprocz>1) {
+
+ /* Allocate memory for send_z */
+ /* Make send_z the size of receive array (max size) */
+ /* (No dynamic reallocation of send_z necessary) */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
+ isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
+
+ if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+
+ ithat_processor=E->parallel.PROCESSORz[lev].pass[ivertical_neighbor];
+
+ /* initialize isend_z and ireceive_z array */
+
+ isend_z[j][ivertical_neighbor]=0;
+ ireceive_z[j][ivertical_neighbor]=0;
+
+ /* sort through receive array and check radius */
+
+ it=0;
+ num_tracers=irec[j];
+ for (kk=1;kk<=num_tracers;kk++) {
+
+ it++;
+
+ ireceive_position=((it-1)*E->trace.number_of_tracer_quantities);
+
+ irad=ireceive_position+2;
+
+ rad=REC[j][irad];
+
+ ival=icheck_that_processor_shell(E,j,ithat_processor,rad);
+
+
+ /* if tracer is in other shell, take out of receive array and give to send_z*/
+
+ if (ival==1) {
+
+ isend_z[j][ivertical_neighbor]++;
+
+ isend_position=(isend_z[j][ivertical_neighbor]-1)*E->trace.number_of_tracer_quantities;
+
+ ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ ipos=ireceive_position+mm;
+ ipos2=isend_position+mm;
+
+ send_z[j][ivertical_neighbor][ipos2]=REC[j][ipos];
+
+
+ /* eject tracer info from REC array, and replace with last tracer in array */
+
+ ipos3=ilast_receiver_position+mm;
+ REC[j][ipos]=REC[j][ipos3];
+
+ }
+
+ it--;
+ irec[j]--;
+
+ } /* end if ival===1 */
+
+ /* Otherwise, leave tracer */
+
+ } /* end kk (cycling through tracers) */
+
+ } /* end ivertical_neighbor */
+
+ } /* end j */
+
+
+ /* Send arrays are now filled. */
+ /* Now send send information to vertical processor neighbor */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+
+ MPI_Sendrecv(&isend_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag,
+ &ireceive_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag,E->parallel.world,&status1);
+
+ /* for testing - remove */
+ /*
+ fprintf(E->trace.fpt,"PROC: %d IVN: %d (P: %d) SEND: %d REC: %d\n",
+ E->parallel.me,ivertical_neighbor,E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ isend_z[j][ivertical_neighbor],ireceive_z[j][ivertical_neighbor]);
+ fflush(E->trace.fpt);
+ */
+
+ } /* end ivertical_neighbor */
+
+ } /* end j */
+
+
+ /* Allocate memory to receive_z arrays */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
+ isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(isize[j],1);
+
+ if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end j */
+
+ /* Send Tracers */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+ isize_send=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ isize_receive=ireceive_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+
+ MPI_Sendrecv(send_z[j][ivertical_neighbor],isize_send,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag+1,
+ receive_z[j][ivertical_neighbor],isize_receive,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag+1,E->parallel.world,&status2);
+
+ }
+ }
+
+ /* Put tracers into REC array */
+
+ /* First, reallocate memory to REC */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ isum[j]=0;
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+ isum[j]=isum[j]+ireceive_z[j][ivertical_neighbor];
+ }
+
+ isum[j]=isum[j]+irec[j];
+
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+
+ if (isize[j]>0) {
+ if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (i981)\n");
+ fprintf(E->trace.fpt,"isize: %d\n",isize[j]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+
+ for (kk=1;kk<=ireceive_z[j][ivertical_neighbor];kk++) {
+ irec[j]++;
+
+ irec_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
+ }
+ }
+
+ }
+ }
+
+ /* Free Vertical Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
+ free(send_z[j][ivertical_neighbor]);
+ free(receive_z[j][ivertical_neighbor]);
+ }
+ }
+
+ } /* endif nprocz>1 */
+
+ /* END OF VERTICAL TRANSPORT */
+
+ /* Put away tracers */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ for (kk=1;kk<=irec[j];kk++) {
+ E->trace.ntracers[j]++;
+
+ if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
+
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_basic_quantities-1);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-1);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",x,y,z,theta,phi,rad);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.ielement[j][E->trace.ntracers[j]]=iel;
+
+ }
+ }
+
+ fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
+ fflush(E->trace.fpt);
+ parallel_process_sync(E);
+
+ /* Free Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ ithiscap=0;
+
+ free(REC[j]);
+
+ free(send[j][ithiscap]);
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++) {
+ ithatcap=kk;
+
+ free(send[j][ithatcap]);
+ free(receive[j][ithatcap]);
+
+ }
+
+ }
+ fprintf(E->trace.fpt,"Leaving lost_souls()\n");
+ fflush(E->trace.fpt);
+
+ 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 generate_random_tracers(struct All_variables *E,
+ int tracers_cap, int j)
+{
+ void cart_to_sphere();
+ void keep_in_sphere();
+ void initialize_tracer_arrays();
+
+ int kk;
+ int ival;
+ int number_of_tries=0;
+ int max_tries;
+
+ double x,y,z;
+ double theta,phi,rad;
+ double dmin,dmax;
+ double random1,random2,random3;
+
+
+ initialize_tracer_arrays(E,j,tracers_cap);
+
+
+ /* Tracers are placed randomly in cap */
+ /* (intentionally using rand() instead of srand() )*/
+
+ dmin=-1.0*E->sphere.ro;
+ dmax=E->sphere.ro;
+
+ while (E->trace.ntracers[j]<tracers_cap) {
+
+ number_of_tries++;
+ max_tries=500*tracers_cap;
+
+ if (number_of_tries>max_tries) {
+ fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
+ fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ random1=(1.0*rand())/(1.0*RAND_MAX);
+ random2=(1.0*rand())/(1.0*RAND_MAX);
+ random3=(1.0*rand())/(1.0*RAND_MAX);
+
+ x=dmin+random1*(dmax-dmin);
+ y=dmin+random2*(dmax-dmin);
+ z=dmin+random3*(dmax-dmin);
+
+ /* 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 */
+
+ keep_in_sphere(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. */
+
+void read_tracer_file(E)
+ 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();
+ int isum_tracers();
+ void initialize_tracer_arrays();
+ void keep_in_sphere();
+ void sphere_to_cart();
+ void cart_to_sphere();
+ void expand_tracer_arrays();
+
+ double x,y,z;
+ double theta,phi,rad;
+ double extra[100];
+
+ FILE *fptracer;
+
+ fptracer=fopen(E->trace.tracer_file,"r");
+ fprintf(E->trace.fpt,"Opening %s\n",E->trace.tracer_file);
+
+ fgets(input_s,200,fptracer);
+ sscanf(input_s,"%d %d",&number_of_tracers,&ncolumns);
+ 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++) {
+
+ initialize_tracer_arrays(E,j,iestimate);
+
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ fgets(input_s,200,fptracer);
+ if (E->trace.number_of_extra_quantities==0) {
+ sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
+ }
+ else if (E->trace.number_of_extra_quantities==1) {
+ sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&extra[0]);
+ }
+ /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
+ /* this part has to be changed... */
+ else {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+
+
+ /* make sure theta, phi is in range, and radius is within bounds */
+
+ keep_in_sphere(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]]=extra[i];
+
+ } /* end kk, number of tracers */
+
+ fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
+ E->trace.ntracers[j]);
+
+ } /* end j */
+
+ fclose(fptracer);
+
+ icheck=isum_tracers(E);
+
+ if (icheck!=number_of_tracers) {
+ fprintf(E->trace.fpt,"ERROR(read_tracer_file) - tracers != number in file\n");
+ fprintf(E->trace.fpt,"Tracers in system: %d\n", icheck);
+ fprintf(E->trace.fpt,"Tracers in file: %d\n", number_of_tracers);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ return;
+}
+
+
+/************** RESTART TRACERS ******************************************/
+/* */
+/* This function restarts tracers written from previous calculation */
+/* and the tracers are read as seperate files for each processor domain. */
+
+void restart_tracers(E)
+ struct All_variables *E;
+{
+
+ char output_file[200];
+ char input_s[1000];
+
+ int i,j,kk;
+ int idum1,ncolumns;
+ int numtracers;
+
+ double rdum1;
+ double theta,phi,rad;
+ double extra[100];
+ double x,y,z;
+
+ void initialize_tracer_arrays();
+ void sphere_to_cart();
+ void keep_in_sphere();
+
+ FILE *fp1;
+
+ if (E->trace.number_of_extra_quantities>99) {
+ fprintf(E->trace.fpt,"ERROR(restart_tracers)-increase size of extra[]\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
+
+ if ( (fp1=fopen(output_file,"r"))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ fprintf(stderr,"Restarting Tracers from %s\n",output_file);
+ fflush(stderr);
+
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fgets(input_s,200,fp1);
+ sscanf(input_s,"%d %d %d %lf",
+ &idum1, &numtracers, &ncolumns, &rdum1);
+
+ /* some error control */
+ if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-wrong # of columns\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ /* allocate memory for tracer arrays */
+
+ initialize_tracer_arrays(E,j,numtracers);
+ E->trace.ntracers[j]=numtracers;
+
+ for (kk=1;kk<=numtracers;kk++) {
+ fgets(input_s,200,fp1);
+ if (E->trace.number_of_extra_quantities==0) {
+ sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
+ }
+ else if (E->trace.number_of_extra_quantities==1) {
+ sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&extra[0]);
+ }
+ /* XXX: if E->trace.number_of_extra_quantities is greater than 1 */
+ /* this part has to be changed... */
+ else {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ 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 */
+
+ keep_in_sphere(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]=extra[i];
+
+ }
+
+ fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
+ fflush(E->trace.fpt);
+
+ }
+ fclose(fp1);
+
+
+ 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);
+ 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. */
+
+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;
+ double myatan();
+
+ 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;
+}
+
+
+
+void init_tracer_flavors(struct All_variables *E)
+{
+ int j, kk, number_of_tracers;
+ double rad;
+
+
+ /* ic_method_for_flavors == 0 (layered structure) */
+ /* any tracer above z_interface is of flavor 0 */
+ /* any tracer below z_interface is of flavor 1 */
+ if (E->trace.ic_method_for_flavors == 0) {
+ 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];
+
+ if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
+ if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
+ }
+ }
+ }
+
+ 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*3], rr[12][ncorners*3];
+ int nox,noy,noz,dims;
+ double x,y,z;
+ double theta,phi,rad;
+
+ dims=E->mesh.nsd;
+ 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<dims; 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++) {
+ for (i=1; i<=ncorners; i++) {
+ theta = rr[kk][(i-1)*dims];
+ phi = rr[kk][(i-1)*dims+1];
+ rad = rr[kk][(i-1)*dims+2];
+
+ 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++) {
+ for (i=1; i<=ncorners; i++) {
+ fprintf(E->trace.fpt, "pass=%d corner=%d sx=(%g, %g, %g)\n",
+ kk, 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 initialize_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 *****************************************/
+
+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 */
+
+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.ilater[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 ************************************************/
+
+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/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -702,6 +702,8 @@
static void low_viscosity_channel_factor(struct All_variables *E, float *F)
{
+ //TODO: fix this
+#if 0
int i, n;
float THETA_LOCAL_ELEM_T, THETA_LOCAL_ELEM, FI_LOCAL_ELEM_T, FI_LOCAL_ELEM, R_LOCAL_ELEM_T, R_LOCAL_ELEM;
@@ -732,13 +734,15 @@
}
}
}
-
+#endif
}
static void low_viscosity_wedge_factor(struct All_variables *E, float *F)
{
+ //TODO: fix this
+#if 0
int i, n;
float THETA_LOCAL_ELEM_T, THETA_LOCAL_ELEM, FI_LOCAL_ELEM_T, FI_LOCAL_ELEM, R_LOCAL_ELEM_T, R_LOCAL_ELEM;
@@ -766,4 +770,5 @@
}
}
}
+#endif
}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-17 00:21:01 UTC (rev 6273)
@@ -525,8 +525,6 @@
int coor;
char coor_file[100];
- char tracer_file[100];
-
int lith_age;
int lith_age_time;
int lith_age_old_cycles;
@@ -717,10 +715,6 @@
struct SBC sbc;
struct Output output;
- /* for regional tracer*/
- struct Tracer Tracer;
-
- /* for global tracer*/
struct TRACE trace;
/* for chemical convection & composition rheology */
@@ -810,7 +804,6 @@
void (* problem_update_node_positions)(void*);
void (* problem_initial_fields)(void*);
void (* problem_tracer_setup)(void*);
- void (* problem_tracer_advection)(void*);
void (* problem_tracer_output)(void*, int);
void (* problem_update_bcs)(void*);
void (* special_process_new_velocity)(void*);
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-17 00:21:01 UTC (rev 6273)
@@ -25,29 +25,12 @@
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
-struct Tracer {
- float *tracer_x;
- float *tracer_y;
- float *tracer_z;
- float *itcolor;
- float *x_space, *z_space, *y_space;
- int NUM_TRACERS;
- int LOCAL_NUM_TRACERS;
- int *LOCAL_ELEMENT;
+/* forward declaration */
+struct All_variables;
- float *THETA_LOC_ELEM;
- float *FI_LOC_ELEM;
- float *R_LOC_ELEM;
- float *THETA_LOC_ELEM_T;
- float *FI_LOC_ELEM_T;
- float *R_LOC_ELEM_T;
-
-};
-
-
struct TRACE{
FILE *fpt;
@@ -86,21 +69,6 @@
int ic_method_for_flavors;
double z_interface;
- /* regular mesh parameters */
- int numtheta[13];
- int numphi[13];
- unsigned int numregel[13];
- unsigned int numregnodes[13];
- double deltheta[13];
- double delphi[13];
- double thetamax[13];
- double thetamin[13];
- double phimax[13];
- double phimin[13];
- int *regnodetoel[13];
- int *regtoel[13][5];
-
-
/* statistical parameters */
int istat_ichoice[13][5];
int istat_isend;
@@ -123,10 +91,52 @@
double sin_phi[13][5];
+
+ /*********************/
+ /* for globall model */
+ /*********************/
+
+ /* regular mesh parameters */
+ int numtheta[13];
+ int numphi[13];
+ unsigned int numregel[13];
+ unsigned int numregnodes[13];
+ double deltheta[13];
+ double delphi[13];
+ double thetamax[13];
+ double thetamin[13];
+ double phimax[13];
+ double phimin[13];
+ int *regnodetoel[13];
+ int *regtoel[13][5];
+
/* gnomonic shape functions */
double *UV[13][3];
double cos_theta_f;
double sin_theta_f;
double *shape_coefs[13][3][10];
+
+
+
+ /**********************/
+ /* for regional model */
+ /**********************/
+
+ double *x_space;
+ double *y_space;
+ double *z_space;
+
+
+
+
+ /*********************/
+ /* function pointers */
+ /*********************/
+
+ int (* iget_element)(struct All_variables*, int, int,
+ double, double, double, double, double, double);
+
+ void (* get_velocity)(struct All_variables*, int, int,
+ double, double, double, double*);
};
Modified: mc/3D/CitcomS/trunk/module/misc.c
===================================================================
--- mc/3D/CitcomS/trunk/module/misc.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/module/misc.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -281,7 +281,7 @@
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
if(E->control.tracer==1)
- E->problem_tracer_advection(E);
+ tracer_advection(E);
Py_INCREF(Py_None);
return Py_None;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-16 20:27:35 UTC (rev 6272)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-03-17 00:21:01 UTC (rev 6273)
@@ -572,42 +572,48 @@
getIntProperty(properties, "tracer", E->control.tracer, fp);
- if(E->parallel.nprocxy == 1) {
- // TODO: change to E->tracer.tracer_file..
- getStringProperty(properties, "tracer_file", E->control.tracer_file, fp);
- }
- else if(E->parallel.nprocxy == 12) {
+ if(E->control.restart)
+ E->trace.ic_method = 2;
+ else {
- if(E->control.restart)
- E->trace.ic_method = 2;
+ getIntProperty(properties, "tracer_ic_method",
+ E->trace.ic_method, fp);
+
+ if (E->trace.ic_method==0) {
+ getIntProperty(properties, "tracers_per_element",
+ E->trace.itperel, fp);
+ }
+ else if (E->trace.ic_method==1) {
+ getStringProperty(properties, "tracer_file",
+ E->trace.tracer_file, fp);
+ }
+ else if (E->trace.ic_method==2) {
+ }
else {
+ fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+ }
- getIntProperty(properties, "tracer_ic_method",
- E->trace.ic_method, fp);
+ getIntProperty(properties, "tracer_flavors", E->trace.nflavors, fp);
- if (E->trace.ic_method==0) {
- getIntProperty(properties, "tracers_per_element",
- E->trace.itperel, fp);
- }
- else if (E->trace.ic_method==1) {
- getStringProperty(properties, "tracer_file",
- E->trace.tracer_file, fp);
- }
- else if (E->trace.ic_method==2) {
- }
- else {
- fprintf(stderr,"Sorry, tracer_ic_method only 0, 1 and 2 available\n");
- fflush(stderr);
- parallel_process_termination();
- }
- }
+ getIntProperty(properties, "ic_method_for_flavors", E->trace.ic_method_for_flavors, fp);
+ if (E->trace.ic_method_for_flavors == 0)
+ getDoubleProperty(properties, "z_interface", E->trace.z_interface, fp);
- getIntProperty(properties, "tracer_flavors", E->trace.nflavors, fp);
+ getIntProperty(properties, "chemical_buoyancy",
+ E->composition.ichemical_buoyancy, fp);
- getIntProperty(properties, "ic_method_for_flavors", E->trace.ic_method_for_flavors, fp);
- if (E->trace.ic_method_for_flavors == 0)
- getDoubleProperty(properties, "z_interface", E->trace.z_interface, fp);
+ if (E->composition.ichemical_buoyancy==1) {
+ getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
+ getDoubleProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, fp);
+ getIntProperty(properties, "reset_initial_composition", E->composition.ireset_initial_composition, fp);
+ }
+
+ if(E->parallel.nprocxy == 12) {
+
getIntProperty(properties, "tracer_advection_scheme", E->trace.itracer_advection_scheme, fp);
getIntProperty(properties, "tracer_interpolation_scheme", E->trace.itracer_interpolation_scheme, fp);
@@ -618,15 +624,6 @@
getIntProperty(properties, "analytical_tracer_test", E->trace.ianalytical_tracer_test, fp);
- getIntProperty(properties, "chemical_buoyancy",
- E->composition.ichemical_buoyancy, fp);
-
- if (E->composition.ichemical_buoyancy==1) {
- getIntProperty(properties, "buoy_type", E->composition.ibuoy_type, fp);
- getDoubleProperty(properties, "buoyancy_ratio", E->composition.buoyancy_ratio, fp);
- getIntProperty(properties, "reset_initial_composition", E->composition.ireset_initial_composition, fp);
- }
-
/*
getIntProperty(properties, "compositional_rheology", E->composition.icompositional_rheology, fp);
More information about the cig-commits
mailing list