diff --git a/src/comm.cpp b/src/comm.cpp index 3c79b0922d..324c9d0546 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -1006,345 +1006,6 @@ void Comm::reverse_comm_compute(Compute *compute) } } -/* ---------------------------------------------------------------------- - communicate atoms to new owning procs via irregular communication - can be used in place of exchange() - unlike exchange(), allows atoms to have moved arbitrarily long distances - first setup irregular comm pattern, then invoke it - atoms must be remapped to be inside simulation box before this is called - for triclinic: - atoms must be in lamda coords (0-1) before this is called -------------------------------------------------------------------------- */ - -void Comm::irregular() -{ - // clear global->local map since atoms move to new procs - // zero out ghosts so map_set() at end will operate only on local atoms - // exchange() doesn't need to zero ghosts b/c borders() - // is called right after and it zeroes ghosts and calls map_set() - - if (map_style) atom->map_clear(); - atom->nghost = 0; - - // subbox bounds for orthogonal or triclinic - - double *sublo,*subhi; - if (triclinic == 0) { - sublo = domain->sublo; - subhi = domain->subhi; - } else { - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } - - // loop over atoms, flag any that are not in my sub-box - // fill buffer with atoms leaving my box, using < and >= - // assign which proc it belongs to via irregular_lookup() - // if irregular_lookup() returns me, due to round-off - // in triclinic x2lamda(), then keep atom and don't send - // when atom is deleted, fill it in with last atom - - AtomVec *avec = atom->avec; - double **x = atom->x; - int nlocal = atom->nlocal; - - int nsend = 0; - int nsendatom = 0; - int *sizes = new int[nlocal]; - int *proclist = new int[nlocal]; - - int i = 0; - while (i < nlocal) { - if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] || - x[i][1] < sublo[1] || x[i][1] >= subhi[1] || - x[i][2] < sublo[2] || x[i][2] >= subhi[2]) { - proclist[nsendatom] = irregular_lookup(x[i]); - if (proclist[nsendatom] != me) { - if (nsend > maxsend) grow_send(nsend,1); - sizes[nsendatom] = avec->pack_exchange(i,&buf_send[nsend]); - nsend += sizes[nsendatom]; - nsendatom++; - avec->copy(nlocal-1,i); - nlocal--; - } else i++; - } else i++; - } - atom->nlocal = nlocal; - - // create irregular communication plan, perform comm, destroy plan - // returned nrecv = size of buffer needed for incoming atoms - - int nrecv; - Plan *plan = irregular_create(nsendatom,sizes,proclist,&nrecv); - if (nrecv > maxrecv) grow_recv(nrecv); - irregular_perform(plan,buf_send,sizes,buf_recv); - irregular_destroy(plan); - - delete [] sizes; - delete [] proclist; - - // add received atoms to my list - - int m = 0; - while (m < nrecv) m += avec->unpack_exchange(&buf_recv[m]); - - // reset global->local map - - if (map_style) atom->map_set(); -} - -/* ---------------------------------------------------------------------- - create an irregular communication plan - n = # of atoms to send - sizes = # of doubles for each atom - proclist = proc to send each atom to (none to me) - return nrecvsize = total # of doubles I will recv -------------------------------------------------------------------------- */ - -Comm::Plan *Comm::irregular_create(int n, int *sizes, int *proclist, - int *nrecvsize) -{ - int i; - - // allocate plan and work vectors - - Plan *plan = (struct Plan *) memory->smalloc(sizeof(Plan),"comm:plan"); - int *list = new int[nprocs]; - int *count = new int[nprocs]; - - // nrecv = # of messages I receive - - for (i = 0; i < nprocs; i++) { - list[i] = 0; - count[i] = 1; - } - for (i = 0; i < n; i++) list[proclist[i]] = 1; - - int nrecv; - MPI_Reduce_scatter(list,&nrecv,count,MPI_INT,MPI_SUM,world); - - // allocate receive arrays - - int *proc_recv = new int[nrecv]; - int *length_recv = new int[nrecv]; - MPI_Request *request = new MPI_Request[nrecv]; - MPI_Status *status = new MPI_Status[nrecv]; - - // nsend = # of messages I send - - for (i = 0; i < nprocs; i++) list[i] = 0; - for (i = 0; i < n; i++) list[proclist[i]] += sizes[i]; - - int nsend = 0; - for (i = 0; i < nprocs; i++) - if (list[i]) nsend++; - - // allocate send arrays - - int *proc_send = new int[nsend]; - int *length_send = new int[nsend]; - int *num_send = new int[nsend]; - int *index_send = new int[n]; - int *offset_send = new int[n]; - - // list still stores size of message for procs I send to - // proc_send = procs I send to - // length_send = total size of message I send to each proc - // to balance pattern of send messages: - // each proc begins with iproc > me, continues until iproc = me - // reset list to store which send message each proc corresponds to - - int iproc = me; - int isend = 0; - for (i = 0; i < nprocs; i++) { - iproc++; - if (iproc == nprocs) iproc = 0; - if (list[iproc] > 0) { - proc_send[isend] = iproc; - length_send[isend] = list[iproc]; - list[iproc] = isend; - isend++; - } - } - - // num_send = # of datums I send to each proc - - for (i = 0; i < nsend; i++) num_send[i] = 0; - for (i = 0; i < n; i++) { - isend = list[proclist[i]]; - num_send[isend]++; - } - - // count = offsets into n-length index_send for each proc I send to - // index_send = list of which datums to send to each proc - // 1st N1 values are datum indices for 1st proc, - // next N2 values are datum indices for 2nd proc, etc - // offset_send = where each datum starts in send buffer - - count[0] = 0; - for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1]; - - for (i = 0; i < n; i++) { - isend = list[proclist[i]]; - index_send[count[isend]++] = i; - if (i) offset_send[i] = offset_send[i-1] + sizes[i-1]; - else offset_send[i] = 0; - } - - // tell receivers how much data I send - // sendmax = largest # of datums I send in a single message - - int sendmax = 0; - for (i = 0; i < nsend; i++) { - MPI_Send(&length_send[i],1,MPI_INT,proc_send[i],0,world); - sendmax = MAX(sendmax,length_send[i]); - } - - // receive incoming messages - // proc_recv = procs I recv from - // length_recv = total size of message each proc sends me - // nrecvsize = total size of data I recv - - *nrecvsize = 0; - for (i = 0; i < nrecv; i++) { - MPI_Recv(&length_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status); - proc_recv[i] = status->MPI_SOURCE; - *nrecvsize += length_recv[i]; - } - - // barrier to insure all MPI_ANY_SOURCE messages are received - // else another proc could proceed to irregular_perform() and send to me - - MPI_Barrier(world); - - // free work vectors - - delete [] count; - delete [] list; - - // initialize plan and return it - - plan->nsend = nsend; - plan->nrecv = nrecv; - plan->sendmax = sendmax; - - plan->proc_send = proc_send; - plan->length_send = length_send; - plan->num_send = num_send; - plan->index_send = index_send; - plan->offset_send = offset_send; - plan->proc_recv = proc_recv; - plan->length_recv = length_recv; - plan->request = request; - plan->status = status; - - return plan; -} - -/* ---------------------------------------------------------------------- - perform irregular communication - plan = previouly computed plan via irregular_create() - sendbuf = list of atoms to send - sizes = # of doubles for each atom - recvbuf = received atoms -------------------------------------------------------------------------- */ - -void Comm::irregular_perform(Plan *plan, double *sendbuf, int *sizes, - double *recvbuf) -{ - int i,m,n,offset; - - // post all receives - - offset = 0; - for (int irecv = 0; irecv < plan->nrecv; irecv++) { - MPI_Irecv(&recvbuf[offset],plan->length_recv[irecv],MPI_DOUBLE, - plan->proc_recv[irecv],0,world,&plan->request[irecv]); - offset += plan->length_recv[irecv]; - } - - // allocate buf for largest send - - double *buf = (double *) memory->smalloc(plan->sendmax*sizeof(double), - "comm::irregular"); - - // send each message - // pack buf with list of datums (datum = one atom) - // m = index of datum in sendbuf - - n = 0; - for (int isend = 0; isend < plan->nsend; isend++) { - offset = 0; - for (i = 0; i < plan->num_send[isend]; i++) { - m = plan->index_send[n++]; - memcpy(&buf[offset],&sendbuf[plan->offset_send[m]], - sizes[m]*sizeof(double)); - offset += sizes[m]; - } - MPI_Send(buf,plan->length_send[isend],MPI_DOUBLE, - plan->proc_send[isend],0,world); - } - - // free temporary send buffer - - memory->sfree(buf); - - // wait on all incoming messages - - if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status); -} - -/* ---------------------------------------------------------------------- - destroy an irregular communication plan -------------------------------------------------------------------------- */ - -void Comm::irregular_destroy(Plan *plan) -{ - delete [] plan->proc_send; - delete [] plan->length_send; - delete [] plan->num_send; - delete [] plan->index_send; - delete [] plan->offset_send; - delete [] plan->proc_recv; - delete [] plan->length_recv; - delete [] plan->request; - delete [] plan->status; - memory->sfree(plan); -} - -/* ---------------------------------------------------------------------- - determine which proc owns atom with x coord - x will be in box (orthogonal) or lamda coords (triclinic) -------------------------------------------------------------------------- */ - -int Comm::irregular_lookup(double *x) -{ - int loc[3]; - if (triclinic == 0) { - double *boxlo = domain->boxlo; - double *boxhi = domain->boxhi; - loc[0] = static_cast - (procgrid[0] * (x[0]-boxlo[0]) / (boxhi[0]-boxlo[0])); - loc[1] = static_cast - (procgrid[1] * (x[1]-boxlo[1]) / (boxhi[1]-boxlo[1])); - loc[2] = static_cast - (procgrid[2] * (x[2]-boxlo[2]) / (boxhi[2]-boxlo[2])); - } else { - loc[0] = static_cast (procgrid[0] * x[0]); - loc[1] = static_cast (procgrid[1] * x[1]); - loc[2] = static_cast (procgrid[2] * x[2]); - } - - if (loc[0] < 0) loc[0] = 0; - if (loc[0] >= procgrid[0]) loc[0] = procgrid[0] - 1; - if (loc[1] < 0) loc[1] = 0; - if (loc[1] >= procgrid[1]) loc[1] = procgrid[1] - 1; - if (loc[2] < 0) loc[2] = 0; - if (loc[2] >= procgrid[2]) loc[2] = procgrid[2] - 1; - - return grid2proc[loc[0]][loc[1]][loc[2]]; -} - /* ---------------------------------------------------------------------- assign nprocs to 3d xprd,yprd,zprd box so as to minimize surface area area = surface area of each of 3 faces of simulation box diff --git a/src/comm.h b/src/comm.h index 5c9bad339d..330088b848 100644 --- a/src/comm.h +++ b/src/comm.h @@ -21,21 +21,17 @@ namespace LAMMPS_NS { class Comm : protected Pointers { public: int me,nprocs; // proc info - int style; // single vs multi-type comm int procgrid[3]; // assigned # of procs in each dim int user_procgrid[3]; // user request for procs in each dim int myloc[3]; // which proc I am in each dim int procneigh[3][2]; // my 6 neighboring procs - int nswap; // # of swaps to perform - int need[3]; // procs I need atoms from in each dim int ghost_velocity; // 1 if ghost atoms have velocity, 0 if not int maxforward_fix; // comm sizes called from Fix,Pair int maxreverse_fix; int maxforward_pair; int maxreverse_pair; double cutghost[3]; // cutoffs used for acquiring ghost atoms - double cutghostuser; // user-specified ghost cutoff - + int ***grid2proc; // which proc owns i,j,k loc in 3d grid Comm(class LAMMPS *); ~Comm(); @@ -55,12 +51,13 @@ class Comm : protected Pointers { void forward_comm_compute(class Compute *); // forward comm from a Compute void reverse_comm_compute(class Compute *); // reverse comm from a Compute - void irregular(); // irregular communication across all procs - void set(int, char **); // set communication style double memory_usage(); private: + int style; // single vs multi-type comm + int nswap; // # of swaps to perform + int need[3]; // procs I need atoms from in each dim int triclinic; // 0 if domain is orthog, 1 if triclinic int maxswap; // max # of swaps memory is allocated for int size_forward; // # of per-atom datums in forward comm @@ -74,11 +71,11 @@ class Comm : protected Pointers { double *slablo,*slabhi; // bounds of slab to send at each swap double **multilo,**multihi; // bounds of slabs for multi-type swap double **cutghostmulti; // cutghost on a per-type basis + double cutghostuser; // user-specified ghost cutoff int *pbc_flag; // general flag for sending atoms thru PBC int **pbc; // dimension flags for PBC adjustments int comm_x_only,comm_f_only; // 1 if only exchange x,f in for/rev comm int map_style; // non-0 if global->local mapping is done - int ***grid2proc; // which proc owns i,j,k loc in 3d grid int bordergroup; // only communicate this group in borders int *firstrecv; // where to put 1st recv atom in each swap @@ -90,21 +87,6 @@ class Comm : protected Pointers { int maxsend,maxrecv; // current size of send/recv buffer int maxforward,maxreverse; // max # of datums in forward/reverse comm - struct Plan { // plan for irregular communication - int nsend; // # of messages to send - int nrecv; // # of messages to recv - int sendmax; // # of doubles in largest send message - int *proc_send; // procs to send to - int *length_send; // # of doubles to send to each proc - int *num_send; // # of datums to send to each proc - int *index_send; // list of which datums to send to each proc - int *offset_send; // where each datum starts in send buffer - int *proc_recv; // procs to recv from - int *length_recv; // # of doubles to recv from each proc - MPI_Request *request; // MPI requests for posted recvs - MPI_Status *status; // MPI statuses for WaitAll - }; - void procs2box(); // map procs to 3d box void cross(double, double, double, double, double, double, @@ -117,11 +99,6 @@ class Comm : protected Pointers { void allocate_multi(int); // allocate multi arrays void free_swap(); // free swap arrays void free_multi(); // free multi arrays - - struct Plan *irregular_create(int, int *, int *, int *); - void irregular_perform(Plan *, double *, int *, double *); - void irregular_destroy(Plan *); - int irregular_lookup(double *); }; } diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index 7ed3bb9e09..29ffe56c4f 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -16,9 +16,11 @@ #include "string.h" #include "displace_atoms.h" #include "atom.h" +#include "modify.h" #include "domain.h" #include "lattice.h" #include "comm.h" +#include "irregular.h" #include "group.h" #include "random_park.h" #include "error.h" @@ -43,13 +45,9 @@ void DisplaceAtoms::command(int narg, char **arg) if (domain->box_exist == 0) error->all("Displace_atoms command before simulation box is defined"); if (narg < 2) error->all("Illegal displace_atoms command"); - - // init entire system since comm->irregular is done - // comm::init needs neighbor::init needs pair::init needs kspace::init, etc - - if (comm->me == 0 && screen) - fprintf(screen,"System init for displace_atoms ...\n"); - lmp->init(); + if (modify->nfix_restart_peratom) + error->all("Cannot displace_atoms after " + "reading restart file with per-atom info"); if (comm->me == 0 && screen) fprintf(screen,"Displacing atoms ...\n"); @@ -168,7 +166,6 @@ void DisplaceAtoms::command(int narg, char **arg) // move atoms randomly // makes atom result independent of what proc owns it via random->reset() - if (style == RANDOM) { RanPark *random = new RanPark(lmp,1); @@ -197,7 +194,7 @@ void DisplaceAtoms::command(int narg, char **arg) // move atoms back inside simulation box and to new processors // use remap() instead of pbc() in case atoms moved a long distance - // use irregular() instead of exchange() in case atoms moved a long distance + // use irregular() in case atoms moved a long distance double **x = atom->x; int *image = atom->image; @@ -206,8 +203,9 @@ void DisplaceAtoms::command(int narg, char **arg) if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->reset_box(); - comm->setup(); - comm->irregular(); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(); + delete irregular; if (domain->triclinic) domain->lamda2x(atom->nlocal); // check if any atoms were lost diff --git a/src/displace_box.cpp b/src/displace_box.cpp index cf2e80f2e8..58962b795d 100644 --- a/src/displace_box.cpp +++ b/src/displace_box.cpp @@ -17,9 +17,11 @@ #include "string.h" #include "displace_box.h" #include "atom.h" +#include "modify.h" #include "domain.h" #include "lattice.h" #include "comm.h" +#include "irregular.h" #include "group.h" #include "error.h" @@ -42,13 +44,9 @@ void DisplaceBox::command(int narg, char **arg) if (domain->box_exist == 0) error->all("Displace_box command before simulation box is defined"); if (narg < 2) error->all("Illegal displace_box command"); - - // init entire system since comm->irregular is done - // comm::init needs neighbor::init needs pair::init needs kspace::init, etc - - if (comm->me == 0 && screen) - fprintf(screen,"System init for displace_box ...\n"); - lmp->init(); + if (modify->nfix_restart_peratom) + error->all("Cannot displace_box after " + "reading restart file with per-atom info"); if (comm->me == 0 && screen) fprintf(screen,"Displacing box ...\n"); @@ -356,10 +354,9 @@ void DisplaceBox::command(int narg, char **arg) } // move atoms back inside simulation box and to new processors - // use remap() instead of pbc() in case box - // moved a long distance relative to atoms - // use irregular() instead of exchange() in case box - // moved a long distance relative to atoms + // use remap() instead of pbc() + // in case box moved a long distance relative to atoms + // use irregular() in case box moved a long distance relative to atoms double **x = atom->x; int *image = atom->image; @@ -368,8 +365,9 @@ void DisplaceBox::command(int narg, char **arg) if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->reset_box(); - comm->setup(); - comm->irregular(); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(); + delete irregular; if (domain->triclinic) domain->lamda2x(atom->nlocal); // clean up diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 4544ac2c3d..ac8a64ba01 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -22,6 +22,7 @@ #include "atom.h" #include "update.h" #include "comm.h" +#include "irregular.h" #include "domain.h" #include "lattice.h" #include "force.h" @@ -301,6 +302,9 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) rfix = NULL; flip = 0; + if (force_reneighbor) irregular = new Irregular(lmp); + else irregular = NULL; + TWOPI = 8.0*atan(1.0); } @@ -311,6 +315,8 @@ FixDeform::~FixDeform() delete [] set; delete [] rfix; + delete irregular; + // reset domain's h_rate = 0.0, since this fix may have made it non-zero double *h_rate = domain->h_rate; @@ -326,7 +332,7 @@ FixDeform::~FixDeform() int FixDeform::setmask() { int mask = 0; - mask |= PRE_EXCHANGE; + if (force_reneighbor) mask |= PRE_EXCHANGE; mask |= END_OF_STEP; return mask; } @@ -570,7 +576,7 @@ void FixDeform::pre_exchange() for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); domain->x2lamda(atom->nlocal); - comm->irregular(); + irregular->migrate_atoms(); domain->lamda2x(atom->nlocal); flip = 0; diff --git a/src/fix_deform.h b/src/fix_deform.h index 90eb880f50..2c9e3a43f1 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -42,6 +42,7 @@ class FixDeform : public Fix { int kspace_flag; // 1 if KSpace invoked, 0 if not int nrigid; // number of rigid fixes int *rfix; // indices of rigid fixes + class Irregular *irregular; // for migrating atoms after box flips double TWOPI; diff --git a/src/modify.cpp b/src/modify.cpp index 139a87161e..2b39865854 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -91,7 +91,6 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) nfix_restart_peratom = 0; id_restart_peratom = style_restart_peratom = NULL; index_restart_peratom = NULL; - nfix_restart_save = 0; ncompute = maxcompute = 0; compute = NULL; @@ -147,9 +146,8 @@ void Modify::init() int i,j; // delete storage of restart info since it is not valid after 1st run - // read_restart sets nfix_restart_save so will persist to actual 1st run - if (!nfix_restart_save) restart_deallocate(); + restart_deallocate(); // create lists of fixes to call at each stage of run diff --git a/src/modify.h b/src/modify.h index 0cdf25600f..5de137b02c 100644 --- a/src/modify.h +++ b/src/modify.h @@ -32,7 +32,6 @@ class Modify : protected Pointers { int restart_pbc_any; // 1 if any fix sets restart_pbc int nfix_restart_global; // stored fix global info from restart file int nfix_restart_peratom; // stored fix peratom info from restart file - int nfix_restart_save; // 1 if init() should not whack restart info class Fix **fix; // list of fixes int *fmask; // bit mask for when each fix is applied diff --git a/src/read_restart.cpp b/src/read_restart.cpp index c9356d0edb..2b8f59e3bd 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -21,6 +21,7 @@ #include "atom_vec.h" #include "domain.h" #include "comm.h" +#include "irregular.h" #include "update.h" #include "modify.h" #include "fix.h" @@ -236,12 +237,10 @@ void ReadRestart::command(int narg, char **arg) delete [] perproc; - // save modify->nfix_restart values, so Modify::init() won't whack them - // create a temporary fix to hold extra atom info - // necessary b/c Atom::init() will destroy atom->extra + // create a temporary fix to hold and migrate extra atom info + // necessary b/c irregular will migrate atoms if (nextra) { - modify->nfix_restart_save = 1; char cextra[8],fixextra[8]; sprintf(cextra,"%d",nextra); sprintf(fixextra,"%d",modify->nfix_restart_peratom); @@ -255,28 +254,22 @@ void ReadRestart::command(int narg, char **arg) delete [] newarg; } - // init entire system before comm->irregular - // comm::init needs neighbor::init needs pair::init needs kspace::init, etc // move atoms to new processors via irregular() - // in case read by totally different proc than wrote restart file + // in case read by different proc than wrote restart file - if (me == 0 && screen) - fprintf(screen," system init for parallel read_restart ...\n"); - lmp->init(); if (domain->triclinic) domain->x2lamda(atom->nlocal); - domain->reset_box(); - comm->setup(); - comm->irregular(); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(); + delete irregular; if (domain->triclinic) domain->lamda2x(atom->nlocal); - // unsave modify->nfix_restart values and atom->extra + // put extra atom info held by fix back into atom->extra // destroy temporary fix if (nextra) { - modify->nfix_restart_save = 0; - atom->nextra_store = nextra; + memory->destroy_2d_double_array(atom->extra); atom->extra = memory->create_2d_double_array(atom->nmax,nextra, - "atom:extra"); + "atom:extra"); int ifix = modify->find_fix("_read_restart"); FixReadRestart *fix = (FixReadRestart *) modify->fix[ifix]; int *count = fix->count;