[pull] develop from lammps:develop

doc/src/fix_wall_gran.rst

@@ -44,7 +44,7 @@ Syntax
   .. parsed-literal::
 
        *xplane* or *yplane* or *zplane* args = lo hi
-         lo,hi = position of lower and upper plane (distance units), either can be NULL)
+         lo,hi = position of lower and upper plane (distance units), each can be NULL or an equal-style variable as v_name (see below)
 
 * zero or more keyword/value pairs may be appended to args
 * keyword = *wiggle* or *shear* or *contacts* or *temperature*
@@ -76,6 +76,9 @@ Examples
    fix 5 all wall/gran granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0
    fix 6 all wall/gran hooke  200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts
 
+   variable zhi equal ramp(20.0,15.0)
+   fix 7 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 v_zhi
+
 Description
 """""""""""
 
@@ -170,6 +173,28 @@ is desired.
 The *zcylinder* wallstyle has been removed.  Pleas use :doc:`fix
 wall/gran/region <fix_wall_gran_region>` instead.
 
+.. versionadded:: TBD
+
+The *lo* and *hi* wall positions can also be set by an equal-style
+:doc:`variable <variable>`, specified as v_name, where "name" is the
+variable name.  The variable is evaluated at every timestep, so the wall
+position can change during a run, for example to model a piston
+compressing a bed of granular particles.  Since the damping and friction
+forces of granular models depend on the relative velocity between
+particle and wall, the velocity of such a wall is inferred from the
+change of the wall position between consecutive timesteps, in the same
+manner as for moving regions with :doc:`fix wall/gran/region
+<fix_wall_gran_region>`: there is no analytic formula for the velocity
+of a wall position defined by a variable.  The wall position should
+therefore be a continuous function of the elapsed time; see the
+discussion of the *move* keyword of the :doc:`region <region>` command
+for examples, including how to make the motion span consecutive runs in
+a continuous fashion.  The wall velocity is zero when the variable is
+evaluated for the first time.  LAMMPS stops with an error if the *lo*
+wall position does not remain below the *hi* wall position during a run.
+A wall position set by a variable cannot be combined with the *wiggle*
+or *shear* keyword.
+
 Optionally, the wall can be moving, if the *wiggle* or *shear*
 keywords are appended.  Both keywords cannot be used together.
 
@@ -309,6 +334,9 @@ LAMMPS was built with that package.  See the :doc:`Build package <Build_package>
 
 Any dimension (xyz) that has a granular wall must be non-periodic.
 
+The *wall/gran/kk* style does not support wall positions set by an
+equal-style variable.
+
 Related commands
 """"""""""""""""
 

src/GRANULAR/fix_wall_gran.cpp

@@ -133,54 +133,44 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
   numwalls = 0;
   if (iarg >= narg) error->all(FLERR, "Illegal fix wall/gran command");
 
-  if (strcmp(arg[iarg],"xplane") == 0) {
+  xstyle[0] = xstyle[1] = NONE;
+  xvar[0] = xvar[1] = -1;
+  xstr[0] = xstr[1] = nullptr;
+  velwall[0] = velwall[1] = prevwall[0] = prevwall[1] = 0.0;
+  velstep = -1;
+  velflag = varflag = 0;
+
+  if ((strcmp(arg[iarg],"xplane") == 0) || (strcmp(arg[iarg],"yplane") == 0) ||
+      (strcmp(arg[iarg],"zplane") == 0)) {
     if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");
-    wallstyle = XPLANE;
+    if (strcmp(arg[iarg],"xplane") == 0) wallstyle = XPLANE;
+    else if (strcmp(arg[iarg],"yplane") == 0) wallstyle = YPLANE;
+    else wallstyle = ZPLANE;
     numwalls = 0;
     if (strcmp(arg[iarg+1],"NULL") == 0) {
       lo = -BIG;
-    } else {
-      lo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+    } else if (utils::strmatch(arg[iarg+1],"^v_")) {
+      delete[] xstr[0];
+      xstr[0] = utils::strdup(arg[iarg+1]+2);
+      xstyle[0] = EQUAL;
+      lo = 0.0;
       ++numwalls;
-    }
-    if (strcmp(arg[iarg+2],"NULL") == 0) {
-      hi = BIG;
-    } else {
-      hi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
-      ++numwalls;
-    }
-    iarg += 3;
-  } else if (strcmp(arg[iarg],"yplane") == 0) {
-    if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");
-    wallstyle = YPLANE;
-    numwalls = 0;
-    if (strcmp(arg[iarg+1],"NULL") == 0) {
-      lo = -BIG;
     } else {
       lo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      xstyle[0] = CONSTANT;
       ++numwalls;
     }
     if (strcmp(arg[iarg+2],"NULL") == 0) {
       hi = BIG;
-    } else {
-      hi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+    } else if (utils::strmatch(arg[iarg+2],"^v_")) {
+      delete[] xstr[1];
+      xstr[1] = utils::strdup(arg[iarg+2]+2);
+      xstyle[1] = EQUAL;
+      hi = 0.0;
       ++numwalls;
-    }
-    iarg += 3;
-  } else if (strcmp(arg[iarg],"zplane") == 0) {
-    if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");
-    wallstyle = ZPLANE;
-    numwalls = 0;
-    if (strcmp(arg[iarg+1],"NULL") == 0) {
-      lo = -BIG;
-    } else {
-      lo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
-      ++numwalls;
-    }
-    if (strcmp(arg[iarg+2],"NULL") == 0) {
-      hi = BIG;
     } else {
       hi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+      xstyle[1] = CONSTANT;
       ++numwalls;
     }
     iarg += 3;
@@ -268,6 +258,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
   if ((wiggle || wshear) && wallstyle == REGION)
     error->all(FLERR,"Cannot wiggle or shear with fix wall/gran/region");
 
+  // walls with a position variable cannot be combined with prescribed wall motion
+
+  if ((xstyle[0] == EQUAL) || (xstyle[1] == EQUAL)) velflag = 1;
+  if ((wiggle || wshear) && velflag)
+    error->all(FLERR,"Cannot wiggle or shear a fix wall/gran wall with a position variable");
+
   // setup oscillations
 
   if (wiggle) omega = 2.0*MY_PI / period;
@@ -338,6 +334,8 @@ FixWallGran::~FixWallGran()
   delete model;
   delete[] tstr;
   delete[] idregion;
+  delete[] xstr[0];
+  delete[] xstr[1];
   memory->destroy(history_one);
   memory->destroy(mass_rigid);
 
@@ -398,6 +396,22 @@ void FixWallGran::init()
     if (! input->variable->equalstyle(tvar))
       error->all(FLERR, "Variable {} for fix wall/gran must be an equal style variable", tstr);
   }
+
+  // locate and check variables for wall positions and velocities
+
+  varflag = 0;
+  if (tstr) varflag = 1;
+  for (int m = 0; m < 2; m++) {
+    if (xstyle[m] == EQUAL) {
+      xvar[m] = input->variable->find(xstr[m]);
+      if (xvar[m] < 0)
+        error->all(FLERR, "Variable {} for fix wall/gran wall position does not exist", xstr[m]);
+      if (! input->variable->equalstyle(xvar[m]))
+        error->all(FLERR, "Variable {} for fix wall/gran wall position must be "
+                   "an equal style variable", xstr[m]);
+      varflag = 1;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -417,7 +431,7 @@ void FixWallGran::setup(int vflag)
 
 void FixWallGran::post_force(int /*vflag*/)
 {
-  int i,j,n;
+  int i,j,n,mwall;
   double dx,dy,dz,del1,del2,rwall,meff;
   double *forces, *torquesi;
   double vwall[3];
@@ -452,11 +466,34 @@ void FixWallGran::post_force(int /*vflag*/)
   }
 
   // set position of wall to initial settings and velocity to 0.0
+  // evaluate variables for wall position and velocity, if defined
   // if wiggle or shear, set wall position and velocity accordingly
 
   double wlo = lo;
   double whi = hi;
   vwall[0] = vwall[1] = vwall[2] = 0.0;
+  if (varflag) modify->clearstep_compute();
+  if (xstyle[0] == EQUAL) wlo = input->variable->compute_equal(xvar[0]);
+  if (xstyle[1] == EQUAL) whi = input->variable->compute_equal(xvar[1]);
+  if ((xstyle[0] != NONE) && (xstyle[1] != NONE) && (wlo >= whi))
+    error->all(FLERR, Error::NOLASTLINE,
+               "Fix wall/gran lo wall position {} must remain below hi wall position {}", wlo, whi);
+
+  // infer velocity of walls with a position variable from the change of the wall
+  // position since the previous evaluation, as done for moving regions: there is
+  // no analytic formula for the velocity of a variable-defined wall position.
+  // the velocity is zero at the very first evaluation; re-evaluations on the
+  // same timestep (e.g. during setup of a continued run) keep the velocity.
+
+  if (velflag && (update->ntimestep != velstep)) {
+    if (velstep >= 0) {
+      if (xstyle[0] == EQUAL) velwall[0] = (wlo - prevwall[0]) / update->dt;
+      if (xstyle[1] == EQUAL) velwall[1] = (whi - prevwall[1]) / update->dt;
+    }
+    prevwall[0] = wlo;
+    prevwall[1] = whi;
+    velstep = update->ntimestep;
+  }
   if (wiggle) {
     double arg = omega * (update->ntimestep - time_origin) * dt;
     if (wallstyle == axis) {
@@ -506,29 +543,36 @@ void FixWallGran::post_force(int /*vflag*/)
       Twall = input->variable->compute_equal(tvar);
     model->Tj = Twall;
   }
+  if (varflag) modify->addstep_compute(update->ntimestep + 1);
 
   for (int i = 0; i < nlocal; i++) {
     if (!(mask[i] & groupbit)) continue;
 
     dx = dy = dz = 0.0;
+    mwall = 0;
 
     if (wallstyle == XPLANE) {
       del1 = x[i][0] - wlo;
       del2 = whi - x[i][0];
       if (del1 < del2) dx = del1;
-      else dx = -del2;
+      else { dx = -del2; mwall = 1; }
     } else if (wallstyle == YPLANE) {
       del1 = x[i][1] - wlo;
       del2 = whi - x[i][1];
       if (del1 < del2) dy = del1;
-      else dy = -del2;
+      else { dy = -del2; mwall = 1; }
     } else if (wallstyle == ZPLANE) {
       del1 = x[i][2] - wlo;
       del2 = whi - x[i][2];
       if (del1 < del2) dz = del1;
-      else dz = -del2;
+      else { dz = -del2; mwall = 1; }
     }
 
+    // for walls with a position variable, use the velocity of the nearer wall
+    // XPLANE,YPLANE,ZPLANE are 0,1,2, i.e. the index of the wall normal direction
+
+    if (velflag) vwall[wallstyle] = velwall[mwall];
+
     // Reset model and copy initial geometric data
     model->dx[0] = dx;
     model->dx[1] = dy;

src/GRANULAR/fix_wall_gran.h

@@ -66,6 +66,16 @@ class FixWallGran : public Fix {
   double Twall;
   char *idregion;
 
+  // wall positions set by equal-style variables
+
+  int xstyle[2], xvar[2];    // style and variable index for lo/hi wall position
+  char *xstr[2];             // variable names for lo/hi wall position
+  double velwall[2];         // current velocity of lo/hi wall
+  double prevwall[2];        // lo/hi wall position at previous evaluation
+  bigint velstep;            // timestep of last wall velocity update
+  int velflag;               // 1 if any wall position is set by a variable
+  int varflag;               // 1 if any wall attribute is set by a variable
+
   int use_history;       // if particle/wall interaction stores history
   int history_update;    // flag for whether shear history is updated
   int size_history;      // # of shear history values per contact

src/GRANULAR/fix_wall_gran_region.cpp

@@ -24,6 +24,7 @@
 #include "error.h"
 #include "input.h"
 #include "memory.h"
+#include "modify.h"
 #include "neighbor.h"
 #include "math_extra.h"
 #include "region.h"
@@ -186,8 +187,11 @@ void FixWallGranRegion::post_force(int /*vflag*/)
   if (heat_flag) {
     temperature = atom->temperature;
     heatflow = atom->heatflow;
-    if (tstr)
+    if (tstr) {
+      modify->clearstep_compute();
       Twall = input->variable->compute_equal(tvar);
+      modify->addstep_compute(update->ntimestep + 1);
+    }
     model->Tj = Twall;
   }
 

src/KOKKOS/fix_wall_gran_old.cpp

@@ -333,9 +333,17 @@ FixWallGranOld::FixWallGranOld(LAMMPS *lmp, int narg, char **arg) :
   }
 
   // wallstyle args
+  // wall positions set by equal-style variables are only supported
+  // by the non-accelerated version of fix wall/gran
 
   idregion = nullptr;
 
+  if (((strcmp(arg[iarg],"xplane") == 0) || (strcmp(arg[iarg],"yplane") == 0) ||
+       (strcmp(arg[iarg],"zplane") == 0)) && (narg > iarg+2)) {
+    if (utils::strmatch(arg[iarg+1],"^v_") || utils::strmatch(arg[iarg+2],"^v_"))
+      error->all(FLERR,"Fix {} does not support wall positions set by a variable", style);
+  }
+
   if (strcmp(arg[iarg],"xplane") == 0) {
     if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");
     wallstyle = XPLANE;