Many minor tweaks, adding self/mass + oxide noshift

doc/src/fix_rheo.rst

@@ -17,7 +17,7 @@ Syntax
 * zmin = minimal number of neighbors for reproducing kernels
 * zero or more keyword/value pairs may be appended to args
 * keyword = *thermal* or *interface/reconstruct* or *surface/detection* or
-            *shift* or *rho/sum* or *density* or *speed/sound*
+            *shift* or *rho/sum* or *density* or *self/mass* or *speed/sound*
 
   .. parsed-literal::
 
@@ -29,6 +29,7 @@ Syntax
          *limit/splash* = threshold for splash particles
        *shift* values = none, turns on velocity shifting
        *rho/sum* values = none, uses the kernel to compute the density of particles
+       *self/mass* values = none, a particle uses its own mass in a rho summation
        *density* values = *rho01*, ... *rho0N* (density)
        *speed/sound* values = *cs0*, ... *csN* (velocity)
 
@@ -106,24 +107,38 @@ threshold for this classification is set by the numerical value of
 
 By default, RHEO integrates particles' densities using a mass diffusion
 equation. Alternatively, one can update densities every timestep by performing
-a kernel summation of the masses of neighboring particles by specifying the *rho/sum* keyword.
+a kernel summation of the masses of neighboring particles by specifying the *rho/sum*
+keyword.
 
-The *density* is used to specify the equilbrium density of each of the N
+The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*.
+Typically, the density :math:`\rho` of a particle is calculated as the sum
+
+.. math::
+   \rho_i = \Sum_{j} W_{ij} M_j
+
+where the summation is over neighbors, :math:`W_{ij}` is the kernel, and :math:`M_j`
+is the mass of particle :math:`j`. The *self/mass* keyword augments this expression
+by replacing :math:`M_j` with :math:`M_i`. This may be useful in simulations of
+multiple fluid phases with large differences in density, :ref:`(Hu) <fix_rheo_hu>`.
+
+The *density* keyword is used to specify the equilbrium density of each of the N
 particle types. It must be followed by N numerical values specifying each
 type's equilibrium density *rho0*.
 
-The *density* is used to specify the speed of sound of each of the N particle
+The *speed/sound* keyword is used to specify the speed of sound of each of the
-types. It must be followed by N numerical values specifying each type's speed
+N particle types. It must be followed by N numerical values specifying each
-of sound *cs*.
+type's speed of sound *cs*.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
+No information about this fix is written to :doc:`binary restart files <restart>`.
+None of the :doc:`fix_modify <fix_modify>` options
 are relevant to this fix.  No global or per-atom quantities are stored
 by this fix for access by various :doc:`output commands <Howto_output>`.
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
@@ -138,7 +153,8 @@ set to all. Only one instance of fix rheo may be defined and it
 must be defined prior to all other RHEO fixes.
 
 This fix is part of the RHEO package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+LAMMPS was built with that package.  See the
+:doc:`Build package <Build_package>` page for more info.
 
 Related commands
 """"""""""""""""
@@ -156,6 +172,10 @@ Default
 
 ----------
 
-.. _howto-howto_rheo_palermo:
+.. _howto_rheo_palermo:
 
 **(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation.
+
+.. _fix_rheo_hu:
+
+**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006).

doc/src/fix_rheo_viscosity.rst

@@ -14,20 +14,26 @@ Syntax
 * rheo/viscosity = style name of this fix command
 * one or more types and viscosity styles must be appended
 * types = lists of types (see below)
-* vstyle = *constant*
+* vstyle = *constant* or *power*
 
   .. parsed-literal::
 
        *constant* args = *eta*
          *eta* = viscosity
 
+       *power* args = *eta*, *gd0*, *K*, *n*
+         *eta* = viscosity
+         *gd0* = critical strain rate
+         *K* = consistency index
+         *n* = power-law exponent
+
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
    fix 1 all rheo/viscosity * constant 1.0
-   fix 1 all rheo/viscosity 1 constant 1.0 2 constant 2.0
+   fix 1 all rheo/viscosity 1 constant 1.0 2 power 0.1 5e-4 0.001 0.5
 
 Description
 """""""""""
@@ -47,18 +53,38 @@ means all types from 1 to n (inclusive).  A trailing asterisk means all types
 from m to :math:`N` (inclusive).  A middle asterisk means all types from m to n
 (inclusive).
 
-The *types* definition is followed by the viscosity style, *vstyle*. Currently,
+The *types* definition is followed by the viscosity style, *vstyle*. Two
-the only option is *constant*. Style *constant* simply applies a constant value
+options are available, *constant* and *power*. Style *constant* simply
-of the viscosity *eta* to each particle of the assigned type.
+applies a constant value of the viscosity *eta* to each particle of the
+assigned type. Style *power* is a Hershchel-Bulkley constitutive equation
+for the stress :math:`\tau`
+
+.. math::
+
+   \tau = \left(\frac{\tau_0}{\dot{\gamma}} + K \dot{\gamma}^{n - 1}\right) \dot{\gamma}, \tau \ge \tau_0
+
+where :math:`\dot{\gamma}` is the strain rate and :math:`tau_0` is the critical
+yield stress, below which :math:`\dot{\gamma} = 0.0`. To avoid divergences, this
+expression is regularized by defining a critical strain rate *gd0*. If the local
+strain rate on a particle falls below this limit, a constant viscosity of *eta*
+is assigned. This implies a value of
+
+.. math::
+   \tau_0 = \eta * \dot{\gamma}_0 - K \dot{\gamma}_0^N
+
+as further discussed in :ref:`(Palermo) <howto_rheo_palermo>`.
+
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
+No information about this fix is written to :doc:`binary restart files <restart>`.
+None of the :doc:`fix_modify <fix_modify>` options
 are relevant to this fix.  No global or per-atom quantities are stored
 by this fix for access by various :doc:`output commands <Howto_output>`.
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
@@ -69,7 +95,8 @@ conjuction with :doc:`fix rheo <fix_rheo>`. The fix group must be
 set to all. Only one instance of fix rheo/viscosity can be defined.
 
 This fix is part of the RHEO package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+LAMMPS was built with that package.  See the
+:doc:`Build package <Build_package>` page for more info.
 
 Related commands
 """"""""""""""""
@@ -82,3 +109,9 @@ Default
 """""""
 
 none
+
+----------
+
+.. _howto_rheo_palermo:
+
+**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation.

src/RHEO/compute_rheo_rho_sum.cpp

@@ -35,7 +35,9 @@ using namespace LAMMPS_NS;
 ComputeRHEORhoSum::ComputeRHEORhoSum(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute RHEO/rho command");
+  if (narg != 4) error->all(FLERR,"Illegal compute RHEO/rho command");
+
+  self_mass_flag = utils::bnumeric(FLERR, arg[3], false, lmp);
 
   comm_forward = 1;
   comm_reverse = 1;
@@ -117,9 +119,15 @@ void ComputeRHEORhoSum::compute_peratom()
       rsq = delx * delx + dely * dely + delz * delz;
       if (rsq < cutsq) {
         w = compute_kernel->calc_w(i, j, delx, dely, delz, sqrt(rsq));
-        rho[i] += w * mass[type[i]];
+
-        if (newton || j < nlocal) {
+        if (self_mass_flag) {
-          rho[j] += w * mass[type[j]];
+          rho[i] += w * mass[type[i]];
+          if (newton || j < nlocal)
+            rho[j] += w * mass[type[j]];
+        } else {
+          rho[i] += w * mass[type[j]];
+          if (newton || j < nlocal)
+            rho[j] += w * mass[type[i]];
         }
       }
     }

src/RHEO/compute_rheo_rho_sum.h

@@ -39,6 +39,7 @@ class ComputeRHEORhoSum : public Compute {
   class FixRHEO *fix_rheo;
 
  private:
+  int self_mass_flag;
   double cut, cutsq;
 
   class NeighList *list;

src/RHEO/fix_rheo.cpp

@@ -60,6 +60,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   interface_flag = 0;
   surface_flag = 0;
   oxidation_flag = 0;
+  self_mass_flag = 0;
 
   int i;
   int n = atom->ntypes;
@@ -124,6 +125,8 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       interface_flag = 1;
     } else if (strcmp(arg[iarg], "rho/sum") == 0) {
       rhosum_flag = 1;
+    } else if (strcmp(arg[iarg], "self/mass")) {
+      self_mass_flag = 1;
     } else if (strcmp(arg[iarg], "density") == 0) {
       if (iarg + n >= narg) error->all(FLERR, "Illegal rho0 option in fix rheo");
       for (i = 1; i <= n; i++)
@@ -142,6 +145,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     iarg += 1;
   }
 
+  if (self_mass_flag && !rhosum_flag)
+    error->all(FLERR, "Cannot use self/mass setting without rho/sum");
+
   if (lmp->citeme) lmp->citeme->add(cite_rheo);
 }
 
@@ -178,7 +184,7 @@ void FixRHEO::post_constructor()
 
   if (rhosum_flag) {
     compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(
-      "rheo_rhosum all RHEO/RHO/SUM"));
+      fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag)));
     compute_rhosum->fix_rheo = this;
   }
 

src/RHEO/fix_rheo.h

@@ -41,6 +41,7 @@ class FixRHEO : public Fix {
   // Model parameters
   double h, cut;
   double *rho0, *csq;
+  int self_mass_flag;
   int zmin_kernel, zmin_surface, zmin_splash;
   int kernel_style, surface_style;
   double divr_surface;

src/RHEO/fix_rheo_oxidation.cpp

@@ -141,6 +141,10 @@ void FixRHEOOxidation::setup_pre_force(int /*vflag*/)
 
 void FixRHEOOxidation::pre_force(int /*vflag*/)
 {
+  int *status = atom->status;
+  for (int i = 0; i < atom->nlocal; i++)
+    if (num_bond[i] != 0)
+      status[i] |= STATUS_NO_SHIFT;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -215,7 +219,7 @@ void FixRHEOOxidation::post_integrate()
       if (bflag) continue;
 
       // Add bonds to owned atoms
-      // If newton bond, add to both, otherwise add to whichever has a smaller tag
+      // If newton bond off, add to both, otherwise add to whichever has a smaller tag
 
       if (!newton_bond || tagi < tagj) {
         if (num_bond[i] == atom->bond_per_atom)

src/RHEO/fix_rheo_thermal.cpp

@@ -674,7 +674,7 @@ void FixRHEOThermal::create_bonds()
       if (rsq > cutsq_bond) continue;
 
       // Add bonds to owned atoms
-      // If newton bond, add to both, otherwise add to whichever has a smaller tag
+      // If newton bond off, add to both, otherwise add to whichever has a smaller tag
       if (i < nlocal && (!newton_bond || tag[i] < tag[j])) {
         if (num_bond[i] == atom->bond_per_atom)
           error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal");

src/RHEO/fix_rheo_viscosity.cpp

@@ -56,6 +56,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
   memory->create(gd0, n + 1, "rheo:gd0");
   memory->create(K, n + 1, "rheo:K");
   memory->create(npow, n + 1, "rheo:npow");
+  memory->create(tau0, n + 1, "rheo:tau0");
   for (i = 1; i <= n; i++) viscosity_style[i] = NONE;
 
   int iarg = 3;
@@ -76,7 +77,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
         viscosity_style[i] = CONSTANT;
         eta[i] = eta_one;
       }
-    } else if (strcmp(arg[iarg], "power") == 0) {
+    } else if (strcmp(arg[iarg + 1], "power") == 0) {
       if (iarg + 5 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/viscosity power", error);
 
       comm_forward = 1;