jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15274 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-07-09 20:41:07 +00:00
parent 18b0426c99
commit 6bc507cd6a
3 changed files with 182 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

136
src/neigh_gran.cpp
View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <string.h>
#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;
}
}

55
src/neigh_request.cpp
View File

@ -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;
}

87
src/neighbor.cpp
View File

@ -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;
}