lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2013-07-03 15:50:00 +00:00
parent d469f8f067
commit 5b5798769f
12 changed files with 147 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/GRANULAR/fix_pour.h
View File

@ -30,6 +30,7 @@ class FixPour : public Fix {
friend class PairGranHooke;
friend class PairGranHookeOMP;
friend class PairGranHookeHistory;
friend class PairGranHookeHistory2;
friend class PairGranHookeHistoryOMP;
friend class PairGranHookeCuda;
@ -61,7 +62,6 @@ class FixPour : public Fix {
int *recvcounts,*displs;
int nfreq,nfirst,ninserted,nper;
double lo_current,hi_current;
class FixShearHistory *fix_history;
class RanPark *random;
int overlap(int);

8
src/comm.cpp
View File

@ -1024,8 +1024,7 @@ void Comm::borders()
// pack up list of border atoms
if (nsend*size_border > maxsend)
grow_send(nsend*size_border,0);
if (nsend*size_border > maxsend) grow_send(nsend*size_border,0);
if (ghost_velocity)
n = avec->pack_border_vel(nsend,sendlist[iswap],buf_send,
pbc_flag[iswap],pbc[iswap]);
@ -1466,10 +1465,7 @@ void Comm::forward_comm_array(int n, double **array)
MPI_Request request;
MPI_Status status;
// check that buf_send and buf_recv are big enough
// NOTE: check that buf_send and buf_recv are big enough
for (iswap = 0; iswap < nswap; iswap++) {

133
src/fix_shear_history.cpp
View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "fix_shear_history.h"
@ -27,8 +28,6 @@
using namespace LAMMPS_NS;
using namespace FixConst;
#define MAXTOUCH 15
/* ---------------------------------------------------------------------- */
FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
@ -47,10 +46,15 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
atom->add_callback(1);
ipage = NULL;
dpage = NULL;
pgsize = oneatom = 0;
// initialize npartner to 0 so neighbor list creation is OK the 1st time
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) npartner[i] = 0;
maxtouch = 0;
}
/* ---------------------------------------------------------------------- */
@ -65,8 +69,10 @@ FixShearHistory::~FixShearHistory()
// delete locally stored arrays
memory->destroy(npartner);
memory->destroy(partner);
memory->destroy(shearpartner);
memory->sfree(partner);
memory->sfree(shearpartner);
delete ipage;
delete dpage;
}
/* ---------------------------------------------------------------------- */
@ -89,6 +95,21 @@ void FixShearHistory::init()
int dim;
computeflag = (int *) pair->extract("computeflag",dim);
// create pages if first time or if neighbor pgsize/oneatom has changed
// note that latter could cause shear history info to be discarded
int create = 0;
if (ipage == NULL) create = 1;
if (pgsize != neighbor->pgsize) create = 1;
if (oneatom != neighbor->oneatom) create = 1;
if (create) {
pgsize = neighbor->pgsize;
oneatom = neighbor->oneatom;
ipage = new MyPage<int>(oneatom,pgsize);
dpage = new MyPage<double[3]>(oneatom,pgsize);
}
}
/* ----------------------------------------------------------------------
@ -119,17 +140,22 @@ void FixShearHistory::setup_pre_exchange()
void FixShearHistory::pre_exchange()
{
int i,j,ii,jj,m,inum,jnum;
int i,j,ii,jj,m,n,inum,jnum;
int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch;
double *shear,*allshear,**firstshear;
// zero npartner for all current atoms
// clear 2 page data structures
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) npartner[i] = 0;
// copy shear info from neighbor list atoms to atom arrays
ipage->reset();
dpage->reset();
// 1st loop over neighbor list
// calculate npartner for each owned atom
int *tag = atom->tag;
NeighList *list = pair->list;
@ -140,6 +166,39 @@ void FixShearHistory::pre_exchange()
firsttouch = list->listgranhistory->firstneigh;
firstshear = list->listgranhistory->firstdouble;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
touch = firsttouch[i];
for (jj = 0; jj < jnum; jj++) {
if (touch[jj]) {
npartner[i]++;
j = jlist[jj];
j &= NEIGHMASK;
if (j < nlocal) npartner[j]++;
}
}
}
// get page chunks to store atom IDs and shear history for my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
n = npartner[i];
partner[i] = ipage->get(n);
shearpartner[i] = dpage->get(n);
if (partner[i] == NULL || shearpartner[i] == NULL)
error->one(FLERR,"Shear history overflow, boost neigh_modify one");
}
// 2nd loop over neighbor list
// store atom IDs and shear history for my atoms
// re-zero npartner to use as counter for all my atoms
for (i = 0; i < nlocal; i++) npartner[i] = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
@ -152,37 +211,28 @@ void FixShearHistory::pre_exchange()
shear = &allshear[3*jj];
j = jlist[jj];
j &= NEIGHMASK;
if (npartner[i] < MAXTOUCH) {
m = npartner[i];
partner[i][m] = tag[j];
shearpartner[i][m][0] = shear[0];
shearpartner[i][m][1] = shear[1];
shearpartner[i][m][2] = shear[2];
}
npartner[i]++;
if (j < nlocal) {
if (npartner[j] < MAXTOUCH) {
m = npartner[j];
partner[j][m] = tag[i];
shearpartner[j][m][0] = -shear[0];
shearpartner[j][m][1] = -shear[1];
shearpartner[j][m][2] = -shear[2];
}
npartner[j]++;
}
}
}
}
// test for too many touching neighbors
// set maxtouch = max # of partners of any owned atom
int flag = 0;
for (i = 0; i < nlocal; i++)
if (npartner[i] >= MAXTOUCH) flag = 1;
int flag_all;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
if (flag_all)
error->all(FLERR,"Too many touching neighbors - boost MAXTOUCH");
maxtouch = 0;
for (i = 0; i < nlocal; i++) maxtouch = MAX(maxtouch,npartner[i]);
}
/* ---------------------------------------------------------------------- */
@ -208,8 +258,10 @@ double FixShearHistory::memory_usage()
{
int nmax = atom->nmax;
double bytes = nmax * sizeof(int);
bytes += nmax*MAXTOUCH * sizeof(int);
bytes += nmax*MAXTOUCH*3 * sizeof(double);
bytes = nmax * sizeof(int *);
bytes = nmax * sizeof(double *);
bytes += ipage->ndatum * sizeof(int);
bytes += dpage->ndatum * sizeof(double[3]);
return bytes;
}
@ -220,8 +272,12 @@ double FixShearHistory::memory_usage()
void FixShearHistory::grow_arrays(int nmax)
{
memory->grow(npartner,nmax,"shear_history:npartner");
memory->grow(partner,nmax,MAXTOUCH,"shear_history:partner");
memory->grow(shearpartner,nmax,MAXTOUCH,3,"shear_history:shearpartner");
partner = (int **) memory->srealloc(partner,nmax*sizeof(int *),
"shear_history:partner");
typedef double (*sptype)[3];
shearpartner = (sptype *)
memory->srealloc(shearpartner,nmax*sizeof(sptype),
"shear_history:shearpartner");
}
/* ----------------------------------------------------------------------
@ -230,13 +286,15 @@ void FixShearHistory::grow_arrays(int nmax)
void FixShearHistory::copy_arrays(int i, int j, int delflag)
{
// just copy pointers for partner and shearpartner
// b/c can't overwrite chunk allocation inside ipage,dpage
// incoming atoms in unpack_exchange just grab new chunks
// so are orphaning chunks for migrating atoms
// OK, b/c will reset ipage,dpage on next reneighboring
npartner[j] = npartner[i];
for (int m = 0; m < npartner[j]; m++) {
partner[j][m] = partner[i][m];
shearpartner[j][m][0] = shearpartner[i][m][0];
shearpartner[j][m][1] = shearpartner[i][m][1];
shearpartner[j][m][2] = shearpartner[i][m][2];
}
partner[j] = partner[i];
shearpartner[j] = shearpartner[i];
}
/* ----------------------------------------------------------------------
@ -254,6 +312,9 @@ void FixShearHistory::set_arrays(int i)
int FixShearHistory::pack_exchange(int i, double *buf)
{
// NOTE: how do I know comm buf is big enough if extreme # of touching neighs
// Comm::BUFEXTRA may need to be increased
int m = 0;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
@ -271,8 +332,13 @@ int FixShearHistory::pack_exchange(int i, double *buf)
int FixShearHistory::unpack_exchange(int nlocal, double *buf)
{
// allocate new chunks from ipage,dpage for incoming values
int m = 0;
npartner[nlocal] = static_cast<int> (buf[m++]);
maxtouch = MAX(maxtouch,npartner[nlocal]);
partner[nlocal] = ipage->get(npartner[nlocal]);
shearpartner[nlocal] = dpage->get(npartner[nlocal]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<int> (buf[m++]);
shearpartner[nlocal][n][0] = buf[m++];
@ -314,7 +380,12 @@ void FixShearHistory::unpack_restart(int nlocal, int nth)
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
// allocate new chunks from ipage,dpage for incoming values
npartner[nlocal] = static_cast<int> (extra[nlocal][m++]);
maxtouch = MAX(maxtouch,npartner[nlocal]);
partner[nlocal] = ipage->get(npartner[nlocal]);
shearpartner[nlocal] = dpage->get(npartner[nlocal]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<int> (extra[nlocal][m++]);
shearpartner[nlocal][n][0] = extra[nlocal][m++];
@ -329,7 +400,11 @@ void FixShearHistory::unpack_restart(int nlocal, int nth)
int FixShearHistory::maxsize_restart()
{
return 4*MAXTOUCH + 2;
// maxtouch_all = max touching partners across all procs
int maxtouch_all;
MPI_Allreduce(&maxtouch,&maxtouch_all,1,MPI_INT,MPI_MAX,world);
return 4*maxtouch_all + 2;
}
/* ----------------------------------------------------------------------

9
src/fix_shear_history.h
View File

@ -21,13 +21,13 @@ FixStyle(SHEAR_HISTORY,FixShearHistory)
#define LMP_FIX_SHEAR_HISTORY_H
#include "fix.h"
#include "my_page.h"
namespace LAMMPS_NS {
class FixShearHistory : public Fix {
friend class Neighbor;
friend class PairGranHookeHistory;
friend class FixPour;
public:
FixShearHistory(class LAMMPS *, int, char **);
@ -53,10 +53,15 @@ class FixShearHistory : public Fix {
protected:
int *npartner; // # of touching partners of each atom
int **partner; // tags for the partners
double ***shearpartner; // 3 shear values with the partner
double (**shearpartner)[3]; // 3 shear values with the partner
int maxtouch; // max # of touching partners for my atoms
class Pair *pair;
int *computeflag; // computeflag in PairGranHookeHistory
int pgsize,oneatom; // copy of settings in Neighbor
MyPage<int> *ipage; // pages of partner atom IDs
MyPage<double[3]> *dpage; // pages of shear history with partners
};
}

3
src/modify.cpp
View File

@ -863,7 +863,8 @@ void Modify::modify_compute(int narg, char **arg)
int icompute;
for (icompute = 0; icompute < ncompute; icompute++)
if (strcmp(arg[0],compute[icompute]->id) == 0) break;
if (icompute == ncompute) error->all(FLERR,"Could not find compute_modify ID");
if (icompute == ncompute)
error->all(FLERR,"Could not find compute_modify ID");
compute[icompute]->modify_params(narg-1,&arg[1]);
}

4
src/neigh_gran.cpp
View File

@ -38,7 +38,7 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
NeighList *listgranhistory;
int *npartner,**partner;
double ***shearpartner;
double (**shearpartner)[3];
int **firsttouch;
double **firstshear;
int **pages_touch;
@ -280,7 +280,7 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
NeighList *listgranhistory;
int *npartner,**partner;
double ***shearpartner;
double (**shearpartner)[3];
int **firsttouch;
double **firstshear;
int **pages_touch;

3
src/neigh_list.cpp
View File

@ -28,10 +28,11 @@ enum{NSQ,BIN,MULTI}; // also in neighbor.cpp
/* ---------------------------------------------------------------------- */
NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp)
NeighList::NeighList(LAMMPS *lmp, int size, int onesize) : Pointers(lmp)
{
maxatoms = 0;
pgsize = size;
oneatom = onesize;
inum = gnum = 0;
ilist = NULL;

6
src/neigh_list.h
View File

@ -15,6 +15,7 @@
#define LMP_NEIGH_LIST_H
#include "pointers.h"
#include "my_page.h"
namespace LAMMPS_NS {
@ -39,11 +40,14 @@ class NeighList : protected Pointers {
double **firstdouble; // ptr to 1st J double value of each I atom
int pgsize; // size of each page
int oneatom; // max size for one atom
int maxpage; // # of pages currently allocated
int **pages; // neighbor list pages for ints
double **dpages; // neighbor list pages for doubles
int dnum; // # of doubles for each pair (0 if none)
MyPage<int> *page;
// atom types to skip when building list
// iskip,ijskip are just ptrs to corresponding request
@ -76,7 +80,7 @@ class NeighList : protected Pointers {
class CudaNeighList *cuda_list; // CUDA neighbor list
NeighList(class LAMMPS *, int);
NeighList(class LAMMPS *, int, int);
~NeighList();
void grow(int); // grow maxlocal
void stencil_allocate(int, int); // allocate stencil arrays

3
src/neigh_stencil.cpp
View File

@ -23,8 +23,7 @@ using namespace LAMMPS_NS;
sx,sy,sz = bin bounds = furthest the stencil could possibly extend
3d creates xyz stencil, 2d creates xy stencil
for half list with newton off:
stencil is all surrounding bins
stencil includes self
stencil is all surrounding bins including self
regardless of triclinic
for half list with newton on:
stencil is bins to the "upper right" of central bin

10
src/neighbor.cpp
View File

@ -119,6 +119,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
old_style = BIN;
old_triclinic = 0;
old_pgsize = pgsize;
old_oneatom = oneatom;
old_nrequest = 0;
old_requests = NULL;
@ -416,12 +418,14 @@ void Neighbor::init()
// test if pairwise lists need to be re-created
// no need to re-create if:
// neigh style and triclinic has not changed and
// neigh style, triclinic, pgsize, oneatom have not changed
// current requests = old requests
int same = 1;
if (style != old_style) same = 0;
if (triclinic != old_triclinic) same = 0;
if (pgsize != old_pgsize) same = 0;
if (oneatom != old_oneatom) same = 0;
if (nrequest != old_nrequest) same = 0;
else
for (i = 0; i < nrequest; i++)
@ -452,7 +456,7 @@ void Neighbor::init()
// pass list ptr back to requestor (except for Command class)
for (i = 0; i < nlist; i++) {
lists[i] = new NeighList(lmp,pgsize);
lists[i] = new NeighList(lmp,pgsize,oneatom);
lists[i]->index = i;
lists[i]->dnum = requests[i]->dnum;
@ -1676,10 +1680,12 @@ void Neighbor::modify_params(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"page") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
old_pgsize = pgsize;
pgsize = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"one") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
old_oneatom = oneatom;
oneatom = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"binsize") == 0) {

2
src/neighbor.h
View File

@ -49,6 +49,8 @@ class Neighbor : protected Pointers {
int old_style; // previous run info to avoid
int old_nrequest; // re-creation of pairwise neighbor lists
int old_triclinic;
int old_pgsize;
int old_oneatom;
class NeighRequest **old_requests;
int nlist; // pairwise neighbor lists

1
src/random_mars.cpp
View File

@ -12,6 +12,7 @@
------------------------------------------------------------------------- */
// Marsaglia random number generator
// see RANMAR in F James, Comp Phys Comm, 60, 329 (1990)
#include "math.h"
#include "random_mars.h"