From 5b5798769f1d5420d9e2e5cb98626713dc90f26f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 3 Jul 2013 15:50:00 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10215 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_pour.h | 2 +- src/comm.cpp | 8 +- src/fix_shear_history.cpp | 153 ++++++++++++++++++++++++++++---------- src/fix_shear_history.h | 9 ++- src/modify.cpp | 3 +- src/neigh_gran.cpp | 4 +- src/neigh_list.cpp | 3 +- src/neigh_list.h | 6 +- src/neigh_stencil.cpp | 3 +- src/neighbor.cpp | 10 ++- src/neighbor.h | 2 + src/random_mars.cpp | 1 + 12 files changed, 147 insertions(+), 57 deletions(-) diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h index 5c596715d7..96fd44bcde 100644 --- a/src/GRANULAR/fix_pour.h +++ b/src/GRANULAR/fix_pour.h @@ -30,6 +30,7 @@ class FixPour : public Fix { friend class PairGranHooke; friend class PairGranHookeOMP; friend class PairGranHookeHistory; + friend class PairGranHookeHistory2; friend class PairGranHookeHistoryOMP; friend class PairGranHookeCuda; @@ -61,7 +62,6 @@ class FixPour : public Fix { int *recvcounts,*displs; int nfreq,nfirst,ninserted,nper; double lo_current,hi_current; - class FixShearHistory *fix_history; class RanPark *random; int overlap(int); diff --git a/src/comm.cpp b/src/comm.cpp index 4090fc2cea..ef7ed1afda 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -1024,8 +1024,7 @@ void Comm::borders() // pack up list of border atoms - if (nsend*size_border > maxsend) - grow_send(nsend*size_border,0); + if (nsend*size_border > maxsend) grow_send(nsend*size_border,0); if (ghost_velocity) n = avec->pack_border_vel(nsend,sendlist[iswap],buf_send, pbc_flag[iswap],pbc[iswap]); @@ -1466,10 +1465,7 @@ void Comm::forward_comm_array(int n, double **array) MPI_Request request; MPI_Status status; - // check that buf_send and buf_recv are big enough - - - + // NOTE: check that buf_send and buf_recv are big enough for (iswap = 0; iswap < nswap; iswap++) { diff --git a/src/fix_shear_history.cpp b/src/fix_shear_history.cpp index f9e321a235..da0def908b 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_shear_history.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "mpi.h" #include "string.h" #include "stdio.h" #include "fix_shear_history.h" @@ -27,8 +28,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -#define MAXTOUCH 15 - /* ---------------------------------------------------------------------- */ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : @@ -47,10 +46,15 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : atom->add_callback(0); atom->add_callback(1); + ipage = NULL; + dpage = NULL; + pgsize = oneatom = 0; + // initialize npartner to 0 so neighbor list creation is OK the 1st time int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) npartner[i] = 0; + maxtouch = 0; } /* ---------------------------------------------------------------------- */ @@ -65,8 +69,10 @@ FixShearHistory::~FixShearHistory() // delete locally stored arrays memory->destroy(npartner); - memory->destroy(partner); - memory->destroy(shearpartner); + memory->sfree(partner); + memory->sfree(shearpartner); + delete ipage; + delete dpage; } /* ---------------------------------------------------------------------- */ @@ -89,6 +95,21 @@ void FixShearHistory::init() int dim; computeflag = (int *) pair->extract("computeflag",dim); + + // create pages if first time or if neighbor pgsize/oneatom has changed + // note that latter could cause shear history info to be discarded + + int create = 0; + if (ipage == NULL) create = 1; + if (pgsize != neighbor->pgsize) create = 1; + if (oneatom != neighbor->oneatom) create = 1; + + if (create) { + pgsize = neighbor->pgsize; + oneatom = neighbor->oneatom; + ipage = new MyPage(oneatom,pgsize); + dpage = new MyPage(oneatom,pgsize); + } } /* ---------------------------------------------------------------------- @@ -119,17 +140,22 @@ void FixShearHistory::setup_pre_exchange() void FixShearHistory::pre_exchange() { - int i,j,ii,jj,m,inum,jnum; + int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *shear,*allshear,**firstshear; // zero npartner for all current atoms + // clear 2 page data structures int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) npartner[i] = 0; - // copy shear info from neighbor list atoms to atom arrays + ipage->reset(); + dpage->reset(); + + // 1st loop over neighbor list + // calculate npartner for each owned atom int *tag = atom->tag; NeighList *list = pair->list; @@ -140,6 +166,39 @@ void FixShearHistory::pre_exchange() firsttouch = list->listgranhistory->firstneigh; firstshear = list->listgranhistory->firstdouble; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + touch = firsttouch[i]; + + for (jj = 0; jj < jnum; jj++) { + if (touch[jj]) { + npartner[i]++; + j = jlist[jj]; + j &= NEIGHMASK; + if (j < nlocal) npartner[j]++; + } + } + } + + // get page chunks to store atom IDs and shear history for my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + n = npartner[i]; + partner[i] = ipage->get(n); + shearpartner[i] = dpage->get(n); + if (partner[i] == NULL || shearpartner[i] == NULL) + error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + } + + // 2nd loop over neighbor list + // store atom IDs and shear history for my atoms + // re-zero npartner to use as counter for all my atoms + + for (i = 0; i < nlocal; i++) npartner[i] = 0; + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; @@ -152,37 +211,28 @@ void FixShearHistory::pre_exchange() shear = &allshear[3*jj]; j = jlist[jj]; j &= NEIGHMASK; - if (npartner[i] < MAXTOUCH) { - m = npartner[i]; - partner[i][m] = tag[j]; - shearpartner[i][m][0] = shear[0]; - shearpartner[i][m][1] = shear[1]; - shearpartner[i][m][2] = shear[2]; - } + m = npartner[i]; + partner[i][m] = tag[j]; + shearpartner[i][m][0] = shear[0]; + shearpartner[i][m][1] = shear[1]; + shearpartner[i][m][2] = shear[2]; npartner[i]++; if (j < nlocal) { - if (npartner[j] < MAXTOUCH) { - m = npartner[j]; - partner[j][m] = tag[i]; - shearpartner[j][m][0] = -shear[0]; - shearpartner[j][m][1] = -shear[1]; - shearpartner[j][m][2] = -shear[2]; - } + m = npartner[j]; + partner[j][m] = tag[i]; + shearpartner[j][m][0] = -shear[0]; + shearpartner[j][m][1] = -shear[1]; + shearpartner[j][m][2] = -shear[2]; npartner[j]++; } } } } - // test for too many touching neighbors + // set maxtouch = max # of partners of any owned atom - int flag = 0; - for (i = 0; i < nlocal; i++) - if (npartner[i] >= MAXTOUCH) flag = 1; - int flag_all; - MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); - if (flag_all) - error->all(FLERR,"Too many touching neighbors - boost MAXTOUCH"); + maxtouch = 0; + for (i = 0; i < nlocal; i++) maxtouch = MAX(maxtouch,npartner[i]); } /* ---------------------------------------------------------------------- */ @@ -208,8 +258,10 @@ double FixShearHistory::memory_usage() { int nmax = atom->nmax; double bytes = nmax * sizeof(int); - bytes += nmax*MAXTOUCH * sizeof(int); - bytes += nmax*MAXTOUCH*3 * sizeof(double); + bytes = nmax * sizeof(int *); + bytes = nmax * sizeof(double *); + bytes += ipage->ndatum * sizeof(int); + bytes += dpage->ndatum * sizeof(double[3]); return bytes; } @@ -220,8 +272,12 @@ double FixShearHistory::memory_usage() void FixShearHistory::grow_arrays(int nmax) { memory->grow(npartner,nmax,"shear_history:npartner"); - memory->grow(partner,nmax,MAXTOUCH,"shear_history:partner"); - memory->grow(shearpartner,nmax,MAXTOUCH,3,"shear_history:shearpartner"); + partner = (int **) memory->srealloc(partner,nmax*sizeof(int *), + "shear_history:partner"); + typedef double (*sptype)[3]; + shearpartner = (sptype *) + memory->srealloc(shearpartner,nmax*sizeof(sptype), + "shear_history:shearpartner"); } /* ---------------------------------------------------------------------- @@ -230,13 +286,15 @@ void FixShearHistory::grow_arrays(int nmax) void FixShearHistory::copy_arrays(int i, int j, int delflag) { + // just copy pointers for partner and shearpartner + // b/c can't overwrite chunk allocation inside ipage,dpage + // incoming atoms in unpack_exchange just grab new chunks + // so are orphaning chunks for migrating atoms + // OK, b/c will reset ipage,dpage on next reneighboring + npartner[j] = npartner[i]; - for (int m = 0; m < npartner[j]; m++) { - partner[j][m] = partner[i][m]; - shearpartner[j][m][0] = shearpartner[i][m][0]; - shearpartner[j][m][1] = shearpartner[i][m][1]; - shearpartner[j][m][2] = shearpartner[i][m][2]; - } + partner[j] = partner[i]; + shearpartner[j] = shearpartner[i]; } /* ---------------------------------------------------------------------- @@ -254,6 +312,9 @@ void FixShearHistory::set_arrays(int i) int FixShearHistory::pack_exchange(int i, double *buf) { + // NOTE: how do I know comm buf is big enough if extreme # of touching neighs + // Comm::BUFEXTRA may need to be increased + int m = 0; buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { @@ -271,8 +332,13 @@ int FixShearHistory::pack_exchange(int i, double *buf) int FixShearHistory::unpack_exchange(int nlocal, double *buf) { + // allocate new chunks from ipage,dpage for incoming values + int m = 0; npartner[nlocal] = static_cast (buf[m++]); + maxtouch = MAX(maxtouch,npartner[nlocal]); + partner[nlocal] = ipage->get(npartner[nlocal]); + shearpartner[nlocal] = dpage->get(npartner[nlocal]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (buf[m++]); shearpartner[nlocal][n][0] = buf[m++]; @@ -314,7 +380,12 @@ void FixShearHistory::unpack_restart(int nlocal, int nth) for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; + // allocate new chunks from ipage,dpage for incoming values + npartner[nlocal] = static_cast (extra[nlocal][m++]); + maxtouch = MAX(maxtouch,npartner[nlocal]); + partner[nlocal] = ipage->get(npartner[nlocal]); + shearpartner[nlocal] = dpage->get(npartner[nlocal]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (extra[nlocal][m++]); shearpartner[nlocal][n][0] = extra[nlocal][m++]; @@ -329,7 +400,11 @@ void FixShearHistory::unpack_restart(int nlocal, int nth) int FixShearHistory::maxsize_restart() { - return 4*MAXTOUCH + 2; + // maxtouch_all = max touching partners across all procs + + int maxtouch_all; + MPI_Allreduce(&maxtouch,&maxtouch_all,1,MPI_INT,MPI_MAX,world); + return 4*maxtouch_all + 2; } /* ---------------------------------------------------------------------- diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index 7dbf953c42..ed73338473 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -21,13 +21,13 @@ FixStyle(SHEAR_HISTORY,FixShearHistory) #define LMP_FIX_SHEAR_HISTORY_H #include "fix.h" +#include "my_page.h" namespace LAMMPS_NS { class FixShearHistory : public Fix { friend class Neighbor; friend class PairGranHookeHistory; - friend class FixPour; public: FixShearHistory(class LAMMPS *, int, char **); @@ -53,10 +53,15 @@ class FixShearHistory : public Fix { protected: int *npartner; // # of touching partners of each atom int **partner; // tags for the partners - double ***shearpartner; // 3 shear values with the partner + double (**shearpartner)[3]; // 3 shear values with the partner + int maxtouch; // max # of touching partners for my atoms class Pair *pair; int *computeflag; // computeflag in PairGranHookeHistory + + int pgsize,oneatom; // copy of settings in Neighbor + MyPage *ipage; // pages of partner atom IDs + MyPage *dpage; // pages of shear history with partners }; } diff --git a/src/modify.cpp b/src/modify.cpp index 8c7e296cd0..1bc1753e33 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -863,7 +863,8 @@ void Modify::modify_compute(int narg, char **arg) int icompute; for (icompute = 0; icompute < ncompute; icompute++) if (strcmp(arg[0],compute[icompute]->id) == 0) break; - if (icompute == ncompute) error->all(FLERR,"Could not find compute_modify ID"); + if (icompute == ncompute) + error->all(FLERR,"Could not find compute_modify ID"); compute[icompute]->modify_params(narg-1,&arg[1]); } diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp index a5a70b2b40..b9fd604f61 100644 --- a/src/neigh_gran.cpp +++ b/src/neigh_gran.cpp @@ -38,7 +38,7 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) NeighList *listgranhistory; int *npartner,**partner; - double ***shearpartner; + double (**shearpartner)[3]; int **firsttouch; double **firstshear; int **pages_touch; @@ -280,7 +280,7 @@ void Neighbor::granular_bin_no_newton(NeighList *list) NeighList *listgranhistory; int *npartner,**partner; - double ***shearpartner; + double (**shearpartner)[3]; int **firsttouch; double **firstshear; int **pages_touch; diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp index e294dcaf48..316b39ed77 100644 --- a/src/neigh_list.cpp +++ b/src/neigh_list.cpp @@ -28,10 +28,11 @@ enum{NSQ,BIN,MULTI}; // also in neighbor.cpp /* ---------------------------------------------------------------------- */ -NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp) +NeighList::NeighList(LAMMPS *lmp, int size, int onesize) : Pointers(lmp) { maxatoms = 0; pgsize = size; + oneatom = onesize; inum = gnum = 0; ilist = NULL; diff --git a/src/neigh_list.h b/src/neigh_list.h index 536b5c2ae2..d29f3a7db6 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -15,6 +15,7 @@ #define LMP_NEIGH_LIST_H #include "pointers.h" +#include "my_page.h" namespace LAMMPS_NS { @@ -39,11 +40,14 @@ class NeighList : protected Pointers { double **firstdouble; // ptr to 1st J double value of each I atom int pgsize; // size of each page + int oneatom; // max size for one atom int maxpage; // # of pages currently allocated int **pages; // neighbor list pages for ints double **dpages; // neighbor list pages for doubles int dnum; // # of doubles for each pair (0 if none) + MyPage *page; + // atom types to skip when building list // iskip,ijskip are just ptrs to corresponding request @@ -76,7 +80,7 @@ class NeighList : protected Pointers { class CudaNeighList *cuda_list; // CUDA neighbor list - NeighList(class LAMMPS *, int); + NeighList(class LAMMPS *, int, int); ~NeighList(); void grow(int); // grow maxlocal void stencil_allocate(int, int); // allocate stencil arrays diff --git a/src/neigh_stencil.cpp b/src/neigh_stencil.cpp index 0474e3348b..a39ebf66bb 100644 --- a/src/neigh_stencil.cpp +++ b/src/neigh_stencil.cpp @@ -23,8 +23,7 @@ using namespace LAMMPS_NS; sx,sy,sz = bin bounds = furthest the stencil could possibly extend 3d creates xyz stencil, 2d creates xy stencil for half list with newton off: - stencil is all surrounding bins - stencil includes self + stencil is all surrounding bins including self regardless of triclinic for half list with newton on: stencil is bins to the "upper right" of central bin diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 2afb7de38f..0f597b64d7 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -119,6 +119,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) old_style = BIN; old_triclinic = 0; + old_pgsize = pgsize; + old_oneatom = oneatom; old_nrequest = 0; old_requests = NULL; @@ -416,12 +418,14 @@ void Neighbor::init() // test if pairwise lists need to be re-created // no need to re-create if: - // neigh style and triclinic has not changed and + // neigh style, triclinic, pgsize, oneatom have not changed // current requests = old requests int same = 1; if (style != old_style) same = 0; if (triclinic != old_triclinic) same = 0; + if (pgsize != old_pgsize) same = 0; + if (oneatom != old_oneatom) same = 0; if (nrequest != old_nrequest) same = 0; else for (i = 0; i < nrequest; i++) @@ -452,7 +456,7 @@ void Neighbor::init() // pass list ptr back to requestor (except for Command class) for (i = 0; i < nlist; i++) { - lists[i] = new NeighList(lmp,pgsize); + lists[i] = new NeighList(lmp,pgsize,oneatom); lists[i]->index = i; lists[i]->dnum = requests[i]->dnum; @@ -1676,10 +1680,12 @@ void Neighbor::modify_params(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"page") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command"); + old_pgsize = pgsize; pgsize = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"one") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command"); + old_oneatom = oneatom; oneatom = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"binsize") == 0) { diff --git a/src/neighbor.h b/src/neighbor.h index 2c58eb18d5..13590dd473 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -49,6 +49,8 @@ class Neighbor : protected Pointers { int old_style; // previous run info to avoid int old_nrequest; // re-creation of pairwise neighbor lists int old_triclinic; + int old_pgsize; + int old_oneatom; class NeighRequest **old_requests; int nlist; // pairwise neighbor lists diff --git a/src/random_mars.cpp b/src/random_mars.cpp index 387cf6a781..fb5725b8c3 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -12,6 +12,7 @@ ------------------------------------------------------------------------- */ // Marsaglia random number generator +// see RANMAR in F James, Comp Phys Comm, 60, 329 (1990) #include "math.h" #include "random_mars.h"