diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 4b60f96714..4234f91a9b 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -182,6 +182,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : shear[i][0] = shear[i][1] = shear[i][2] = 0.0; time_origin = update->ntimestep; + laststep = -1; } /* ---------------------------------------------------------------------- */ @@ -282,6 +283,9 @@ void FixWallGran::post_force(int vflag) int *mask = atom->mask; int nlocal = atom->nlocal; + if (update->ntimestep > laststep) shearupdate = 1; + else shearupdate = 0; + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -338,6 +342,8 @@ void FixWallGran::post_force(int vflag) } } } + + laststep = update->ntimestep; } /* ---------------------------------------------------------------------- */ @@ -489,18 +495,22 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, // shear history effects - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; + if (shearupdate) { + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; + } shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); // rotate shear displacements rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; rsht = rsht*rsqinv; - shear[0] -= rsht*dx; - shear[1] -= rsht*dy; - shear[2] -= rsht*dz; + if (shearupdate) { + shear[0] -= rsht*dx; + shear[1] -= rsht*dy; + shear[2] -= rsht*dz; + } // tangential forces = shear + tangential velocity damping diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 0f6bb0145f..ae260e478e 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -55,8 +55,10 @@ class FixWallGran : public Fix { int nlevels_respa; int time_origin; + bigint laststep; int *touch; double **shear; + int shearupdate; void hooke(double, double, double, double, double *, double *, double *, double *, double *, double, double); diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 2184abd41d..87e341927a 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -21,6 +21,7 @@ #include "string.h" #include "pair_gran_hertz_history.h" #include "atom.h" +#include "update.h" #include "force.h" #include "neigh_list.h" #include "error.h" @@ -55,6 +56,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -168,9 +172,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag) touch[jj] = 1; shear = &allshear[3*jj]; - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; + if (shearupdate) { + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; + } shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); @@ -178,9 +184,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag) rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; rsht *= rsqinv; - shear[0] -= rsht*delx; - shear[1] -= rsht*dely; - shear[2] -= rsht*delz; + if (shearupdate) { + shear[0] -= rsht*delx; + shear[1] -= rsht*dely; + shear[2] -= rsht*delz; + } // tangential forces = shear + tangential velocity damping @@ -237,6 +245,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag) } } } + + laststep = update->ntimestep; } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index fbf7043098..087ffb4e74 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -49,6 +49,8 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) no_virial_compute = 1; history = 1; fix_history = NULL; + + laststep = -1; } /* ---------------------------------------------------------------------- */ @@ -88,6 +90,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; + int shearupdate = 0; + if (update->ntimestep > laststep) shearupdate = 1; + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -199,9 +204,12 @@ void PairGranHookeHistory::compute(int eflag, int vflag) touch[jj] = 1; shear = &allshear[3*jj]; - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; + + if (shearupdate) { + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; + } shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); @@ -209,9 +217,11 @@ void PairGranHookeHistory::compute(int eflag, int vflag) rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; rsht *= rsqinv; - shear[0] -= rsht*delx; - shear[1] -= rsht*dely; - shear[2] -= rsht*delz; + if (shearupdate) { + shear[0] -= rsht*delx; + shear[1] -= rsht*dely; + shear[2] -= rsht*delz; + } // tangential forces = shear + tangential velocity damping @@ -268,6 +278,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag) } } } + + laststep = update->ntimestep; } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index f005aa4e7c..a27b6ce916 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -46,7 +46,10 @@ class PairGranHookeHistory : public Pair { double dt; int freeze_group_bit; int history; + + bigint laststep; class FixShearHistory *fix_history; + int shearupdate; double *onerad_dynamic,*onerad_frozen; double *maxrad_dynamic,*maxrad_frozen; diff --git a/src/fix.h b/src/fix.h index 1d5910969c..ca7008a8b2 100644 --- a/src/fix.h +++ b/src/fix.h @@ -92,6 +92,7 @@ class Fix : protected Pointers { virtual void init() {} virtual void init_list(int, class NeighList *) {} virtual void setup(int) {} + virtual void setup_pre_exchange() {} virtual void setup_pre_force(int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} diff --git a/src/fix_shear_history.cpp b/src/fix_shear_history.cpp index 5763516bbf..328756e49e 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_shear_history.cpp @@ -33,8 +33,12 @@ using namespace LAMMPS_NS; FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { + // set time_depend so that history will be preserved correctly + // across multiple runs via laststep setting in granular pair styles + restart_peratom = 1; create_attribute = 1; + time_depend = 1; // perform initial allocation of atom-based arrays // register with atom class @@ -85,6 +89,13 @@ void FixShearHistory::init() error->all("Pair style granular with history requires atoms have IDs"); } +/* ---------------------------------------------------------------------- */ + +void FixShearHistory::setup_pre_exchange() +{ + pre_exchange(); +} + /* ---------------------------------------------------------------------- copy shear partner info from neighbor lists to atom arrays so can be exchanged with atoms diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index b6dffa0a4d..8bb3a9814f 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -34,6 +34,7 @@ class FixShearHistory : public Fix { ~FixShearHistory(); int setmask(); void init(); + void setup_pre_exchange(); void pre_exchange(); double memory_usage(); diff --git a/src/modify.cpp b/src/modify.cpp index aa88dd3a07..271e463229 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -255,7 +255,17 @@ void Modify::setup(int vflag) } /* ---------------------------------------------------------------------- - setup pre_force call, only for relevant fixes + setup pre_exchange call, only for fixes that define pre_exchange +------------------------------------------------------------------------- */ + +void Modify::setup_pre_exchange() +{ + for (int i = 0; i < n_pre_exchange; i++) + fix[list_pre_exchange[i]]->setup_pre_exchange(); +} + +/* ---------------------------------------------------------------------- + setup pre_force call, only for fixes that define pre_force ------------------------------------------------------------------------- */ void Modify::setup_pre_force(int vflag) diff --git a/src/modify.h b/src/modify.h index 01a32431ce..81d9401e69 100644 --- a/src/modify.h +++ b/src/modify.h @@ -43,6 +43,7 @@ class Modify : protected Pointers { ~Modify(); void init(); void setup(int); + void setup_pre_exchange(); void setup_pre_force(int); void initial_integrate(int); void post_integrate(); diff --git a/src/respa.cpp b/src/respa.cpp index 806c7bea44..8aa18c0eee 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -327,6 +327,7 @@ void Respa::setup() // build neighbor lists atom->setup(); + modify->setup_pre_exchange(); if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); diff --git a/src/verlet.cpp b/src/verlet.cpp index 0546e55576..2afa3a89ed 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -90,6 +90,7 @@ void Verlet::setup() // build neighbor lists atom->setup(); + modify->setup_pre_exchange(); if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box();