From d413eb9eeed2e00b6bfdfa2d49d5dd789a5b4a75 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 11 Mar 2018 17:10:42 -0600 Subject: [PATCH] faster, simpler, 'more completely random' way to do this --- src/create_atoms.cpp | 177 ++++++++++++++++++++++--------------------- src/create_atoms.h | 4 +- 2 files changed, 93 insertions(+), 88 deletions(-) diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 8a012e6a67..96faab6b0d 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -56,7 +56,7 @@ void CreateAtoms::command(int narg, char **arg) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); - + if (domain->box_exist == 0) error->all(FLERR,"Create_atoms command before simulation box is defined"); if (modify->nfix_restart_peratom) @@ -110,13 +110,13 @@ void CreateAtoms::command(int narg, char **arg) remapflag = 0; mode = ATOM; int molseed; - int nboxseed; + int insertseed; varflag = 0; vstr = xstr = ystr = zstr = NULL; quatone[0] = quatone[1] = quatone[2] = 0.0; - nlattpts = created_Nmask = nbox = nboxflag = 0; + nlattpts = ninsert = insertflag = 0; Nmask = NULL; - + nbasis = domain->lattice->nbasis; basistype = new int[nbasis]; for (int i = 0; i < nbasis; i++) basistype[i] = ntype; @@ -200,9 +200,9 @@ void CreateAtoms::command(int narg, char **arg) iarg += 5; } else if (strcmp(arg[iarg],"insert") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); - nbox = force->inumeric(FLERR,arg[iarg+1]); - nboxseed = force->inumeric(FLERR,arg[iarg+2]); - nboxflag = 1; + ninsert = force->inumeric(FLERR,arg[iarg+1]); + insertseed = force->inumeric(FLERR,arg[iarg+2]); + insertflag = 1; iarg += 3; } else error->all(FLERR,"Illegal create_atoms command"); } @@ -241,10 +241,10 @@ void CreateAtoms::command(int narg, char **arg) ranmol = new RanMars(lmp,molseed+comm->me); } - if (nboxflag) { - ranbox = new RanMars(lmp,nboxseed+comm->me); + if (insertflag) { + ranbox = new RanMars(lmp,insertseed+comm->me); } - + // error check and further setup for variable test if (!vstr && (xstr || ystr || zstr)) @@ -381,7 +381,7 @@ void CreateAtoms::command(int narg, char **arg) else add_lattice(); // init per-atom fix/compute/variable values for created atoms - + atom->data_fix_compute_variable(nlocal_previous,atom->nlocal); // set new total # of atoms and error check @@ -542,7 +542,7 @@ void CreateAtoms::command(int narg, char **arg) delete [] ystr; delete [] zstr; delete [] Nmask; - + // print status if (comm->me == 0) { @@ -764,7 +764,7 @@ void CreateAtoms::add_lattice() double *coord; int i,j,k,m; - + // first pass: count how many particles will be inserted // second pass: filter to N number of particles (and insert) int maskcntr = 0; @@ -774,35 +774,35 @@ void CreateAtoms::add_lattice() for (j = jlo; j <= jhi; j++) for (i = ilo; i <= ihi; i++) for (m = 0; m < nbasis; m++) { - + x[0] = i + basis[m][0]; x[1] = j + basis[m][1]; x[2] = k + basis[m][2]; - + // convert from lattice coords to box coords - + domain->lattice->lattice2box(x[0],x[1],x[2]); - + // if a region was specified, test if atom is in it - + if (style == REGION) if (!domain->regions[nregion]->match(x[0],x[1],x[2])) continue; - + // if variable test specified, eval variable - + if (varflag && vartest(x) == 0) continue; - + // test if atom/molecule position is in my subbox - + if (triclinic) { domain->x2lamda(x,lamda); coord = lamda; } else coord = x; - + if (coord[0] < sublo[0] || coord[0] >= subhi[0] || coord[1] < sublo[1] || coord[1] >= subhi[1] || coord[2] < sublo[2] || coord[2] >= subhi[2]) continue; - + // add the atom or entire molecule to my list of atoms if (pass == 0) nlattpts++; else { @@ -821,72 +821,79 @@ void CreateAtoms::add_lattice() void CreateAtoms::lattice_mask() { - // Nmask will be used to choose which lattice points to insert - Nmask = new int[nlattpts]; - for (int i = 0; i < nlattpts; i++) - Nmask[i] = 1; - - if (nbox == 0 && nboxflag && me == 0) error->warning(FLERR,"Specifying an 'insert' value of '0' is equivalent to no 'insert' keyword"); - - int nboxme = 0; - if (nbox > 0) { + if (ninsert == 0 && insertflag && me == 0) + error->warning(FLERR,"Specifying an 'insert' value of '0' is equivalent to no 'insert' keyword"); + + if (ninsert > 0) { + int nboxme; if (nprocs > 1) { int *allnlattpts = new int[nprocs](); - + int allnboxmes[nprocs]; + int cumsum_lattpts[nprocs]; + for (size_t i = 0; i < nprocs; i++) + allnboxmes[i] = 0; + MPI_Allgather(&nlattpts, 1, MPI_INT, allnlattpts, 1, MPI_INT, world); - - int total_lattpts = 0; - for (int i = 0; i < nprocs; i++) - total_lattpts += allnlattpts[i]; - - if (nbox > total_lattpts) error->all(FLERR,"Attempting to insert more particles than available lattice points"); - if (nbox < 0) error->all(FLERR,"Cannot insert a negative number of particles (in this universe)"); - - // adjust nboxme based on personal number of eligible lattice points - - nboxme = round(nbox*nlattpts/total_lattpts); - - // limit nboxme to available lattice points on this processor - - if (nboxme > nlattpts) nboxme = nlattpts; - - // readjust so that all nboxme's add to exactly nbox - - int sum_nboxme; - MPI_Allreduce(&nboxme,&sum_nboxme,1,MPI_INT,MPI_SUM,world); - - if (nbox != sum_nboxme) { - - int diff = nbox - sum_nboxme; - int diff_sign = 1; - if (diff < 0) { - diff_sign = -1; - diff = -diff; + + if (me == 0) { + int total_lattpts = 0; + for (int i = 0; i < nprocs; i++) { + total_lattpts += allnlattpts[i]; + cumsum_lattpts[i] = total_lattpts; } - - // go around adding or subtracting one from all processors (if there's space) until equal - int which_proc = 0; - int success_me = 0, success_all = 0; - int failed_me = 0, failed_all = 0; - while (success_all < diff && failed_all < nprocs) { - if (me == which_proc++ && !failed_me) - if (nboxme < nlattpts) { - nboxme += diff_sign; - success_me++; - } else { - failed_me == 1; + + + if (ninsert > total_lattpts) + error->all(FLERR,"Attempting to insert more particles than available lattice points"); + if (ninsert < 0) + error->all(FLERR,"Cannot insert a negative number of particles (in this universe)"); + + // using proc 0, let's insert N particles onto available lattice points by instead + // poking [nlattpts - N] holes into the lattice, randomly + int allNmask[total_lattpts]; + for (int i = 0; i < total_lattpts; i++) + allNmask[i] = 1; + int nholes = total_lattpts - ninsert; + int noptions = total_lattpts; + for (int i = 0; i < nholes; i++) { + int hindex = ceil(ranbox->uniform()*noptions); + int optcount = 0; + for (int j = 0; j < total_lattpts; j++) { + if (allNmask[j] == 1) { + optcount++; + if (optcount == hindex) { + allNmask[j] = 0; + noptions--; + break; + } } - MPI_Allreduce(&success_me,&success_all,1,MPI_INT,MPI_SUM,world); - MPI_Allreduce(&failed_me,&failed_all,1,MPI_INT,MPI_SUM,world); - if (which_proc > nprocs) - which_proc = 0; + } + } + + // tell individual procs their nboxme (personal # to insert) + int iproc = 0; + while (allnlattpts[iproc] == 0 && iproc < nprocs) + iproc++; + for (int i = 0; i < total_lattpts; i++) { + if (i == cumsum_lattpts[iproc]) { + iproc++; + while (allnlattpts[iproc] == 0 && iproc < nprocs) + iproc++; + } + allnboxmes[iproc] += allNmask[i]; } - // should have guaranteed success at this point - if (failed_all == nprocs) error->warning(FLERR,"This is an uncaught error. Ask developer"); } + MPI_Scatter(&allnboxmes, 1, MPI_INT, &nboxme, 1, MPI_INT, 0, world); delete [] allnlattpts; - } else nboxme = nbox; - + } else nboxme = ninsert; + + // probably faster to have individual processors 're-choose' their random points + // Nmask will be used to indicate which lattice points to insert + if (nlattpts < nboxme) printf("WHOAA\n"); //debug + Nmask = new int[nlattpts]; + for (int i = 0; i < nlattpts; i++) + Nmask[i] = 1; + // let's insert N particles onto available lattice points by instead // poking [nlattpts - N] holes into the lattice, randomly int nholes = nlattpts - nboxme; @@ -901,14 +908,12 @@ void CreateAtoms::lattice_mask() Nmask[j] = 0; noptions--; break; - } + } } } } } - - created_Nmask = 1; - + } /* ---------------------------------------------------------------------- diff --git a/src/create_atoms.h b/src/create_atoms.h index fa5ee77f8d..7ed4293bec 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -43,7 +43,7 @@ class CreateAtoms : protected Pointers { class Molecule *onemol; class RanMars *ranmol; class RanMars *ranbox; - + int triclinic; double sublo[3],subhi[3]; // epsilon-extended proc sub-box for adding atoms @@ -54,7 +54,7 @@ class CreateAtoms : protected Pointers { void add_molecule(double *, double * = NULL); int nlattpts; // number of eligible lattice points int *Nmask; // used to insert N number of particles on lattice - int created_Nmask,nbox,nboxflag; + int ninsert,insertflag; int vartest(double *); // evaluate a variable with new atom position };