Merge pull request #4407 from jtclemm/small-patches

Bug fix for domain class, minor edits to other classes

doc/src/bond_bpm_rotational.rst

@@ -132,8 +132,8 @@ or :doc:`read_restart <read_restart>` commands:
 * :math:`k_b`           (force*distance/radians units)
 * :math:`f_{r,c}`       (force units)
 * :math:`f_{s,c}`       (force units)
-* :math:`\tau_{b,c}`    (force*distance units)
 * :math:`\tau_{t,c}`    (force*distance units)
+* :math:`\tau_{b,c}`    (force*distance units)
 * :math:`\gamma_n`      (force/velocity units)
 * :math:`\gamma_s`      (force/velocity units)
 * :math:`\gamma_r`      (force*distance/velocity units)
@@ -154,8 +154,11 @@ page on BPMs.
 
 If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
 during a simulation run. This will prevent some unnecessary calculation.
-However, if a bond reaches a damage criterion greater than one,
-it will trigger an error.
+The recommended bond communication distance no longer depends on bond failure
+coefficients (which are ignored) but instead corresponds to the typical heurestic
+maximum strain used by typical non-bpm bond styles. Similar behavior to *break no*
+can also be attained by setting arbitrarily high values for all four failure
+coefficients. One cannot use *break no* with *smooth yes*.
 
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed

doc/src/bond_bpm_spring.rst

@@ -117,8 +117,11 @@ page on BPMs.
 
 If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
 during a simulation run. This will prevent some unnecessary calculation.
-However, if a bond reaches a strain greater than :math:`\epsilon_c`,
-it will trigger an error.
+The recommended bond communication distance no longer depends on the value of
+:math:`\epsilon_c` (which is ignored) but instead corresponds to the typical
+heurestic maximum strain used by typical non-bpm bond styles. Similar behavior
+to *break no* can also be attained by setting an arbitrarily high value of
+:math:`\epsilon_c`. One cannot use *break no* with *smooth yes*.
 
 .. versionadded:: TBD
 

doc/src/pair_granular.rst

@@ -235,12 +235,15 @@ given by:
 
 .. math::
 
-   \eta_n = \alpha (m_{eff}k_n)^{1/2}
+   \eta_n = \alpha (m_{eff}k_{nd})^{1/2}
 
-For normal contact models based on material parameters, :math:`k_n = 4/3Ea`. This
-damping model is not compatible with cohesive normal models such as *JKR* or *DMT*.
-The parameter :math:`\alpha` is related to the restitution coefficient *e*
-according to:
+where :math:`k_{nd}` is an effective harmonic stiffness equal to the ratio of
+the normal force to the overlap. For example, :math:`k_{nd} = 4/3Ea` for a
+Hertz contact model based on material parameters with :math:`a` being
+the contact radius of :math:`\sqrt{\delta R}`. For Hooke, :math:`k_{nd}`
+is simply the spring constant or :math:`k_{n}`. This damping model is not
+compatible with cohesive normal models such as *JKR* or *DMT*. The parameter
+:math:`\alpha` is related to the restitution coefficient *e* according to:
 
 .. math::
 
@@ -251,25 +254,26 @@ of the normal contact model parameters should be between 0 and 1, but
 no error check is performed on this.
 
 The *coeff_restitution* model is useful when a specific normal coefficient of
-restitution :math:`e` is required. In these models, the normal coefficient of
-restitution :math:`e` is specified as an input. Following the approach of
-:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model,
-*coeff_restitution* calculates the damping coefficient as:
+restitution :math:`e` is required. It operates much like the *Tsuji* model
+but, the normal coefficient of restitution :math:`e` is specified as an input
+in place of the usual :math:`\eta_{n0}` value in the normal model. Following
+the approach of :ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke*
+normal model, *coeff_restitution* then calculates the damping coefficient as:
 
 .. math::
 
-   \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
+   \eta_n = \sqrt{\frac{4m_{eff}k_{nd}}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
 
+where :math:`k_{nd}` is the same stiffness defined in the above *Tsuji* model.
 For any other normal model, e.g. the *hertz* and *hertz/material* models, the damping
 coefficient is:
 
 .. math::
 
-   \eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}(R_{eff} \delta_{ij})^{\frac{1}{4}}\sqrt{\frac{3}{2}k_n m_{eff}} ,
+   \eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}\sqrt{\frac{3}{2}k_{nd} m_{eff}} ,
 
-where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
-*coeff_restitution* accounts for the effective mass, effective radius, and
-pairwise overlaps (except when used with the *hooke* normal model) when calculating
+Since *coeff_restitution* accounts for the effective mass, effective radius,
+and pairwise overlaps (except when used with the *hooke* normal model) when calculating
 the damping coefficient, it accurately reproduces the specified coefficient of
 restitution for both monodisperse and polydisperse particle pairs.  This damping
 model is not compatible with cohesive normal models such as *JKR* or *DMT*.

examples/bpm/pour/in.bpm.pour

@@ -13,9 +13,9 @@ region          wall_cyl cylinder z 0.0 0.0 10.0 EDGE EDGE side in
 region          dropzone cylinder z 0.0 0.0 10.0 40.0 50.0 side in
 
 pair_style      gran/hertz/history 1.0 NULL 0.5 NULL 0.1 1
-bond_style      bpm/rotational
+bond_style      bpm/rotational break no smooth no
 pair_coeff      1 1
-bond_coeff      1 1.0 0.2 0.01 0.01 2.0 0.4 0.02 0.02 0.2 0.04 0.002 0.002
+bond_coeff      1 1.0 0.2 0.01 0.01 0.0 0.0 0.0 0.0 0.2 0.04 0.002 0.002
 
 compute         nbond all nbond/atom
 compute         tbond all reduce sum c_nbond

src/BPM/bond_bpm.cpp

@@ -377,8 +377,13 @@ double BondBPM::equilibrium_distance(int /*i*/)
     r0_max_estimate = temp;
   }
 
-  // Divide out heuristic prefactor added in comm class
-  return max_stretch * r0_max_estimate / 1.5;
+  double dist = r0_max_estimate;
+
+  // Add stretch and remove heuristic prefactor (added in comm class)
+  if (break_flag)
+    dist *= max_stretch / 1.5;
+
+  return dist;
 }
 
 /* ----------------------------------------------------------------------

src/BPM/bond_bpm_rotational.cpp

@@ -523,7 +523,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
     breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
                               torque1on2, torque2on1);
 
-    if (breaking >= 1.0) {
+    if ((breaking >= 1.0) && break_flag) {
       bondlist[n][2] = 0;
       process_broken(i1, i2);
       continue;
@@ -684,6 +684,9 @@ void BondBPMRotational::settings(int narg, char **arg)
       error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
     }
   }
+
+  if (smooth_flag && !break_flag)
+    error->all(FLERR, "Illegal bond bpm command, must turn off smoothing with break no option");
 }
 
 /* ----------------------------------------------------------------------

src/BPM/bond_bpm_spring.cpp

@@ -245,7 +245,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
     r = sqrt(rsq);
     e = (r - r0) / r0;
 
-    if (fabs(e) > ecrit[type]) {
+    if ((fabs(e) > ecrit[type]) && break_flag) {
       bondlist[n][2] = 0;
       process_broken(i1, i2);
 
@@ -485,6 +485,9 @@ void BondBPMSpring::settings(int narg, char **arg)
       error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
     }
   }
+
+  if (smooth_flag && !break_flag)
+    error->all(FLERR, "Illegal bond bpm command, must turn off smoothing with break no option");
 }
 
 /* ----------------------------------------------------------------------

src/EXTRA-FIX/fix_deform_pressure.cpp

@@ -407,6 +407,16 @@ void FixDeformPressure::init()
 
   set_box.vol_start = domain->xprd * domain->yprd * domain->zprd;
 
+  // reset cumulative counters to match resetting "start" variables in parent
+
+  for (int i = 0; i < 7; i++) {
+    set_extra[i].cumulative_remap = 0.0;
+    set_extra[i].cumulative_shift = 0.0;
+    set_extra[i].cumulative_vshift[0] = 0.0;
+    set_extra[i].cumulative_vshift[1] = 0.0;
+    set_extra[i].cumulative_vshift[2] = 0.0;
+  }
+
   // check optional variables for PRESSURE or PMEAN style
 
   for (int i = 0; i < 7; i++) {

src/GRANULAR/gran_sub_mod_damping.cpp

@@ -152,9 +152,8 @@ double GranSubModDampingTsuji::calculate_forces()
 ------------------------------------------------------------------------- */
 
 GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) :
-    GranSubModDamping(gm, lmp)
+    GranSubModDampingTsuji(gm, lmp)
 {
-  allow_cohesion = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -171,17 +170,3 @@ void GranSubModDampingCoeffRestitution::init()
     damp /= sqrt(MY_PI * MY_PI + logcor * logcor);
   }
 }
-
-/* ---------------------------------------------------------------------- */
-
-double GranSubModDampingCoeffRestitution::calculate_forces()
-{
-  // in case argument <= 0 due to precision issues
-  double sqrt1;
-  if (gm->delta > 0.0)
-    sqrt1 = MAX(0.0, gm->meff * gm->Fnormal / gm->delta);
-  else
-    sqrt1 = 0.0;
-  damp_prefactor = damp * sqrt(sqrt1);
-  return -damp_prefactor * gm->vnnr;
-}

src/GRANULAR/gran_sub_mod_damping.h

@@ -87,11 +87,10 @@ namespace Granular_NS {
 
   /* ---------------------------------------------------------------------- */
 
-  class GranSubModDampingCoeffRestitution : public GranSubModDamping {
+  class GranSubModDampingCoeffRestitution : public GranSubModDampingTsuji {
    public:
     GranSubModDampingCoeffRestitution(class GranularModel *, class LAMMPS *);
     void init() override;
-    double calculate_forces() override;
   };
 
   /* ---------------------------------------------------------------------- */

src/domain.cpp

@@ -1668,7 +1668,7 @@ void Domain::remap(double *x)
    adjust 3 image flags encoded in image accordingly
    resulting coord must satisfy lo <= coord < hi
    MAX is important since coord - prd < lo can happen when coord = hi
-   for triclinic, point is converted to lamda coords (0-1) before doing remap
+   for triclinic, point is converted to lamda coords (0-1) within remap()
    image = 10 bits for each dimension
    increment/decrement in wrap-around fashion
 ------------------------------------------------------------------------- */
@@ -1679,9 +1679,7 @@ void Domain::remap_all()
   imageint *image = atom->image;
   int nlocal = atom->nlocal;
 
-  if (triclinic) x2lamda(nlocal);
   for (int i = 0; i < nlocal; i++) remap(x[i],image[i]);
-  if (triclinic) lamda2x(nlocal);
 }
 
 /* ----------------------------------------------------------------------