Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F122174674
npair_half_bin_newton_ssa.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jul 16, 08:42
Size
10 KB
Mime Type
text/x-c
Expires
Fri, Jul 18, 08:42 (2 d)
Engine
blob
Format
Raw Data
Handle
27444678
Attached To
rLAMMPS lammps
npair_half_bin_newton_ssa.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
James Larentzos and Timothy I. Mattox (Engility Corporation)
------------------------------------------------------------------------- */
#include "npair_half_bin_newton_ssa.h"
#include "neighbor.h"
#include "nstencil_ssa.h"
#include "nbin_ssa.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "group.h"
#include "memory.h"
#include "my_page.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp)
{
ssa_maxPhaseCt = 0;
ssa_maxPhaseLen = 0;
ssa_phaseCt = 0;
ssa_phaseLen = NULL;
ssa_itemLoc = NULL;
ssa_itemLen = NULL;
}
/* ---------------------------------------------------------------------- */
NPairHalfBinNewtonSSA::~NPairHalfBinNewtonSSA()
{
ssa_maxPhaseCt = 0;
ssa_maxPhaseLen = 0;
ssa_phaseCt = 0;
memory->destroy(ssa_phaseLen);
memory->destroy(ssa_itemLoc);
memory->destroy(ssa_itemLen);
}
/* ----------------------------------------------------------------------
binned neighbor list construction with full Newton's 3rd law
for use by Shardlow Spliting Algorithm
each owned atom i checks its own bin and other bins in Newton stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void NPairHalfBinNewtonSSA::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *tag = atom->tag;
tagint *molecule = atom->molecule;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NStencilSSA *ns_ssa = dynamic_cast<NStencilSSA*>(ns);
if (!ns_ssa) error->one(FLERR, "NStencil wasn't a NStencilSSA object");
int *nstencil_ssa = &(ns_ssa->nstencil_ssa[0]);
int nstencil_full = ns_ssa->nstencil;
NBinSSA *nb_ssa = dynamic_cast<NBinSSA*>(nb);
if (!nb_ssa) error->one(FLERR, "NBin wasn't a NBinSSA object");
int *bins = nb_ssa->bins;
int *binhead = nb_ssa->binhead;
int *gairhead_ssa = &(nb_ssa->gairhead_ssa[0]);
int inum = 0;
int gnum = 0;
int xbin,ybin,zbin,xbin2,ybin2,zbin2;
int **stencilxyz = ns_ssa->stencilxyz;
int lbinxlo = nb_ssa->lbinxlo;
int lbinxhi = nb_ssa->lbinxhi;
int lbinylo = nb_ssa->lbinylo;
int lbinyhi = nb_ssa->lbinyhi;
int lbinzlo = nb_ssa->lbinzlo;
int lbinzhi = nb_ssa->lbinzhi;
int sx1 = ns_ssa->sx + 1;
int sy1 = ns_ssa->sy + 1;
int sz1 = ns_ssa->sz + 1;
ssa_phaseCt = sz1*sy1*sx1;
xbin = (lbinxhi - lbinxlo + sx1 - 1) / sx1 + 1;
ybin = (lbinyhi - lbinylo + sy1 - 1) / sy1 + 1;
zbin = (lbinzhi - lbinzlo + sz1 - 1) / sz1 + 1;
int phaseLenEstimate = xbin*ybin*zbin;
if (ssa_phaseCt > ssa_maxPhaseCt) {
ssa_maxPhaseCt = ssa_phaseCt;
ssa_maxPhaseLen = 0;
memory->destroy(ssa_phaseLen);
memory->destroy(ssa_itemLoc);
memory->destroy(ssa_itemLen);
memory->create(ssa_phaseLen,ssa_maxPhaseCt,"NPairHalfBinNewtonSSA:ssa_phaseLen");
}
if (phaseLenEstimate > ssa_maxPhaseLen) {
ssa_maxPhaseLen = phaseLenEstimate;
memory->destroy(ssa_itemLoc);
memory->destroy(ssa_itemLen);
memory->create(ssa_itemLoc,ssa_maxPhaseCt,ssa_maxPhaseLen,"NPairHalfBinNewtonSSA:ssa_itemLoc");
memory->create(ssa_itemLen,ssa_maxPhaseCt,ssa_maxPhaseLen,"NPairHalfBinNewtonSSA:ssa_itemLen");
}
ipage->reset();
int workPhase = 0;
// loop over bins with local atoms, storing half of the neighbors
for (int zoff = ns_ssa->sz; zoff >= 0; --zoff) {
for (int yoff = ns_ssa->sy; yoff >= 0; --yoff) {
for (int xoff = ns_ssa->sx; xoff >= 0; --xoff) {
int workItem = 0;
for (zbin = lbinzlo + zoff; zbin < lbinzhi; zbin += sz1) {
for (ybin = lbinylo + yoff - ns_ssa->sy; ybin < lbinyhi; ybin += sy1) {
for (xbin = lbinxlo + xoff - ns_ssa->sx; xbin < lbinxhi; xbin += sx1) {
if (workItem >= phaseLenEstimate) error->one(FLERR,"phaseLenEstimate was too small");
ssa_itemLoc[workPhase][workItem] = inum; // record where workItem starts in ilist
for (int subphase = 0; subphase < 4; subphase++) {
int s_ybin = ybin + ((subphase & 0x2) ? ns_ssa->sy : 0);
int s_xbin = xbin + ((subphase & 0x1) ? ns_ssa->sx : 0);
int ibin, ct;
if ((s_ybin < lbinylo) || (s_ybin >= lbinyhi)) continue;
if ((s_xbin < lbinxlo) || (s_xbin >= lbinxhi)) continue;
ibin = zbin*nb_ssa->mbiny*nb_ssa->mbinx
+ s_ybin*nb_ssa->mbinx
+ s_xbin;
for (i = binhead[ibin]; i >= 0; i = bins[i]) {
n = 0;
neighptr = ipage->vget();
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all local atoms in the current stencil "subphase"
for (k = nstencil_ssa[subphase]; k < nstencil_ssa[subphase+1]; k++) {
const int jbin = ibin+stencil[k];
if (jbin != ibin) j = binhead[jbin];
else j = bins[i]; // same bin as i, so start just past i in the bin
for (; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (n > 0) {
firstneigh[inum] = neighptr;
numneigh[inum] = n;
ilist[inum++] = i;
}
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
}
// record where workItem ends in ilist
ssa_itemLen[workPhase][workItem] = inum - ssa_itemLoc[workPhase][workItem];
if (ssa_itemLen[workPhase][workItem] > 0) workItem++;
}
}
}
// record where workPhase ends
ssa_phaseLen[workPhase++] = workItem;
}
}
}
if (ssa_phaseCt != workPhase) error->one(FLERR,"ssa_phaseCt was wrong");
list->AIRct_ssa[0] = list->inum = inum;
// loop over AIR ghost atoms, storing their local neighbors
// since these are ghosts, must check if stencil bin is out of bounds
for (int airnum = 1; airnum <= 7; airnum++) {
int locAIRct = 0;
for (i = gairhead_ssa[airnum]; i >= 0; i = bins[i]) {
n = 0;
neighptr = ipage->vget();
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
ibin = coord2bin(x[i],xbin,ybin,zbin);
// loop over AIR ghost atoms in all bins in "full" stencil
// Note: the non-AIR ghost atoms have already been filtered out
for (k = 0; k < nstencil_full; k++) {
xbin2 = xbin + stencilxyz[k][0];
ybin2 = ybin + stencilxyz[k][1];
zbin2 = zbin + stencilxyz[k][2];
// Skip it if this bin is outside the extent of local bins
if (xbin2 < lbinxlo || xbin2 >= lbinxhi ||
ybin2 < lbinylo || ybin2 >= lbinyhi ||
zbin2 < lbinzlo || zbin2 >= lbinzhi) continue;
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[j],nspecial[j],tag[i]);
else {
int jmol = molindex[j];
if (jmol >= 0) {
int jatom = molatom[j];
which = find_special(onemols[jmol]->special[jatom],
onemols[jmol]->nspecial[jatom],
tag[i] - (tag[j] - jatom - 1));
} else which = 0;
}
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (n > 0) {
firstneigh[inum + gnum] = neighptr;
numneigh[inum + gnum] = n;
ilist[inum + (gnum++)] = i;
++locAIRct;
}
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor (ghost) list overflow, boost neigh_modify one");
}
list->AIRct_ssa[airnum] = locAIRct;
}
list->gnum = gnum;
}
Event Timeline
Log In to Comment