From 6bc507cd6aed552d5662ae3aad1a112c73701e03 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Sat, 9 Jul 2016 20:41:07 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15274 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/neigh_gran.cpp | 136 ++++++++++++++++++++++++------------------ src/neigh_request.cpp | 55 ++++++++++++----- src/neighbor.cpp | 87 ++++++++++++++++++++------- 3 files changed, 182 insertions(+), 96 deletions(-) diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp index 3549212c51..a537fbae09 100644 --- a/src/neigh_gran.cpp +++ b/src/neigh_gran.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include #include "neighbor.h" #include "neigh_list.h" #include "atom.h" @@ -30,7 +31,7 @@ using namespace LAMMPS_NS; void Neighbor::granular_nsq_no_newton(NeighList *list) { - int i,j,m,n,nn,bitmask; + int i,j,m,n,nn,bitmask,dnum,dnumbytes; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; @@ -75,6 +76,8 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) firstshear = listgranhistory->firstdouble; ipage_touch = listgranhistory->ipage; dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); } int inum = 0; @@ -120,20 +123,17 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } @@ -170,7 +170,7 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) void Neighbor::granular_nsq_newton(NeighList *list) { - int i,j,m,n,nn,itag,jtag,bitmask; + int i,j,m,n,nn,itag,jtag,bitmask,dnum,dnumbytes; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; @@ -215,6 +215,8 @@ void Neighbor::granular_nsq_newton(NeighList *list) firstshear = listgranhistory->firstdouble; ipage_touch = listgranhistory->ipage; dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); } int inum = 0; @@ -277,20 +279,17 @@ void Neighbor::granular_nsq_newton(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } @@ -316,6 +315,18 @@ void Neighbor::granular_nsq_newton(NeighList *list) list->inum = inum; } +/* ---------------------------------------------------------------------- + granular particles + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + shear history must be accounted for when a neighbor pair is added + pair added to list if atoms i and j are both owned and i < j + pair added if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::granular_nsq_newton_onesided(NeighList *list) +{ +} + /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with partial Newton's 3rd law @@ -327,12 +338,13 @@ void Neighbor::granular_nsq_newton(NeighList *list) void Neighbor::granular_bin_no_newton(NeighList *list) { - int i,j,k,m,n,nn,ibin; + int i,j,k,m,n,nn,ibin,dnum,dnumbytes; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; double *shearptr; + NeighList *listgranhistory; int *npartner; tagint **partner; @@ -376,6 +388,8 @@ void Neighbor::granular_bin_no_newton(NeighList *list) firstshear = listgranhistory->firstdouble; ipage_touch = listgranhistory->ipage; dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); } int inum = 0; @@ -426,20 +440,17 @@ void Neighbor::granular_bin_no_newton(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } @@ -476,7 +487,7 @@ void Neighbor::granular_bin_no_newton(NeighList *list) void Neighbor::granular_bin_newton(NeighList *list) { - int i,j,k,m,n,nn,ibin; + int i,j,k,m,n,nn,ibin,dnum,dnumbytes; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; @@ -525,6 +536,8 @@ void Neighbor::granular_bin_newton(NeighList *list) firstshear = listgranhistory->firstdouble; ipage_touch = listgranhistory->ipage; dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); } int inum = 0; @@ -579,20 +592,17 @@ void Neighbor::granular_bin_newton(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } @@ -623,20 +633,17 @@ void Neighbor::granular_bin_newton(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } @@ -663,6 +670,18 @@ void Neighbor::granular_bin_newton(NeighList *list) list->inum = inum; } +/* ---------------------------------------------------------------------- + granular particles + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + shear history must be accounted for when a neighbor pair is added + pair added to list if atoms i and j are both owned and i < j + pair added if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::granular_bin_newton_onesided(NeighList *list) +{ +} + /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with Newton's 3rd law for triclinic @@ -673,7 +692,7 @@ void Neighbor::granular_bin_newton(NeighList *list) void Neighbor::granular_bin_newton_tri(NeighList *list) { - int i,j,k,m,n,nn,ibin; + int i,j,k,m,n,nn,ibin,dnum,dnumbytes; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; @@ -722,6 +741,8 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) firstshear = listgranhistory->firstdouble; ipage_touch = listgranhistory->ipage; dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); } int inum = 0; @@ -781,20 +802,17 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; - shearptr[nn++] = shearpartner[i][m][0]; - shearptr[nn++] = shearpartner[i][m][1]; - shearptr[nn++] = shearpartner[i][m][2]; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } else { touchptr[n] = 0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; - shearptr[nn++] = 0.0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; } } diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index a9c0924387..cc8c0d7c99 100644 --- a/src/neigh_request.cpp +++ b/src/neigh_request.cpp @@ -26,23 +26,28 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) id = 0; unprocessed = 1; - // default is pair request + // class user of list: default is pair request + // only one is set to 1 pair = 1; fix = compute = command = 0; - // default is half neighbor list + // kind of list: default is half neighbor list + // only one is set to 1 half = 1; full = 0; full_cluster = 0; - gran = granhistory = 0; + gran = 0; respainner = respamiddle = respaouter = 0; half_from_full = 0; + // combination of settings, mutiple can be set to 1 // default is every reneighboring // default is use newton_pair setting in force // default is encode special bond flags + // default is no granular history (when gran = 1) + // default is no one-sided sphere/surface interactions (when gran = 1) // default is no auxiliary floating point values // default is no neighbors of ghosts // default is no multi-threaded neighbor list build @@ -52,6 +57,8 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) occasional = 0; newton = 0; //special = 1; + granhistory = 0; + granonesided = 0; dnum = 0; ghost = 0; omp = 0; @@ -59,7 +66,8 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) kokkos_host = kokkos_device = 0; ssa = 0; - // default is no copy or skip + // copy/skip/derive info, default is no copy or skip + // none or only one option is set copy = 0; skip = 0; @@ -119,7 +127,6 @@ int NeighRequest::identical(NeighRequest *other) if (half != other->half_original) same = 0; if (full != other->full) same = 0; if (gran != other->gran) same = 0; - if (granhistory != other->granhistory) same = 0; if (respainner != other->respainner) same = 0; if (respamiddle != other->respamiddle) same = 0; if (respaouter != other->respaouter) same = 0; @@ -128,6 +135,8 @@ int NeighRequest::identical(NeighRequest *other) if (newton != other->newton) same = 0; if (occasional != other->occasional) same = 0; //if (special != other->special) same = 0; + if (granhistory != other->granhistory) same = 0; + if (granonesided != other->granonesided) same = 0; if (dnum != other->dnum) same = 0; if (ghost != other->ghost) same = 0; if (omp != other->omp) same = 0; @@ -150,6 +159,10 @@ int NeighRequest::same_kind(NeighRequest *other) { int same = 1; + // class caller will be same (pair), b/c called by pair hybrid styles + + // kind must be the the same + if (half != other->half) same = 0; if (full != other->full) same = 0; if (gran != other->gran) same = 0; @@ -158,12 +171,20 @@ int NeighRequest::same_kind(NeighRequest *other) if (respamiddle != other->respamiddle) same = 0; if (respaouter != other->respaouter) same = 0; if (half_from_full != other->half_from_full) same = 0; + + // settings that must be the same + // other settings do not need to be the same (e.g. granonesided) + if (newton != other->newton) same = 0; if (ghost != other->ghost) same = 0; if (omp != other->omp) same = 0; if (intel != other->intel) same = 0; + if (kokkos_host != other->kokkos_host) same = 0; + if (kokkos_device != other->kokkos_device) same = 0; if (ssa != other->ssa) same = 0; + // copy/skip/derive info does not need to be the same + return same; } @@ -192,26 +213,28 @@ int NeighRequest::same_skip(NeighRequest *other) } /* ---------------------------------------------------------------------- - set kind and other values of this request to that of other request + set kind and optional values in this request to those of other request + not copy/skip/derive values since this may be non-skip parent of child ------------------------------------------------------------------------- */ void NeighRequest::copy_request(NeighRequest *other) { - half = 0; - - if (other->half) half = 1; - if (other->full) full = 1; - if (other->gran) gran = 1; - if (other->granhistory) granhistory = 1; - if (other->respainner) respainner = 1; - if (other->respamiddle) respamiddle = 1; - if (other->respaouter) respaouter = 1; - if (other->half_from_full) half_from_full = 1; + half = other->half; + full = other->full; + gran = other->gran; + granhistory = other->granhistory; + respainner = other->respainner; + respamiddle = other->respamiddle; + respaouter = other->respaouter; + half_from_full = other->half_from_full; newton = other->newton; + granonesided = other->granonesided; dnum = other->dnum; ghost = other->ghost; omp = other->omp; intel = other->intel; + kokkos_host = other->kokkos_host; + kokkos_device = other->kokkos_device; ssa = other->ssa; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index f82a20acd9..307ea10dd8 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -153,6 +153,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) old_check = dist_check; old_cutoff = cutneighmax; + zeroes = NULL; + // bond lists maxbond = 0; @@ -213,6 +215,8 @@ Neighbor::~Neighbor() for (int i = 0; i < old_nrequest; i++) delete old_requests[i]; memory->sfree(old_requests); + delete [] zeroes; + memory->destroy(bondlist); memory->destroy(anglelist); memory->destroy(dihedrallist); @@ -582,13 +586,11 @@ void Neighbor::init() } // granhistory: set preceeding list's listgranhistory to this list - // also set preceeding list's ptr to FixShearHistory + // also set FH ptr in preceeding list to FSH class created by pair if (requests[i]->granhistory) { lists[i-1]->listgranhistory = lists[i]; - for (int ifix = 0; ifix < modify->nfix; ifix++) - if (strcmp(modify->fix[ifix]->style,"SHEAR_HISTORY") == 0) - lists[i-1]->fix_history = (FixShearHistory *) modify->fix[ifix]; + lists[i-1]->fix_history = requests[i]->fix_history; processed = 1; // respaouter: point this list at preceeding 1/2 inner/middle lists @@ -669,11 +671,21 @@ void Neighbor::init() } // allocate initial pages for each list, except if listcopy set + // allocate dnum vector of zeroes if set + int dnummax = 0; for (i = 0; i < nrequest; i++) { if (!lists[i]) continue; - if (!lists[i]->listcopy) + if (!lists[i]->listcopy) { lists[i]->setup_pages(pgsize,oneatom,requests[i]->dnum); + dnummax = MAX(dnummax,requests[i]->dnum); + } + } + + if (dnummax) { + delete [] zeroes; + zeroes = new double[dnummax]; + for (i = 0; i < dnummax; i++) zeroes[i] = 0.0; } // set ptrs to pair_build and stencil_create functions for each list @@ -1028,17 +1040,17 @@ int Neighbor::request(void *requestor, int instance) determine which pair_build function each neigh list needs based on settings of neigh request copy -> copy_from function - skip -> granular function if gran with granhistory, + skip -> granular function if gran, several options respa function if respaouter, skip_from function for everything else ssa -> special case for USER-DPD pair styles half_from_full, half, full, gran, respaouter -> - choose by newton and rq->newton and tri settings + choose by newton and rq->newton and triclinic settings style NSQ options = newton off, newton on style BIN options = newton off, newton on and not tri, newton on and tri stlye MULTI options = same options as BIN if none of these, ptr = NULL since pair_build is not invoked for this list - use "else if" b/c skip,copy can be set in addition to half,full,etc + use "else if" logic b/c skip,copy can be set in addition to half,full,etc ------------------------------------------------------------------------- */ void Neighbor::choose_build(int index, NeighRequest *rq) @@ -1050,8 +1062,19 @@ void Neighbor::choose_build(int index, NeighRequest *rq) if (rq->copy) pb = &Neighbor::copy_from; else if (rq->skip) { - if (rq->gran && lists[index]->listgranhistory) - pb = &Neighbor::skip_from_granular; + if (rq->gran) { + NeighRequest *otherrq = requests[rq->otherlist]; + if (otherrq->newton == 0) { + pb = &Neighbor::skip_from_granular; + } else if (otherrq->newton == 1) { + error->all(FLERR,"Neighbor build method not supported"); + } else if (otherrq->newton == 2) { + if (rq->granonesided == 0) + pb = &Neighbor::skip_from_granular_off2on; + else if (rq->granonesided == 1) + pb = &Neighbor::skip_from_granular_off2on_onesided; + } + } else if (rq->respaouter) pb = &Neighbor::skip_from_respa; else pb = &Neighbor::skip_from; @@ -1145,16 +1168,39 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } } else if (rq->gran) { - if (style == NSQ) { - if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton; - else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton; - } else if (style == BIN) { - if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton; - else if (triclinic == 0) pb = &Neighbor::granular_bin_newton; - else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri; - } else if (style == MULTI) - error->all(FLERR,"Neighbor multi not yet enabled for granular"); - + if (rq->newton == 0) { + if (style == NSQ) { + if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton; + else if (newton_pair == 1) { + if (rq->granonesided == 0) pb = &Neighbor::granular_nsq_newton; + else pb = &Neighbor::granular_nsq_newton_onesided; + } + } else if (style == BIN) { + if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton; + else if (newton_pair == 1) { + if (triclinic == 0) { + if (rq->granonesided == 0) pb = &Neighbor::granular_bin_newton; + else pb = &Neighbor::granular_bin_newton_onesided; + } else if (triclinic == 1) { + if (rq->granonesided == 0) + pb = &Neighbor::granular_bin_newton_tri; + else error->all(FLERR,"Neighbor build method not supported"); + } + } + } else if (style == MULTI) + error->all(FLERR,"Neighbor multi not yet enabled for granular"); + } else if (rq->newton == 1) { + error->all(FLERR,"Neighbor build method not yet supported"); + } else if (rq->newton == 2) { + if (style == NSQ) pb = &Neighbor::granular_nsq_no_newton; + else if (style == BIN) { + if (triclinic == 0) pb = &Neighbor::granular_bin_no_newton; + else if (triclinic == 1) + error->all(FLERR,"Neighbor build method not yet supported"); + } else if (style == MULTI) + error->all(FLERR,"Neighbor multi not yet enabled for granular"); + } + } else if (rq->respaouter) { if (style == NSQ) { if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton; @@ -2232,4 +2278,3 @@ int Neighbor::exclude_setting() { return exclude; } -