[cig-commits] r18794 - mc/3D/CitcomS/branches/eheien_dev/lib
emheien at geodynamics.org
emheien at geodynamics.org
Thu Jul 21 18:55:47 PDT 2011
Author: emheien
Date: 2011-07-21 18:55:47 -0700 (Thu, 21 Jul 2011)
New Revision: 18794
Modified:
mc/3D/CitcomS/branches/eheien_dev/lib/Checkpoints.c
mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c
mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Output.c
mc/3D/CitcomS/branches/eheien_dev/lib/Output_gzdir.c
mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
Log:
Changed tracer code to use Tracer object and initial work on Spherical/Cartesian coordinate objects
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Checkpoints.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Checkpoints.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -247,32 +247,32 @@
static void tracer_checkpoint(struct All_variables *E, FILE *fp)
{
- int m, i;
+ int m, i, ntracers, tracer_size;
+ double *tracer_data;
+ TracerList::iterator tr;
+ Tracer temp;
write_sentinel(fp);
- fwrite(&(E->trace.number_of_basic_quantities), sizeof(int), 1, fp);
- fwrite(&(E->trace.number_of_extra_quantities), sizeof(int), 1, fp);
+ tracer_size = temp.size();
+ tracer_data = (double *)malloc(tracer_size*sizeof(double));
+
+ fwrite(&tracer_size, sizeof(int), 1, fp);
fwrite(&(E->trace.nflavors), sizeof(int), 1, fp);
fwrite(&(E->trace.ilast_tracer_count), sizeof(int), 1, fp);
- for(m=1; m<=E->sphere.caps_per_proc; m++)
- fwrite(&(E->trace.ntracers[m]), sizeof(int), 1, fp);
-
- /* the 0-th element of basicq/extraq/ielement is not init'd
- * and won't be used when read it. */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- for(i=0; i<6; i++) {
- fwrite(E->trace.basicq[m][i], sizeof(double),
- E->trace.ntracers[m]+1, fp);
- }
- for(i=0; i<E->trace.number_of_extra_quantities; i++) {
- fwrite(E->trace.extraq[m][i], sizeof(double),
- E->trace.ntracers[m]+1, fp);
- }
- fwrite(E->trace.ielement[m], sizeof(int),
- E->trace.ntracers[m]+1, fp);
+ ntracers = E->trace.tracers[m].size();
+ fwrite(&ntracers, sizeof(int), 1, fp);
+
+ // Write each tracer to the memory data structure, then copy that to the file
+ for (tr=E->trace.tracers[m].begin();tr!=E->trace.tracers[m].end();++tr) {
+ tr->writeToMem(tracer_data);
+ fwrite(tracer_data, tracer_size*sizeof(double), 1, fp);
+ }
}
+
+ free(tracer_data);
return;
}
@@ -280,38 +280,30 @@
static void read_tracer_checkpoint(struct All_variables *E, FILE *fp)
{
- void count_tracers_of_flavors(struct All_variables *E);
- void allocate_tracer_arrays();
+ int m, i, itmp, ntracers, tracer_size;
+ double *tracer_data;
+ Tracer new_tracer;
- int m, i, itmp;
-
read_sentinel(fp, E->parallel.me);
- fread(&itmp, sizeof(int), 1, fp);
- if (itmp != E->trace.number_of_basic_quantities) {
- fprintf(stderr, "Error in reading checkpoint file: tracer basicq, me=%d\n",
+ fread(&tracer_size, sizeof(int), 1, fp);
+ if (tracer_size != new_tracer.size()) {
+ fprintf(stderr, "Error in reading checkpoint file: tracer size, me=%d\n",
E->parallel.me);
- fprintf(stderr, "%d\n", itmp);
+ fprintf(stderr, "%d\n", tracer_size);
exit(-1);
}
+
+ // Allocate memory for reading tracer data in
+ tracer_data = (double *)malloc(tracer_size*sizeof(double));
fread(&itmp, sizeof(int), 1, fp);
- if (itmp != E->trace.number_of_extra_quantities) {
- fprintf(stderr, "Error in reading checkpoint file: tracer extraq, me=%d\n",
- E->parallel.me);
- fprintf(stderr, "%d\n", itmp);
- exit(-1);
-
- }
-
- fread(&itmp, sizeof(int), 1, fp);
if (itmp != E->trace.nflavors) {
fprintf(stderr, "Error in reading checkpoint file: tracer nflavors, me=%d\n",
E->parallel.me);
fprintf(stderr, "%d\n", itmp);
exit(-1);
-
}
fread(&itmp, sizeof(int), 1, fp);
@@ -319,25 +311,19 @@
/* # of tracers, allocate memory */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- fread(&itmp, sizeof(int), 1, fp);
- allocate_tracer_arrays(E, m, itmp);
- E->trace.ntracers[m] = itmp;
- }
+ fread(&ntracers, sizeof(int), 1, fp);
+ allocate_tracer_arrays(E, m, ntracers);
- /* read tracer data */
- for(m=1; m<=E->sphere.caps_per_proc; m++) {
- for(i=0; i<6; i++) {
- fread(E->trace.basicq[m][i], sizeof(double),
- E->trace.ntracers[m]+1, fp);
- }
- for(i=0; i<E->trace.number_of_extra_quantities; i++) {
- fread(E->trace.extraq[m][i], sizeof(double),
- E->trace.ntracers[m]+1, fp);
- }
- fread(E->trace.ielement[m], sizeof(int),
- E->trace.ntracers[m]+1, fp);
+ /* read tracer data */
+ for (i=0;i<ntracers;++i) {
+ fread(tracer_data, tracer_size*sizeof(double), 1, fp);
+ new_tracer.readFromMem(tracer_data);
+ E->trace.tracers[m].push_back(new_tracer);
+ }
}
+ free(tracer_data);
+
/* init E->trace.ntracer_flavor */
count_tracers_of_flavors(E);
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -318,6 +318,7 @@
flavor = i + 1;
E->composition.comp_el[j][i][e] =
E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
+ fprintf(E->trace.fpt, "%d %d %d %f\n", j, i, e, E->composition.comp_el[j][i][e]);
}
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -161,19 +161,19 @@
/* Determine number of tracer quantities */
/* advection_quantites - those needed for advection */
- E->trace.number_of_basic_quantities=12;
+ // ERIC: E->trace.number_of_basic_quantities=12;
/* extra_quantities - used for flavors, composition, etc. */
/* (can be increased for additional science i.e. tracing chemistry */
- E->trace.number_of_extra_quantities = 0;
- if (E->trace.nflavors > 0)
- E->trace.number_of_extra_quantities += 1;
+ // ERIC: E->trace.number_of_extra_quantities = 0;
+ // if (E->trace.nflavors > 0)
+ // E->trace.number_of_extra_quantities += 1;
- E->trace.number_of_tracer_quantities =
- E->trace.number_of_basic_quantities +
- E->trace.number_of_extra_quantities;
+ // ERIC: E->trace.number_of_tracer_quantities =
+ // E->trace.number_of_basic_quantities +
+ // E->trace.number_of_extra_quantities;
/* Fixed positions in tracer array */
@@ -181,25 +181,6 @@
/* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on science being done */
-
- /* Some error control regarding size of pointer arrays */
-
- if (E->trace.number_of_basic_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_tracer_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
write_trace_instructions(E);
@@ -288,7 +269,6 @@
double *receive_z[13][3];
double *REC[13];
- void expand_tracer_arrays();
int icheck_that_processor_shell();
double CPU_time0();
@@ -297,6 +277,8 @@
int number_of_caps=12;
int lev=E->mesh.levmax;
int num_ngb = E->parallel.TNUM_PASS[lev][j];
+
+ Tracer temp_tracer;
/* Note, if for some reason, the number of neighbors exceeds */
/* 50, which is unlikely, the MPI arrays must be increased. */
@@ -312,7 +294,7 @@
fprintf(E->trace.fpt, "Entering lost_souls()\n");
- E->trace.istat_isend=E->trace.ilater[j];
+ E->trace.istat_isend=E->trace.escaped_tracers[j].size();
/** debug **
for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
@@ -327,7 +309,7 @@
/* initialize isend and ireceive */
/* # of neighbors in the horizontal plane */
- isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=E->trace.escaped_tracers[j].size()*temp_tracer.size();
for (kk=0;kk<=num_ngb;kk++) isend[j][kk]=0;
for (kk=0;kk<=num_ngb;kk++) ireceive[j][kk]=0;
@@ -407,7 +389,7 @@
/* Allocate memory in receive arrays */
for (ithatcap=0;ithatcap<=num_ngb;ithatcap++) {
- isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ isize[j]=ireceive[j][ithatcap]*temp_tracer.size();
itemp_size=citmax(1,isize[j]);
@@ -428,7 +410,7 @@
if (E->parallel.nprocz>1) {
ithatcap=ithiscap;
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ isize[j]=isend[j][ithatcap]*temp_tracer.size();
for (mm=0;mm<isize[j];mm++) {
receive[j][ithatcap][mm]=send[j][ithatcap][mm];
}
@@ -440,12 +422,12 @@
for (kk=1;kk<=num_ngb;kk++) {
idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- isize[j]=isend[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=isend[j][kk]*temp_tracer.size();
MPI_Isend(send[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
- isize[j]=ireceive[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=ireceive[j][kk]*temp_tracer.size();
MPI_Irecv(receive[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
@@ -475,7 +457,7 @@
/* Allocate Memory for REC array */
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=isum[j]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
@@ -495,9 +477,9 @@
for (pp=0;pp<ireceive[j][ithatcap];pp++) {
irec[j]++;
- ipos=pp*E->trace.number_of_tracer_quantities;
+ ipos=pp*temp_tracer.size();
- for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
+ for (mm=0;mm<temp_tracer.size();mm++) {
ipos2=ipos+mm;
REC[j][irec_position]=receive[j][ithatcap][ipos2];
@@ -520,7 +502,7 @@
/* (No dynamic reallocation of send_z necessary) */
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=itracers_subject_to_vertical_transport[j]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
@@ -546,7 +528,7 @@
num_tracers=irec[j];
for (kk=1;kk<=num_tracers;kk++) {
- ireceive_position=it*E->trace.number_of_tracer_quantities;
+ ireceive_position=it*temp_tracer.size();
it++;
irad=ireceive_position+2;
@@ -561,12 +543,12 @@
if (ival==1) {
- isend_position=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ isend_position=isend_z[j][ivertical_neighbor]*temp_tracer.size();
isend_z[j][ivertical_neighbor]++;
- ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ilast_receiver_position=(irec[j]-1)*temp_tracer.size();
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ for (mm=0;mm<=(temp_tracer.size()-1);mm++) {
ipos=ireceive_position+mm;
ipos2=isend_position+mm;
@@ -626,7 +608,7 @@
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=ireceive_z[j][kk]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
@@ -643,12 +625,12 @@
idestination_proc = E->parallel.PROCESSORz[lev].pass[kk];
- isize_send=isend_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize_send=isend_z[j][kk]*temp_tracer.size();
MPI_Isend(send_z[j][kk],isize_send,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
- isize_receive=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize_receive=ireceive_z[j][kk]*temp_tracer.size();
MPI_Irecv(receive_z[j][kk],isize_receive,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
@@ -670,7 +652,7 @@
isum[j]=isum[j]+irec[j];
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=isum[j]*temp_tracer.size();
if (isize[j]>0) {
if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL) {
@@ -681,16 +663,15 @@
}
}
-
for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
for (kk=0;kk<ireceive_z[j][ivertical_neighbor];kk++) {
- irec_position=irec[j]*E->trace.number_of_tracer_quantities;
+ irec_position=irec[j]*temp_tracer.size();
irec[j]++;
- ireceive_position=kk*E->trace.number_of_tracer_quantities;
+ ireceive_position=kk*temp_tracer.size();
- for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
+ for (mm=0;mm<temp_tracer.size();mm++) {
REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
}
}
@@ -711,31 +692,19 @@
for (kk=0;kk<irec[j];kk++) {
- E->trace.ntracers[j]++;
+ Tracer new_tracer;
- if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
+ ireceive_position=kk*new_tracer.size();
+
+ new_tracer.readFromMem(&REC[j][ireceive_position]);
- ireceive_position=kk*E->trace.number_of_tracer_quantities;
+ theta=new_tracer.theta();
+ phi=new_tracer.phi();
+ rad=new_tracer.rad();
+ x=new_tracer.x();
+ y=new_tracer.y();
+ z=new_tracer.z();
- for (mm=0;mm<E->trace.number_of_basic_quantities;mm++) {
- ipos=ireceive_position+mm;
-
- E->trace.basicq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
- }
- for (mm=0;mm<E->trace.number_of_extra_quantities;mm++) {
- ipos=ireceive_position+E->trace.number_of_basic_quantities+mm;
-
- E->trace.extraq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
- }
-
- theta=E->trace.basicq[j][0][E->trace.ntracers[j]];
- phi=E->trace.basicq[j][1][E->trace.ntracers[j]];
- rad=E->trace.basicq[j][2][E->trace.ntracers[j]];
- x=E->trace.basicq[j][3][E->trace.ntracers[j]];
- y=E->trace.basicq[j][4][E->trace.ntracers[j]];
- z=E->trace.basicq[j][5][E->trace.ntracers[j]];
-
-
iel=(E->trace.iget_element)(E,j,-99,x,y,z,theta,phi,rad);
if (iel<1) {
@@ -745,8 +714,8 @@
exit(10);
}
- E->trace.ielement[j][E->trace.ntracers[j]]=iel;
-
+ new_tracer.set_ielement(iel);
+ E->trace.tracers[j].push_back(new_tracer);
}
if(E->control.verbose){
fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
@@ -778,22 +747,21 @@
{
const int j = 1;
int kk, pp;
- int numtracers, ithatcap, icheck;
+ int ithatcap, icheck;
int isend_position, ipos;
int lev = E->mesh.levmax;
double theta, phi, rad;
double x, y, z;
+ TracerList::iterator tr;
/* transfer tracers from rlater to send */
- numtracers=E->trace.ilater[j];
+ for (tr=E->trace.escaped_tracers[j].begin();tr!=E->trace.escaped_tracers[j].end();++tr) {
+ rad=tr->rad();
+ x=tr->x();
+ y=tr->y();
+ z=tr->z();
- 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) {
@@ -829,15 +797,13 @@
/* assign tracer to send */
- isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
+ isend_position=(isend[j][ithatcap]-1)*tr->size();
+ tr->writeToMem(&send[j][ithatcap][isend_position]);
- for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++) {
- ipos=isend_position+pp;
- send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
- }
-
} /* end kk, assigning tracers */
+ E->trace.escaped_tracers[j].clear();
+
return;
}
@@ -1971,11 +1937,11 @@
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
- E->trace.number_of_basic_quantities);
+ 12); // ERIC: E->trace.number_of_basic_quantities);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
- E->trace.number_of_extra_quantities);
+ 1); // ERIC: E->trace.number_of_extra_quantities);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
- E->trace.number_of_tracer_quantities);
+ 13); // ERIC: E->trace.number_of_tracer_quantities);
/* analytical test */
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Output.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Output.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -667,38 +667,34 @@
void output_tracer(struct All_variables *E, int cycles)
{
- int i, j, n, ncolumns;
- char output_file[255];
- FILE *fp1;
-
- sprintf(output_file,"%s.tracer.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
- fp1 = output_open(output_file, "w");
-
- ncolumns = 3 + E->trace.number_of_extra_quantities;
-
- 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");
- }
-
- }
-
- fclose(fp1);
- return;
+ int i, j, n, ncolumns;
+ char output_file[255];
+ FILE *fp1;
+ TracerList::iterator tr;
+
+ sprintf(output_file,"%s.tracer.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file, "w");
+
+ ncolumns = 3 + 1; // ERIC: E->trace.number_of_extra_quantities;
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.tracers[j].size(),
+ ncolumns, E->monitor.elapsed_time);
+
+ for(tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ /* write basic quantities (coordinate) */
+ fprintf(fp1,"%.12e %.12e %.12e", tr->theta(), tr->phi(), tr->rad());
+
+ /* write extra quantities */
+ fprintf(fp1," %.12e", tr->flavor());
+ fprintf(fp1, "\n");
+ }
+
+ }
+
+ fclose(fp1);
+ return;
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Output_gzdir.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Output_gzdir.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -1075,39 +1075,35 @@
void gzdir_output_tracer(struct All_variables *E, int cycles)
{
- int i, j, n, ncolumns;
- char output_file[255];
- gzFile *fp1;
-
- snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz",
- E->control.data_dir,cycles,
- E->parallel.me, cycles);
- fp1 = gzdir_output_open(output_file,"w");
-
- ncolumns = 3 + E->trace.number_of_extra_quantities;
-
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
- ncolumns, E->monitor.elapsed_time);
-
- for(n=1;n<=E->trace.ntracers[j];n++) {
- /* write basic quantities (coordinate) */
- gzprintf(fp1,"%9.5e %9.5e %9.5e",
- E->trace.basicq[j][0][n],
- E->trace.basicq[j][1][n],
- E->trace.basicq[j][2][n]);
-
- /* write extra quantities */
- for (i=0; i<E->trace.number_of_extra_quantities; i++) {
- gzprintf(fp1," %9.5e", E->trace.extraq[j][i][n]);
- }
- gzprintf(fp1, "\n");
- }
-
- }
-
- gzclose(fp1);
- return;
+ int i, j, n, ncolumns;
+ char output_file[255];
+ gzFile *fp1;
+ TracerList::iterator tr;
+
+ snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz",
+ E->control.data_dir,cycles,
+ E->parallel.me, cycles);
+ fp1 = gzdir_output_open(output_file,"w");
+
+ ncolumns = 3 + 1;
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ gzprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.tracers[j].size(),
+ ncolumns, E->monitor.elapsed_time);
+
+ for(tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ /* write basic quantities (coordinate) */
+ gzprintf(fp1,"%9.5e %9.5e %9.5e", tr->x(), tr->y(), tr->z());
+
+ /* write extra quantities */
+ gzprintf(fp1," %9.5e", tr->flavor());
+ gzprintf(fp1, "\n");
+ }
+
+ }
+
+ gzclose(fp1);
+ return;
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -47,9 +47,6 @@
static void write_trace_instructions(struct All_variables *E);
static void make_mesh_ijk(struct All_variables *E);
-static void put_lost_tracers(struct All_variables *E,
- int *send_size, double *send,
- int kk, int j);
static void put_found_tracers(struct All_variables *E,
int recv_size, double *recv,
int j);
@@ -92,48 +89,11 @@
/* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
- /* Determine number of tracer quantities */
-
- /* advection_quantites - those needed for advection */
- E->trace.number_of_basic_quantities=12;
-
- /* extra_quantities - used for flavors, composition, etc. */
- /* (can be increased for additional science i.e. tracing chemistry */
-
- E->trace.number_of_extra_quantities = 0;
- if (E->trace.nflavors > 0)
- E->trace.number_of_extra_quantities += 1;
-
-
- E->trace.number_of_tracer_quantities =
- E->trace.number_of_basic_quantities +
- E->trace.number_of_extra_quantities;
-
-
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on science being done */
-
- /* Some error control regarding size of pointer arrays */
-
- if (E->trace.number_of_basic_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_tracer_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
write_trace_instructions(E);
/* The bounding box of neiboring processors */
@@ -210,11 +170,11 @@
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
- E->trace.number_of_basic_quantities);
+ 12);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
- E->trace.number_of_extra_quantities);
+ 1);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
- E->trace.number_of_tracer_quantities);
+ 13);
@@ -634,18 +594,18 @@
double *send[2];
double *recv[2];
- void expand_tracer_arrays();
- int icheck_that_processor_shell();
-
int ipass;
MPI_Status status[4];
MPI_Request request[4];
+ Tracer temp_tracer;
+ TracerList::iterator tr;
+
double CPU_time0();
double begin_time = CPU_time0();
- E->trace.istat_isend = E->trace.ilater[j];
+ E->trace.istat_isend = E->trace.escaped_tracers[j].size();
/* the bounding box */
for (d=0; d<E->mesh.nsd; d++) {
@@ -689,8 +649,8 @@
/* Allocate Maximum Memory to Send Arrays */
- max_send_size = citmax(2*E->trace.ilater[j], E->trace.ntracers[j]/100);
- itemp_size = max_send_size * E->trace.number_of_tracer_quantities;
+ max_send_size = citmax(2*E->trace.escaped_tracers[j].size(), E->trace.tracers[j].size()/100);
+ itemp_size = max_send_size * temp_tracer.size();
if ((send[0] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
@@ -707,29 +667,42 @@
for (d=0; d<E->mesh.nsd; d++) {
- int original_size = E->trace.ilater[j];
+ int original_size = E->trace.escaped_tracers[j].size();
int idb;
- int kk = 1;
int isend[2], irecv[2];
isend[0] = isend[1] = 0;
/* move out-of-bound tracers to send array */
- while (kk<=E->trace.ilater[j]) {
+ for (tr=E->trace.escaped_tracers[j].begin();tr!=E->trace.escaped_tracers[j].end();) {
double coord;
/* Is the tracer within the bounds in the d-th dimension */
- coord = E->trace.rlater[j][d][kk];
+ switch (d) {
+ case 0:
+ coord = tr->theta();
+ break;
+ case 1:
+ coord = tr->phi();
+ break;
+ case 2:
+ coord = tr->rad();
+ break;
+ }
if (coord < bounds[d][0]) {
- put_lost_tracers(E, &(isend[0]), send[0], kk, j);
+ tr->writeToMem(&send[0][isend[0] * temp_tracer.size()]);
+ isend[0]++;
+ tr = E->trace.escaped_tracers[j].erase(tr);
}
else if (coord >= bounds[d][1]) {
- put_lost_tracers(E, &(isend[1]), send[1], kk, j);
+ tr->writeToMem(&send[1][isend[1] * temp_tracer.size()]);
+ isend[1]++;
+ tr = E->trace.escaped_tracers[j].erase(tr);
}
else {
/* check next tracer */
- kk++;
+ tr++;
}
/* reallocate send if size too small */
@@ -737,7 +710,7 @@
(isend[1] > max_send_size - 5)) {
isize = max_send_size + max_send_size/4 + 10;
- itemp_size = isize * E->trace.number_of_tracer_quantities;
+ itemp_size = isize * temp_tracer.size();
if ((send[0] = (double *)realloc(send[0],
itemp_size*sizeof(double)))
@@ -766,10 +739,10 @@
/* check the total # of tracers is conserved */
- if ((isend[0] + isend[1] + E->trace.ilater[j]) != original_size) {
+ if ((isend[0] + isend[1] + E->trace.escaped_tracers[j].size()) != original_size) {
fprintf(E->trace.fpt, "original_size: %d, rlater_size: %d, "
"send_size: %d\n",
- original_size, E->trace.ilater[j], kk);
+ original_size, E->trace.escaped_tracers[j].size(), kk);
}
@@ -827,7 +800,7 @@
/* Allocate memory in receive arrays */
for (i=0; i<2; i++) {
- isize = irecv[i] * E->trace.number_of_tracer_quantities;
+ isize = irecv[i] * temp_tracer.size();
itemp_size = citmax(1, isize);
if ((recv[i] = (double *)malloc(itemp_size*sizeof(double)))
@@ -846,11 +819,11 @@
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
- isize = isend[i] * E->trace.number_of_tracer_quantities;
+ isize = isend[i] * temp_tracer.size();
MPI_Isend(send[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
- isize = irecv[i] * E->trace.number_of_tracer_quantities;
+ isize = irecv[i] * temp_tracer.size();
MPI_Irecv(recv[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
@@ -889,7 +862,8 @@
/* rlater should be empty by now */
- if (E->trace.ilater[j] > 0) {
+ /* ERIC: restore this later
+ if (E->trace.ilater[j] > 0) {
fprintf(E->trace.fpt, "Error(regional_lost_souls) lost tracers\n");
for (kk=1; kk<=E->trace.ilater[j]; kk++) {
fprintf(E->trace.fpt, "lost #%d xx=(%e, %e, %e)\n", kk,
@@ -899,7 +873,7 @@
}
fflush(E->trace.fpt);
exit(10);
- }
+ }*/
/* Free Arrays */
@@ -912,33 +886,6 @@
}
-static void put_lost_tracers(struct All_variables *E,
- int *send_size, double *send,
- int kk, int j)
-{
- int ilast_tracer, isend_position, ipos;
- int pp;
-
- /* move the tracer from rlater to send */
- isend_position = (*send_size) * E->trace.number_of_tracer_quantities;
-
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
- ipos = isend_position + pp;
- send[ipos] = E->trace.rlater[j][pp][kk];
- }
- (*send_size)++;
-
- /* eject the tracer from rlater */
- ilast_tracer = E->trace.ilater[j];
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
- E->trace.rlater[j][pp][kk] = E->trace.rlater[j][pp][ilast_tracer];
- }
- E->trace.ilater[j]--;
-
- return;
-}
-
-
/****************************************************************/
/* Put the received tracers in basiq & extraq, if within bounds */
/* Otherwise, append to rlater for sending to another proc */
@@ -947,16 +894,13 @@
int recv_size, double *recv,
int j)
{
- void expand_tracer_arrays();
- void expand_later_array();
- int icheck_processor_shell();
-
int kk, pp;
int ipos, ilast, inside, iel;
double theta, phi, rad;
+ Tracer new_tracer;
for (kk=0; kk<recv_size; kk++) {
- ipos = kk * E->trace.number_of_tracer_quantities;
+ ipos = kk * new_tracer.size();
theta = recv[ipos];
phi = recv[ipos + 1];
rad = recv[ipos + 2];
@@ -978,20 +922,8 @@
if (inside) {
- E->trace.ntracers[j]++;
- ilast = E->trace.ntracers[j];
+ new_tracer.readFromMem(&recv[ipos]);
- if (E->trace.ntracers[j] > (E->trace.max_ntracers[j]-5))
- expand_tracer_arrays(E, j);
-
- for (pp=0; pp<E->trace.number_of_basic_quantities; pp++)
- E->trace.basicq[j][pp][ilast] = recv[ipos+pp];
-
- ipos += E->trace.number_of_basic_quantities;
- for (pp=0; pp<E->trace.number_of_extra_quantities; pp++)
- E->trace.extraq[j][pp][ilast] = recv[ipos+pp];
-
-
/* found the element */
iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
@@ -1004,31 +936,14 @@
exit(10);
}
- E->trace.ielement[j][ilast] = iel;
-
+ new_tracer.set_ielement(iel);
+ E->trace.tracers[j].push_back(new_tracer);
}
else {
- if (E->trace.ilatersize[j]==0) {
+ new_tracer.readFromMem(&recv[ipos]);
+
+ E->trace.escaped_tracers[j].push_back(new_tracer);
- E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
-
- for (kk=0;kk<E->trace.number_of_tracer_quantities;kk++) {
- if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(put_found_tracers)-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end first particle initiating memory allocation */
-
- E->trace.ilater[j]++;
- ilast = E->trace.ilater[j];
-
- if (E->trace.ilater[j] > (E->trace.ilatersize[j]-5))
- expand_later_array(E, j);
-
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++)
- E->trace.rlater[j][pp][ilast] = recv[ipos+pp];
} /* end of if-else */
/** debug **
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-07-22 01:55:47 UTC (rev 18794)
@@ -56,8 +56,6 @@
int icheck_that_processor_shell(struct All_variables *E,
int j, int nprocessor, double rad);
-void expand_later_array(struct All_variables *E, int j);
-void expand_tracer_arrays(struct All_variables *E, int j);
void tracer_post_processing(struct All_variables *E);
void allocate_tracer_arrays(struct All_variables *E,
int j, int number_of_tracers);
@@ -79,9 +77,6 @@
static void check_sum(struct All_variables *E);
static int isum_tracers(struct All_variables *E);
static void init_tracer_flavors(struct All_variables *E);
-static void reduce_tracer_arrays(struct All_variables *E);
-static void put_away_later(struct All_variables *E, int j, int it);
-static void eject_tracer(struct All_variables *E, int j, int it);
int read_double_vector(FILE *, int , double *);
void cart_to_sphere(struct All_variables *,
double , double , double ,
@@ -94,6 +89,60 @@
+void SphericalCoord::writeToMem(double *mem) const {
+ mem[0] = _theta;
+ mem[1] = _phi;
+ mem[2] = _rad;
+}
+
+void SphericalCoord::readFromMem(const double *mem) {
+ _theta = mem[0];
+ _phi = mem[1];
+ _rad = mem[2];
+}
+
+void CartesianCoord::writeToMem(double *mem) const {
+ mem[0] = _x;
+ mem[1] = _y;
+ mem[2] = _z;
+}
+
+void CartesianCoord::readFromMem(const double *mem) {
+ _x = mem[0];
+ _y = mem[1];
+ _z = mem[2];
+}
+
+
+size_t Tracer::size(void) {
+ return _sc.size() // spherical coords
+ + _cc.size() // Cartesian coords
+ + _cc0.size() // Original coords
+ + _Vc.size() // Velocity
+ + 1 // flavor
+ + 1; // ielement
+}
+
+void Tracer::writeToMem(double *mem) const {
+ _sc.writeToMem(&mem[0]);
+ _cc.writeToMem(&mem[3]);
+ _cc0.writeToMem(&mem[6]);
+ _Vc.writeToMem(&mem[9]);
+ mem[12] = _flavor;
+ mem[13] = _ielement;
+}
+
+void Tracer::readFromMem(const double *mem) {
+ _sc.readFromMem(&mem[0]);
+ _cc.readFromMem(&mem[3]);
+ _cc0.readFromMem(&mem[6]);
+ _Vc.readFromMem(&mem[9]);
+ _flavor = mem[12];
+ _ielement = mem[13];
+}
+
+
+
void tracer_input(struct All_variables *E)
{
void full_tracer_input();
@@ -286,7 +335,6 @@
int i;
/* reset statistical counters */
-
E->trace.istat_isend=0;
E->trace.istat_elements_checked=0;
E->trace.istat1=0;
@@ -353,13 +401,11 @@
/* 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;
@@ -368,27 +414,25 @@
double theta_pred,phi_pred,rad_pred;
double x_pred,y_pred,z_pred;
double velocity_vector[4];
+
+ TracerList::iterator tr;
- void cart_to_sphere();
-
dt=E->advection.timestep;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- numtracers=E->trace.ntracers[j];
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- for (kk=1;kk<=numtracers;kk++) {
+ theta0=tr->theta();
+ phi0=tr->phi();
+ rad0=tr->rad();
+ x0=tr->x();
+ y0=tr->y();
+ z0=tr->z();
- 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];
+ nelem=tr->ielement();
(E->trace.get_velocity)(E,j,nelem,theta0,phi0,rad0,velocity_vector);
x_pred=x0+velocity_vector[1]*dt;
@@ -403,26 +447,14 @@
/* 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;
+ tr->setCoords(SphericalCoord(theta_pred, phi_pred, rad_pred), CartesianCoord(x_pred, y_pred, 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 */
+ tr->setOrigVals(CartesianCoord(x0, y0, z0),
+ CartesianCoord(velocity_vector[1], velocity_vector[2], velocity_vector[3]));
-
} /* end kk, predicting tracers */
} /* end caps */
@@ -452,7 +484,6 @@
{
int j;
- int kk;
int nelem;
@@ -465,32 +496,32 @@
double velocity_vector[4];
double Vx0,Vy0,Vz0;
double Vx_pred,Vy_pred,Vz_pred;
+
+ TracerList::iterator tr;
- void cart_to_sphere();
-
dt=E->advection.timestep;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- for (kk=1;kk<=E->trace.ntracers[j];kk++) {
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
- 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];
+ theta_pred=tr->theta();
+ phi_pred=tr->phi();
+ rad_pred=tr->rad();
+ x_pred=tr->x();
+ y_pred=tr->y();
+ z_pred=tr->z();
- x0=E->trace.basicq[j][6][kk];
- y0=E->trace.basicq[j][7][kk];
- z0=E->trace.basicq[j][8][kk];
+ x0=tr->getOrigCartesianPos()._x;
+ y0=tr->getOrigCartesianPos()._y;
+ z0=tr->getOrigCartesianPos()._z;
- Vx0=E->trace.basicq[j][9][kk];
- Vy0=E->trace.basicq[j][10][kk];
- Vz0=E->trace.basicq[j][11][kk];
+ Vx0=tr->getCartesianVel()._x;
+ Vy0=tr->getCartesianVel()._y;
+ Vz0=tr->getCartesianVel()._z;
- nelem=E->trace.ielement[j][kk];
+ nelem=tr->ielement();
(E->trace.get_velocity)(E,j,nelem,theta_pred,phi_pred,rad_pred,velocity_vector);
@@ -507,12 +538,7 @@
/* 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;
+ tr->setCoords(SphericalCoord(theta_cor, phi_cor, rad_cor), CartesianCoord(x_cor, y_cor, z_cor));
} /* end kk, correcting tracers */
} /* end caps */
@@ -537,22 +563,15 @@
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;
+
+ TracerList::iterator tr;
- void put_away_later();
- void eject_tracer();
- void reduce_tracer_arrays();
- void sphere_to_cart();
- void full_lost_souls();
- void regional_lost_souls();
-
double CPU_time0();
double begin_time = CPU_time0();
@@ -562,45 +581,35 @@
/* initialize arrays and statistical counters */
- E->trace.ilater[j]=E->trace.ilatersize[j]=0;
-
E->trace.istat1=0;
for (kk=0;kk<=4;kk++) {
E->trace.istat_ichoice[j][kk]=0;
}
- //TODO: use while-loop instead of for-loop
- /* important to index by it, not kk */
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();) {
- it=0;
- num_tracers=E->trace.ntracers[j];
+ theta=tr->theta();
+ phi=tr->phi();
+ rad=tr->rad();
+ x=tr->x();
+ y=tr->y();
+ z=tr->z();
- for (kk=1;kk<=num_tracers;kk++) {
+ iprevious_element=tr->ielement();
- 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;
+ tr->set_ielement(iel);
if (iel == -99) {
/* tracer is inside other processors */
- put_away_later(E,j,it);
- eject_tracer(E,j,it);
- it--;
+ E->trace.escaped_tracers[j].push_back(*tr);
+ tr=E->trace.tracers[j].erase(tr);
+ //fprintf(stderr, "ejected!\n" );
} else if (iel == -1) {
/* tracer is inside this processor,
* but cannot find its element.
@@ -608,11 +617,11 @@
if (E->trace.itracer_warnings) exit(10);
+ tr=E->trace.tracers[j].erase(tr);
+ } else {
+ tr++;
+ }
- eject_tracer(E,j,it);
- it--;
- }
-
} /* end tracers */
} /* end j */
@@ -630,21 +639,7 @@
else
regional_lost_souls(E);
- /* Free later arrays */
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
- if (E->trace.ilatersize[j]>0) {
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++) {
- free(E->trace.rlater[j][kk]);
- }
- }
- } /* end j */
-
-
- /* Adjust Array Sizes */
-
- reduce_tracer_arrays(E);
-
E->trace.find_tracers_time += CPU_time0() - begin_time;
return;
@@ -660,7 +655,7 @@
{
int j, flavor, e, kk;
- int numtracers;
+ TracerList::iterator tr;
for (j=1; j<=E->sphere.caps_per_proc; j++) {
@@ -669,12 +664,10 @@
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];
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ e = tr->ielement();
+ flavor = tr->flavor();
E->trace.ntracer_flavor[j][flavor][e]++;
}
}
@@ -702,6 +695,9 @@
void initialize_tracers(struct All_variables *E)
{
+ E->trace.tracers = new TracerList[13];
+ E->trace.escaped_tracers = new TracerList[13];
+
if (E->trace.ic_method==0)
make_tracer_array(E);
else if (E->trace.ic_method==1)
@@ -721,7 +717,6 @@
if(E->parallel.me==0)
fprintf(stderr, "Sum of Tracers: %d\n", E->trace.ilast_tracer_count);
-
/* find elements */
find_tracers(E);
@@ -731,7 +726,7 @@
if (E->trace.nflavors > 0)
count_tracers_of_flavors(E);
-
+
return;
}
@@ -747,9 +742,6 @@
int j;
double processor_fraction;
- void generate_random_tracers();
- void init_tracer_flavors();
-
if (E->parallel.me==0) fprintf(stderr,"Making Tracer Array\n");
for (j=1;j<=E->sphere.caps_per_proc;j++) {
@@ -780,7 +772,6 @@
static void generate_random_tracers(struct All_variables *E,
int tracers_cap, int j)
{
- void cart_to_sphere();
int kk;
int ival;
int number_of_tries=0;
@@ -791,7 +782,6 @@
double xmin,xmax,ymin,ymax,zmin,zmax;
double random1,random2,random3;
-
allocate_tracer_arrays(E,j,tracers_cap);
/* Finding the min/max of the cartesian coordinates. */
@@ -814,7 +804,7 @@
/* Tracers are placed randomly in cap */
/* (intentionally using rand() instead of srand() )*/
- while (E->trace.ntracers[j]<tracers_cap) {
+ while (E->trace.tracers[j].size()<tracers_cap) {
number_of_tries++;
max_tries=100*tracers_cap;
@@ -847,7 +837,6 @@
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);
@@ -860,16 +849,10 @@
(E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
- E->trace.ntracers[j]++;
- kk=E->trace.ntracers[j];
-
- E->trace.basicq[j][0][kk]=theta;
- E->trace.basicq[j][1][kk]=phi;
- E->trace.basicq[j][2][kk]=rad;
- E->trace.basicq[j][3][kk]=x;
- E->trace.basicq[j][4][kk]=y;
- E->trace.basicq[j][5][kk]=z;
-
+ Tracer new_tracer;
+
+ new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x, y, z));
+ E->trace.tracers[j].push_back(new_tracer);
} /* end while */
return;
@@ -894,12 +877,6 @@
int icushion;
int i, j;
-
- int icheck_processor_shell();
- void sphere_to_cart();
- void cart_to_sphere();
- void expand_tracer_arrays();
-
double x,y,z;
double theta,phi,rad;
double buffer[100];
@@ -917,7 +894,7 @@
number_of_tracers, ncolumns);
/* some error control */
- if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ if (1+3 != ncolumns) {
fprintf(E->trace.fpt,"ERROR(read tracer file)-wrong # of columns\n");
fflush(E->trace.fpt);
exit(10);
@@ -939,7 +916,7 @@
for (kk=1;kk<=number_of_tracers;kk++) {
int len, ncol;
- ncol = 3 + E->trace.number_of_extra_quantities;
+ ncol = 3 + 1; // ERIC: E->trace.number_of_extra_quantities;
len = read_double_vector(fptracer, ncol, buffer);
if (len != ncol) {
@@ -974,25 +951,16 @@
/* if still here, tracer is in processor domain */
+ Tracer new_tracer;
+
+ new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x,y,z));
+ new_tracer.set_flavor(buffer[3]);
+ E->trace.tracers[j].push_back(new_tracer);
- E->trace.ntracers[j]++;
-
- if (E->trace.ntracers[j]>=(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
-
- E->trace.basicq[j][0][E->trace.ntracers[j]]=theta;
- E->trace.basicq[j][1][E->trace.ntracers[j]]=phi;
- E->trace.basicq[j][2][E->trace.ntracers[j]]=rad;
- E->trace.basicq[j][3][E->trace.ntracers[j]]=x;
- E->trace.basicq[j][4][E->trace.ntracers[j]]=y;
- E->trace.basicq[j][5][E->trace.ntracers[j]]=z;
-
- for (i=0; i<E->trace.number_of_extra_quantities; i++)
- E->trace.extraq[j][i][E->trace.ntracers[j]]=buffer[i+3];
-
} /* end kk, number of tracers */
fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
- E->trace.ntracers[j]);
+ E->trace.tracers[j].size());
/** debug **
for (kk=1; kk<=E->trace.ntracers[j]; kk++) {
@@ -1050,14 +1018,8 @@
FILE *fp1;
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-increase size of extra[]\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
/* deal with different output formats */
#ifdef USE_GZDIR
if(strcmp(E->output.format, "ascii-gz") == 0){
@@ -1094,7 +1056,8 @@
/* some error control */
- if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ // ERIC: if (E->trace.number_of_extra_quantities+3 != ncolumns) {
+ if (1+3 != ncolumns) {
fprintf(E->trace.fpt,"ERROR(read_old_tracer_file)-wrong # of columns\n");
fflush(E->trace.fpt);
exit(10);
@@ -1103,11 +1066,10 @@
/* allocate memory for tracer arrays */
allocate_tracer_arrays(E,j,numtracers);
- E->trace.ntracers[j]=numtracers;
for (kk=1;kk<=numtracers;kk++) {
int len, ncol;
- ncol = 3 + E->trace.number_of_extra_quantities;
+ ncol = 3 + 1;
len = read_double_vector(fp1, ncol, buffer);
if (len != ncol) {
@@ -1126,16 +1088,12 @@
(E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
- E->trace.basicq[j][0][kk]=theta;
- E->trace.basicq[j][1][kk]=phi;
- E->trace.basicq[j][2][kk]=rad;
- E->trace.basicq[j][3][kk]=x;
- E->trace.basicq[j][4][kk]=y;
- E->trace.basicq[j][5][kk]=z;
+ Tracer new_tracer;
+
+ new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x, y, z));
+ new_tracer.set_flavor(buffer[3]);
+ E->trace.tracers[j].push_back(new_tracer);
- for (i=0; i<E->trace.number_of_extra_quantities; i++)
- E->trace.extraq[j][i][kk]=buffer[i+3];
-
}
/** debug **
@@ -1211,7 +1169,7 @@
imycount = 0;
for (j=1; j<=E->sphere.caps_per_proc; j++)
- imycount = imycount + E->trace.ntracers[j];
+ imycount = imycount + E->trace.tracers[j].size();
MPI_Allreduce(&imycount,&iallcount,1,MPI_INT,MPI_SUM,E->parallel.world);
@@ -1265,53 +1223,51 @@
static void init_tracer_flavors(struct All_variables *E)
{
- int j, kk, number_of_tracers;
- int i;
- double flavor;
- double rad;
-
+ int j, kk, i;
+ double flavor, rad;
+ TracerList::iterator tr;
+
switch(E->trace.ic_method_for_flavors){
- case 0:
- /* ic_method_for_flavors == 0 (layered structure) */
- /* any tracer above z_interface[i] is of flavor i */
- /* any tracer below z_interface is of flavor (nflavors-1) */
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
- number_of_tracers = E->trace.ntracers[j];
- for (kk=1;kk<=number_of_tracers;kk++) {
- rad = E->trace.basicq[j][2][kk];
-
- flavor = E->trace.nflavors - 1;
- for (i=0; i<E->trace.nflavors-1; i++) {
- if (rad > E->trace.z_interface[i]) {
- flavor = i;
- break;
- }
- }
- E->trace.extraq[j][0][kk] = flavor;
- }
- }
- break;
-
- case 1: /* from grd in top n layers */
- case 99: /* (will override restart) */
+ case 0:
+ /* ic_method_for_flavors == 0 (layered structure) */
+ /* any tracer above z_interface[i] is of flavor i */
+ /* any tracer below z_interface is of flavor (nflavors-1) */
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ rad = tr->rad();
+
+ flavor = E->trace.nflavors - 1;
+ for (i=0; i<E->trace.nflavors-1; i++) {
+ if (rad > E->trace.z_interface[i]) {
+ flavor = i;
+ break;
+ }
+ }
+ tr->set_flavor(flavor);
+ }
+ }
+ break;
+
+ case 1: /* from grd in top n layers */
+ case 99: /* (will override restart) */
#ifndef USE_GGRD
- fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
- E->trace.ic_method_for_flavors);
- parallel_process_termination();
+ fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
+ E->trace.ic_method_for_flavors);
+ parallel_process_termination();
#else
- ggrd_init_tracer_flavors(E);
+ ggrd_init_tracer_flavors(E);
#endif
- break;
-
-
- default:
-
- fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
- parallel_process_termination();
- break;
+ break;
+
+
+ default:
+
+ fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
}
-
+
return;
}
@@ -1444,39 +1400,7 @@
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++) {
@@ -1489,8 +1413,6 @@
}
- fprintf(E->trace.fpt,"Physical size of tracer arrays (max_ntracers): %d\n",
- E->trace.max_ntracers[j]);
fflush(E->trace.fpt);
return;
@@ -1498,234 +1420,10 @@
-/****** EXPAND TRACER ARRAYS *****************************************/
-void expand_tracer_arrays(struct All_variables *E, int j)
-{
- int inewsize;
- int kk;
- int icushion;
- /* expand basicq and ielement by 20% */
- icushion=100;
-
- inewsize=E->trace.max_ntracers[j]+E->trace.max_ntracers[j]/5+icushion;
-
- if ((E->trace.ielement[j]=(int *)realloc(E->trace.ielement[j],inewsize*sizeof(int)))==NULL) {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (ielement)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
- for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++) {
- if ((E->trace.basicq[j][kk]=(double *)realloc(E->trace.basicq[j][kk],inewsize*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++) {
- if ((E->trace.extraq[j][kk]=(double *)realloc(E->trace.extraq[j][kk],inewsize*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory 78 (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
-
- fprintf(E->trace.fpt,"Expanding physical memory of ielement, basicq, and extraq to %d from %d\n",
- inewsize,E->trace.max_ntracers[j]);
-
- E->trace.max_ntracers[j]=inewsize;
-
- return;
-}
-
-
-
-
-/****** REDUCE TRACER ARRAYS *****************************************/
-
-static void reduce_tracer_arrays(struct All_variables *E)
-{
-
- int inewsize;
- int kk;
- int iempty_space;
- int j;
-
- int icushion=100;
-
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
-
- /* if physical size is double tracer size, reduce it */
-
- iempty_space=(E->trace.max_ntracers[j]-E->trace.ntracers[j]);
-
- if (iempty_space>(E->trace.ntracers[j]+icushion)) {
-
-
- inewsize=E->trace.ntracers[j]+E->trace.ntracers[j]/4+icushion;
-
- if (inewsize<1) {
- fprintf(E->trace.fpt,"Error(reduce tracer arrays)-something up (hdf3)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
-
- if ((E->trace.ielement[j]=(int *)realloc(E->trace.ielement[j],inewsize*sizeof(int)))==NULL) {
- fprintf(E->trace.fpt,"ERROR(reduce tracer arrays )-no memory (ielement)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
-
-
- for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++) {
- if ((E->trace.basicq[j][kk]=(double *)realloc(E->trace.basicq[j][kk],inewsize*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++) {
- if ((E->trace.extraq[j][kk]=(double *)realloc(E->trace.extraq[j][kk],inewsize*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory 783 (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
-
- fprintf(E->trace.fpt,"Reducing physical memory of ielement, basicq, and extraq to %d from %d\n",
- E->trace.max_ntracers[j],inewsize);
-
- E->trace.max_ntracers[j]=inewsize;
-
- } /* end if */
-
- } /* end j */
-
- return;
-}
-
-
-/********** PUT AWAY LATER ************************************/
-/* */
-/* rlater has a similar structure to basicq */
-/* ilatersize is the physical memory and */
-/* ilater is the number of tracers */
-
-static void put_away_later(struct All_variables *E, int j, int it)
-{
- int kk;
- void expand_later_array();
-
-
- /* The first tracer in initiates memory allocation. */
- /* Memory is freed after parallel communications */
-
- if (E->trace.ilatersize[j]==0) {
-
- E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
-
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++) {
- if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(put_away_later)-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end first particle initiating memory allocation */
-
-
- /* Put tracer in later array */
-
- E->trace.ilater[j]++;
-
- if (E->trace.ilater[j] >= (E->trace.ilatersize[j]-5)) expand_later_array(E,j);
-
- /* stack basic and extra quantities together (basic first) */
-
- for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
- E->trace.rlater[j][kk][E->trace.ilater[j]]=E->trace.basicq[j][kk][it];
-
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- E->trace.rlater[j][E->trace.number_of_basic_quantities+kk][E->trace.ilater[j]]=E->trace.extraq[j][kk][it];
-
-
- return;
-}
-
-
-/****** EXPAND LATER ARRAY *****************************************/
-
-void expand_later_array(struct All_variables *E, int j)
-{
-
- int inewsize;
- int kk;
- int icushion;
-
- /* expand rlater by 20% */
-
- icushion=100;
-
- inewsize=E->trace.ilatersize[j]+E->trace.ilatersize[j]/5+icushion;
-
- for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++) {
- if ((E->trace.rlater[j][kk]=(double *)realloc(E->trace.rlater[j][kk],inewsize*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(expand later array )-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
-
-
- fprintf(E->trace.fpt,"Expanding physical memory of rlater to %d from %d\n",
- inewsize,E->trace.ilatersize[j]);
-
- E->trace.ilatersize[j]=inewsize;
-
- return;
-}
-
-
-/***** EJECT TRACER ************************************************/
-
-static void eject_tracer(struct All_variables *E, int j, int it)
-{
-
- int ilast_tracer;
- int kk;
-
-
- ilast_tracer=E->trace.ntracers[j];
-
- /* put last tracer in ejected tracer position */
-
- E->trace.ielement[j][it]=E->trace.ielement[j][ilast_tracer];
-
- for (kk=0;kk<=((E->trace.number_of_basic_quantities)-1);kk++)
- E->trace.basicq[j][kk][it]=E->trace.basicq[j][kk][ilast_tracer];
-
- for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
- E->trace.extraq[j][kk][it]=E->trace.extraq[j][kk][ilast_tracer];
-
-
-
- E->trace.ntracers[j]--;
-
- return;
-}
-
-
-
/********** ICHECK PROCESSOR SHELL *************/
/* returns -99 if rad is below current shell */
/* returns 0 if rad is above current shell */
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-07-21 22:47:49 UTC (rev 18793)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-07-22 01:55:47 UTC (rev 18794)
@@ -93,8 +93,9 @@
int _ielement;
public:
- Tracer(void) : _sc(), _cc(), _cc0(), _Vc() {};
- Tracer(SphericalCoord new_sc, CartesianCoord new_cc) : _sc(new_sc), _cc(new_cc), _cc0(), _Vc() {};
+ Tracer(void) : _sc(), _cc(), _cc0(), _Vc(), _flavor(0), _ielement(-99) {};
+ Tracer(SphericalCoord new_sc, CartesianCoord new_cc) :
+ _sc(new_sc), _cc(new_cc), _cc0(), _Vc(), _flavor(0), _ielement(-99) {};
CartesianCoord getCartesianPos(void) const { return _cc; };
SphericalCoord getSphericalPos(void) const { return _sc; };
@@ -148,22 +149,27 @@
double box_cushion;
/* tracer arrays */
- int number_of_basic_quantities;
- int number_of_extra_quantities;
- int number_of_tracer_quantities;
+
+ TracerList *tracers;
+
+ TracerList *escaped_tracers;
+
+ //int number_of_basic_quantities;
+ //int number_of_extra_quantities;
+ //int number_of_tracer_quantities;
- double *basicq[13][100];
- double *extraq[13][100];
+ //double *basicq[13][100];
+ //double *extraq[13][100];
- int ntracers[13];
- int max_ntracers[13];
- int *ielement[13];
+ //int ntracers[13];
+ //int max_ntracers[13];
+ //int *ielement[13];
int number_of_tracers;
- int ilatersize[13];
- int ilater[13];
- double *rlater[13][100];
+ //int ilatersize[13];
+ //int ilater[13];
+ //double *rlater[13][100];
/* tracer flavors */
int nflavors;
More information about the CIG-COMMITS
mailing list