From 951e7c916afe83e4df2bb8cd143327ce29ce7c54 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 17 Jun 2016 17:39:55 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15187 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/pair_gran_hertz_history.cpp | 19 +- src/GRANULAR/pair_gran_hooke.cpp | 15 +- src/GRANULAR/pair_gran_hooke_history.cpp | 28 ++- src/comm.h | 1 + src/comm_brick.cpp | 49 ++++- src/comm_brick.h | 2 + src/comm_tiled.cpp | 12 ++ src/comm_tiled.h | 2 + src/fix.h | 1 + src/fix_shear_history.cpp | 262 +++++++++++++++++++++-- src/fix_shear_history.h | 12 ++ src/neigh_gran.cpp | 257 ++++++++++++++++++++-- 12 files changed, 613 insertions(+), 47 deletions(-) diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e14dc7110f..0ba4113c99 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -89,6 +89,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; @@ -254,7 +255,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; - if (j < nlocal) { + if (newton_pair || j < nlocal) { f[j][0] -= fx; f[j][1] -= fy; f[j][2] -= fz; @@ -268,6 +269,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag) } } } + + if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- @@ -321,7 +324,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, if (rsq >= radsum*radsum) { fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; return 0.0; } @@ -441,12 +444,22 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, } else fs1 = fs2 = fs3 = fs = 0.0; } - // set all forces and return no energy + // set force and return no energy fforce = ccel; + + // set single_extra quantities + svector[0] = fs1; svector[1] = fs2; svector[2] = fs3; svector[3] = fs; + svector[4] = vn1; + svector[5] = vn2; + svector[6] = vn3; + svector[7] = vt1; + svector[8] = vt2; + svector[9] = vt3; + return 0.0; } diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 9ff23a7553..66e1c3dd7a 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -239,7 +239,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, if (rsq >= radsum*radsum) { fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; return 0.0; } @@ -318,14 +318,25 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; - // set all forces and return no energy + // set force and return no energy fforce = ccel; + + + // set single_extra quantities + svector[0] = -ft*vtr1; svector[1] = -ft*vtr2; svector[2] = -ft*vtr3; svector[3] = sqrt(svector[0]*svector[0] + svector[1]*svector[1] + svector[2]*svector[2]); + svector[4] = vn1; + svector[5] = vn2; + svector[6] = vn3; + svector[7] = vt1; + svector[8] = vt2; + svector[9] = vt3; + return 0.0; } diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index b8b0381d58..c1a96e56b1 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -46,8 +46,8 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) history = 1; fix_history = NULL; - single_extra = 4; - svector = new double[4]; + single_extra = 10; + svector = new double[10]; neighprev = 0; @@ -131,6 +131,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; @@ -295,7 +296,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; - if (j < nlocal) { + if (newton_pair || j < nlocal) { f[j][0] -= fx; f[j][1] -= fy; f[j][2] -= fz; @@ -309,6 +310,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag) } } } + + if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- @@ -413,13 +416,8 @@ void PairGranHookeHistory::init_style() dt = update->dt; // if shear history is stored: - // check if newton flag is valid // if first init, create Fix needed for storing shear history - if (history && force->newton_pair == 1) - error->all(FLERR, - "Pair granular with shear history requires newton pair off"); - if (history && fix_history == NULL) { char **fixarg = new char*[3]; fixarg[0] = (char *) "SHEAR_HISTORY"; @@ -624,7 +622,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, if (rsq >= radsum*radsum) { fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; return 0.0; } @@ -741,13 +739,23 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, } else fs1 = fs2 = fs3 = fs = 0.0; } - // set all forces and return no energy + // set force and return no energy fforce = ccel; + + // set single_extra quantities + svector[0] = fs1; svector[1] = fs2; svector[2] = fs3; svector[3] = fs; + svector[4] = vn1; + svector[5] = vn2; + svector[6] = vn3; + svector[7] = vt1; + svector[8] = vt2; + svector[9] = vt3; + return 0.0; } diff --git a/src/comm.h b/src/comm.h index cac8b2fc50..15b42111d8 100644 --- a/src/comm.h +++ b/src/comm.h @@ -80,6 +80,7 @@ class Comm : protected Pointers { virtual void reverse_comm_pair(class Pair *) = 0; virtual void forward_comm_fix(class Fix *, int size=0) = 0; virtual void reverse_comm_fix(class Fix *, int size=0) = 0; + virtual void reverse_comm_fix_variable(class Fix *) = 0; virtual void forward_comm_compute(class Compute *) = 0; virtual void reverse_comm_compute(class Compute *) = 0; virtual void forward_comm_dump(class Dump *) = 0; diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index 0aea55b82c..c44b41928f 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -1037,6 +1037,50 @@ void CommBrick::reverse_comm_fix(Fix *fix, int size) } } +/* ---------------------------------------------------------------------- + reverse communication invoked by a Fix with variable size data + query fix for pack size to insure buf_send is big enough + handshake sizes before each Irecv/Send to insure buf_recv is big enough +------------------------------------------------------------------------- */ + +void CommBrick::reverse_comm_fix_variable(Fix *fix) +{ + int iswap,nsend,nrecv; + double *buf; + MPI_Request request; + + for (iswap = nswap-1; iswap >= 0; iswap--) { + + // pack buffer + + nsend = fix->pack_reverse_comm_size(recvnum[iswap],firstrecv[iswap]); + if (nsend > maxsend) grow_send(nsend,0); + nsend = fix->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send); + + // exchange with another proc + // if self, set recv buffer to send buffer + + if (sendproc[iswap] != me) { + MPI_Sendrecv(&nsend,1,MPI_INT,recvproc[iswap],0, + &nrecv,1,MPI_INT,sendproc[iswap],0,world, + MPI_STATUS_IGNORE); + if (sendnum[iswap]) { + if (nrecv > maxrecv) grow_recv(nrecv); + MPI_Irecv(buf_recv,maxrecv,MPI_DOUBLE,sendproc[iswap],0, + world,&request); + } + if (recvnum[iswap]) + MPI_Send(buf_send,nsend,MPI_DOUBLE,recvproc[iswap],0,world); + if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE); + buf = buf_recv; + } else buf = buf_send; + + // unpack buffer + + fix->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf); + } +} + /* ---------------------------------------------------------------------- forward communication invoked by a Compute nsize used only to set recv buffer limit @@ -1055,7 +1099,7 @@ void CommBrick::forward_comm_compute(Compute *compute) // pack buffer n = compute->pack_forward_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flag[iswap],pbc[iswap]); + buf_send,pbc_flag[iswap],pbc[iswap]); // exchange with another proc // if self, set recv buffer to send buffer @@ -1276,7 +1320,8 @@ int CommBrick::exchange_variable(int n, double *inbuf, double *&outbuf) nrecv += nrecv1; if (procgrid[dim] > 2) { MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0, - &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,MPI_STATUS_IGNORE); + &nrecv2,1,MPI_INT,procneigh[dim][0],0,world, + MPI_STATUS_IGNORE); nrecv += nrecv2; } else nrecv2 = 0; diff --git a/src/comm_brick.h b/src/comm_brick.h index ddcf8e6fbd..461ee84b40 100644 --- a/src/comm_brick.h +++ b/src/comm_brick.h @@ -37,6 +37,8 @@ class CommBrick : public Comm { // forward comm from a Fix virtual void reverse_comm_fix(class Fix *, int size=0); // reverse comm from a Fix + virtual void reverse_comm_fix_variable(class Fix *); + // variable size reverse comm from a Fix virtual void forward_comm_compute(class Compute *); // forward from a Compute virtual void reverse_comm_compute(class Compute *); // reverse from a Compute virtual void forward_comm_dump(class Dump *); // forward comm from a Dump diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index 058de915ba..c07a8e5fba 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -1180,6 +1180,18 @@ void CommTiled::reverse_comm_fix(Fix *fix, int size) } } +/* ---------------------------------------------------------------------- + reverse communication invoked by a Fix with variable size data + query fix for all pack sizes to insure buf_send is big enough + handshake sizes before irregular comm to insure buf_recv is big enough + NOTE: how to setup one big buf recv with correct offsets ?? +------------------------------------------------------------------------- */ + +void CommTiled::reverse_comm_fix_variable(Fix *fix) +{ + error->all(FLERR,"Reverse comm fix variable not yet supported by CommTiled"); +} + /* ---------------------------------------------------------------------- forward communication invoked by a Compute nsize used only to set recv buffer limit diff --git a/src/comm_tiled.h b/src/comm_tiled.h index dd86ecfdba..5f08c8a703 100644 --- a/src/comm_tiled.h +++ b/src/comm_tiled.h @@ -37,6 +37,8 @@ class CommTiled : public Comm { // forward comm from a Fix virtual void reverse_comm_fix(class Fix *, int size=0); // reverse comm from a Fix + virtual void reverse_comm_fix_variable(class Fix *); + // variable size reverse comm from a Fix void forward_comm_compute(class Compute *); // forward from a Compute void reverse_comm_compute(class Compute *); // reverse from a Compute void forward_comm_dump(class Dump *); // forward comm from a Dump diff --git a/src/fix.h b/src/fix.h index 2eddd9bd35..500b45a90f 100644 --- a/src/fix.h +++ b/src/fix.h @@ -170,6 +170,7 @@ class Fix : protected Pointers { virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} + virtual int pack_reverse_comm_size(int, int) {return 0;} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} diff --git a/src/fix_shear_history.cpp b/src/fix_shear_history.cpp index da69a95603..5951029328 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_shear_history.cpp @@ -29,6 +29,8 @@ using namespace LAMMPS_NS; using namespace FixConst; +enum{NPARTNER,PERPARTNER}; + /* ---------------------------------------------------------------------- */ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : @@ -37,6 +39,12 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : restart_peratom = 1; create_attribute = 1; + newton_pair = force->newton_pair; + + if (newton_pair) comm_reverse = 1; // just for single npartner value + // variable-size history communicated via + // reverse_comm_fix_variable() + // perform initial allocation of atom-based arrays // register with atom class @@ -56,6 +64,8 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) npartner[i] = 0; maxtouch = 0; + + nlocal_neigh = nall_neigh = 0; } /* ---------------------------------------------------------------------- */ @@ -131,40 +141,51 @@ void FixShearHistory::allocate_pages() /* ---------------------------------------------------------------------- copy shear partner info from neighbor lists to atom arrays should be called whenever neighbor list stores current history info - and need to have atoms store the info + and need to store the info with owned atoms e.g. so atoms can migrate to new procs or between runs when atoms may be added or deleted (neighbor list becomes out-of-date) the next granular neigh list build will put this info back into neigh list - called during run before atom exchanges + called during run before atom exchanges, including for restart files called at end of run via post_run() do not call during setup of run (setup_pre_exchange) b/c there is no guarantee of a current neigh list (even on continued run) if run command does a 2nd run with pre = no, then no neigh list will be built, but old neigh list will still have the info + newton ON and newton OFF versions + newton ON does reverse comm to acquire shear partner info from ghost atoms + newton OFF works with smaller vectors that don't include ghost info ------------------------------------------------------------------------- */ void FixShearHistory::pre_exchange() +{ + if (newton_pair) pre_exchange_newton(); + else pre_exchange_no_newton(); +} + +/* ---------------------------------------------------------------------- */ + +void FixShearHistory::pre_exchange_newton() { int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *shear,*allshear,**firstshear; - // nlocal may include atoms added since last neigh build + // NOTE: all operations until very end are with + // nlocal_neigh <= current nlocal and nall_neigh + // b/c previous neigh list was built with nlocal_neigh,nghost_neigh + // nlocal can be larger if other fixes added atoms at this pre_exchange() - int nlocal = atom->nlocal; - - // zero npartner for all current atoms + // zero npartner for owned+ghost atoms // clear 2 page data structures - for (i = 0; i < nlocal; i++) npartner[i] = 0; + for (i = 0; i < nall_neigh; i++) npartner[i] = 0; ipage->reset(); dpage->reset(); // 1st loop over neighbor list - // calculate npartner for each owned atom - // nlocal_neigh = nlocal when neigh list was built, may be smaller than nlocal + // calculate npartner for owned+ghost atoms tagint *tag = atom->tag; NeighList *list = pair->list; @@ -175,8 +196,133 @@ void FixShearHistory::pre_exchange() firsttouch = list->listgranhistory->firstneigh; firstshear = list->listgranhistory->firstdouble; - int nlocal_neigh = 0; - if (inum) nlocal_neigh = ilist[inum-1] + 1; + 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; + npartner[j]++; + } + } + } + + // perform reverse comm to augment owned npartner counts with ghost counts + + commflag = NPARTNER; + comm->reverse_comm_fix(this,0); + + // get page chunks to store atom IDs and shear history for owned+ghost 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"); + } + } + + for (i = nlocal_neigh; i < nall_neigh; i++) { + 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 owned+ghost atoms + // re-zero npartner to use as counter + + for (i = 0; i < nall_neigh; i++) npartner[i] = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + allshear = firstshear[i]; + jnum = numneigh[i]; + touch = firsttouch[i]; + + for (jj = 0; jj < jnum; jj++) { + if (touch[jj]) { + shear = &allshear[3*jj]; + j = jlist[jj]; + j &= NEIGHMASK; + 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]; + 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]; + } + } + } + + // perform reverse comm to augment + // owned atom partner/shearpartner with ghost info + // use variable variant b/c size of packed data can be arbitrarily large + // if many touching neighbors for large particle + + commflag = PERPARTNER; + comm->reverse_comm_fix_variable(this); + + // set maxtouch = max # of partners of any owned atom + // bump up comm->maxexchange_fix if necessary + + maxtouch = 0; + for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); + comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxtouch+1); + + // zero npartner values from previous nlocal_neigh to current nlocal + + int nlocal = atom->nlocal; + for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixShearHistory::pre_exchange_no_newton() +{ + int i,j,ii,jj,m,n,inum,jnum; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *shear,*allshear,**firstshear; + + // NOTE: all operations until very end are with nlocal_neigh <= current nlocal + // b/c previous neigh list was built with nlocal_neigh + // nlocal can be larger if other fixes added atoms at this pre_exchange() + + // zero npartner for owned atoms + // clear 2 page data structures + + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; + + ipage->reset(); + dpage->reset(); + + // 1st loop over neighbor list + // calculate npartner for owned atoms + + tagint *tag = atom->tag; + NeighList *list = pair->list; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + firsttouch = list->listgranhistory->firstneigh; + firstshear = list->listgranhistory->firstdouble; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -194,7 +340,7 @@ void FixShearHistory::pre_exchange() } } - // get page chunks to store atom IDs and shear history for my atoms + // get page chunks to store atom IDs and shear history for owned atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -206,10 +352,10 @@ void FixShearHistory::pre_exchange() } // 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 + // store atom IDs and shear history for owned atoms + // re-zero npartner to use as counter - for (i = 0; i < nlocal; i++) npartner[i] = 0; + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -243,10 +389,15 @@ void FixShearHistory::pre_exchange() // set maxtouch = max # of partners of any owned atom // bump up comm->maxexchange_fix if necessary - + maxtouch = 0; - for (i = 0; i < nlocal; i++) maxtouch = MAX(maxtouch,npartner[i]); + for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxtouch+1); + + // zero npartner values from previous nlocal_neigh to current nlocal + + int nlocal = atom->nlocal; + for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0; } /* ---------------------------------------------------------------------- */ @@ -324,6 +475,83 @@ void FixShearHistory::set_arrays(int i) npartner[i] = 0; } +/* ---------------------------------------------------------------------- + only called by Comm::reverse_comm_fix_variable for PERPARTNER mode +------------------------------------------------------------------------- */ + +int FixShearHistory::pack_reverse_comm_size(int n, int first) +{ + int i,j,k,last; + + int m = 0; + last = first + n; + + for (i = first; i < last; i++) + m += 1 + 4*npartner[i]; + + return m; +} + +/* ---------------------------------------------------------------------- + two modes: NPARTNER and PERPARTNER +------------------------------------------------------------------------- */ + +int FixShearHistory::pack_reverse_comm(int n, int first, double *buf) +{ + int i,j,k,last; + + int m = 0; + last = first + n; + + if (commflag == NPARTNER) { + for (i = first; i < last; i++) { + buf[m++] = npartner[i]; + } + } else if (commflag == PERPARTNER) { + for (i = first; i < last; i++) { + buf[m++] = npartner[i]; + for (int k = 0; k < npartner[i]; k++) { + buf[m++] = partner[i][k]; + buf[m++] = shearpartner[i][k][0]; + buf[m++] = shearpartner[i][k][1]; + buf[m++] = shearpartner[i][k][2]; + } + } + } + + return m; +} + +/* ---------------------------------------------------------------------- + two modes: NPARTNER and PERPARTNER +------------------------------------------------------------------------- */ + +void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,k,kk,ncount; + + int m = 0; + + if (commflag == NPARTNER) { + for (i = 0; i < n; i++) { + j = list[i]; + npartner[j] += static_cast (buf[m++]); + } + } else if (commflag == PERPARTNER) { + for (i = 0; i < n; i++) { + j = list[i]; + ncount = static_cast (buf[m++]); + for (int k = 0; k < ncount; k++) { + kk = npartner[j]++; + partner[j][kk] = static_cast (buf[m++]); + shearpartner[j][kk][0] = buf[m++]; + shearpartner[j][kk][1] = buf[m++]; + shearpartner[j][kk][2] = buf[m++]; + } + } + } +} + /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index beb9683bca..d624c9efef 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -44,6 +44,10 @@ class FixShearHistory : public Fix { void grow_arrays(int); void copy_arrays(int, int, int); void set_arrays(int); + + int pack_reverse_comm_size(int, int); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); int pack_exchange(int, double *); int unpack_exchange(int, double *); int pack_restart(int, double *); @@ -52,11 +56,17 @@ class FixShearHistory : public Fix { int maxsize_restart(); protected: + int newton_pair; + int nlocal_neigh; // nlocal at last time neigh list was built + int nall_neigh; // ditto for nlocal+nghost + int *npartner; // # of touching partners of each atom tagint **partner; // global atom IDs for the partners double (**shearpartner)[3]; // 3 shear values with the partner int maxtouch; // max # of touching partners for my atoms + int commflag; // mode of reverse comm to get ghost info + class Pair *pair; int *computeflag; // computeflag in PairGranHookeHistory @@ -64,6 +74,8 @@ class FixShearHistory : public Fix { MyPage *ipage; // pages of partner atom IDs MyPage *dpage; // pages of shear history with partners + void pre_exchange_newton(); + void pre_exchange_no_newton(); void allocate_pages(); }; diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp index e6ef93c5eb..0475a52278 100644 --- a/src/neigh_gran.cpp +++ b/src/neigh_gran.cpp @@ -65,6 +65,8 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) FixShearHistory *fix_history = list->fix_history; if (fix_history) { + fix_history->nlocal_neigh = nlocal; + fix_history->nall_neigh = nall; npartner = fix_history->npartner; partner = fix_history->partner; shearpartner = fix_history->shearpartner; @@ -160,7 +162,6 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) /* ---------------------------------------------------------------------- granular particles N^2 / 2 search for neighbor pairs with full Newton's 3rd law - no shear history is allowed for this option pair added to list if atoms i and j are both owned and i < j if j is ghost only me or other proc adds pair decision based on itag,jtag tests @@ -168,10 +169,20 @@ void Neighbor::granular_nsq_no_newton(NeighList *list) void Neighbor::granular_nsq_newton(NeighList *list) { - int i,j,n,itag,jtag,bitmask; + int i,j,m,n,nn,itag,jtag,bitmask; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr; + int *neighptr,*touchptr; + double *shearptr; + + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage *ipage_touch; + MyPage *dpage_shear; double **x = atom->x; double *radius = atom->radius; @@ -191,12 +202,35 @@ void Neighbor::granular_nsq_newton(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; + FixShearHistory *fix_history = list->fix_history; + if (fix_history) { + fix_history->nlocal_neigh = nlocal; + fix_history->nall_neigh = nall; + 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; + } + int inum = 0; ipage->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); + if (fix_history) { + nn = 0; + touchptr = ipage_touch->vget(); + shearptr = dpage_shear->vget(); + } itag = tag[i]; xtmp = x[i][0]; @@ -233,7 +267,34 @@ void Neighbor::granular_nsq_newton(NeighList *list) radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutsq) neighptr[n++] = j; + if (rsq <= cutsq) { + neighptr[n++] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + 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]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } } ilist[inum++] = i; @@ -242,6 +303,13 @@ void Neighbor::granular_nsq_newton(NeighList *list) 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; @@ -297,6 +365,8 @@ void Neighbor::granular_bin_no_newton(NeighList *list) 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; @@ -398,17 +468,26 @@ void Neighbor::granular_bin_no_newton(NeighList *list) /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with full Newton's 3rd law - no shear history is allowed for this option each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::granular_bin_newton(NeighList *list) { - int i,j,k,n,ibin; + int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr; + int *neighptr,*touchptr; + double *shearptr; + + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage *ipage_touch; + MyPage *dpage_shear; // bin local & ghost atoms @@ -418,6 +497,7 @@ void Neighbor::granular_bin_newton(NeighList *list) double **x = atom->x; double *radius = atom->radius; + tagint *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; tagint *molecule = atom->molecule; @@ -431,12 +511,35 @@ void Neighbor::granular_bin_newton(NeighList *list) int *stencil = list->stencil; MyPage *ipage = list->ipage; + 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; + } + int inum = 0; ipage->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); + if (fix_history) { + nn = 0; + touchptr = ipage_touch->vget(); + shearptr = dpage_shear->vget(); + } xtmp = x[i][0]; ytmp = x[i][1]; @@ -465,7 +568,34 @@ void Neighbor::granular_bin_newton(NeighList *list) radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutsq) neighptr[n++] = j; + if (rsq <= cutsq) { + neighptr[n] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + 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]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } } // loop over all atoms in other bins in stencil, store every pair @@ -482,7 +612,34 @@ void Neighbor::granular_bin_newton(NeighList *list) radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutsq) neighptr[n++] = j; + if (rsq <= cutsq) { + neighptr[n] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + 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]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } } } @@ -492,6 +649,13 @@ void Neighbor::granular_bin_newton(NeighList *list) 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; @@ -500,17 +664,26 @@ void Neighbor::granular_bin_newton(NeighList *list) /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with Newton's 3rd law for triclinic - no shear history is allowed for this option each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::granular_bin_newton_tri(NeighList *list) { - int i,j,k,n,ibin; + int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr; + int *neighptr,*touchptr; + double *shearptr; + + NeighList *listgranhistory; + int *npartner; + tagint **partner; + double (**shearpartner)[3]; + int **firsttouch; + double **firstshear; + MyPage *ipage_touch; + MyPage *dpage_shear; // bin local & ghost atoms @@ -520,6 +693,7 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) double **x = atom->x; double *radius = atom->radius; + tagint *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; tagint *molecule = atom->molecule; @@ -533,12 +707,35 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) int *stencil = list->stencil; MyPage *ipage = list->ipage; + 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; + } + int inum = 0; ipage->reset(); + if (fix_history) { + ipage_touch->reset(); + dpage_shear->reset(); + } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); + if (fix_history) { + nn = 0; + touchptr = ipage_touch->vget(); + shearptr = dpage_shear->vget(); + } xtmp = x[i][0]; ytmp = x[i][1]; @@ -572,7 +769,34 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutsq) neighptr[n++] = j; + if (rsq <= cutsq) { + neighptr[n++] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + 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]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } } } @@ -582,6 +806,13 @@ void Neighbor::granular_bin_newton_tri(NeighList *list) 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;