From ee6f9b96211ae8d0ae713e9f4e5bd4550f4a4841 Mon Sep 17 00:00:00 2001 From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa> Date: Sat, 9 Jul 2016 20:41:38 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15276 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/neigh_derive.cpp | 436 ++++++++++++++++++++++++++++++++++++++++--- src/neighbor.h | 20 +- src/pair.h | 2 +- src/pair_hybrid.cpp | 77 ++++++-- 4 files changed, 485 insertions(+), 50 deletions(-) diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp index 5de3750968..7161421449 100644 --- a/src/neigh_derive.cpp +++ b/src/neigh_derive.cpp @@ -11,12 +11,23 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include <string.h> #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" +#include "fix_shear_history.h" #include "my_page.h" #include "error.h" + + + +// DEBUG +#include "update.h" +#include "comm.h" + + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- @@ -222,16 +233,28 @@ void Neighbor::skip_from(NeighList *list) /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip - this is for granular lists with history, copy the history values from parent + if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void Neighbor::skip_from_granular(NeighList *list) { - int i,j,ii,jj,n,nn,itype,jnum,joriginal; + int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes; + tagint jtag; int *neighptr,*jlist,*touchptr,*touchptr_skip; double *shearptr,*shearptr_skip; + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage<int> *ipage_touch; + MyPage<double> *dpage_shear; + + tagint *tag = atom->tag; int *type = atom->type; + int nlocal = atom->nlocal; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -241,23 +264,33 @@ void Neighbor::skip_from_granular(NeighList *list) int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; int **firstneigh_skip = list->listskip->firstneigh; - int **firsttouch_skip = list->listskip->listgranhistory->firstneigh; - double **firstshear_skip = list->listskip->listgranhistory->firstdouble; int inum_skip = list->listskip->inum; int *iskip = list->iskip; int **ijskip = list->ijskip; - NeighList *listgranhistory = list->listgranhistory; - int **firsttouch = listgranhistory->firstneigh; - double **firstshear = listgranhistory->firstdouble; - MyPage<int> *ipage_touch = listgranhistory->ipage; - MyPage<double> *dpage_shear = listgranhistory->dpage; + FixShearHistory *fix_history = list->fix_history; + if (fix_history) { + fix_history->nlocal_neigh = nlocal; + fix_history->nall_neigh = nlocal + atom->nghost; + npartner = fix_history->npartner; + partner = fix_history->partner; + shearpartner = fix_history->shearpartner; + listgranhistory = list->listgranhistory; + firsttouch = listgranhistory->firstneigh; + firstshear = listgranhistory->firstdouble; + ipage_touch = listgranhistory->ipage; + dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); + } int inum = 0; ipage->reset(); - ipage_touch->reset(); - dpage_shear->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] @@ -268,15 +301,16 @@ void Neighbor::skip_from_granular(NeighList *list) itype = type[i]; if (iskip[itype]) continue; - n = nn = 0; + n = 0; neighptr = ipage->vget(); - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); + if (fix_history) { + nn = 0; + touchptr = ipage_touch->vget(); + shearptr = dpage_shear->vget(); + } - // loop over parent non-skip granular list and its history info + // loop over parent non-skip granular list and optionally its history info - touchptr_skip = firsttouch_skip[i]; - shearptr_skip = firstshear_skip[i]; jlist = firstneigh_skip[i]; jnum = numneigh_skip[i]; @@ -285,10 +319,166 @@ void Neighbor::skip_from_granular(NeighList *list) j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr[n] = joriginal; - touchptr[n++] = touchptr_skip[jj]; - shearptr[nn++] = shearptr_skip[3*jj]; - shearptr[nn++] = shearptr_skip[3*jj+1]; - shearptr[nn++] = shearptr_skip[3*jj+2]; + + // no numeric test for current touch + // just use FSH partner list to infer it + // would require distance calculation for spheres + // more complex calculation for surfs + + if (fix_history) { + jtag = tag[j]; + for (m = 0; m < npartner[i]; m++) + if (partner[i][m] == jtag) break; + if (m < npartner[i]) { + touchptr[n] = 1; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; + } else { + touchptr[n] = 0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; + } + } + + n++; + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + + if (fix_history) { + firsttouch[i] = touchptr; + firstshear[i] = shearptr; + ipage_touch->vgot(n); + dpage_shear->vgot(nn); + } + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + build skip list for subset of types from parent list + iskip and ijskip flag which atom types and type pairs to skip + parent non-skip list used newton off, this skip list is newton on + if list requests it, preserve shear history via fix shear/history +------------------------------------------------------------------------- */ + +void Neighbor::skip_from_granular_off2on(NeighList *list) +{ + int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes; + tagint itag,jtag; + int *neighptr,*jlist,*touchptr,*touchptr_skip; + double *shearptr,*shearptr_skip; + + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage<int> *ipage_touch; + MyPage<double> *dpage_shear; + + tagint *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage<int> *ipage = list->ipage; + + int *ilist_skip = list->listskip->ilist; + int *numneigh_skip = list->listskip->numneigh; + int **firstneigh_skip = list->listskip->firstneigh; + int inum_skip = list->listskip->inum; + + int *iskip = list->iskip; + int **ijskip = list->ijskip; + + FixShearHistory *fix_history = list->fix_history; + if (fix_history) { + fix_history->nlocal_neigh = nlocal; + fix_history->nall_neigh = nlocal + atom->nghost; + npartner = fix_history->npartner; + partner = fix_history->partner; + shearpartner = fix_history->shearpartner; + listgranhistory = list->listgranhistory; + firsttouch = listgranhistory->firstneigh; + firstshear = listgranhistory->firstdouble; + ipage_touch = listgranhistory->ipage; + dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); + } + + int inum = 0; + ipage->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } + + // loop over atoms in other list + // skip I atom entirely if iskip is set for type[I] + // skip I,J pair if ijskip is set for type[I],type[J] + + for (ii = 0; ii < inum_skip; ii++) { + i = ilist_skip[ii]; + itype = type[i]; + if (iskip[itype]) continue; + itag = tag[i]; + + n = 0; + neighptr = ipage->vget(); + if (fix_history) { + nn = 0; + touchptr = ipage_touch->vget(); + shearptr = dpage_shear->vget(); + } + + // loop over parent non-skip granular list and optionally its history info + + jlist = firstneigh_skip[i]; + jnum = numneigh_skip[i]; + + for (jj = 0; jj < jnum; jj++) { + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (ijskip[itype][type[j]]) continue; + + // only keep I,J when J = ghost if Itag < Jtag + + jtag = tag[j]; + if (j >= nlocal && jtag < itag) continue; + + neighptr[n] = joriginal; + + // no numeric test for current touch + // just use FSH partner list to infer it + // would require distance calculation for spheres + // more complex calculation for surfs + + if (fix_history) { + for (m = 0; m < npartner[i]; m++) + if (partner[i][m] == jtag) break; + if (m < npartner[i]) { + touchptr[n] = 1; + memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes); + nn += dnum; + } else { + touchptr[n] = 0; + memcpy(&shearptr[nn],zeroes,dnumbytes); + nn += dnum; + } + } + + n++; } ilist[inum++] = i; @@ -298,10 +488,206 @@ void Neighbor::skip_from_granular(NeighList *list) if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); + if (fix_history) { + firsttouch[i] = touchptr; + firstshear[i] = shearptr; + ipage_touch->vgot(n); + dpage_shear->vgot(nn); + } + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + build skip list for subset of types from parent list + iskip and ijskip flag which atom types and type pairs to skip + parent non-skip list used newton off and was not onesided, + this skip list is newton on and onesided + if list requests it, preserve shear history via fix shear/history +------------------------------------------------------------------------- */ + +void Neighbor::skip_from_granular_off2on_onesided(NeighList *list) +{ + int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,flip,dnum,dnumbytes,tmp; + tagint itag,jtag; + int *surf,*neighptr,*jlist; + + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage<int> *ipage_touch; + MyPage<double> *dpage_shear; + + tagint *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage<int> *ipage = list->ipage; + + int *ilist_skip = list->listskip->ilist; + int *numneigh_skip = list->listskip->numneigh; + int **firstneigh_skip = list->listskip->firstneigh; + int inum_skip = list->listskip->inum; + + int *iskip = list->iskip; + int **ijskip = list->ijskip; + + if (domain->dimension == 2) surf = atom->line; + else surf = atom->tri; + + FixShearHistory *fix_history = list->fix_history; + if (fix_history) { + fix_history->nlocal_neigh = nlocal; + fix_history->nall_neigh = nlocal + atom->nghost; + npartner = fix_history->npartner; + partner = fix_history->partner; + shearpartner = fix_history->shearpartner; + listgranhistory = list->listgranhistory; + firsttouch = listgranhistory->firstneigh; + firstshear = listgranhistory->firstdouble; + ipage_touch = listgranhistory->ipage; + dpage_shear = listgranhistory->dpage; + dnum = listgranhistory->dnum; + dnumbytes = dnum * sizeof(double); + } + + int inum = 0; + ipage->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } + + // two loops over parent list required, one to count, one to store + // because onesided constraint means pair I,J may be stored with I or J + // so don't know in advance how much space to alloc for each atom's neighs + + // first loop over atoms in other list to count neighbors + // skip I atom entirely if iskip is set for type[I] + // skip I,J pair if ijskip is set for type[I],type[J] + + for (i = 0; i < nlocal; i++) numneigh[i] = 0; + + for (ii = 0; ii < inum_skip; ii++) { + i = ilist_skip[ii]; + itype = type[i]; + if (iskip[itype]) continue; + itag = tag[i]; + + n = 0; + + // loop over parent non-skip granular list + + jlist = firstneigh_skip[i]; + jnum = numneigh_skip[i]; + + for (jj = 0; jj < jnum; jj++) { + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (ijskip[itype][type[j]]) continue; + + // flip I,J if necessary to satisfy onesided constraint + // do not keep if I is now ghost + + if (surf[i] >= 0) { + if (j >= nlocal) continue; + tmp = i; + i = j; + j = tmp; + flip = 1; + } else flip = 0; + + numneigh[i]++; + if (flip) i = j; + } + } + + // allocate all per-atom neigh list chunks, including history + + for (i = 0; i < nlocal; i++) { + if (numneigh[i] == 0) continue; + n = numneigh[i]; + firstneigh[i] = ipage->get(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + if (fix_history) { + firsttouch[i] = ipage_touch->get(n); + firstshear[i] = dpage_shear->get(dnum*n); + } + } + + // second loop over atoms in other list to store neighbors + // skip I atom entirely if iskip is set for type[I] + // skip I,J pair if ijskip is set for type[I],type[J] + + for (i = 0; i < nlocal; i++) numneigh[i] = 0; + + for (ii = 0; ii < inum_skip; ii++) { + i = ilist_skip[ii]; + itype = type[i]; + if (iskip[itype]) continue; + itag = tag[i]; + + // loop over parent non-skip granular list and optionally its history info + + jlist = firstneigh_skip[i]; + jnum = numneigh_skip[i]; + + for (jj = 0; jj < jnum; jj++) { + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (ijskip[itype][type[j]]) continue; + + // flip I,J if necessary to satisfy onesided constraint + // do not keep if I is now ghost + + if (surf[i] >= 0) { + if (j >= nlocal) continue; + tmp = i; + i = j; + j = tmp; + flip = 1; + } else flip = 0; + + // store j in neigh list, not joriginal, like other neigh methods + // OK, b/c there is no special list flagging for surfs + + firstneigh[i][numneigh[i]] = j; + + // no numeric test for current touch + // just use FSH partner list to infer it + // would require complex calculation for surfs + + if (fix_history) { + jtag = tag[j]; + n = numneigh[i]; + nn = dnum*n; + for (m = 0; m < npartner[i]; m++) + if (partner[i][m] == jtag) break; + if (m < npartner[i]) { + firsttouch[i][n] = 1; + memcpy(&firstshear[i][nn],shearpartner[i][m],dnumbytes); + } else { + firsttouch[i][n] = 0; + memcpy(&firstshear[i][nn],zeroes,dnumbytes); + } + } + + numneigh[i]++; + if (flip) i = j; + } + + // only add atom I to ilist if it has neighbors + // fix shear/history allows for this in pre_exchange_onesided() + + if (numneigh[i]) ilist[inum++] = i; } list->inum = inum; diff --git a/src/neighbor.h b/src/neighbor.h index b44e0fde00..ceec3249ef 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -174,6 +174,8 @@ class Neighbor : protected Pointers { int *glist; // lists to grow atom arrays every reneigh int *slist; // lists to grow stencil arrays every reneigh + double *zeroes; // vector of zeroes for shear history init + // USER-DPD package int *bins_ssa; // ptr to next atom in each bin used by SSA @@ -236,17 +238,12 @@ class Neighbor : protected Pointers { void full_bin_ghost(class NeighList *); void full_multi(class NeighList *); - void half_from_full_no_newton(class NeighList *); - void half_from_full_newton(class NeighList *); - void skip_from(class NeighList *); - void skip_from_granular(class NeighList *); - void skip_from_respa(class NeighList *); - void copy_from(class NeighList *); - void granular_nsq_no_newton(class NeighList *); void granular_nsq_newton(class NeighList *); + void granular_nsq_newton_onesided(class NeighList *); void granular_bin_no_newton(class NeighList *); void granular_bin_newton(class NeighList *); + void granular_bin_newton_onesided(class NeighList *); void granular_bin_newton_tri(class NeighList *); void respa_nsq_no_newton(class NeighList *); @@ -255,6 +252,15 @@ class Neighbor : protected Pointers { void respa_bin_newton(class NeighList *); void respa_bin_newton_tri(class NeighList *); + void half_from_full_no_newton(class NeighList *); + void half_from_full_newton(class NeighList *); + void skip_from(class NeighList *); + void skip_from_granular(class NeighList *); + void skip_from_granular_off2on(class NeighList *); + void skip_from_granular_off2on_onesided(class NeighList *); + void skip_from_respa(class NeighList *); + void copy_from(class NeighList *); + // include prototypes for multi-threaded neighbor lists // builds or their corresponding dummy versions diff --git a/src/pair.h b/src/pair.h index f70bca127a..cba2ee9bb8 100644 --- a/src/pair.h +++ b/src/pair.h @@ -170,7 +170,7 @@ class Pair : protected Pointers { virtual void write_data(FILE *) {} virtual void write_data_all(FILE *) {} - virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} + virtual int pack_forward_comm(int, int *, double *, int, int *) {printf("WHAT\n"); return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_forward_comm_kokkos(int, DAT::tdual_int_2d, int, DAT::tdual_xfloat_1d &, diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index ced2f6f34a..cf7207d7f1 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -491,7 +491,7 @@ void PairHybrid::init_style() for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style(); - // create skip lists for each pair neigh request + // create skip lists inside each pair neigh request // any kind of list can have its skip flag set at this stage for (i = 0; i < neighbor->nrequest; i++) { @@ -539,6 +539,7 @@ void PairHybrid::init_style() // if any skipping occurs // set request->skip and copy iskip and ijskip into request // else delete iskip and ijskip + // no skipping if pair style assigned to all type pairs skip = 0; for (itype = 1; itype <= ntypes; itype++) @@ -618,7 +619,8 @@ void PairHybrid::setup() } /* ---------------------------------------------------------------------- - combine sub-style neigh list requests and create new ones if needed + examine sub-style neigh list requests + create new parent requests if needed, to derive sub-style requests from ------------------------------------------------------------------------- */ void PairHybrid::modify_requests() @@ -626,29 +628,37 @@ void PairHybrid::modify_requests() int i,j; NeighRequest *irq,*jrq; - // loop over pair requests only - // if list is skip list and not copy, look for non-skip list of same kind - // if one exists, point at that one via otherlist - // else make new non-skip request of same kind and point at that one - // don't bother to set ID for new request, since pair hybrid ignores list - // only exception is half_from_full: - // ignore it, turn off skip, since it will derive from its skip parent - // after possible new request creation, unset skip flag and otherlist - // for these derived lists: granhistory, rRESPA inner/middle - // this prevents neighbor from treating them as skip lists - // copy list check is for pair style = hybrid/overlay - // which invokes this routine + // loop over pair requests only, including those added during looping + + int nrequest_original = neighbor->nrequest; for (i = 0; i < neighbor->nrequest; i++) { if (!neighbor->requests[i]->pair) continue; + // nothing more to do if this request: + // is not a skip list + // is a copy or half_from_full or granhistory list + // copy list setup is from pair style = hybrid/overlay + // which invokes this method at end of its modify_requests() + // if granhistory, turn off skip, since each gran sub-style + // its own history list, parent gran list does not have history + // if half_from_full, turn off skip, since it will derive + // from its full parent and its skip status + irq = neighbor->requests[i]; - if (irq->skip == 0 || irq->copy) continue; - if (irq->half_from_full) { + if (irq->skip == 0) continue; + if (irq->copy) continue; + if (irq->granhistory || irq->half_from_full) { irq->skip = 0; continue; } + // look for another list that matches via same_kind() and is not a skip list + // if one exists, point at that one via otherlist + // else make new parent request via copy_request() and point at that one + // new parent list is not a skip list + // parent does not need its ID set, since pair hybrid does not use it + for (j = 0; j < neighbor->nrequest; j++) { if (!neighbor->requests[j]->pair) continue; jrq = neighbor->requests[j]; @@ -662,11 +672,44 @@ void PairHybrid::modify_requests() irq->otherlist = newrequest; } - if (irq->granhistory || irq->respainner || irq->respamiddle) { + // for rRESPA inner/middle lists, + // which just created or set otherlist to parent: + // unset skip flag and otherlist + // this prevents neighbor from treating them as skip lists + + if (irq->respainner || irq->respamiddle) { irq->skip = 0; irq->otherlist = -1; } } + + // adjustments to newly added granular parent requests (gran = 1) + // parent newton = 2 if has children with granonesided = 0 and 1 + // else newton = 0 = setting of children + // parent gran onesided = 0 if has children with granonesided = 0 and 1 + // else onesided = setting of children + + for (i = nrequest_original; i < neighbor->nrequest; i++) { + if (!neighbor->requests[i]->pair) continue; + if (!neighbor->requests[i]->gran) continue; + irq = neighbor->requests[i]; + + int onesided = -1; + for (j = 0; j < nrequest_original; j++) { + if (!neighbor->requests[j]->pair) continue; + if (!neighbor->requests[j]->gran) continue; + jrq = neighbor->requests[j]; + + if (onesided < 0) onesided = jrq->granonesided; + else if (onesided != jrq->granonesided) onesided = 2; + if (onesided == 2) break; + } + + if (onesided == 2) { + irq->newton = 2; + irq->granonesided = 0; + } + } } /* ---------------------------------------------------------------------- -- GitLab