diff --git a/doc/src/reset_ids.rst b/doc/src/reset_ids.rst index 13a374fc29..bc7ad31927 100644 --- a/doc/src/reset_ids.rst +++ b/doc/src/reset_ids.rst @@ -8,7 +8,14 @@ Syntax .. code-block:: LAMMPS - reset_ids + reset_ids keyword values ... + + * zero or more keyword/value pairs may be appended + * keyword = *sort* + +.. parsed-literal:: + + *sort* value = *yes* or *no* Examples """""""" @@ -16,6 +23,7 @@ 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 between 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/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp index 6bc0ab587b..23cff55ca2 100644 --- a/src/MANYBODY/pair_adp.cpp +++ b/src/MANYBODY/pair_adp.cpp @@ -554,7 +554,7 @@ void PairADP::read_file(char *filename) ValueTokenizer values = reader.next_values(1); file->nelements = values.next_int(); - if (values.count() != file->nelements + 1) + if ((int)values.count() != file->nelements + 1) error->one(FLERR,"Incorrect element names in ADP potential file"); file->elements = new char*[file->nelements]; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index df90609e7d..73af9e1acf 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -3488,7 +3488,7 @@ void PairAIREBO::read_file(char *filename) // global parameters current_section = "global parameters"; - for(int i = 0; i < params.size(); i++) { + for(int i = 0; i < (int)params.size(); i++) { *params[i] = reader.next_double(); } diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp index a71edccab0..5b039c8eb5 100644 --- a/src/MANYBODY/pair_comb3.cpp +++ b/src/MANYBODY/pair_comb3.cpp @@ -312,7 +312,7 @@ double PairComb3::init_one(int i, int j) void PairComb3::read_lib() { - int i,j,k,l,m; + int i,j,k,l; int ii,jj,kk,ll,mm,iii; // open library file on proc 0 diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index ec409d1786..984075d407 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -469,8 +469,6 @@ void PairEAM::read_file(char *filename) PotentialFileReader reader(lmp, filename, "EAM"); try { - char * line = nullptr; - reader.skip_line(); ValueTokenizer values = reader.next_values(2); diff --git a/src/MANYBODY/pair_eam_alloy.cpp b/src/MANYBODY/pair_eam_alloy.cpp index d164682e9e..42587e9434 100644 --- a/src/MANYBODY/pair_eam_alloy.cpp +++ b/src/MANYBODY/pair_eam_alloy.cpp @@ -130,7 +130,7 @@ void PairEAMAlloy::read_file(char *filename) ValueTokenizer values = reader.next_values(1); file->nelements = values.next_int(); - if (values.count() != file->nelements + 1) + if ((int)values.count() != file->nelements + 1) error->one(FLERR,"Incorrect element names in EAM potential file"); file->elements = new char*[file->nelements]; diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp index 076f6d5e7d..eb8fc4e768 100644 --- a/src/MANYBODY/pair_eam_cd.cpp +++ b/src/MANYBODY/pair_eam_cd.cpp @@ -518,7 +518,7 @@ void PairEAMCD::read_h_coeff(char *filename) int degree = values.next_int(); nhcoeff = degree+1; - if (values.count() != nhcoeff + 1 || nhcoeff < 1) + if ((int)values.count() != nhcoeff + 1 || nhcoeff < 1) error->one(FLERR, "Failed to read h(x) function coefficients in EAM file."); hcoeff = new double[nhcoeff]; diff --git a/src/MANYBODY/pair_eam_fs.cpp b/src/MANYBODY/pair_eam_fs.cpp index 5d2339c4b5..448a1d2701 100644 --- a/src/MANYBODY/pair_eam_fs.cpp +++ b/src/MANYBODY/pair_eam_fs.cpp @@ -122,8 +122,6 @@ void PairEAMFS::read_file(char *filename) PotentialFileReader reader(lmp, filename, "EAM"); try { - char * line = nullptr; - reader.skip_line(); reader.skip_line(); reader.skip_line(); @@ -132,7 +130,7 @@ void PairEAMFS::read_file(char *filename) ValueTokenizer values = reader.next_values(1); file->nelements = values.next_int(); - if (values.count() != file->nelements + 1) + if ((int)values.count() != file->nelements + 1) error->one(FLERR,"Incorrect element names in EAM potential file"); file->elements = new char*[file->nelements]; diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 4fa35cf5fd..56c4773b2b 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -72,6 +72,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : create_attribute = 1; dof_flag = 1; enforce2d_flag = 1; + stores_ids = 1; MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index b1ce975005..32d49fd7b7 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -64,6 +64,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : thermo_virial = 1; create_attribute = 1; dof_flag = 1; + stores_ids = 1; // error check diff --git a/src/USER-OMP/pair_eam_alloy_omp.cpp b/src/USER-OMP/pair_eam_alloy_omp.cpp index ea2130721b..e8f56a9700 100644 --- a/src/USER-OMP/pair_eam_alloy_omp.cpp +++ b/src/USER-OMP/pair_eam_alloy_omp.cpp @@ -130,7 +130,7 @@ void PairEAMAlloyOMP::read_file(char *filename) ValueTokenizer values = reader.next_values(1); file->nelements = values.next_int(); - if (values.count() != file->nelements + 1) + if ((int)values.count() != file->nelements + 1) error->one(FLERR,"Incorrect element names in EAM potential file"); file->elements = new char*[file->nelements]; diff --git a/src/USER-OMP/pair_eam_fs_omp.cpp b/src/USER-OMP/pair_eam_fs_omp.cpp index 71d58469c8..c68d5f6019 100644 --- a/src/USER-OMP/pair_eam_fs_omp.cpp +++ b/src/USER-OMP/pair_eam_fs_omp.cpp @@ -130,7 +130,7 @@ void PairEAMFSOMP::read_file(char *filename) ValueTokenizer values = reader.next_values(1); file->nelements = values.next_int(); - if (values.count() != file->nelements + 1) + if ((int)values.count() != file->nelements + 1) error->one(FLERR,"Incorrect element names in EAM potential file"); file->elements = new char*[file->nelements]; diff --git a/src/atom_vec.h b/src/atom_vec.h index eaff9b888a..8ccf922c4b 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -209,25 +209,6 @@ class AtomVec : protected Pointers { bool *threads; - // union data struct for packing 32-bit and 64-bit ints into double bufs - // this avoids aliasing issues by having 2 pointers (double,int) - // to same buf memory - // constructor for 32-bit int prevents compiler - // from possibly calling the double constructor when passed an int - // copy to a double *buf: - // buf[m++] = ubuf(foo).d, where foo is a 32-bit or 64-bit int - // copy from a double *buf: - // foo = (int) ubuf(buf[m++]).i;, where (int) or (tagint) match foo - // the cast prevents compiler warnings about possible truncation - - union ubuf { - double d; - int64_t i; - ubuf(double arg) : d(arg) {} - ubuf(int64_t arg) : i(arg) {} - ubuf(int arg) : i(arg) {} - }; - // local methods void grow_nmax(); diff --git a/src/comm.cpp b/src/comm.cpp index 9102e2991d..b00308a6e5 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -256,7 +256,6 @@ void Comm::init_exchange() int nfix = modify->nfix; Fix **fix = modify->fix; - int onefix; maxexchange_fix = 0; for (int i = 0; i < nfix; i++) maxexchange_fix += fix[i]->maxexchange; diff --git a/src/compute.h b/src/compute.h index 66974bf106..b6d053dd0e 100644 --- a/src/compute.h +++ b/src/compute.h @@ -161,17 +161,6 @@ class Compute : protected Pointers { return j >> SBBITS & 3; } - // 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) {} - }; - // private methods void adjust_dof_fix(); diff --git a/src/fix.cpp b/src/fix.cpp index 8a82e35d94..c15f5638e2 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -80,6 +80,7 @@ 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; diff --git a/src/fix.h b/src/fix.h index 1cb1cc3ca4..ab8dddf1d9 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 @@ -241,17 +242,6 @@ class Fix : protected Pointers { void v_tally(int, int *, double, double *); void v_tally(int, double *); void v_tally(int, int, double); - - // 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) {} - }; }; namespace FixConst { diff --git a/src/library.cpp b/src/library.cpp index 7840567866..1cacf9e7d8 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -1667,7 +1667,7 @@ int lammps_style_name(void* ptr, char * category, int index, char * buffer, int Info info(lmp); auto styles = info.get_available_styles(category); - if (index < styles.size()) { + if (index < (int)styles.size()) { strncpy(buffer, styles[index].c_str(), max_size); return true; } diff --git a/src/lmptype.h b/src/lmptype.h index 1513fe5782..68af28af61 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -171,6 +171,46 @@ typedef int bigint; #endif + /// Data structe for packing 32-bit and 64-bit integers + /// into double (communication) buffers + /// + /// Using this union avoids aliasing issues by having member types + /// (double, int) referencing the same buffer memory location. + /// + /// The explicit constructor for 32-bit integers prevents compilers + /// from (incorrectly) calling the double constructor when storing + /// an int into a double buffer. + /* +\verbatim embed:rst +**Usage:** + +To copy an integer into a double buffer: + +.. code-block: c++ + + double buf[2]; + int foo = 1; + tagint bar = 2<<40; + buf[1] = ubuf(foo).d; + buf[2] = ubuf(bar).d; + +To copy from a double buffer back to an int: + +.. code-block: c++ + + foo = (int) ubuf(buf[1]).i; + bar = (tagint) ubuf(buf[2]).i; + +The typecast prevents compiler warnings about possible truncations. +\endverbatim + */ + union ubuf { + double d; + int64_t i; + ubuf(const double &arg) : d(arg) {} + ubuf(const int64_t &arg) : i(arg) {} + ubuf(const int &arg) : i(arg) {} + }; } // preprocessor macros for compiler specific settings diff --git a/src/pair.h b/src/pair.h index d578908b20..d0aee8cd7f 100644 --- a/src/pair.h +++ b/src/pair.h @@ -258,17 +258,6 @@ class Pair : protected Pointers { double, double, double, double, double, double); void virial_fdotr_compute(); - // 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) {} - }; - inline int sbmask(int j) const { return j >> SBBITS & 3; } diff --git a/src/reset_ids.cpp b/src/reset_ids.cpp index 20759d2ef6..b8255cd273 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]; } @@ -228,7 +263,7 @@ void ResetIDs::command(int narg, char ** /* arg */) 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..902fe51519 100644 --- a/src/reset_ids.h +++ b/src/reset_ids.h @@ -26,10 +26,33 @@ 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(); }; } diff --git a/src/table_file_reader.cpp b/src/table_file_reader.cpp index 194db8f8eb..2f3317cb2a 100644 --- a/src/table_file_reader.cpp +++ b/src/table_file_reader.cpp @@ -38,7 +38,7 @@ TableFileReader::~TableFileReader() { char * TableFileReader::find_section_start(const std::string & keyword) { char * line = nullptr; - while (line = reader->next_line()) { + while ((line = reader->next_line())) { ValueTokenizer values(line); std::string word = values.next_string(); diff --git a/src/utils.cpp b/src/utils.cpp index 88677bf542..f3a5770f5d 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -547,7 +547,7 @@ std::string utils::get_potential_date(const std::string & path, const std::strin reader.ignore_comments = false; char * line = nullptr; - while (line = reader.next_line()) { + while ((line = reader.next_line())) { ValueTokenizer values(line); while (values.has_next()) { std::string word = values.next_string(); diff --git a/unittest/utils/test_fmtlib.cpp b/unittest/utils/test_fmtlib.cpp index 6d354cf69f..0b1ec16368 100644 --- a/unittest/utils/test_fmtlib.cpp +++ b/unittest/utils/test_fmtlib.cpp @@ -42,45 +42,63 @@ TEST(FmtLib, insert_neg_int) { } TEST(FmtLib, insert_bigint) { - if (sizeof(bigint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG) const bigint val = 9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word 9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_neg_bigint) { - if (sizeof(bigint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG) const bigint val = -9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word -9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_tagint) { - if (sizeof(tagint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) const tagint val = 9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word 9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_neg_tagint) { - if (sizeof(tagint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) const tagint val = -9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word -9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_imageint) { - if (sizeof(imageint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) const imageint val = 9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word 9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_neg_imageint) { - if (sizeof(imageint) == 4) GTEST_SKIP(); +#if defined(LAMMPS_BIGBIG) const imageint val = -9945234592L; auto text = fmt::format("word {}",val); ASSERT_THAT(text, Eq("word -9945234592")); +#else + GTEST_SKIP(); +#endif } TEST(FmtLib, insert_double) {