Merge pull request lammps#4366 from jtclemm/rheo

Updates for RHEO package, including new optional comm features

doc/src/Developer_comm_ops.rst

@@ -79,19 +79,19 @@ containing ``double`` values.  To correctly store integers that may be
 64-bit (bigint, tagint, imageint) in the buffer, you need to use the
 :ref:`ubuf union <communication_buffer_coding_with_ubuf>` construct.
 
-The *Fix*, *Compute*, and *Dump* classes can also invoke the same kind
-of forward and reverse communication operations using the same *Comm*
-class methods.  Likewise, the same pack/unpack methods and
+The *Fix*, *Bond*, *Compute*, and *Dump* classes can also invoke the
+same kind of forward and reverse communication operations using the
+same *Comm* class methods.  Likewise, the same pack/unpack methods and
 comm_forward/comm_reverse variables must be defined by the calling
-*Fix*, *Compute*, or *Dump* class.
+*Fix*, *Bond*, *Compute*, or *Dump* class.
 
-For *Fix* classes, there is an optional second argument to the
+For all of these classes, there is an optional second argument to the
 *forward_comm()* and *reverse_comm()* call which can be used when the
-fix performs multiple modes of communication, with different numbers
-of values per atom.  The fix should set the *comm_forward* and
+class performs multiple modes of communication, with different numbers
+of values per atom.  The class should set the *comm_forward* and
 *comm_reverse* variables to the maximum value, but can invoke the
 communication for a particular mode with a smaller value.  For this
-to work, the *pack_forward_comm()*, etc methods typically use a class
+to work, the *pack_forward_comm()*, etc. methods typically use a class
 member variable to choose which values to pack/unpack into/from the
 buffer.
 

doc/src/Howto_rheo.rst

@@ -15,8 +15,9 @@ details of the system, or develop new capabilities. For instance, the numerics
 associated with calculating gradients, reproducing kernels, etc. are separated
 into distinct classes to simplify the development of new integration schemes
 which can call these calculations. Additional numerical details can be found in
-:ref:`(Clemmer) <howto_rheo_clemmer>`. Example movies illustrating some of these
-capabilities are found at https://www.lammps.org/movies.html#rheopackage.
+:ref:`(Palermo) <howto_rheo_palermo>` and :ref:`(Clemmer) <howto_rheo_clemmer>`.
+Example movies illustrating some of these capabilities are found at
+https://www.lammps.org/movies.html#rheopackage.
 
 Note, if you simply want to run a traditional SPH simulation, the :ref:`SPH package
 <PKG-SPH>` package is likely better suited for your application. It has fewer advanced
@@ -70,7 +71,7 @@ particles to solid (e.g. with the :doc:`set <set>` command), (b) create bpm
 bonds between the particles (see the :doc:`bpm howto <Howto_bpm>` page for
 more details), and (c) use :doc:`pair rheo/solid <pair_rheo_solid>` to
 apply repulsive contact forces between distinct solid bodies. Akin to pair rheo,
-pair rheo/solid considers a particles fluid/solid phase to determine whether to
+pair rheo/solid considers a particle's fluid/solid phase to determine whether to
 apply forces. However, unlike pair rheo, pair rheo/solid does obey special bond
 settings such that contact forces do not have to be calculated between two bonded
 solid particles in the same elastic body.
@@ -79,10 +80,10 @@ In systems with thermal evolution, fix rheo/thermal can optionally set a
 melting/solidification temperature allowing particles to dynamically swap their
 state between fluid and solid when the temperature exceeds or drops below the
 critical temperature, respectively. Using the *react* option, one can specify a maximum
-bond length and a bond type. Then, when solidifying, particles will search their
+bond length and a bond type. Then, when solidifying, particles search their
 local neighbors and automatically create bonds with any neighboring solid particles
-in range. For BPM bond styles, bonds will then use the immediate position of the two
-particles to calculate a reference state. When melting, particles will delete any
+in range. For BPM bond styles, bonds then use the immediate position of the two
+particles to calculate a reference state. When melting, particles delete any
 bonds of the specified type when reverting to a fluid state. Special bonds are updated
 as bonds are created/broken.
 
@@ -107,6 +108,10 @@ criteria for creating/deleting a bond or altering force calculations).
 
 ----------
 
+.. _howto_rheo_palermo:
+
+**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024).
+
 .. _howto_rheo_clemmer:
 
 **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

doc/src/fix_rheo.rst

@@ -16,21 +16,36 @@ Syntax
 * kstyle = *quintic* or *RK0* or *RK1* or *RK2*
 * 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 *self/mass* or *speed/sound*
+* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *speed/sound*
 
   .. parsed-literal::
 
-       *thermal* values = none, turns on thermal evolution
-       *interface/reconstruct* values = none, reconstructs interfaces with solid particles
-       *surface/detection* values = *sdstyle* *limit* *limit/splash*
-         *sdstyle* = *coordination* or *divergence*
-         *limit* = threshold for surface particles
-         *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)
+       *thermal* turns on thermal evolution
+         values = none
+       *interface/reconstruct* reconstructs interfaces with solid particles
+         values = none
+       *surface/detection* detects free-surfaces with an absence of particles
+         values = *sdstyle* *limit* *limit/splash*
+           *sdstyle* = *coordination* or *divergence*
+           *limit* = threshold for surface particles
+           *limit/splash* = threshold for splash particles (unitless)
+       *shift* turns on velocity shifting
+         values = none
+         optional args = *exclude/type* or *scale/cross/type*
+           *exclude/type* values = *types*
+             *types* = list of types
+           *scale/cross/type* values = *shiftscale* *cmin* *wmin*
+             *shiftscale* = fraction of shifting in normal direction to preserve (unitless)
+             *cmin* = minimum color function value required for scaling (unitless)
+             *wmin* = minimum local same-type support required for any shifting (unitless)
+       *rho/sum* density evolution performed by a kernel summation
+         values = none
+         optional args = *self/mass*
+           *self/mass* values = none, a particle uses its own mass in summation
+       *density* specify equilibrium densities for each atom type
+         values = *rho01*, ... *rho0N* (density)
+       *speed/sound* specify speeds of sound for each atom type
+         values = *cs0*, ... *csN* (velocity)
 
 Examples
 """"""""
@@ -39,15 +54,19 @@ Examples
 
    fix 1 all rheo 3.0 quintic 0 thermal density 0.1 0.1 speed/sound 10.0 1.0
    fix 1 all rheo 3.0 RK1 10 shift surface/detection coordination 40
+   fix 1 all rheo 3.0 RK1 10 shift exclude/type 2*4 scale/cross/type 0.05 0.02 0.5
+   fix 1 all rheo 3.0 RK1 10 rhosum self/mass
 
 Description
 """""""""""
 
 .. versionadded:: 29Aug2024
 
 Perform time integration for RHEO particles, updating positions, velocities,
-and densities. For an overview of other features available in the RHEO package,
-see :doc:`the RHEO howto <Howto_rheo>`.
+and densities. For a detailed breakdown of the integration timestep and
+numerical details, see :ref:`(Palermo) <rheo_palermo>`. For an overview
+and list of other features available in the RHEO package, see
+:doc:`the RHEO howto <Howto_rheo>`.
 
 The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four
 kernels are currently available. The *quintic* kernel is a standard quintic
@@ -70,16 +89,51 @@ and velocity of solid particles are alternatively reconstructed for every
 fluid-solid interaction to ensure no-slip and pressure-balanced boundaries.
 This is done by estimating the location of the fluid-solid interface and
 extrapolating fluid particle properties across the interface to calculate a
-temporary apparent density and velocity for a solid particle.
+temporary apparent density and velocity for a solid particle. The numerical
+details are the same as those described in
+:ref:`(Palermo) <fix_rheo_palermo>` except there is an additional
+restriction that the reconstructed solid density cannot be less than the
+equilibrium density. This prevents fluid particles from sticking to solid
+surfaces.
 
 A modified form of Fickian particle shifting can be enabled with the
 *shift* keyword. This effectively shifts particle positions to generate a
-more uniform spatial distribution. Shifting currently does not consider the
+more uniform spatial distribution. By default, shifting does not consider the
 type of a particle and therefore may be inappropriate in systems consisting
-of multiple fluid phases.
+of multiple atom types representing multiple fluid phases. However, two
+optional subarguments can follow the *shift* keyword, *exclude/type* and
+*scale/cross/type* to adjust shifting at fluid interfaces.
+
+The *exclude/type* option lets the user specify a list of atom types which
+are not shifted, *types*. A wild-card asterisk can be used in place
+of or in conjunction with the *types* argument to toggle shifting for
+multiple atom types.  This takes the form "\*" or "\*n" or "m\*"
+or "m\*n".  If :math:`N` is the number of atom types, then an asterisk with
+no numeric values means all types from 1 to :math:`N`.  A leading asterisk
+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 *scale/cross/type* option is designed to handle interfaces between fluids
+made up of different atom types. Similar to the method by
+:ref:`(Yang) <fix_rheo_yang>`, a color function is calculated and used to
+estimate a local interfacial normal vector. Shifting along this normal direction
+is rescaled by a factor of *scaleshift*, such that a value of *scaleshift* of
+zero implies there is no shifting in the normal direction and a value of
+*scaleshift* of one implies no change in behavior. This scaling is only applied
+to atoms with a color function value greater than *cmin*. To handle scenarios
+of a small inclusion of one fluid type (e.g. a single atom) inside another,
+the degree of same-type support is calculated
 
-In systems with free surfaces, the *surface/detection* keyword can be used
-to classify the location of particles as being within the bulk fluid, on a
+.. math::
+   W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij}
+
+where :math:`\delta_{ij}` is zero if atoms :math:`i` and :math:`j` have different
+types but unity otherwise. If :math:`W_{i,\mathrm{same}}` is ever less than the
+specified value of *wmin*, shifting is turned off for particle :math:`i`
+
+In systems with free surfaces (atom-vacuum), the *surface/detection* keyword
+can classify the location of particles as being within the bulk fluid, on a
 free surface, or isolated from other particles in a splash or droplet.
 Shifting is then disabled in the normal direction away from the free surface
 to prevent particles from diffusing away. Surface detection can also be used
@@ -101,10 +155,9 @@ 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.
-
-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 over neighbors
+keyword. Following this keyword, one may include the optional *self/mass* subargument
+which modifies the behavior of the density summation. Typically, the density
+:math:`\rho` of a particle is calculated as the sum over neighbors
 
 .. math::
    \rho_i = \sum_{j} W_{ij} M_j
@@ -120,7 +173,9 @@ equilibrium density *rho0*.
 
 The *speed/sound* keyword is used to specify the speed of sound of each of the
 N particle types. It must be followed by N numerical values specifying each type's
-speed of sound *cs*.
+speed of sound *cs*. These values may be ignored if the pressure equation of
+state has a non-constant speed of sound, as discussed further in
+:doc:`fix rheo/pressure <fix_rheo_pressure>`.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -163,6 +218,14 @@ Default
 
 ----------
 
+.. _rheo_palermo:
+
+**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024).
+
+.. _rheo_yang:
+
+**(Yang)** Yang, Rakhsha, Hu, Negrut, J. Comp. Physics, 458, 111079 (2022).
+
 .. _fix_rheo_hu:
 
-**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006).
+**(Hu)** Hu, and Adams, J. Comp. Physics, 213, 844-861 (2006).

doc/src/fix_rheo_pressure.rst

@@ -14,13 +14,16 @@ Syntax
 * rheo/pressure = style name of this fix command
 * one or more types and pressure styles must be appended
 * types = lists of types (see below)
-* pstyle = *linear* or *taitwater* or *cubic*
+* pstyle = *linear* or *tait/water* or *tait/general* or *cubic* or *ideal/gas* or *background*
 
   .. parsed-literal::
 
        *linear* args = none
-       *taitwater* args = none
+       *tait/water* args = none
+       *tait/general* args = exponent :math:`gamma` (unitless)
        *cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2)
+       *ideal/gas* args = heat capacity ratio :math:`gamma` (unitless)
+       *background* args = background pressure :math:`P[b]` (pressure)
 
 Examples
 """"""""
@@ -29,6 +32,7 @@ Examples
 
    fix 1 all rheo/pressure * linear
    fix 1 all rheo/pressure 1 linear 2 cubic 10.0
+   fix 1 all rheo/pressure * linear * background 0.1
 
 Description
 """""""""""
@@ -40,21 +44,20 @@ define different equations of state for different atom types. An equation
 must be specified for every atom type.
 
 One first defines the atom *types*. A wild-card asterisk can be used in place
-of or in conjunction with the *types* argument to set the coefficients for
-multiple pairs of atom types.  This takes the form "\*" or "\*n" or "m\*"
-or "m\*n".  If :math:`N` is the number of atom types, then an asterisk with
-no numeric values means all types from 1 to :math:`N`.  A leading asterisk
-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).
+of or in conjunction with the *types* argument to set values for multiple atom
+types.  This takes the form "\*" or "\*n" or "m\*" or "m\*n".  If :math:`N` is
+the number of atom types, then an asterisk with no numeric values means all types
+from 1 to :math:`N`.  A leading asterisk 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 pressure style, *pstyle*. Current
 options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear
 equation of state with a particle pressure :math:`P` calculated as
 
 .. math::
 
-   P = c (\rho - \rho_0)
+   P = c^2 (\rho - \rho_0)
 
 where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density,
 and :math:`\rho` is the current density of a particle. The numerical values of
@@ -63,14 +66,39 @@ is a cubic equation of state which has an extra argument :math:`A_3`,
 
 .. math::
 
-   P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) .
+   P = c^2 ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) .
 
-Style *taitwater* is Tait's equation of state:
+Style *tait/water* is Tait's equation of state:
 
 .. math::
 
    P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr].
 
+Style *tait/general* generalizes this equation of state
+
+.. math::
+
+   P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr].
+
+where :math:`\gamma` is an exponent.
+
+Style *ideal/gas* is the ideal gas equation of state
+
+.. math::
+
+   P = (\gamma - 1) \rho e
+
+where :math:`\gamma` is the heat capacity ratio and :math:`e` is the internal energy of
+a particle per unit mass. This style is only compatible with atom style rheo/thermal.
+Note that when using this style, the speed of sound is no longer constant such that the
+value of :math:`c` specified in :doc:`fix rheo <fix_rheo>` is not used.
+
+The *background* style acts differently than the rest as it
+only adds a constant background pressure shift :math:`P[b]`
+to all atoms of the designated types. Therefore, this style
+must be used in conjunction with another style that specifies
+an equation of state.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

doc/src/fix_rheo_thermal.rst

@@ -70,13 +70,13 @@ of the energy is used to shift energies. This may be inappropriate in systems
 with multiple atom types with different specific heats.
 
 For each property, one must first define a list of atom types. A wild-card
-asterisk can be used in place of or in conjunction with the *types* argument
-to set the coefficients for multiple pairs of atom types.  This takes the
-form "\*" or "\*n" or "m\*" or "m\*n".  If :math:`N` is the number of atom
-types, then an asterisk with no numeric values means all types from 1 to
-:math:`N`.  A leading asterisk 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).
+asterisk can be used in place of or in conjunction with the *types* argument to
+set values for multiple atom types.  This takes the form "\*" or "\*n" or "m\*"
+or "m\*n".  If :math:`N` is the number of atom types, then an asterisk with no
+numeric values means all types from 1 to :math:`N`.  A leading asterisk 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 for each property is followed by the style. Currently,
 the only option is *constant*. Style *constant* simply applies a constant value

doc/src/fix_rheo_viscosity.rst

@@ -45,13 +45,12 @@ viscosities for different atom types, but a viscosity must be specified for
 every atom type.
 
 One first defines the atom *types*. A wild-card asterisk can be used in place
-of or in conjunction with the *types* argument to set the coefficients for
-multiple pairs of atom types.  This takes the form "\*" or "\*n" or "m\*"
-or "m\*n".  If :math:`N` is the number of atom types, then an asterisk with
-no numeric values means all types from 1 to :math:`N`.  A leading asterisk
-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).
+of or in conjunction with the *types* argument to set values for multiple atom
+types.  This takes the form "\*" or "\*n" or "m\*" or "m\*n".  If :math:`N` is
+the number of atom types, then an asterisk with no numeric values means all types
+from 1 to :math:`N`.  A leading asterisk 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*. Two
 options are available, *constant* and *power*. Style *constant* simply