diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 92d5aed075d..6eccda3e585 100644 --- a/doc/src/fix_wall_gran.rst +++ b/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 ` instead. +.. versionadded:: TBD + +The *lo* and *hi* wall positions can also be set by an equal-style +:doc:`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 +`: 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 ` 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 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 """""""""""""""" diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 186f19c1d6f..030682ad5b4 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/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; diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index abb769ee5bb..45dd3e66061 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/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 diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index bb2069516f3..6cae95fa934 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/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; } diff --git a/src/KOKKOS/fix_wall_gran_old.cpp b/src/KOKKOS/fix_wall_gran_old.cpp index 1356ebbd1d3..1a4d8b783fe 100644 --- a/src/KOKKOS/fix_wall_gran_old.cpp +++ b/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;