[cig-commits] r6304 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Mar 19 16:49:52 PDT 2007
Author: tan2
Date: 2007-03-19 16:49:51 -0700 (Mon, 19 Mar 2007)
New Revision: 6304
Modified:
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/tracer_defs.h
Log:
lost_souls() cannot be shared by both versions. The parallel communications of
both versions are too different. So, the original lost_souls() is renamed to
full_lost_souls(). A new regional_lost_souls() is implemented. This code has
been tested in 3x3x3 processors, running over 100 timesteps. The number of
total tracers is conserved.
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-19 23:49:51 UTC (rev 6304)
@@ -90,7 +90,9 @@
int *ntheta, int *nphi);
static void define_uv_space(struct All_variables *E);
static void determine_shape_coefficients(struct All_variables *E);
-static void pdebug(struct All_variables *E, int i);
+static void full_put_lost_tracers(struct All_variables *E,
+ int isend[13][13], double *send[13][13]);
+void pdebug(struct All_variables *E, int i);
@@ -261,9 +263,667 @@
}
-void full_put_lost_tracers(struct All_variables *E,
- int isend[13][13], double *send[13][13])
+/************** 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 receiving from cap n */
+
+
+void full_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 icheck_that_processor_shell();
+
+ 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 */
+ full_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;
+}
+
+
+static 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;
@@ -2535,19 +3195,14 @@
}
-/*********** KEEP IN SPHERE *********************************************/
+/*********** KEEP WITHIN BOUNDS *****************************************/
/* */
/* This function makes sure the particle is within the sphere, and */
/* phi and theta are within the proper degree range. */
-void keep_in_sphere(E,x,y,z,theta,phi,rad)
- struct All_variables *E;
- double *x;
- double *y;
- double *z;
- double *theta;
- double *phi;
- double *rad;
+void full_keep_within_bounds(struct All_variables *E,
+ double *x, double *y, double *z,
+ double *theta, double *phi, double *rad)
{
fix_theta_phi(theta, phi);
fix_radius(E,rad,theta,phi,x,y,z);
@@ -2988,7 +3643,7 @@
/**** PDEBUG ***********************************************************/
-static void pdebug(struct All_variables *E, int i)
+void pdebug(struct All_variables *E, int i)
{
fprintf(E->trace.fpt,"HERE (Before Sync): %d\n",i);
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-19 23:49:51 UTC (rev 6304)
@@ -382,11 +382,11 @@
for(j=1;j<=vpts;j++)
for(i=1;i<=ends;i++) {
n = E->ien[m][el].node[i];
- volume1 += E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j-1];
- integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j-1];
- }
+ volume1 += E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
+ integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
+ }
- }
+ }
MPI_Allreduce(&volume1 ,&volume ,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
@@ -395,7 +395,6 @@
if(average && volume != 0.0)
integral /= volume;
-
return((double)integral);
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2007-03-19 23:49:51 UTC (rev 6304)
@@ -411,16 +411,17 @@
} /* end for m */
} /* end for level */
-/* ii=0; */
-/* for (m=1;m<=E->sphere.caps_per_proc;m++) */
-/* for (node=1;node<=E->lmesh.nno;node++) */
-/* if(E->node[m][node] & SKIPS) */
-/* ii+=1; */
+ /* count # of global nodes, ignoring overlapping nodes */
+ ii=0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for (node=1;node<=E->lmesh.nno;node++)
+ if(E->node[m][node] & SKIPS)
+ ++ii;
-/* MPI_Allreduce(&ii, &node ,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); */
+ MPI_Allreduce(&ii, &node, 1, MPI_INT, MPI_SUM, E->parallel.world);
-/* E->mesh.nno = E->lmesh.nno*E->parallel.nproc - node - 2*E->mesh.noz; */
-/* E->mesh.neq = E->mesh.nno*3; */
+ E->mesh.nno = E->lmesh.nno*E->parallel.nproc - node - 2*E->mesh.noz;
+ E->mesh.neq = E->mesh.nno*3;
if (E->control.verbose) {
@@ -601,6 +602,26 @@
} /* end for level */
+ if(E->control.verbose) {
+ for(lev=E->mesh.gridmax;lev>=E->mesh.gridmin;lev--) {
+ fprintf(E->fp_out,"output_communication route surface for lev=%d \n",lev);
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ fprintf(E->fp_out," me= %d cap=%d pass %d \n",E->parallel.me,E->sphere.capid[m],E->parallel.TNUM_PASS[lev][m]);
+ for (k=1;k<=E->parallel.TNUM_PASS[lev][m];k++) {
+ fprintf(E->fp_out,"proc %d and pass %d to proc %d with %d eqn and %d node\n",E->parallel.me,k,E->parallel.PROCESSOR[lev][m].pass[k],E->parallel.NUM_NEQ[lev][m].pass[k],E->parallel.NUM_NODE[lev][m].pass[k]);
+/* fprintf(E->fp_out,"Eqn:\n"); */
+/* for (ii=1;ii<=E->parallel.NUM_NEQ[lev][m].pass[k];ii++) */
+/* fprintf(E->fp_out,"%d %d\n",ii,E->parallel.EXCHANGE_ID[lev][m][ii].pass[k]); */
+/* fprintf(E->fp_out,"Node:\n"); */
+/* for (ii=1;ii<=E->parallel.NUM_NODE[lev][m].pass[k];ii++) */
+/* fprintf(E->fp_out,"%d %d\n",ii,E->parallel.EXCHANGE_NODE[lev][m][ii].pass[k]); */
+ }
+ }
+
+ }
+ fflush(E->fp_out);
+ }
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2007-03-19 23:49:51 UTC (rev 6304)
@@ -41,8 +41,17 @@
#include "element_definitions.h"
#include "global_defs.h"
+#include "composition_related.h"
+
+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);
void regional_tracer_setup(struct All_variables *E)
@@ -120,6 +129,8 @@
parallel_process_termination();
}
+ write_trace_instructions(E);
+
/* The bounding box of neiboring processors */
get_neighboring_caps(E);
@@ -147,6 +158,66 @@
}
+/**** WRITE TRACE INSTRUCTIONS ***************/
+static void write_trace_instructions(struct All_variables *E)
+{
+ fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
+ fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
+
+ if (E->trace.ic_method==0) {
+ fprintf(E->trace.fpt,"Generating New Tracer Array\n");
+ fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
+ }
+ if (E->trace.ic_method==1) {
+ fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
+ }
+ if (E->trace.ic_method==2) {
+ fprintf(E->trace.fpt,"Restarting Tracers\n");
+ }
+
+ fprintf(E->trace.fpt,"Number of tracer flavors: %d\n", E->trace.nflavors);
+
+ if (E->trace.nflavors && E->trace.ic_method==0) {
+ fprintf(E->trace.fpt,"Initialized tracer flavors by: %d\n", E->trace.ic_method_for_flavors);
+ if (E->trace.ic_method_for_flavors == 0) {
+ fprintf(E->trace.fpt,"Layered tracer flavors\n");
+ fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
+ }
+ else {
+ fprintf(E->trace.fpt,"Sorry-This IC methods for Flavors are Unavailable %d\n",E->trace.ic_method_for_flavors);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ }
+
+
+
+ /* more obscure stuff */
+
+ fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
+ fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
+ E->trace.number_of_basic_quantities);
+ fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
+ E->trace.number_of_extra_quantities);
+ fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
+ E->trace.number_of_tracer_quantities);
+
+
+
+ if (E->trace.itracer_warnings==0) {
+ fprintf(E->trace.fpt,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fprintf(stderr,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ }
+
+ write_composition_instructions(E);
+
+
+ return;
+}
+
+
static void make_mesh_ijk(struct All_variables *E)
{
int m,i,j,k,node;
@@ -226,7 +297,7 @@
int regional_iget_element(struct All_variables *E,
int m, int iprevious_element,
- double x, double y, double z,
+ double dummy1, double dummy2, double dummy3,
double theta, double phi, double rad)
{
int e, i, j, k;
@@ -237,7 +308,7 @@
ely = E->lmesh.ely;
elz = E->lmesh.elz;
- //TODO: take care of upper bound
+ //TODO: take care of south west bound
/* Search neighboring elements if the previous element is known */
@@ -329,8 +400,8 @@
double theta_min, theta_max;
double phi_min, phi_max;
- /* corner 2 is the lower-left corner */
- /* corner 4 is the upper-right corner */
+ /* corner 2 is the north-west corner */
+ /* corner 4 is the south-east corner */
theta_min = E->trace.theta_cap[icap][2];
theta_max = E->trace.theta_cap[icap][4];
@@ -342,7 +413,7 @@
(phi >= phi_min) && (phi < phi_max))
return 1;
- //TODO: deal with upper right bounds
+ //TODO: deal with south west bounds
return 0;
}
@@ -478,75 +549,467 @@
}
-void regional_put_lost_tracers(struct All_variables *E,
- int isend[13][13], double *send[13][13])
+void regional_keep_within_bounds(struct All_variables *E,
+ double *x, double *y, double *z,
+ double *theta, double *phi, double *rad)
{
- int j, kk, pp;
- int numtracers, ithatcap, icheck;
- int isend_position, ipos;
+ void sphere_to_cart();
+ int changed = 0;
+
+ if (*theta > E->control.theta_max - E->trace.box_cushion) {
+ *theta = E->control.theta_max - E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (*theta < E->control.theta_min + E->trace.box_cushion) {
+ *theta = E->control.theta_min + E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (*phi > E->control.fi_max - E->trace.box_cushion) {
+ *phi = E->control.fi_max - E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (*phi < E->control.fi_min + E->trace.box_cushion) {
+ *phi = E->control.fi_min + E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (*rad > E->sphere.ro - E->trace.box_cushion) {
+ *rad = E->sphere.ro - E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (*rad < E->sphere.ri + E->trace.box_cushion) {
+ *rad = E->sphere.ri + E->trace.box_cushion;
+ changed = 1;
+ }
+
+ if (changed)
+ sphere_to_cart(E, *theta, *phi, *rad, x, y, z);
+
+
+ return;
+}
+
+
+void regional_lost_souls(struct All_variables *E)
+{
+ /* This part only works if E->sphere.caps_per_proc==1 */
+ const int j = 1;
int lev = E->mesh.levmax;
- double theta, phi, rad;
+ int i, d, kk;
+ int max_send_size, isize, itemp_size;
- for (j=1; j<=E->sphere.caps_per_proc; j++) {
+ int ngbr_rank[6+1];
- /* transfer tracers from rlater to send */
+ double bounds[3][2];
+ double *send[2];
+ double *recv[2];
- numtracers = E->trace.ilater[j];
+ void expand_tracer_arrays();
+ int icheck_that_processor_shell();
- 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];
+ int ipass;
+ MPI_Status status[4];
+ MPI_Request request[4];
- /* 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;
+ E->trace.istat_isend = E->trace.ilater[j];
+ /* the bounding box */
+ for (d=0; d<E->mesh.nsd; d++) {
+ bounds[d][0] = E->sx[j][d+1][1];
+ bounds[d][1] = E->sx[j][d+1][E->lmesh.nno];
+ }
+
+ /* set up ranks for neighboring procs */
+ /* if ngbr_rank is -1, there is no neighbor on this side */
+ ipass = 1;
+ for (kk=1; kk<=6; kk++) {
+ if (E->parallel.NUM_PASS[lev][j].bound[kk] == 1) {
+ ngbr_rank[kk] = E->parallel.PROCESSOR[lev][j].pass[ipass];
+ ipass++;
+ }
+ else {
+ ngbr_rank[kk] = -1;
+ }
+ }
+
+ /* debug *
+ 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]);
+ }
+
+ for (d=0; d<E->mesh.nsd; d++) {
+ fprintf(E->trace.fpt, "bounds(dim=%d) = (%e, %e)\n",
+ d, bounds[d][0], bounds[d][1]);
+ }
+
+ for (kk=1; kk<=6; kk++) {
+ fprintf(E->trace.fpt, "pass=%d neighbor_rank=%d\n",
+ kk, ngbr_rank[kk]);
+ }
+ fflush(E->trace.fpt);
+ parallel_process_sync(E);
+ /**/
+
+
+ /* Allocate Maximum Memory to Send Arrays */
+ max_send_size = max(2*E->trace.ilater[j], E->trace.ntracers[j]/100);
+ itemp_size = max_send_size * E->trace.number_of_tracer_quantities;
+
+ if ((send[0] = (double *)malloc(itemp_size*sizeof(double)))
+ == NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (u388)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((send[1] = (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 (d=0; d<E->mesh.nsd; d++) {
+ int original_size = E->trace.ilater[j];
+ int idb;
+ int kk = 1;
+ int isend[2], irecv[2];
+ isend[0] = isend[1] = 0;
+
+
+ /* move out-of-bound tracers to send array */
+ while (kk<=E->trace.ilater[j]) {
+ double coord;
+
+ /* Is the tracer within the bounds in the d-th dimension */
+ coord = E->trace.rlater[j][d][kk];
+
+ if (coord < bounds[d][0]) {
+ put_lost_tracers(E, &(isend[0]), send[0], kk, j);
}
+ else if (coord >= bounds[d][1]) {
+ put_lost_tracers(E, &(isend[1]), send[1], kk, j);
+ }
+ else {
+ /* check next tracer */
+ kk++;
+ }
- /* check neighboring caps */
+ /* reallocate send if size too small */
+ if ((isend[0] > max_send_size - 5) ||
+ (isend[1] > max_send_size - 5)) {
- 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;
+ isize = max_send_size + max_send_size/4 + 10;
+ itemp_size = isize * E->trace.number_of_tracer_quantities;
+
+ if ((send[0] = (double *)realloc(send[0],
+ itemp_size*sizeof(double)))
+ == NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (s4)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((send[1] = (double *)realloc(send[1],
+ itemp_size*sizeof(double)))
+ == NULL) {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (s5)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ fprintf(E->trace.fpt,"Expanding physical memory of send to "
+ "%d from %d\n",
+ isize, max_send_size);
+
+ max_send_size = isize;
}
- /* 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");
+ } /* end of while kk */
+
+
+ /* check the total # of tracers is conserved */
+ if ((isend[0] + isend[1] + E->trace.ilater[j]) != original_size) {
+ fprintf(E->trace.fpt, "original_size: %d, rlater_size: %d, "
+ "send_size: %d\n",
+ original_size, E->trace.ilater[j], kk);
+ }
+
+
+ /** debug **
+ for (i=0; i<2; i++) {
+ for (kk=0; kk<isend[i]; kk++) {
+ fprintf(E->trace.fpt, "dim:%d side:%d kk=%d coord[kk]=%e\n",
+ d, i, kk,
+ send[i][kk*E->trace.number_of_tracer_quantities+d]);
+ }
+ }
+ fflush(E->trace.fpt);
+ /**/
+
+
+ /* Send info to other processors regarding number of send tracers */
+
+ /* check whether there is a neighbor in this pass*/
+ idb = 0;
+ for (i=0; i<2; i++) {
+ int target_rank;
+ kk = d*2 + i + 1;
+ target_rank = ngbr_rank[kk];
+ if (target_rank >= 0) {
+ MPI_Isend(&isend[i], 1, MPI_INT, target_rank,
+ 11, E->parallel.world, &request[idb++]);
+
+ MPI_Irecv(&irecv[i], 1, MPI_INT, target_rank,
+ 11, E->parallel.world, &request[idb++]);
+ }
+ else {
+ irecv[i] = 0;
+ }
+ } /* end of for i */
+
+
+ /* Wait for non-blocking calls to complete */
+ MPI_Waitall(idb, request, status);
+
+
+ /** debug **
+ for (i=0; i<2; i++) {
+ int target_rank;
+ kk = d*2 + i + 1;
+ target_rank = ngbr_rank[kk];
+ if (target_rank >= 0) {
+ fprintf(E->trace.fpt, "%d: %d send %d to proc %d\n",
+ d, i, isend[i], target_rank);
+ fprintf(E->trace.fpt, "%d: %d recv %d from proc %d\n",
+ d, i, irecv[i], target_rank);
+ }
+ }
+ parallel_process_sync(E);
+ /**/
+
+ /* Allocate memory in receive arrays */
+ for (i=0; i<2; i++) {
+ isize = irecv[i] * E->trace.number_of_tracer_quantities;
+ itemp_size = max(1, isize);
+
+ if ((recv[i] = (double *)malloc(itemp_size*sizeof(double)))
+ == NULL) {
+ fprintf(E->trace.fpt, "Error(lost souls)-no memory (c721)\n");
fflush(E->trace.fpt);
exit(10);
}
+ }
- foundit:
- isend[j][ithatcap]++;
+ /* Now, send the tracers to proper procs */
+ idb = 0;
+ for (i=0; i<2; i++) {
+ int target_rank;
+ kk = d*2 + i + 1;
+ target_rank = ngbr_rank[kk];
+ if (target_rank >= 0) {
+ isize = isend[i] * E->trace.number_of_tracer_quantities;
+ MPI_Isend(send[i], isize, MPI_DOUBLE, target_rank,
+ 12, E->parallel.world, &request[idb++]);
- /* assign tracer to send */
+ isize = irecv[i] * E->trace.number_of_tracer_quantities;
+ MPI_Irecv(recv[i], isize, MPI_DOUBLE, target_rank,
+ 12, E->parallel.world, &request[idb++]);
- 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];
+
+ /* Wait for non-blocking calls to complete */
+ MPI_Waitall(idb, request, status);
+
+
+ /** debug **
+ for (i=0; i<2; i++) {
+ for (kk=1; kk<=irecv[i]; kk++) {
+ fprintf(E->trace.fpt, "recv: %d %e %e %e\n",
+ kk,
+ recv[i][(kk-1)*E->trace.number_of_tracer_quantities],
+ recv[i][(kk-1)*E->trace.number_of_tracer_quantities+1],
+ recv[i][(kk-1)*E->trace.number_of_tracer_quantities+2]);
}
+ }
+ fflush(E->trace.fpt);
+ parallel_process_sync(E);
+ /**/
- } /* end kk, assigning tracers */
+ /* put the received tracers */
+ for (i=0; i<2; i++) {
+ put_found_tracers(E, irecv[i], recv[i], j);
+ }
- } /* end j */
+ free(recv[0]);
+ free(recv[1]);
+ } /* end of for d */
+
+
+ /* rlater should be empty by now */
+ if (E->trace.ilater[j] > 0) {
+ fprintf(E->trace.fpt, "Error(regional_lost_souls) lost tracers\n");
+ for (kk=1; kk<=E->trace.ilater[j]; kk++) {
+ fprintf(E->trace.fpt, "lost #%d xx=(%e, %e, %e)\n", kk,
+ E->trace.rlater[j][0][kk],
+ E->trace.rlater[j][1][kk],
+ E->trace.rlater[j][2][kk]);
+ }
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ /* Free Arrays */
+
+ free(send[0]);
+ free(send[1]);
+
return;
}
+
+
+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 */
+
+static void put_found_tracers(struct All_variables *E,
+ 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;
+
+ for (kk=0; kk<recv_size; kk++) {
+ ipos = kk * E->trace.number_of_tracer_quantities;
+ theta = recv[ipos];
+ phi = recv[ipos + 1];
+ rad = recv[ipos + 2];
+
+ /* check whether this tracer is inside this proc */
+ /* check radius first, since it is cheaper */
+ inside = icheck_processor_shell(E, j, rad);
+ if (inside == 1)
+ inside = regional_icheck_cap(E, 0, theta, phi, rad, rad);
+ else
+ inside = 0;
+
+ /** debug **
+ fprintf(E->trace.fpt, "kk=%d, inside=%d, xx=(%e, %e, %e)\n",
+ kk, inside, theta, phi, rad);
+ fprintf(E->trace.fpt, "before: %d %d\n",
+ E->trace.ilater[j], E->trace.ntracers[j]);
+ /**/
+
+ if (inside) {
+
+ E->trace.ntracers[j]++;
+ ilast = E->trace.ntracers[j];
+
+ if (E->trace.ntracers[j] > (E->trace.max_ntracers[j]-5))
+ expand_tracer_arrays(E, j);
+
+ for (pp=0; pp<E->trace.number_of_basic_quantities; pp++)
+ E->trace.basicq[j][pp][ilast] = recv[ipos+pp];
+
+ ipos += E->trace.number_of_basic_quantities;
+ for (pp=0; pp<E->trace.number_of_extra_quantities; pp++)
+ E->trace.extraq[j][pp][ilast] = recv[ipos+pp];
+
+
+ /* found the element */
+ iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
+
+ if (iel<1) {
+ fprintf(E->trace.fpt, "Error(regional lost souls) - "
+ "element not here?\n");
+ fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n",
+ theta, phi, rad);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.ielement[j][ilast] = iel;
+
+ }
+ else {
+ if (E->trace.ilatersize[j]==0) {
+
+ E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
+
+ for (kk=0;kk<E->trace.number_of_tracer_quantities;kk++) {
+ if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
+ fprintf(E->trace.fpt,"AKM(put_found_tracers)-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end first particle initiating memory allocation */
+
+ E->trace.ilater[j]++;
+ ilast = E->trace.ilater[j];
+
+ if (E->trace.ilater[j] > (E->trace.ilatersize[j]-5))
+ expand_later_array(E, j);
+
+ for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++)
+ E->trace.rlater[j][pp][ilast] = recv[ipos+pp];
+ } /* end of if-else */
+
+ /** debug **
+ fprintf(E->trace.fpt, "after: %d %d\n",
+ E->trace.ilater[j], E->trace.ntracers[j]);
+ fflush(E->trace.fpt);
+ /**/
+
+ } /* end of for kk */
+
+ return;
+}
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-03-19 23:49:51 UTC (rev 6304)
@@ -42,14 +42,20 @@
#include <math.h>
#include "global_defs.h"
#include "parsing.h"
+#include "parallel_related.h"
+#include "composition_related.h"
+
+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 count_tracers_of_flavors(struct All_variables *E);
static void predict_tracers(struct All_variables *E);
static void correct_tracers(struct All_variables *E);
-static void lost_souls(struct All_variables *E);
static void make_tracer_array(struct All_variables *E);
static void generate_random_tracers(struct All_variables *E,
int tracers_cap, int j);
@@ -60,13 +66,9 @@
static void init_tracer_flavors(struct All_variables *E);
static void initialize_tracer_arrays(struct All_variables *E,
int j, int number_of_tracers);
-static void expand_tracer_arrays(struct All_variables *E, int j);
static void reduce_tracer_arrays(struct All_variables *E);
static void put_away_later(struct All_variables *E, int j, int it);
-static void expand_later_array(struct All_variables *E, int j);
static void eject_tracer(struct All_variables *E, int j, int it);
-static int icheck_that_processor_shell(struct All_variables *E,
- int j, int nprocessor, double rad);
void tracer_input(struct All_variables *E)
@@ -149,9 +151,11 @@
void tracer_initial_settings(struct All_variables *E)
{
+ void full_keep_within_bounds();
void full_tracer_setup();
void full_get_velocity();
int full_iget_element();
+ void regional_keep_within_bounds();
void regional_tracer_setup();
void regional_get_velocity();
int regional_iget_element();
@@ -159,12 +163,14 @@
if(E->parallel.nprocxy == 1) {
E->problem_tracer_setup = regional_tracer_setup;
+ E->trace.keep_within_bounds = regional_keep_within_bounds;
E->trace.get_velocity = regional_get_velocity;
E->trace.iget_element = regional_iget_element;
}
else {
E->problem_tracer_setup = full_tracer_setup;
+ E->trace.keep_within_bounds = full_keep_within_bounds;
E->trace.get_velocity = full_get_velocity;
E->trace.iget_element = full_iget_element;
}
@@ -210,7 +216,6 @@
void tracer_post_processing(struct All_variables *E)
{
- void get_bulk_composition();
char output_file[200];
double convection_time,tracer_time;
@@ -223,23 +228,25 @@
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) {
+ if (E->composition.on) {
+ fprintf(E->trace.fpt,"Empty elements filled with old compositional "
+ "values: %d (%f percent)\n", E->trace.istat_iempty,
+ (100.0*E->trace.istat_iempty)/E->lmesh.nel);
+ E->trace.istat_iempty=0;
+
+
get_bulk_composition(E);
+
if (E->parallel.me==0) {
fprintf(E->fp,"composition: %e %e\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);
@@ -283,7 +290,6 @@
double velocity_vector[4];
void cart_to_sphere();
- void keep_in_sphere();
void find_tracers();
@@ -314,7 +320,7 @@
/* 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);
+ (E->trace.keep_within_bounds)(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
/* Current Coordinates are always kept in positions 0-5. */
@@ -383,7 +389,6 @@
double Vx_pred,Vy_pred,Vz_pred;
void cart_to_sphere();
- void keep_in_sphere();
void find_tracers();
@@ -421,7 +426,7 @@
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);
+ (E->trace.keep_within_bounds)(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
/* Fill in Current Positions (other positions are no longer important) */
@@ -468,6 +473,8 @@
void eject_tracer();
void reduce_tracer_arrays();
void sphere_to_cart();
+ void full_lost_souls();
+ void regional_lost_souls();
time_stat1=CPU_time0();
@@ -477,7 +484,7 @@
/* initialize arrays and statistical counters */
- E->trace.ilater[j]=0;
+ E->trace.ilater[j]=E->trace.ilatersize[j]=0;
E->trace.istat1=0;
for (kk=0;kk<=4;kk++) {
@@ -529,12 +536,15 @@
parallel_process_termination();
*/
- lost_souls(E);
+ if (E->parallel.nprocxy == 12)
+ full_lost_souls(E);
+ else
+ regional_lost_souls(E);
/* Free later arrays */
for (j=1;j<=E->sphere.caps_per_proc;j++) {
- if (E->trace.ilater[j]>0) {
+ if (E->trace.ilatersize[j]>0) {
for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++) {
free(E->trace.rlater[j][kk]);
}
@@ -553,666 +563,7 @@
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 */
-
-static 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 */
@@ -1336,7 +687,6 @@
int tracers_cap, int j)
{
void cart_to_sphere();
- void keep_in_sphere();
void initialize_tracer_arrays();
int kk;
@@ -1398,7 +748,7 @@
/* Made it, so record tracer information */
- keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+ (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
E->trace.ntracers[j]++;
kk=E->trace.ntracers[j];
@@ -1437,7 +787,6 @@
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();
@@ -1495,7 +844,7 @@
/* make sure theta, phi is in range, and radius is within bounds */
- keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+ (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
/* check whether tracer is within processor domain */
@@ -1572,7 +921,6 @@
void initialize_tracer_arrays();
void sphere_to_cart();
- void keep_in_sphere();
FILE *fp1;
@@ -1631,7 +979,7 @@
/* 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.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
E->trace.basicq[j][0][kk]=theta;
E->trace.basicq[j][1][kk]=phi;
@@ -1956,7 +1304,7 @@
/****** EXPAND TRACER ARRAYS *****************************************/
-static void expand_tracer_arrays(struct All_variables *E, int j)
+void expand_tracer_arrays(struct All_variables *E, int j)
{
int inewsize;
@@ -2086,7 +1434,7 @@
/* The first tracer in initiates memory allocation. */
/* Memory is freed after parallel communications */
- if (E->trace.ilater[j]==0) {
+ if (E->trace.ilatersize[j]==0) {
E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
@@ -2121,7 +1469,7 @@
/****** EXPAND LATER ARRAY *****************************************/
-static void expand_later_array(struct All_variables *E, int j)
+void expand_later_array(struct All_variables *E, int j)
{
int inewsize;
@@ -2235,8 +1583,8 @@
/* is at the surface (then both boundaries are */
/* included). */
-static int icheck_that_processor_shell(struct All_variables *E,
- int j, int nprocessor, double rad)
+int icheck_that_processor_shell(struct All_variables *E,
+ int j, int nprocessor, double rad)
{
int icheck_processor_shell();
int me = E->parallel.me;
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-19 23:44:12 UTC (rev 6303)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-19 23:49:51 UTC (rev 6304)
@@ -139,4 +139,8 @@
void (* get_velocity)(struct All_variables*, int, int,
double, double, double, double*);
+
+ void (* keep_within_bounds)(struct All_variables*,
+ double*, double*, double*,
+ double*, double*, double*);
};
More information about the cig-commits
mailing list