Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F122206618
neighbor.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, 15:36
Size
11 KB
Mime Type
text/x-c
Expires
Fri, Jul 18, 15:36 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
27449960
Attached To
rLAMMPS lammps
neighbor.cpp
View Options
/*
//@HEADER
// ************************************************************************
//
// Kokkos v. 2.0
// Copyright (2014) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact H. Carter Edwards (hcedwar@sandia.gov)
//
// ************************************************************************
//@HEADER
*/
#include <system.h>
#include <cstdio>
#include <Kokkos_Core.hpp>
#define SMALL 1.0e-6
#define FACTOR 0.999
/* BinningFunctor puts atoms into bins of the simulation box
* Neighborlists are then created by checking only distances of atoms
* in adjacent bins. That makes neighborlist construction a O(N) operation.
*/
struct BinningFunctor {
typedef t_int_2d::execution_space execution_space;
System s;
int atoms_per_bin;
BinningFunctor(System _s): s(_s) {
atoms_per_bin = s.bins.dimension_1();
}
KOKKOS_INLINE_FUNCTION
void operator() (const int &i) const
{
const int ibin = coord2bin(s.d_x(i, 0), s.d_x(i, 1), s.d_x(i, 2));
const int ac = Kokkos::atomic_fetch_add(&s.bincount[ibin], 1);
if(ac < atoms_per_bin) {
s.bins(ibin, ac) = i;
} else if(s.d_resize(0) < ac) {
s.d_resize(0) = ac;
}
}
KOKKOS_INLINE_FUNCTION
int coord2bin(double x, double y, double z) const
{
int ix, iy, iz;
if(x >= s.box.xprd)
ix = (int)((x - s.box.xprd) * s.bininvx) + s.nbinx - s.mbinxlo;
else if(x >= 0.0)
ix = (int)(x * s.bininvx) - s.mbinxlo;
else
ix = (int)(x * s.bininvx) - s.mbinxlo - 1;
if(y >= s.box.yprd)
iy = (int)((y - s.box.yprd) * s.bininvy) + s.nbiny - s.mbinylo;
else if(y >= 0.0)
iy = (int)(y * s.bininvy) - s.mbinylo;
else
iy = (int)(y * s.bininvy) - s.mbinylo - 1;
if(z >= s.box.zprd)
iz = (int)((z - s.box.zprd) * s.bininvz) + s.nbinz - s.mbinzlo;
else if(z >= 0.0)
iz = (int)(z * s.bininvz) - s.mbinzlo;
else
iz = (int)(z * s.bininvz) - s.mbinzlo - 1;
return (iz * s.mbiny * s.mbinx + iy * s.mbinx + ix + 1);
}
};
/* Build the actual neighborlist*/
struct BuildFunctor {
typedef t_int_2d::execution_space execution_space;
System s;
int maxneighs;
BuildFunctor(System _s): s(_s) {
maxneighs = s.neighbors.dimension_1();
}
KOKKOS_INLINE_FUNCTION
void operator() (const int &i) const
{
int n = 0;
const t_int_1d_const_um bincount_c = s.bincount;
const double xtmp = s.d_x(i, 0);
const double ytmp = s.d_x(i, 1);
const double ztmp = s.d_x(i, 2);
const int ibin = coord2bin(xtmp, ytmp, ztmp);
// loop over all bins in neighborhood (includes ibin)
for(int k = 0; k < s.nstencil; k++) {
const int jbin = ibin + s.d_stencil[k];
// get subview of jbin
const t_int_1d_const_um loc_bin =
Kokkos::subview(s.bins,jbin,Kokkos::ALL());
if(ibin == jbin)
for(int m = 0; m < bincount_c[jbin]; m++) {
const int j = loc_bin[m];
//for same bin as atom i skip j if i==j
if (j == i) continue;
const double delx = xtmp - s.d_x(j, 0);
const double dely = ytmp - s.d_x(j, 1);
const double delz = ztmp - s.d_x(j, 2);
const double rsq = delx * delx + dely * dely + delz * delz;
if(rsq <= s.neigh_cutsq && n<maxneighs) s.neighbors(i,n++) = j;
}
else {
for(int m = 0; m < bincount_c[jbin]; m++) {
const int j = loc_bin[m];
const double delx = xtmp - s.d_x(j, 0);
const double dely = ytmp - s.d_x(j, 1);
const double delz = ztmp - s.d_x(j, 2);
const double rsq = delx * delx + dely * dely + delz * delz;
if(rsq <= s.neigh_cutsq && n<maxneighs) s.neighbors(i,n++) = j;
}
}
}
s.numneigh[i] = n;
if(n >= maxneighs) {
if(n >= s.d_resize(0)) s.d_resize(0) = n;
}
}
KOKKOS_INLINE_FUNCTION
int coord2bin(double x, double y, double z) const
{
int ix, iy, iz;
if(x >= s.box.xprd)
ix = (int)((x - s.box.xprd) * s.bininvx) + s.nbinx - s.mbinxlo;
else if(x >= 0.0)
ix = (int)(x * s.bininvx) - s.mbinxlo;
else
ix = (int)(x * s.bininvx) - s.mbinxlo - 1;
if(y >= s.box.yprd)
iy = (int)((y - s.box.yprd) * s.bininvy) + s.nbiny - s.mbinylo;
else if(y >= 0.0)
iy = (int)(y * s.bininvy) - s.mbinylo;
else
iy = (int)(y * s.bininvy) - s.mbinylo - 1;
if(z >= s.box.zprd)
iz = (int)((z - s.box.zprd) * s.bininvz) + s.nbinz - s.mbinzlo;
else if(z >= 0.0)
iz = (int)(z * s.bininvz) - s.mbinzlo;
else
iz = (int)(z * s.bininvz) - s.mbinzlo - 1;
return (iz * s.mbiny * s.mbinx + iy * s.mbinx + ix + 1);
}
};
/* Reset an array to zero */
struct MemsetZeroFunctor {
typedef t_x_array::execution_space execution_space ;
void* ptr;
KOKKOS_INLINE_FUNCTION void operator()(const int i) const {
((int*)ptr)[i] = 0;
}
};
/* Calculate distance of two bins */
double bindist(System &s, int i, int j, int k)
{
double delx, dely, delz;
if(i > 0)
delx = (i - 1) * s.binsizex;
else if(i == 0)
delx = 0.0;
else
delx = (i + 1) * s.binsizex;
if(j > 0)
dely = (j - 1) * s.binsizey;
else if(j == 0)
dely = 0.0;
else
dely = (j + 1) * s.binsizey;
if(k > 0)
delz = (k - 1) * s.binsizez;
else if(k == 0)
delz = 0.0;
else
delz = (k + 1) * s.binsizez;
return (delx * delx + dely * dely + delz * delz);
}
/* Setup the neighborlist construction
* Determine binsizes, a stencil for defining adjacency, etc.
*/
void neigh_setup(System &s) {
s.neigh_cutsq = s.neigh_cut * s.neigh_cut;
/*
c bins must evenly divide into box size,
c becoming larger than cutneigh if necessary
c binsize = 1/2 of cutoff is near optimal
if (flag == 0) {
nbinx = 2.0 * xprd / cutneigh;
nbiny = 2.0 * yprd / cutneigh;
nbinz = 2.0 * zprd / cutneigh;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
}
*/
s.binsizex = s.box.xprd / s.nbinx;
s.binsizey = s.box.yprd / s.nbiny;
s.binsizez = s.box.zprd / s.nbinz;
s.bininvx = 1.0 / s.binsizex;
s.bininvy = 1.0 / s.binsizey;
s.bininvz = 1.0 / s.binsizez;
double coord = s.box.xlo - s.neigh_cut - SMALL * s.box.xprd;
s.mbinxlo = static_cast<int>(coord * s.bininvx);
if(coord < 0.0) s.mbinxlo = s.mbinxlo - 1;
coord = s.box.xhi + s.neigh_cut + SMALL * s.box.xprd;
int mbinxhi = static_cast<int>(coord * s.bininvx);
coord = s.box.ylo - s.neigh_cut - SMALL * s.box.yprd;
s.mbinylo = static_cast<int>(coord * s.bininvy);
if(coord < 0.0) s.mbinylo = s.mbinylo - 1;
coord = s.box.yhi + s.neigh_cut + SMALL * s.box.yprd;
int mbinyhi = static_cast<int>(coord * s.bininvy);
coord = s.box.zlo - s.neigh_cut - SMALL * s.box.zprd;
s.mbinzlo = static_cast<int>(coord * s.bininvz);
if(coord < 0.0) s.mbinzlo = s.mbinzlo - 1;
coord = s.box.zhi + s.neigh_cut + SMALL * s.box.zprd;
int mbinzhi = static_cast<int>(coord * s.bininvz);
/* extend bins by 1 in each direction to insure stencil coverage */
s.mbinxlo = s.mbinxlo - 1;
mbinxhi = mbinxhi + 1;
s.mbinx = mbinxhi - s.mbinxlo + 1;
s.mbinylo = s.mbinylo - 1;
mbinyhi = mbinyhi + 1;
s.mbiny = mbinyhi - s.mbinylo + 1;
s.mbinzlo = s.mbinzlo - 1;
mbinzhi = mbinzhi + 1;
s.mbinz = mbinzhi - s.mbinzlo + 1;
/*
compute bin stencil of all bins whose closest corner to central bin
is within neighbor cutoff
for partial Newton (newton = 0),
stencil is all surrounding bins including self
for full Newton (newton = 1),
stencil is bins to the "upper right" of central bin, does NOT include self
next(xyz) = how far the stencil could possibly extend
factor < 1.0 for special case of LJ benchmark so code will create
correct-size stencil when there are 3 bins for every 5 lattice spacings
*/
int nextx = static_cast<int>(s.neigh_cut * s.bininvx);
if(nextx * s.binsizex < FACTOR * s.neigh_cut) nextx++;
int nexty = static_cast<int>(s.neigh_cut * s.bininvy);
if(nexty * s.binsizey < FACTOR * s.neigh_cut) nexty++;
int nextz = static_cast<int>(s.neigh_cut * s.bininvz);
if(nextz * s.binsizez < FACTOR * s.neigh_cut) nextz++;
int nmax = (2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1);
s.d_stencil = t_int_1d("stencil", nmax);
s.h_stencil = Kokkos::create_mirror_view(s.d_stencil);
s.nstencil = 0;
int kstart = -nextz;
for(int k = kstart; k <= nextz; k++) {
for(int j = -nexty; j <= nexty; j++) {
for(int i = -nextx; i <= nextx; i++) {
if(bindist(s,i, j, k) < s.neigh_cutsq) {
s.h_stencil(s.nstencil++) = k * s.mbiny * s.mbinx + j * s.mbinx + i;
}
}
}
}
/* Allocate neighbor arrays */
Kokkos::deep_copy(s.d_stencil, s.h_stencil);
s.mbins = s.mbinx * s.mbiny * s.mbinz;
s.bincount = t_int_1d("bincount", s.mbins);
s.bins = t_int_2d("bins", s.mbins, 8);
s.neighbors = t_neighbors("neighbors",s.natoms,80);
s.numneigh = t_int_1d("numneigh",s.natoms);
s.d_resize = t_int_scalar("resize");
s.h_resize = Kokkos::create_mirror_view(s.d_resize);
}
/* Build the neighborlist
* This is a try and rerun algorithm for handling the case where the bins array
* and the neighbors array are not big enough. So if one is too small, it will
* reallocate and rerun the binnind algorithm or the neighborlist construction.
*/
void neigh_build(System &s) {
/* Binning of atoms */
s.h_resize(0) = 1;
while(s.h_resize(0) > 0) {
s.h_resize(0) = 0;
Kokkos::deep_copy(s.d_resize, s.h_resize);
MemsetZeroFunctor f_zero;
f_zero.ptr = (void*) s.bincount.ptr_on_device();
Kokkos::parallel_for(s.mbins, f_zero);
execution_space::fence();
BinningFunctor f(s);
Kokkos::parallel_for(s.natoms, f);
execution_space::fence();
/* Check if bins was large enough, if nor reallocated and rerun */
deep_copy(s.h_resize, s.d_resize);
if(s.h_resize(0)) {
int atoms_per_bin = s.h_resize(0)+2;
s.bins = t_int_2d("bins", s.mbins, atoms_per_bin);
}
}
/* Neighborlist construction */
s.h_resize(0) = 1;
while(s.h_resize(0)) {
s.h_resize(0) = 0;
Kokkos::deep_copy(s.d_resize, s.h_resize);
BuildFunctor f(s);
Kokkos::parallel_for(s.nlocal, f);
execution_space::fence();
/* Check if neighbors was large enough, if nor reallocated and rerun */
deep_copy(s.h_resize, s.d_resize);
if(s.h_resize(0)) {
int maxneighs = s.h_resize(0) * 1.2;
s.neighbors = t_neighbors("neighbors", s.natoms, maxneighs);
}
}
}
Event Timeline
Log In to Comment