From d02363b8fbc157e07ab195aeaffab7f34c1b7c11 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 10 Jun 2020 16:33:05 -0600 Subject: [PATCH] adding a reset_ids sort option --- doc/src/reset_ids.rst | 48 +++++- src/fix.cpp | 3 +- src/fix.h | 1 + src/reset_ids.cpp | 367 ++++++++++++++++++++++++++++++++++++++++-- src/reset_ids.h | 34 ++++ 5 files changed, 435 insertions(+), 18 deletions(-) diff --git a/doc/src/reset_ids.rst b/doc/src/reset_ids.rst index 13a374fc29..62fb5d4679 100644 --- a/doc/src/reset_ids.rst +++ b/doc/src/reset_ids.rst @@ -8,14 +8,22 @@ Syntax .. code-block:: LAMMPS - reset_ids + reset_ids keyword values ... -Examples +* zero or more keyword/value pairs may be appended +* keyword = *sort* + + .. parsed-literal:: + + *sort* value = *yes* or *no* + + Examples """""""" .. code-block:: LAMMPS reset_ids + reset_ids sort yes Description """"""""""" @@ -33,11 +41,32 @@ e.g. due to atoms moving outside a simulation box with fixed boundaries (see the "boundary command"), or due to evaporation (see the "fix evaporate" command). -Note that the resetting of IDs is not really a compression, where gaps -in atom IDs are removed by decrementing atom IDs that are larger. -Instead the IDs for all atoms are erased, and new IDs are assigned so -that the atoms owned by an individual processor have consecutive IDs, -as the :doc:`create_atoms ` command explains. +If the *sort* keyword is used with a setting of *yes*, then the +assignment of new atom IDs will be the same no matter how many +processors LAMMPS is running on. This is done by first doing a +spatial sort of all the atoms into bins and sorting them within each +bin. Because the set of bins is independent of the number of +processors, this enables a consistent assignment of new IDs to each +atom. + +This can be useful to do after using the "create_atoms" command and/or +"replicate" command. In general those commands do not guarantee +assignment of the same atom ID to the same physical atom when LAMMPS +is run on different numbers of processors. Enforcing consistent IDs +can be useful for debugging or comparing output from two different +runs. + +Note that the spatial sort requires communication of atom IDs and +coordinates bewteen processors in an all-to-all manner. This is done +efficiently in LAMMPS, but it is more expensive than how atom IDs are +reset without sorting. + +Note that whether sorting or not, the resetting of IDs is not a +compression, where gaps in atom IDs are removed by decrementing atom +IDs that are larger. Instead the IDs for all atoms are erased, and +new IDs are assigned so that the atoms owned by an individual +processor have consecutive IDs, as the :doc:`create_atoms +` command explains. .. note:: @@ -59,4 +88,7 @@ Related commands :doc:`delete_atoms ` -**Default:** none +Default +""""""" + +By default, *sort* is no. diff --git a/src/fix.cpp b/src/fix.cpp index 8a82e35d94..1ac5a5fcdb 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -80,7 +80,8 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : maxexchange = 0; maxexchange_dynamic = 0; pre_exchange_migrate = 0; - + stores_ids = 0; + scalar_flag = vector_flag = array_flag = 0; peratom_flag = local_flag = 0; global_freq = local_freq = peratom_freq = -1; diff --git a/src/fix.h b/src/fix.h index 1cb1cc3ca4..134b1550e9 100644 --- a/src/fix.h +++ b/src/fix.h @@ -64,6 +64,7 @@ class Fix : protected Pointers { int maxexchange; // max # of per-atom values for Comm::exchange() int maxexchange_dynamic; // 1 if fix sets maxexchange dynamically int pre_exchange_migrate; // 1 if fix migrates atoms in pre_exchange() + int stores_ids; // 1 if fix stores atom IDs int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists diff --git a/src/reset_ids.cpp b/src/reset_ids.cpp index 20759d2ef6..2398c6c39a 100644 --- a/src/reset_ids.cpp +++ b/src/reset_ids.cpp @@ -13,10 +13,13 @@ #include "reset_ids.h" #include +#include #include "atom.h" #include "atom_vec.h" #include "domain.h" #include "comm.h" +#include "modify.h" +#include "fix.h" #include "special.h" #include "memory.h" #include "error.h" @@ -24,25 +27,54 @@ using namespace LAMMPS_NS; +#if defined(LMP_QSORT) +// allocate space for static class variable +// prototype for non-class function +ResetIDs::AtomRvous *ResetIDs::sortrvous; +static int compare_coords(const void *, const void *); +#else +#include "mergesort.h" +// prototype for non-class function +static int compare_coords(const int, const int, void *); +#endif + +#define PERBIN 10 +#define BIG 1.0e20 + /* ---------------------------------------------------------------------- */ ResetIDs::ResetIDs(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- */ -void ResetIDs::command(int narg, char ** /* arg */) +void ResetIDs::command(int narg, char **arg) { if (domain->box_exist == 0) error->all(FLERR,"Reset_ids command before simulation box is defined"); - if (narg != 0) error->all(FLERR,"Illegal reset_ids command"); if (atom->tag_enable == 0) error->all(FLERR,"Cannot use reset_ids unless atoms have IDs"); - // NOTE: check if any fixes exist which store atom IDs? - // if so, this operation will mess up the fix + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->stores_ids) + error->all(FLERR,"Cannot use reset_ids when a fix exists that stores atom IDs"); if (comm->me == 0) utils::logmesg(lmp,"Resetting atom IDs ...\n"); + // process args + + int sortflag = 0; + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"sort") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_ids command"); + if (strcmp(arg[iarg+1],"yes") == 0) sortflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) sortflag = 0; + else error->all(FLERR,"Illegal reset_ids command"); + iarg += 2; + } else error->all(FLERR,"Illegal reset_ids command"); + } + // create an atom map if one doesn't exist already int mapflag = 0; @@ -84,9 +116,12 @@ void ResetIDs::command(int narg, char ** /* arg */) tag[i] = 0; } - // assign new contiguous IDs to owned atoms via tag_extend() + // assign new contiguous IDs to owned atoms + // if sortflag = no: use fast tag_extend() + // if sortflag = yes: use slower full spatial sort plus rendezvous comm - atom->tag_extend(); + if (sortflag == 0) atom->tag_extend(); + else sort(); // newIDs = copy of new IDs // restore old IDs, consistent with existing atom map @@ -96,7 +131,7 @@ void ResetIDs::command(int narg, char ** /* arg */) memory->create(newIDs,nall,1,"reset_ids:newIDs"); for (int i = 0; i < nlocal; i++) { - newIDs[i][0] = tag[i]; + newIDs[i][0] = ubuf(tag[i]).d; tag[i] = oldIDs[i]; } @@ -225,10 +260,10 @@ void ResetIDs::command(int narg, char ** /* arg */) } // reset IDs and atom map for owned atoms - + atom->map_clear(); atom->nghost = 0; - for (int i = 0; i < nlocal; i++) tag[i] = static_cast (newIDs[i][0]); + for (int i = 0; i < nlocal; i++) tag[i] = (tagint) ubuf(newIDs[i][0]).i; atom->map_init(); atom->map_set(); @@ -251,3 +286,317 @@ void ResetIDs::command(int narg, char ** /* arg */) memory->destroy(oldIDs); memory->destroy(newIDs); } + +/* ---------------------------------------------------------------------- + spatial sort of atoms followed by rendezvous comm to reset atom IDs +------------------------------------------------------------------------- */ + +void ResetIDs::sort() +{ + double mylo[3],myhi[3],bboxlo[3],bboxhi[3]; + + int me = comm->me; + int nprocs = comm->nprocs; + int dim = domain->dimension; + + // bboxlo,bboxhi = bounding box on all atoms in system + // expanded by 0.01 percent + // bbox should work for orthogonal or triclinic system + + double **x = atom->x; + int nlocal = atom->nlocal; + + mylo[0] = mylo[1] = mylo[2] = BIG; + myhi[0] = myhi[1] = myhi[2] = -BIG; + + for (int i = 0; i < nlocal; i++) { + mylo[0] = MIN(mylo[0],x[i][0]); + mylo[1] = MIN(mylo[1],x[i][1]); + mylo[2] = MIN(mylo[2],x[i][2]); + myhi[0] = MAX(myhi[0],x[i][0]); + myhi[1] = MAX(myhi[1],x[i][1]); + myhi[2] = MAX(myhi[2],x[i][2]); + } + + if (dim == 2) mylo[2] = myhi[2] = 0.0; + + MPI_Allreduce(mylo,bboxlo,3,MPI_DOUBLE,MPI_MIN,world); + MPI_Allreduce(myhi,bboxhi,3,MPI_DOUBLE,MPI_MAX,world); + + bboxlo[0] -= 0.0001 * (bboxhi[0]-bboxlo[0]); + bboxlo[1] -= 0.0001 * (bboxhi[1]-bboxlo[1]); + bboxlo[2] -= 0.0001 * (bboxhi[2]-bboxlo[2]); + bboxhi[0] += 0.0001 * (bboxhi[0]-bboxlo[0]); + bboxhi[1] += 0.0001 * (bboxhi[1]-bboxlo[1]); + bboxhi[2] += 0.0001 * (bboxhi[2]-bboxlo[2]); + + // nbin_estimate = estimate of total number of bins, each with PERBIN atoms + // binsize = edge length of a cubic bin + // nbin xyz = bin count in each dimension + + bigint nbin_estimate = atom->natoms / PERBIN; + double vol; + if (dim == 2) vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]); + else vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]) * (bboxhi[2]-bboxlo[2]); + double binsize = pow(vol/nbin_estimate,1.0/dim); + + int nbinx = static_cast ((bboxhi[0]-bboxlo[0]) / binsize) + 1; + int nbiny = static_cast ((bboxhi[1]-bboxlo[1]) / binsize) + 1; + int nbinz = static_cast ((bboxhi[2]-bboxlo[2]) / binsize) + 1; + + double invx = 1.0 / (bboxhi[0]-bboxlo[0]); + double invy = 1.0 / (bboxhi[1]-bboxlo[1]); + double invz; + if (dim == 2) invz = 0.0; + else invz = 1.0 / (bboxhi[2]-bboxlo[2]); + + // nbins = actual total number of bins + // bins are split evenly across procs + // lo-numbered procs may have own one bin less than rest of procs + // nlo = # of bins/proc for lowest procs + // nhi = nlo+1 = # of bins/proc for highest procs + // nplo = # of procs in low group + // nbinlo = # of bins owned by low group procs + // binlo to binhi-1 = bin indices this proc will own in Rvous decomp + // bins are numbered from 0 to Nbins-1 + + bigint nbins = (bigint) nbinx*nbiny*nbinz; + int nlo = nbins / nprocs; + int nhi = nlo + 1; + int nplo = nprocs - (nbins % nprocs); + bigint nbinlo = nplo*nlo; + + if (me < nplo) { + binlo = me * nlo; + binhi = (me+1) * nlo; + } else { + binlo = nbinlo + (me-nplo) * nhi; + binhi = nbinlo + (me+1-nplo) * nhi; + } + + // fill atombuf with info on my atoms + // ibin = which bin the atom is in + // proclist = proc that owns ibin + + int *proclist; + memory->create(proclist,nlocal,"special:proclist"); + AtomRvous *atombuf = (AtomRvous *) + memory->smalloc((bigint) nlocal*sizeof(AtomRvous),"resetIDs:idbuf"); + + int ibinx,ibiny,ibinz,iproc; + bigint ibin; + + for (int i = 0; i < nlocal; i++) { + ibinx = static_cast ((x[i][0]-bboxlo[0])*invx * nbinx); + ibiny = static_cast ((x[i][1]-bboxlo[1])*invy * nbiny); + ibinz = static_cast ((x[i][2]-bboxlo[2])*invz * nbinz); + ibin = (bigint) ibinz*nbiny*nbinx + (bigint) ibiny*nbinx + ibinx; + + if (ibin < nbinlo) iproc = ibin / nlo; + else iproc = nplo + (ibin-nbinlo) / nhi; + proclist[i] = iproc; + + atombuf[i].ibin = ibin; + atombuf[i].proc = me; + atombuf[i].ilocal = i; + atombuf[i].x[0] = x[i][0]; + atombuf[i].x[1] = x[i][1]; + atombuf[i].x[2] = x[i][2]; + } + + // perform rendezvous operation, send atombuf to other procs + + char *buf; + int nreturn = comm->rendezvous(1,nlocal,(char *) atombuf,sizeof(AtomRvous), + 0,proclist, + sort_bins,0,buf,sizeof(IDRvous),(void *) this); + IDRvous *outbuf = (IDRvous *) buf; + + memory->destroy(proclist); + memory->sfree(atombuf); + + // set new ID for all owned atoms + + int ilocal; + + for (int i = 0; i < nreturn; i++) { + ilocal = outbuf[i].ilocal; + atom->tag[ilocal] = outbuf[i].newID; + } + + memory->sfree(outbuf); +} + +/* ---------------------------------------------------------------------- + process data for atoms assigned to me in rendezvous decomposition + inbuf = list of N AtomRvous datums + outbuf = list of N IDRvous datums, sent back to sending proc +------------------------------------------------------------------------- */ + +int ResetIDs::sort_bins(int n, char *inbuf, + int &flag, int *&proclist, char *&outbuf, + void *ptr) +{ + int i,ibin,index; + + ResetIDs *rptr = (ResetIDs *) ptr; + Memory *memory = rptr->memory; + Error *error = rptr->error; + MPI_Comm world = rptr->world; + bigint binlo = rptr->binlo; + bigint binhi = rptr->binhi; + + // nbins = my subset of bins from binlo to binhi-1 + // initialize linked lists of my Rvous atoms in each bin + + int *next,*head,*last,*count; + + int nbins = binhi - binlo; + memory->create(next,n,"resetIDs:next"); + memory->create(head,nbins,"resetIDs:head"); + memory->create(last,nbins,"resetIDs:last"); + memory->create(count,nbins,"resetIDs:count"); + + for (i = 0; i < n; i++) next[i] = -1; + + for (ibin = 0; ibin < nbins; ibin++) { + head[ibin] = last[ibin] = -1; + count[ibin] = 0; + } + + AtomRvous *in = (AtomRvous *) inbuf; + + for (i = 0; i < n; i++) { + if (in[i].ibin < binlo || in[i].ibin >= binhi) { + //printf("BAD me %d i %d %d ibin %d binlohi %d %d\n", + // rptr->comm->me,i,n,in[i].ibin,binlo,binhi); + error->one(FLERR,"Bad spatial bin assignment in reset_ids sort"); + } + ibin = in[i].ibin - binlo; + if (head[ibin] < 0) head[ibin] = i; + if (last[ibin] >= 0) next[last[ibin]] = i; + last[ibin] = i; + count[ibin]++; + } + + int maxcount = 0; + for (ibin = 0; ibin < nbins; ibin++) + maxcount = MAX(maxcount,count[ibin]); + + int *order; + memory->create(order,maxcount,"resetIDs:order"); + + // sort atoms in each bin by their spatial position + // recreate linked list for each bin based on sorted order + + for (ibin = 0; ibin < nbins; ibin++) { + index = head[ibin]; + for (i = 0; i < count[ibin]; i++) { + order[i] = index; + index = next[index]; + } + + +#if defined(LMP_QSORT) + sortrvous = in; + qsort(order,count[ibin],sizeof(int),compare_coords); +#else + merge_sort(order,count[ibin],(void *) in,compare_coords); +#endif + + head[ibin] = last[ibin] = -1; + for (i = 0; i < count[ibin]; i++) { + if (head[ibin] < 0) head[ibin] = order[i]; + if (last[ibin] >= 0) next[last[ibin]] = order[i]; + last[ibin] = order[i]; + } + } + + // MPI_Scan() to find out where my atoms begin in globally sorted list + + tagint ntag = n; + tagint nprev; + MPI_Scan(&ntag,&nprev,1,MPI_LMP_TAGINT,MPI_SUM,world); + nprev -= n; + + // loop over bins and sorted atoms in each bin, reset ids to be consecutive + // pass outbuf of IDRvous datums back to comm->rendezvous + + int nout = n; + memory->create(proclist,nout,"resetIDs:proclist"); + IDRvous *out = (IDRvous *) + memory->smalloc(nout*sizeof(IDRvous),"resetIDs:out"); + + tagint one = nprev + 1; + for (ibin = 0; ibin < nbins; ibin++) { + index = head[ibin]; + for (i = 0; i < count[ibin]; i++) { + proclist[index] = in[index].proc; + out[index].newID = one++; + out[index].ilocal = in[index].ilocal; + index = next[index]; + } + } + + outbuf = (char *) out; + + // clean up + // Comm::rendezvous will delete proclist and out (outbuf) + + memory->destroy(next); + memory->destroy(head); + memory->destroy(last); + memory->destroy(count); + memory->destroy(order); + + // flag = 2: new outbuf + + flag = 2; + return nout; +} + +#if defined(LMP_QSORT) + +/* ---------------------------------------------------------------------- + comparison function invoked by qsort() + accesses static class member sortrvous, set before call to qsort() +------------------------------------------------------------------------- */ + +int compare_coords(const void *iptr, const void *jptr) +{ + int i = *((int *) iptr); + int j = *((int *) jptr); + ResetIDs::AtomRvous *rvous = ResetIDs::sortrvous; + double *xi = rvous[i].x; + double *xj = rvous[j].x; + if (xi[0] < xj[0]) return -1; + if (xi[0] > xj[0]) return 1; + if (xi[1] < xj[1]) return -1; + if (xi[1] > xj[1]) return 1; + if (xi[2] < xj[2]) return -1; + if (xi[2] > xj[2]) return 1; + return 0; +} + +#else + +/* ---------------------------------------------------------------------- + comparison function invoked by merge_sort() + void pointer contains sortrvous +------------------------------------------------------------------------- */ + +int compare_coords(const int i, const int j, void *ptr) +{ + ResetIDs::AtomRvous *rvous = (ResetIDs::AtomRvous *) ptr; + double *xi = rvous[i].x; + double *xj = rvous[j].x; + if (xi[0] < xj[0]) return -1; + if (xi[0] > xj[0]) return 1; + if (xi[1] < xj[1]) return -1; + if (xi[1] > xj[1]) return 1; + if (xi[2] < xj[2]) return -1; + if (xi[2] > xj[2]) return 1; + return 0; +} + +#endif diff --git a/src/reset_ids.h b/src/reset_ids.h index d146651cc9..69582bd5c3 100644 --- a/src/reset_ids.h +++ b/src/reset_ids.h @@ -26,10 +26,44 @@ namespace LAMMPS_NS { class ResetIDs : protected Pointers { public: + struct AtomRvous { + bigint ibin; + int proc,ilocal; + double x[3]; + }; + + struct IDRvous { + tagint newID; + int ilocal; + }; + + #if defined(LMP_QSORT) + // static variable across all ResetID objects, for qsort callback + static AtomRvous *sortrvous; +#endif + ResetIDs(class LAMMPS *); void command(int, char **); private: + bigint binlo,binhi; + + // callback functions for rendezvous communication + + static int sort_bins(int, char *, int &, int *&, char *&, void *); + + void sort(); + + // union data struct for packing 32-bit and 64-bit ints into double bufs + // see atom_vec.h for documentation + + union ubuf { + double d; + int64_t i; + ubuf(double arg) : d(arg) {} + ubuf(int64_t arg) : i(arg) {} + ubuf(int arg) : i(arg) {} + }; }; }