/******************************************************** % % Written by: % -- % John L. Weatherwax 2006-05-29 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- Note: rather than figure out Cartesian communicators/topologies right now I'll just map each processor into my own block row representation (from left to right, top to bottom) of our 2-d domain, i.e. for 12 processors with Ms=3 rows and Ns=4 columns they will fill a virtual domain like 0=(0,0) 1=(0,1) 2=(0,2) 3=(0,3) 4=(1,0) 5=(1,1) 6=(1,2) 7=(1,3) 8=(2,0) 9=(2,1) 10=(2,2) 11=(2,3) which means that process 0 will have control over particles in a restricted spatial domain Note: everything is indexed starting at 0. */ #include #include #include #include "mpi.h" #include "global_grid.h" #include "sim_consts.h" #include "periodicGridMaps.h" #include "saveBoxes.h" void doInitialGridConstruction(int NP, int my_rank, int Ms, int Ns, int NbPP){ double my_lx,my_ly,my_ux,my_uy,my_Rdiag,my_Ldiag,tx,ty; int bi,my_row,my_col; procInfo.NP = NP; procInfo.Ms = Ms; procInfo.Ns = Ns; procInfo.bx = simConsts.Gx/Ns; if( procInfo.bx < 0.5*simConsts.R ){ /* check simulation assumptions ... we don't have TOO many processors */ fprintf(stderr,"bx < R/2.\n"); fprintf(stderr,"one solution is to DECREASE the number of processors one runs with.\n"); MPI_Abort(MPI_COMM_WORLD,1); } procInfo.by = simConsts.Gy/Ms; if( procInfo.by < 0.5*simConsts.R ){ /* check simulation assumptions ... we don't have TOO many processors */ fprintf(stderr,"by < R/2.\n"); fprintf(stderr,"one solution is to DECREASE the number of processors one runs with.\n"); MPI_Abort(MPI_COMM_WORLD,1); } procInfo.Nb = NP*NbPP; /* assign rank/row/column information */ procInfo.my_rank = my_rank; my_row = my_rank/Ns; procInfo.my_row = my_row; my_col = my_rank - (procInfo.my_row)*Ns; procInfo.my_col = my_col; /* WWX: Currently my_Rdiag,my_Ldiag not used ... */ my_Rdiag = my_row + my_col; procInfo.my_Rdiag = my_Rdiag; my_Ldiag = my_col - my_row; procInfo.my_Ldiag = my_Ldiag; /* assign box center information */ procInfo.my_cx = (procInfo.my_col)*procInfo.bx + 0.5*procInfo.bx; procInfo.my_cy = (Ms-1-procInfo.my_row)*procInfo.by + 0.5*procInfo.by; /* assign lower corner of box information */ my_lx = procInfo.my_cx - 0.5*procInfo.bx; procInfo.my_lx = my_lx; my_ly = procInfo.my_cy - 0.5*procInfo.by; procInfo.my_ly = my_ly; /* assign upper corner of box information */ my_ux = my_lx + procInfo.bx; procInfo.my_ux = my_ux; my_uy = my_ly + procInfo.by; procInfo.my_uy = my_uy; /* fill the array "isFaceExternal" --WWX (2006-10-13): checked this logic against the above global grid cartoon ... no errors found */ for( bi=0; bi < 4; bi++ ) procInfo.isFaceExternal[bi]=0; /* by default we assume this box has no external faces */ for( bi=0; bi < Ns; bi++ ){ /* tag the boxes that have their northern face on a global boundary */ if( procInfo.my_rank==bi ) procInfo.isFaceExternal[0]=1; } for( bi=0; bi < Ms; bi++ ){ /* tag the boxes that have their eastern face on a global boundary */ if( procInfo.my_rank==(Ns*(bi+1)-1) ) procInfo.isFaceExternal[1]=1; } for( bi=0; bi < Ns; bi++ ){ /* tag the boxes that have their southern face on a global boundary */ if( procInfo.my_rank==((Ms-1)*Ns+bi) ) procInfo.isFaceExternal[2]=1; } for( bi=0; bi < Ms; bi++ ){ /* tag the boxes that have their western face on a global boundary */ if( procInfo.my_rank==(bi*Ns) ) procInfo.isFaceExternal[3]=1; } /* fill the global domain corners ... the corners of the portion of the entire \mathbb{R}^2 for which each processing box will oversee the particles contained within */ /* initialize global domain to my box corners */ procInfo.gdNW[0] = procInfo.my_lx; procInfo.gdNW[1] = procInfo.my_uy; procInfo.gdNE[0] = procInfo.my_ux; procInfo.gdNE[1] = procInfo.my_uy; procInfo.gdSE[0] = procInfo.my_ux; procInfo.gdSE[1] = procInfo.my_ly; procInfo.gdSW[0] = procInfo.my_lx; procInfo.gdSW[1] = procInfo.my_ly; /* correct these assignments if a boundary of our box is external */ if( procInfo.isFaceExternal[0]==1 ){ /* the north face is external */ procInfo.gdNE[1] = +DBL_MAX; procInfo.gdNW[1] = +DBL_MAX; } if( procInfo.isFaceExternal[1]==1 ){ /* the east face is external */ procInfo.gdNE[0] = +DBL_MAX; procInfo.gdSE[0] = +DBL_MAX; } if( procInfo.isFaceExternal[2]==1 ){ /* the south face is external */ procInfo.gdSE[1] = -DBL_MAX; procInfo.gdSW[1] = -DBL_MAX; } if( procInfo.isFaceExternal[3]==1 ){ /* the west face is external */ procInfo.gdSW[0] = -DBL_MAX; procInfo.gdNW[0] = -DBL_MAX; } /* compute my local neighbors ... needed for message passing ... a PERIODIC grid is assumed and the index is the global processor index */ procInfo.neighbor[0] = northwestNb(Ms,Ns,my_row,my_col); procInfo.neighbor[1] = northNb(Ms,Ns,my_row,my_col); procInfo.neighbor[2] = northeastNb(Ms,Ns,my_row,my_col); procInfo.neighbor[3] = westNb(Ms,Ns,my_row, my_col); procInfo.neighbor[4] = my_rank; /* myself */ procInfo.neighbor[5] = eastNb(Ms,Ns,my_row,my_col); procInfo.neighbor[6] = southwestNb(Ms,Ns,my_row,my_col); procInfo.neighbor[7] = southNb(Ms,Ns,my_row,my_col); procInfo.neighbor[8] = southeastNb(Ms,Ns,my_row,my_col); /* Fill the processors with NbPP initial particles */ procInfo.NbPP = NbPP; for( bi=0; bi < NbPP; bi++ ){ /* generate two uniform random numbers between 0 <= tx,ty <= 1 */ tx = (double) rand()/RAND_MAX; ty = (double) rand()/RAND_MAX; /* generate random particle positions (inside my box) */ procInfo.px[bi] = my_lx + tx*(my_ux-my_lx); procInfo.py[bi] = my_ly + ty*(my_uy-my_ly); /* velocities are mapped between -1 and +1 */ /*procInfo.vx[bi] = 2*((double) rand()/RAND_MAX)-1; procInfo.vy[bi] = 2*((double) rand()/RAND_MAX)-1; */ procInfo.vx[bi] = 0.0; procInfo.vy[bi] = 0.0; /* initialize forces */ procInfo.fx[bi] = 0.0; procInfo.fy[bi] = 0.0; /* initialize masses */ procInfo.m[bi]=1.0; } /* Set up the underlying charge mesh ... */ /* not implemented yet ... */ saveBoxes(my_rank); }