ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-12-10 09:07:57 -05:00
parent 6f4466f300 a78aee5731
commit 02db7be6e2
44 changed files with 2247 additions and 1561 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
doc/src/Developer_comm_ops.rst
View File

@ -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 64-bit (bigint, tagint, imageint) in the buffer, you need to use the
:ref:`ubuf union <communication_buffer_coding_with_ubuf>` construct. :ref:`ubuf union <communication_buffer_coding_with_ubuf>` construct.
The *Fix*, *Compute*, and *Dump* classes can also invoke the same kind The *Fix*, *Bond*, *Compute*, and *Dump* classes can also invoke the
of forward and reverse communication operations using the same *Comm* same kind of forward and reverse communication operations using the
class methods. Likewise, the same pack/unpack methods and same *Comm* class methods. Likewise, the same pack/unpack methods and
comm_forward/comm_reverse variables must be defined by the calling 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 *forward_comm()* and *reverse_comm()* call which can be used when the
fix performs multiple modes of communication, with different numbers class performs multiple modes of communication, with different numbers
of values per atom. The fix should set the *comm_forward* and of values per atom. The class should set the *comm_forward* and
*comm_reverse* variables to the maximum value, but can invoke the *comm_reverse* variables to the maximum value, but can invoke the
communication for a particular mode with a smaller value. For this 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 member variable to choose which values to pack/unpack into/from the
buffer. buffer.

17
doc/src/Howto_rheo.rst
View File

@ -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 associated with calculating gradients, reproducing kernels, etc. are separated
into distinct classes to simplify the development of new integration schemes into distinct classes to simplify the development of new integration schemes
which can call these calculations. Additional numerical details can be found in which can call these calculations. Additional numerical details can be found in
:ref:`(Clemmer) <howto_rheo_clemmer>`. Example movies illustrating some of these :ref:`(Palermo) <howto_rheo_palermo>` and :ref:`(Clemmer) <howto_rheo_clemmer>`.
capabilities are found at https://www.lammps.org/movies.html#rheopackage. 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 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 <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 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 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, 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 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 settings such that contact forces do not have to be calculated between two bonded
solid particles in the same elastic body. 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 melting/solidification temperature allowing particles to dynamically swap their
state between fluid and solid when the temperature exceeds or drops below the 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 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 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 in range. For BPM bond styles, bonds then use the immediate position of the two
particles to calculate a reference state. When melting, particles will delete any 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 bonds of the specified type when reverting to a fluid state. Special bonds are updated
as bonds are created/broken. 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: .. _howto_rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

9
doc/src/bond_bpm_rotational.rst
View File

@ -132,8 +132,8 @@ or :doc:`read_restart <read_restart>` commands:
* :math:`k_b` (force*distance/radians units) * :math:`k_b` (force*distance/radians units)
* :math:`f_{r,c}` (force units) * :math:`f_{r,c}` (force units)
* :math:`f_{s,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_{t,c}` (force*distance units)
* :math:`\tau_{b,c}` (force*distance units)
* :math:`\gamma_n` (force/velocity units) * :math:`\gamma_n` (force/velocity units)
* :math:`\gamma_s` (force/velocity units) * :math:`\gamma_s` (force/velocity units)
* :math:`\gamma_r` (force*distance/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 If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation. during a simulation run. This will prevent some unnecessary calculation.
However, if a bond reaches a damage criterion greater than one, The recommended bond communication distance no longer depends on bond failure
it will trigger an error. 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 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 break during the simulation. Whenever a bond breaks, data is processed

7
doc/src/bond_bpm_spring.rst
View File

@ -117,8 +117,11 @@ page on BPMs.
If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation. during a simulation run. This will prevent some unnecessary calculation.
However, if a bond reaches a strain greater than :math:`\epsilon_c`, The recommended bond communication distance no longer depends on the value of
it will trigger an error. :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 .. versionadded:: TBD

113
doc/src/fix_rheo.rst
View File

@ -16,21 +16,36 @@ Syntax
* kstyle = *quintic* or *RK0* or *RK1* or *RK2* * kstyle = *quintic* or *RK0* or *RK1* or *RK2*
* zmin = minimal number of neighbors for reproducing kernels * zmin = minimal number of neighbors for reproducing kernels
* zero or more keyword/value pairs may be appended to args * 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:: .. parsed-literal::
*thermal* values = none, turns on thermal evolution *thermal* turns on thermal evolution
*interface/reconstruct* values = none, reconstructs interfaces with solid particles values = none
*surface/detection* values = *sdstyle* *limit* *limit/splash* *interface/reconstruct* reconstructs interfaces with solid particles
*sdstyle* = *coordination* or *divergence* values = none
*limit* = threshold for surface particles *surface/detection* detects free-surfaces with an absence of particles
*limit/splash* = threshold for splash particles values = *sdstyle* *limit* *limit/splash*
*shift* values = none, turns on velocity shifting *sdstyle* = *coordination* or *divergence*
*rho/sum* values = none, uses the kernel to compute the density of particles *limit* = threshold for surface particles
*self/mass* values = none, a particle uses its own mass in a rho summation *limit/splash* = threshold for splash particles (unitless)
*density* values = *rho01*, ... *rho0N* (density) *shift* turns on velocity shifting
*speed/sound* values = *cs0*, ... *csN* (velocity) 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 Examples
"""""""" """"""""
@ -39,6 +54,8 @@ 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 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 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 Description
""""""""""" """""""""""
@ -46,8 +63,10 @@ Description
.. versionadded:: 29Aug2024 .. versionadded:: 29Aug2024
Perform time integration for RHEO particles, updating positions, velocities, Perform time integration for RHEO particles, updating positions, velocities,
and densities. For an overview of other features available in the RHEO package, and densities. For a detailed breakdown of the integration timestep and
see :doc:`the RHEO howto <Howto_rheo>`. 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 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 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. fluid-solid interaction to ensure no-slip and pressure-balanced boundaries.
This is done by estimating the location of the fluid-solid interface and This is done by estimating the location of the fluid-solid interface and
extrapolating fluid particle properties across the interface to calculate a 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 A modified form of Fickian particle shifting can be enabled with the
*shift* keyword. This effectively shifts particle positions to generate a *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 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.
In systems with free surfaces, the *surface/detection* keyword can be used The *exclude/type* option lets the user specify a list of atom types which
to classify the location of particles as being within the bulk fluid, on a 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
.. 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. 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 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 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 By default, RHEO integrates particles' densities using a mass diffusion
equation. Alternatively, one can update densities every timestep by performing equation. Alternatively, one can update densities every timestep by performing
a kernel summation of the masses of neighboring particles by specifying the *rho/sum* a kernel summation of the masses of neighboring particles by specifying the *rho/sum*
keyword. keyword. Following this keyword, one may include the optional *self/mass* subargument
which modifies the behavior of the density summation. Typically, the density
The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*. :math:`\rho` of a particle is calculated as the sum over neighbors
Typically, the density :math:`\rho` of a particle is calculated as the sum over neighbors
.. math:: .. math::
\rho_i = \sum_{j} W_{ij} M_j \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 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 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 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: .. _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).

52
doc/src/fix_rheo_pressure.rst
View File

@ -14,13 +14,16 @@ Syntax
* rheo/pressure = style name of this fix command * rheo/pressure = style name of this fix command
* one or more types and pressure styles must be appended * one or more types and pressure styles must be appended
* types = lists of types (see below) * 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:: .. parsed-literal::
*linear* args = none *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) *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 Examples
"""""""" """"""""
@ -29,6 +32,7 @@ Examples
fix 1 all rheo/pressure * linear fix 1 all rheo/pressure * linear
fix 1 all rheo/pressure 1 linear 2 cubic 10.0 fix 1 all rheo/pressure 1 linear 2 cubic 10.0
fix 1 all rheo/pressure * linear * background 0.1
Description Description
""""""""""" """""""""""
@ -40,13 +44,12 @@ define different equations of state for different atom types. An equation
must be specified for every atom type. must be specified for every atom type.
One first defines the atom *types*. A wild-card asterisk can be used in place 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 of or in conjunction with the *types* argument to set values for multiple atom
multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with the number of atom types, then an asterisk with no numeric values means all types
no numeric values means all types from 1 to :math:`N`. A leading asterisk from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive).
means all types from 1 to n (inclusive). A trailing asterisk means all types A trailing asterisk means all types from m to :math:`N` (inclusive). A middle
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n asterisk means all types from m to n (inclusive).
(inclusive).
The *types* definition is followed by the pressure style, *pstyle*. Current The *types* definition is followed by the pressure style, *pstyle*. Current
options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear
@ -54,7 +57,7 @@ equation of state with a particle pressure :math:`P` calculated as
.. math:: .. 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, 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 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:: .. 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:: .. math::
P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr]. 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 Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

14
doc/src/fix_rheo_thermal.rst
View File

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

13
doc/src/fix_rheo_viscosity.rst
View File

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

32
doc/src/pair_granular.rst
View File

@ -235,12 +235,15 @@ given by:
.. math:: .. 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 where :math:`k_{nd}` is an effective harmonic stiffness equal to the ratio of
damping model is not compatible with cohesive normal models such as *JKR* or *DMT*. the normal force to the overlap. For example, :math:`k_{nd} = 4/3Ea` for a
The parameter :math:`\alpha` is related to the restitution coefficient *e* Hertz contact model based on material parameters with :math:`a` being
according to: 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:: .. 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. no error check is performed on this.
The *coeff_restitution* model is useful when a specific normal coefficient of 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 required. It operates much like the *Tsuji* model
restitution :math:`e` is specified as an input. Following the approach of but, the normal coefficient of restitution :math:`e` is specified as an input
:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model, in place of the usual :math:`\eta_{n0}` value in the normal model. Following
*coeff_restitution* calculates the damping coefficient as: the approach of :ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke*
normal model, *coeff_restitution* then calculates the damping coefficient as:
.. math:: .. 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 For any other normal model, e.g. the *hertz* and *hertz/material* models, the damping
coefficient is: coefficient is:
.. math:: .. 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 Since *coeff_restitution* accounts for the effective mass, effective radius,
*coeff_restitution* accounts for the effective mass, effective radius, and and pairwise overlaps (except when used with the *hooke* normal model) when calculating
pairwise overlaps (except when used with the *hooke* normal model) when calculating
the damping coefficient, it accurately reproduces the specified coefficient of the damping coefficient, it accurately reproduces the specified coefficient of
restitution for both monodisperse and polydisperse particle pairs. This damping restitution for both monodisperse and polydisperse particle pairs. This damping
model is not compatible with cohesive normal models such as *JKR* or *DMT*. model is not compatible with cohesive normal models such as *JKR* or *DMT*.

4
examples/bpm/pour/in.bpm.pour
View File

@ -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 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 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 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 nbond all nbond/atom
compute tbond all reduce sum c_nbond compute tbond all reduce sum c_nbond

1117
examples/bpm/pour/log.29Aug2024.bpm.pour.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff

1096
examples/bpm/pour/log.4May2022.bpm.pour.g++.4
View File

File diff suppressed because it is too large Load Diff

9
src/BPM/bond_bpm.cpp
View File

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

5
src/BPM/bond_bpm_rotational.cpp
View File

@ -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, breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1); torque1on2, torque2on1);
if (breaking >= 1.0) { if ((breaking >= 1.0) && break_flag) {
bondlist[n][2] = 0; bondlist[n][2] = 0;
process_broken(i1, i2); process_broken(i1, i2);
continue; continue;
@ -684,6 +684,9 @@ void BondBPMRotational::settings(int narg, char **arg)
error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]); 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");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

5
src/BPM/bond_bpm_spring.cpp
View File

@ -245,7 +245,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
r = sqrt(rsq); r = sqrt(rsq);
e = (r - r0) / r0; e = (r - r0) / r0;
if (fabs(e) > ecrit[type]) { if ((fabs(e) > ecrit[type]) && break_flag) {
bondlist[n][2] = 0; bondlist[n][2] = 0;
process_broken(i1, i2); 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]); 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");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

10
src/EXTRA-FIX/fix_deform_pressure.cpp
View File

@ -407,6 +407,16 @@ void FixDeformPressure::init()
set_box.vol_start = domain->xprd * domain->yprd * domain->zprd; 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 // check optional variables for PRESSURE or PMEAN style
for (int i = 0; i < 7; i++) { for (int i = 0; i < 7; i++) {

17
src/GRANULAR/gran_sub_mod_damping.cpp
View File

@ -152,9 +152,8 @@ double GranSubModDampingTsuji::calculate_forces()
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) : 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); 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;
}

3
src/GRANULAR/gran_sub_mod_damping.h
View File

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

89
src/KOKKOS/comm_kokkos.cpp
View File

@ -378,10 +378,10 @@ void CommKokkos::forward_comm(Fix *fix, int size)
{ {
if (fix->execution_space == Host || !fix->forward_comm_device || forward_fix_comm_classic) { if (fix->execution_space == Host || !fix->forward_comm_device || forward_fix_comm_classic) {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::forward_comm(fix,size); CommBrick::forward_comm(fix, size);
} else { } else {
k_sendlist.sync<LMPDeviceType>(); k_sendlist.sync<LMPDeviceType>();
forward_comm_device<LMPDeviceType>(fix,size); forward_comm_device<LMPDeviceType>(fix, size);
} }
} }
@ -455,10 +455,10 @@ void CommKokkos::forward_comm_device(Fix *fix, int size)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Fix reverse communication invoked by a Fix
size/nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix size = 0 (default) -> use comm_reverse from Fix
size > 0 -> Fix passes max size per atom size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes, the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::reverse_comm(Fix *fix, int size) void CommKokkos::reverse_comm(Fix *fix, int size)
@ -482,72 +482,94 @@ void CommKokkos::reverse_comm_variable(Fix *fix)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Compute forward communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::forward_comm(Compute *compute) void CommKokkos::forward_comm(Compute *compute, int size)
{ {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::forward_comm(compute); CommBrick::forward_comm(compute, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Bond forward communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::forward_comm(Bond *bond) void CommKokkos::forward_comm(Bond *bond, int size)
{ {
CommBrick::forward_comm(bond); CommBrick::forward_comm(bond, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Bond reverse communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::reverse_comm(Bond *bond) void CommKokkos::reverse_comm(Bond *bond, int size)
{ {
CommBrick::reverse_comm(bond); CommBrick::reverse_comm(bond, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Compute reverse communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::reverse_comm(Compute *compute) void CommKokkos::reverse_comm(Compute *compute, int size)
{ {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::reverse_comm(compute); CommBrick::reverse_comm(compute, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Dump forward communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommKokkos::forward_comm(Pair *pair) void CommKokkos::forward_comm(Pair *pair, int size)
{ {
if (pair->execution_space == Host || forward_pair_comm_classic) { if (pair->execution_space == Host || forward_pair_comm_classic) {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::forward_comm(pair); CommBrick::forward_comm(pair, size);
} else { } else {
k_sendlist.sync<LMPDeviceType>(); k_sendlist.sync<LMPDeviceType>();
forward_comm_device<LMPDeviceType>(pair); forward_comm_device<LMPDeviceType>(pair, size);
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
template<class DeviceType> template<class DeviceType>
void CommKokkos::forward_comm_device(Pair *pair) void CommKokkos::forward_comm_device(Pair *pair, int size)
{ {
int iswap,n; int iswap,n,nsize;
MPI_Request request; MPI_Request request;
DAT::tdual_xfloat_1d k_buf_tmp; DAT::tdual_xfloat_1d k_buf_tmp;
int nsize = pair->comm_forward; if (size) nsize = size;
else nsize = pair->comm_forward;
KokkosBase* pairKKBase = dynamic_cast<KokkosBase*>(pair); KokkosBase* pairKKBase = dynamic_cast<KokkosBase*>(pair);
int nmax = max_buf_pair; int nmax = max_buf_pair;
@ -623,29 +645,30 @@ void CommKokkos::grow_buf_fix(int n) {
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void CommKokkos::reverse_comm(Pair *pair) void CommKokkos::reverse_comm(Pair *pair, int size)
{ {
if (pair->execution_space == Host || !pair->reverse_comm_device || reverse_pair_comm_classic) { if (pair->execution_space == Host || !pair->reverse_comm_device || reverse_pair_comm_classic) {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::reverse_comm(pair); CommBrick::reverse_comm(pair, size);
} else { } else {
k_sendlist.sync<LMPDeviceType>(); k_sendlist.sync<LMPDeviceType>();
reverse_comm_device<LMPDeviceType>(pair); reverse_comm_device<LMPDeviceType>(pair, size);
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
template<class DeviceType> template<class DeviceType>
void CommKokkos::reverse_comm_device(Pair *pair) void CommKokkos::reverse_comm_device(Pair *pair, int size)
{ {
int iswap,n; int iswap,n,nsize;
MPI_Request request; MPI_Request request;
DAT::tdual_xfloat_1d k_buf_tmp; DAT::tdual_xfloat_1d k_buf_tmp;
KokkosBase* pairKKBase = dynamic_cast<KokkosBase*>(pair); KokkosBase* pairKKBase = dynamic_cast<KokkosBase*>(pair);
int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); if (size) nsize = size;
else nsize = MAX(pair->comm_reverse, pair->comm_reverse_off);
int nmax = max_buf_pair; int nmax = max_buf_pair;
for (iswap = 0; iswap < nswap; iswap++) { for (iswap = 0; iswap < nswap; iswap++) {
@ -702,18 +725,18 @@ void CommKokkos::reverse_comm_device(Pair *pair)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void CommKokkos::forward_comm(Dump *dump) void CommKokkos::forward_comm(Dump *dump, int size)
{ {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::forward_comm(dump); CommBrick::forward_comm(dump, size);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void CommKokkos::reverse_comm(Dump *dump) void CommKokkos::reverse_comm(Dump *dump, int size)
{ {
k_sendlist.sync<LMPHostType>(); k_sendlist.sync<LMPHostType>();
CommBrick::reverse_comm(dump); CommBrick::reverse_comm(dump, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

26
src/KOKKOS/comm_kokkos.h
View File

@ -45,24 +45,24 @@ class CommKokkos : public CommBrick {
void exchange() override; // move atoms to new procs void exchange() override; // move atoms to new procs
void borders() override; // setup list of atoms to comm void borders() override; // setup list of atoms to comm
void forward_comm(class Pair *) override; // forward comm from a Pair void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair
void reverse_comm(class Pair *) override; // reverse comm from a Pair void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair
void forward_comm(class Bond *) override; // forward comm from a Bond void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond
void reverse_comm(class Bond *) override; // reverse comm from a Bond void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond
void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix
void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix
void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix
void forward_comm(class Compute *) override; // forward from a Compute void forward_comm(class Compute *, int size = 0) override; // forward from a Compute
void reverse_comm(class Compute *) override; // reverse from a Compute void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute
void forward_comm(class Dump *) override; // forward comm from a Dump void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump
void reverse_comm(class Dump *) override; // reverse comm from a Dump void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump
void forward_comm_array(int, double **) override; // forward comm of array void forward_comm_array(int, double **) override; // forward comm of array
template<class DeviceType> void forward_comm_device(); template<class DeviceType> void forward_comm_device();
template<class DeviceType> void reverse_comm_device(); template<class DeviceType> void reverse_comm_device();
template<class DeviceType> void forward_comm_device(Pair *pair); template<class DeviceType> void forward_comm_device(Pair *pair, int size=0);
template<class DeviceType> void reverse_comm_device(Pair *pair); template<class DeviceType> void reverse_comm_device(Pair *pair, int size=0);
template<class DeviceType> void forward_comm_device(Fix *fix, int size=0); template<class DeviceType> void forward_comm_device(Fix *fix, int size=0);
template<class DeviceType> void exchange_device(); template<class DeviceType> void exchange_device();
template<class DeviceType> void borders_device(); template<class DeviceType> void borders_device();

88
src/KOKKOS/comm_tiled_kokkos.cpp
View File

@ -417,42 +417,58 @@ void CommTiledKokkos::borders()
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Pair forward communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::forward_comm(Pair *pair) void CommTiledKokkos::forward_comm(Pair *pair, int size)
{ {
CommTiled::forward_comm(pair); CommTiled::forward_comm(pair, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Pair reverse communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Pair
size > 0 -> Pair passes max size per atom
the latter is only useful if Pair does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::reverse_comm(Pair *pair) void CommTiledKokkos::reverse_comm(Pair *pair, int size)
{ {
CommTiled::reverse_comm(pair); CommTiled::reverse_comm(pair, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Bond forward communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::forward_comm(Bond *bond) void CommTiledKokkos::forward_comm(Bond *bond, int size)
{ {
CommTiled::forward_comm(bond); CommTiled::forward_comm(bond, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Bond reverse communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::reverse_comm(Bond *bond) void CommTiledKokkos::reverse_comm(Bond *bond, int size)
{ {
CommTiled::reverse_comm(bond); CommTiled::reverse_comm(bond, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -466,21 +482,21 @@ void CommTiledKokkos::reverse_comm(Bond *bond)
void CommTiledKokkos::forward_comm(Fix *fix, int size) void CommTiledKokkos::forward_comm(Fix *fix, int size)
{ {
CommTiled::forward_comm(fix,size); CommTiled::forward_comm(fix, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Fix reverse communication invoked by a Fix
size/nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix size = 0 (default) -> use comm_reverse from Fix
size > 0 -> Fix passes max size per atom size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes, the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::reverse_comm(Fix *fix, int size) void CommTiledKokkos::reverse_comm(Fix *fix, int size)
{ {
CommTiled::reverse_comm(fix,size); CommTiled::reverse_comm(fix, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -497,42 +513,58 @@ void CommTiledKokkos::reverse_comm_variable(Fix *fix)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Compute forward communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::forward_comm(Compute *compute) void CommTiledKokkos::forward_comm(Compute *compute, int size)
{ {
CommTiled::forward_comm(compute); CommTiled::forward_comm(compute, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Compute reverse communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::reverse_comm(Compute *compute) void CommTiledKokkos::reverse_comm(Compute *compute, int size)
{ {
CommTiled::reverse_comm(compute); CommTiled::reverse_comm(compute, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Dump forward communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::forward_comm(Dump *dump) void CommTiledKokkos::forward_comm(Dump *dump, int size)
{ {
CommTiled::forward_comm(dump); CommTiled::forward_comm(dump, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Dump reverse communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiledKokkos::reverse_comm(Dump *dump) void CommTiledKokkos::reverse_comm(Dump *dump, int size)
{ {
CommTiled::reverse_comm(dump); CommTiled::reverse_comm(dump, size);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

22
src/KOKKOS/comm_tiled_kokkos.h
View File

@ -46,17 +46,17 @@ class CommTiledKokkos : public CommTiled {
void exchange() override; // move atoms to new procs void exchange() override; // move atoms to new procs
void borders() override; // setup list of atoms to comm void borders() override; // setup list of atoms to comm
void forward_comm(class Pair *) override; // forward comm from a Pair void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair
void reverse_comm(class Pair *) override; // reverse comm from a Pair void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair
void forward_comm(class Bond *) override; // forward comm from a Bond void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond
void reverse_comm(class Bond *) override; // reverse comm from a Bond void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond
void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix
void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix
void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix
void forward_comm(class Compute *) override; // forward from a Compute void forward_comm(class Compute *, int size = 0) override; // forward from a Compute
void reverse_comm(class Compute *) override; // reverse from a Compute void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute
void forward_comm(class Dump *) override; // forward comm from a Dump void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump
void reverse_comm(class Dump *) override; // reverse comm from a Dump void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump
void forward_comm_array(int, double **) override; // forward comm of array void forward_comm_array(int, double **) override; // forward comm of array

6
src/RHEO/compute_rheo_grad.cpp
View File

@ -290,8 +290,7 @@ void ComputeRHEOGrad::compute_peratom()
void ComputeRHEOGrad::forward_gradients() void ComputeRHEOGrad::forward_gradients()
{ {
comm_stage = COMMGRAD; comm_stage = COMMGRAD;
comm_forward = ncomm_grad; comm->forward_comm(this, ncomm_grad);
comm->forward_comm(this);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -299,8 +298,7 @@ void ComputeRHEOGrad::forward_gradients()
void ComputeRHEOGrad::forward_fields() void ComputeRHEOGrad::forward_fields()
{ {
comm_stage = COMMFIELD; comm_stage = COMMFIELD;
comm_forward = ncomm_field; comm->forward_comm(this, ncomm_field);
comm->forward_comm(this);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

15
src/RHEO/compute_rheo_interface.cpp
View File

@ -97,6 +97,9 @@ void ComputeRHEOInterface::init()
auto fixes = modify->get_fix_by_style("rheo/pressure"); auto fixes = modify->get_fix_by_style("rheo/pressure");
fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]); fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
if (!fix_pressure->invertible_pressure)
error->all(FLERR, "RHEO interface reconstruction incompatible with pressure equation of state");
neighbor->add_request(this, NeighConst::REQ_DEFAULT); neighbor->add_request(this, NeighConst::REQ_DEFAULT);
} }
@ -178,7 +181,7 @@ void ComputeRHEOInterface::compute_peratom()
dot = 0; dot = 0;
for (a = 0; a < 3; a++) dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a]; for (a = 0; a < 3; a++) dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a];
rho[i] += w * (fix_pressure->calc_pressure(rho[j], jtype) - rho[j] * dot); rho[i] += w * (fix_pressure->calc_pressure(rho[j], j) - rho[j] * dot);
normwf[i] += w; normwf[i] += w;
} }
} }
@ -192,7 +195,7 @@ void ComputeRHEOInterface::compute_peratom()
dot = 0; dot = 0;
for (a = 0; a < 3; a++) dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a]; for (a = 0; a < 3; a++) dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a];
rho[j] += w * (fix_pressure->calc_pressure(rho[i], itype) + rho[i] * dot); rho[j] += w * (fix_pressure->calc_pressure(rho[i], i) + rho[i] * dot);
normwf[j] += w; normwf[j] += w;
} }
} }
@ -210,7 +213,7 @@ void ComputeRHEOInterface::compute_peratom()
if (status[i] & PHASECHECK) { if (status[i] & PHASECHECK) {
if (normwf[i] != 0.0) { if (normwf[i] != 0.0) {
// Stores rho for solid particles 1+Pw in Adami Adams 2012 // Stores rho for solid particles 1+Pw in Adami Adams 2012
rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], type[i])); rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], i));
} else { } else {
rho[i] = rho0[itype]; rho[i] = rho0[itype];
} }
@ -218,8 +221,7 @@ void ComputeRHEOInterface::compute_peratom()
} }
comm_stage = 1; comm_stage = 1;
comm_forward = 2; comm->forward_comm(this, 2);
comm->forward_comm(this);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -369,9 +371,8 @@ void ComputeRHEOInterface::store_forces()
} }
// Forward comm forces // Forward comm forces
comm_forward = 3;
comm_stage = 0; comm_stage = 0;
comm->forward_comm(this); comm->forward_comm(this, 3);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

142
src/RHEO/compute_rheo_kernel.cpp
View File

@ -90,7 +90,6 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
comm_forward = ncor * Mdim; comm_forward = ncor * Mdim;
} }
comm_forward_save = comm_forward;
corrections_calculated = 0; corrections_calculated = 0;
lapack_error_flag = 0; lapack_error_flag = 0;
} }
@ -137,6 +136,9 @@ void ComputeRHEOKernel::init()
} }
} }
if (correction_order != -1)
fix_rheo->coordination_flag = 1;
nmax_store = atom->nmax; nmax_store = atom->nmax;
memory->create(coordination, nmax_store, "rheo:coordination"); memory->create(coordination, nmax_store, "rheo:coordination");
if (kernel_style == RK0) { if (kernel_style == RK0) {
@ -186,27 +188,26 @@ double ComputeRHEOKernel::calc_w_self()
double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r) double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r)
{ {
if (kernel_style == WENDLANDC4)
return calc_w_wendlandc4(r);
if (kernel_style == QUINTIC)
return calc_w_quintic(r);
double w = 0.0; double w = 0.0;
int corrections_i, corrections_j, corrections; int corrections_i = check_corrections(i);
int corrections_j = check_corrections(j);
if (kernel_style == WENDLANDC4) return calc_w_wendlandc4(r); int corrections = corrections_i && corrections_j;
if (kernel_style != QUINTIC) {
corrections_i = check_corrections(i);
corrections_j = check_corrections(j);
corrections = corrections_i & corrections_j;
} else {
corrections = 0;
}
if (!corrections) if (!corrections)
w = calc_w_quintic(r); return calc_w_quintic(r);
else if (kernel_style == RK0)
double dx[3] = {delx, dely, delz};
if (kernel_style == RK0)
w = calc_w_rk0(i, j, r); w = calc_w_rk0(i, j, r);
else if (kernel_style == RK1) else if (kernel_style == RK1)
w = calc_w_rk1(i, j, delx, dely, delz, r); w = calc_w_rk1(i, j, dx, r);
else if (kernel_style == RK2) else if (kernel_style == RK2)
w = calc_w_rk2(i, j, delx, dely, delz, r); w = calc_w_rk2(i, j, dx, r);
return w; return w;
} }
@ -215,26 +216,27 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double delz, double r) double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double delz, double r)
{ {
if (kernel_style == WENDLANDC4)
return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji);
if (kernel_style == QUINTIC)
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
double wp; double wp;
int corrections_i, corrections_j; int corrections_i = check_corrections(i);
int corrections_j = check_corrections(j);
if (kernel_style == WENDLANDC4) return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji); wp = calc_dw_scalar_quintic(delx, dely, delz, r);
if (kernel_style != QUINTIC) {
corrections_i = check_corrections(i);
corrections_j = check_corrections(j);
}
// Calc wp and default dW's, a bit inefficient but can redo later
wp = calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
// Overwrite if there are corrections // Overwrite if there are corrections
double dxij[3] = {delx, dely, delz};
double dxji[3] = {-delx, -dely, -delz};
if (kernel_style == RK1) { if (kernel_style == RK1) {
if (corrections_i) calc_dw_rk1(i, delx, dely, delz, r, dWij); if (corrections_i) calc_dw_rk1(i, dxij, r, dWij);
if (corrections_j) calc_dw_rk1(j, -delx, -dely, -delz, r, dWji); if (corrections_j) calc_dw_rk1(j, dxji, r, dWji);
} else if (kernel_style == RK2) { } else if (kernel_style == RK2) {
if (corrections_i) calc_dw_rk2(i, delx, dely, delz, r, dWij); if (corrections_i) calc_dw_rk2(i, dxij, r, dWij);
if (corrections_j) calc_dw_rk2(j, -delx, -dely, -delz, r, dWji); if (corrections_j) calc_dw_rk2(j, dxji, r, dWji);
} }
return wp; return wp;
@ -275,10 +277,9 @@ double ComputeRHEOKernel::calc_w_quintic(double r)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r, double ComputeRHEOKernel::calc_dw_scalar_quintic(double delx, double dely, double delz, double r)
double *dW1, double *dW2)
{ {
double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv; double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s;
s = r * 3.0 * cutinv; s = r * 3.0 * cutinv;
@ -300,14 +301,23 @@ double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz,
} }
wp *= pre_wp; wp *= pre_wp;
wprinv = wp / r;
return wp;
}
/* ---------------------------------------------------------------------- */
double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r,
double *dW1, double *dW2)
{
double wp = calc_dw_scalar_quintic(delx, dely, delz, r);
double wprinv = wp / r;
dW1[0] = delx * wprinv; dW1[0] = delx * wprinv;
dW1[1] = dely * wprinv; dW1[1] = dely * wprinv;
dW1[2] = delz * wprinv; dW1[2] = delz * wprinv;
dW2[0] = -delx * wprinv; scale3(-1.0, dW1, dW2);
dW2[1] = -dely * wprinv;
dW2[2] = -delz * wprinv;
return wp; return wp;
} }
@ -361,9 +371,7 @@ double ComputeRHEOKernel::calc_dw_wendlandc4(double delx, double dely, double de
dW1[1] = dely * wprinv; dW1[1] = dely * wprinv;
dW1[2] = delz * wprinv; dW1[2] = delz * wprinv;
dW2[0] = -delx * wprinv; scale3(-1.0, dW1, dW2);
dW2[1] = -dely * wprinv;
dW2[2] = -delz * wprinv;
return wp; return wp;
} }
@ -384,14 +392,11 @@ double ComputeRHEOKernel::calc_w_rk0(int i, int j, double r)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, double delz, double r) double ComputeRHEOKernel::calc_w_rk1(int i, int j, double *dx, double r)
{ {
int b; int b;
double w, dx[3], H[MAX_MDIM]; double w, H[MAX_MDIM];
dx[0] = delx;
dx[1] = dely;
dx[2] = delz;
w = calc_w_quintic(r); w = calc_w_quintic(r);
if (dim == 2) { if (dim == 2) {
@ -426,13 +431,11 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, double delz, double r) double ComputeRHEOKernel::calc_w_rk2(int i, int j, double *dx, double r)
{ {
int b; int b;
double w, dx[3], H[MAX_MDIM]; double w, H[MAX_MDIM];
dx[0] = delx;
dx[1] = dely;
dx[2] = delz;
w = calc_w_quintic(r); w = calc_w_quintic(r);
if (dim == 2) { if (dim == 2) {
@ -476,14 +479,10 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz, double r, void ComputeRHEOKernel::calc_dw_rk1(int i, double *dx, double r, double *dW)
double *dW)
{ {
int a, b; int a, b;
double w, dx[3], H[MAX_MDIM]; double w, H[MAX_MDIM];
dx[0] = delx;
dx[1] = dely;
dx[2] = delz;
w = calc_w_quintic(r); w = calc_w_quintic(r);
@ -501,8 +500,8 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz
// dWij[] = dWx dWy (dWz) // dWij[] = dWx dWy (dWz)
//compute derivative operators //compute derivative operators
zero3(dW);
for (a = 0; a < dim; a++) { for (a = 0; a < dim; a++) {
dW[a] = 0.0;
for (b = 0; b < Mdim; b++) { for (b = 0; b < Mdim; b++) {
//First derivative kernels //First derivative kernels
dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z)
@ -513,14 +512,10 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz, double r, void ComputeRHEOKernel::calc_dw_rk2(int i, double *dx, double r, double *dW)
double *dW)
{ {
int a, b; int a, b;
double w, dx[3], H[MAX_MDIM]; double w, H[MAX_MDIM];
dx[0] = delx;
dx[1] = dely;
dx[2] = delz;
w = calc_w_quintic(r); w = calc_w_quintic(r);
@ -547,8 +542,8 @@ void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz
// dWij[] = dWx dWy (dWz) // dWij[] = dWx dWy (dWz)
//compute derivative operators //compute derivative operators
zero3(dW);
for (a = 0; a < dim; a++) { for (a = 0; a < dim; a++) {
dW[a] = 0.0;
for (b = 0; b < Mdim; b++) { for (b = 0; b < Mdim; b++) {
//First derivative kernels //First derivative kernels
dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz)
@ -564,7 +559,7 @@ void ComputeRHEOKernel::compute_peratom()
lapack_error_flag = 0; lapack_error_flag = 0;
lapack_error_tags.clear(); lapack_error_tags.clear();
if (kernel_style == QUINTIC) return; if (correction_order == -1) return;
corrections_calculated = 1; corrections_calculated = 1;
int i, j, ii, jj, inum, jnum, a, b, lapack_error; int i, j, ii, jj, inum, jnum, a, b, lapack_error;
@ -666,6 +661,7 @@ void ComputeRHEOKernel::compute_peratom()
w = calc_w_quintic(r); w = calc_w_quintic(r);
rhoj = rho[j]; rhoj = rho[j];
if (interface_flag) if (interface_flag)
if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j); if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j);
@ -781,7 +777,6 @@ void ComputeRHEOKernel::compute_peratom()
// communicate calculated quantities // communicate calculated quantities
comm_stage = 1; comm_stage = 1;
comm_forward = comm_forward_save;
comm->forward_comm(this); comm->forward_comm(this);
} }
@ -829,8 +824,7 @@ void ComputeRHEOKernel::compute_coordination()
// communicate calculated quantities // communicate calculated quantities
comm_stage = 0; comm_stage = 0;
comm_forward = 1; comm->forward_comm(this, 1);
comm->forward_comm(this);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -853,7 +847,9 @@ void ComputeRHEOKernel::grow_arrays(int nmax)
int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
int * /*pbc*/) int * /*pbc*/)
{ {
int a, b;
int m = 0; int m = 0;
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
int j = list[i]; int j = list[i];
if (comm_stage == 0) { if (comm_stage == 0) {
@ -862,8 +858,8 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb
if (kernel_style == RK0) { if (kernel_style == RK0) {
buf[m++] = C0[j]; buf[m++] = C0[j];
} else { } else {
for (int a = 0; a < ncor; a++) for (a = 0; a < ncor; a++)
for (int b = 0; b < Mdim; b++) buf[m++] = C[j][a][b]; for (b = 0; b < Mdim; b++) buf[m++] = C[j][a][b];
} }
} }
} }
@ -874,8 +870,10 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb
void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
{ {
int a, b;
int m = 0; int m = 0;
int last = first + n; int last = first + n;
for (int i = first; i < last; i++) { for (int i = first; i < last; i++) {
if (comm_stage == 0) { if (comm_stage == 0) {
coordination[i] = buf[m++]; coordination[i] = buf[m++];
@ -883,8 +881,8 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
if (kernel_style == RK0) { if (kernel_style == RK0) {
C0[i] = buf[m++]; C0[i] = buf[m++];
} else { } else {
for (int a = 0; a < ncor; a++) for (a = 0; a < ncor; a++)
for (int b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; for (b = 0; b < Mdim; b++) C[i][a][b] = buf[m++];
} }
} }
} }

11
src/RHEO/compute_rheo_kernel.h
View File

@ -40,6 +40,7 @@ class ComputeRHEOKernel : public Compute {
double calc_w(int, int, double, double, double, double); double calc_w(int, int, double, double, double, double);
double calc_dw(int, int, double, double, double, double); double calc_dw(int, int, double, double, double, double);
double calc_w_quintic(double); double calc_w_quintic(double);
double calc_dw_scalar_quintic(double, double, double, double);
double calc_dw_quintic(double, double, double, double, double *, double *); double calc_dw_quintic(double, double, double, double, double *, double *);
double calc_w_wendlandc4(double); double calc_w_wendlandc4(double);
double calc_dw_wendlandc4(double, double, double, double, double *, double *); double calc_dw_wendlandc4(double, double, double, double, double *, double *);
@ -51,7 +52,7 @@ class ComputeRHEOKernel : public Compute {
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
int comm_stage, comm_forward_save; int comm_stage;
int interface_flag; int interface_flag;
int lapack_error_flag; int lapack_error_flag;
std::unordered_set<tagint> lapack_error_tags; std::unordered_set<tagint> lapack_error_tags;
@ -69,10 +70,10 @@ class ComputeRHEOKernel : public Compute {
int check_corrections(int); int check_corrections(int);
double calc_w_rk0(int, int, double); double calc_w_rk0(int, int, double);
double calc_w_rk1(int, int, double, double, double, double); double calc_w_rk1(int, int, double*, double);
double calc_w_rk2(int, int, double, double, double, double); double calc_w_rk2(int, int, double*, double);
void calc_dw_rk1(int, double, double, double, double, double *); void calc_dw_rk1(int, double*, double, double*);
void calc_dw_rk2(int, double, double, double, double, double *); void calc_dw_rk2(int, double*, double, double*);
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
#endif #endif

12
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -79,7 +79,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
size_peratom_cols = nvalues; size_peratom_cols = nvalues;
pressure_flag = thermal_flag = interface_flag = 0; pressure_flag = thermal_flag = interface_flag = 0;
surface_flag = shift_flag = shell_flag = 0; surface_flag = shift_flag = shell_flag = coordination_flag = 0;
// parse input values // parse input values
// customize a new keyword by adding to if statement // customize a new keyword by adding to if statement
@ -109,6 +109,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
surface_flag = 1; surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
} else if (strcmp(arg[iarg], "coordination") == 0) { } else if (strcmp(arg[iarg], "coordination") == 0) {
coordination_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
} else if (strcmp(arg[iarg], "pressure") == 0) { } else if (strcmp(arg[iarg], "pressure") == 0) {
pressure_flag = 1; pressure_flag = 1;
@ -223,6 +224,11 @@ void ComputeRHEOPropertyAtom::compute_peratom()
{ {
invoked_peratom = update->ntimestep; invoked_peratom = update->ntimestep;
// calculate optional values, if needed
if (coordination_flag && !(fix_rheo->coordination_flag))
compute_kernel->compute_coordination();
// grow vector or array if necessary // grow vector or array if necessary
if (atom->nmax > nmax) { if (atom->nmax > nmax) {
@ -434,7 +440,7 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n)
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) if (mask[i] & groupbit)
buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); buf[n] = fix_pressure->calc_pressure(rho[i], i);
else else
buf[n] = 0.0; buf[n] = 0.0;
n += nvalues; n += nvalues;
@ -484,7 +490,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n)
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (index == index_transpose) if (index == index_transpose)
p = fix_pressure->calc_pressure(rho[i], type[i]); p = fix_pressure->calc_pressure(rho[i], i);
else else
p = 0.0; p = 0.0;
buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p; buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p;

2
src/RHEO/compute_rheo_property_atom.h
View File

@ -36,7 +36,7 @@ class ComputeRHEOPropertyAtom : public Compute {
private: private:
int nvalues, nmax; int nvalues, nmax;
int pressure_flag, thermal_flag, interface_flag; int pressure_flag, thermal_flag, interface_flag;
int surface_flag, shift_flag, shell_flag; int surface_flag, shift_flag, shell_flag, coordination_flag;
int *avec_index; int *avec_index;
int *col_index, *col_t_index; int *col_index, *col_t_index;
double *buf; double *buf;

14
src/RHEO/compute_rheo_surface.cpp
View File

@ -81,6 +81,8 @@ void ComputeRHEOSurface::init()
threshold_splash = fix_rheo->zmin_splash; threshold_splash = fix_rheo->zmin_splash;
interface_flag = fix_rheo->interface_flag; interface_flag = fix_rheo->interface_flag;
fix_rheo->coordination_flag = 1;
cutsq = cut * cut; cutsq = cut * cut;
// need an occasional half neighbor list // need an occasional half neighbor list
@ -197,10 +199,8 @@ void ComputeRHEOSurface::compute_peratom()
// reverse gradC and divr, forward divr // reverse gradC and divr, forward divr
comm_stage = 0; comm_stage = 0;
comm_reverse = dim * dim + 1; if (newton) comm->reverse_comm(this, dim * dim + 1);
comm_forward = 1; comm->forward_comm(this, 1);
if (newton) comm->reverse_comm(this);
comm->forward_comm(this);
// calculate nsurface for local atoms // calculate nsurface for local atoms
// Note, this isn't forwarded to ghosts // Note, this isn't forwarded to ghosts
@ -297,10 +297,8 @@ void ComputeRHEOSurface::compute_peratom()
// forward/reverse status and rsurface // forward/reverse status and rsurface
comm_stage = 1; comm_stage = 1;
comm_reverse = 2; if (newton) comm->reverse_comm(this, 2);
comm_forward = 2; comm->forward_comm(this, 2);
if (newton) comm->reverse_comm(this);
comm->forward_comm(this);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

345
src/RHEO/compute_rheo_vshift.cpp
View File

@ -27,6 +27,7 @@
#include "error.h" #include "error.h"
#include "fix_rheo.h" #include "fix_rheo.h"
#include "force.h" #include "force.h"
#include "math_extra.h"
#include "memory.h" #include "memory.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
@ -36,20 +37,22 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS; using namespace RHEO_NS;
using namespace MathExtra;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr), Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(nullptr),
compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) fix_rheo(nullptr), rho0(nullptr), list(nullptr), compute_interface(nullptr),
compute_kernel(nullptr), compute_surface(nullptr)
{ {
if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command"); if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command");
comm_forward = 0;
comm_reverse = 3; comm_reverse = 3;
surface_flag = 0; surface_flag = 0;
nmax_store = atom->nmax; nmax_store = 0;
memory->create(vshift, nmax_store, 3, "rheo:vshift");
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -73,9 +76,19 @@ void ComputeRHEOVShift::init()
compute_surface = fix_rheo->compute_surface; compute_surface = fix_rheo->compute_surface;
rho0 = fix_rheo->rho0; rho0 = fix_rheo->rho0;
shift_type = fix_rheo->shift_type;
cut = fix_rheo->cut; cut = fix_rheo->cut;
cutsq = cut * cut; cutsq = cut * cut;
cutthird = cut / 3.0; cutthird = cut / 3.0;
cross_type_flag = fix_rheo->shift_cross_type_flag;
if (cross_type_flag) {
scale = fix_rheo->shift_scale;
wmin = fix_rheo->shift_wmin;
cmin = fix_rheo->shift_cmin;
comm_forward = 1;
comm_reverse = 4;
}
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -120,6 +133,13 @@ void ComputeRHEOVShift::compute_peratom()
if (nmax_store < atom->nmax) { if (nmax_store < atom->nmax) {
memory->grow(vshift, atom->nmax, 3, "rheo:vshift"); memory->grow(vshift, atom->nmax, 3, "rheo:vshift");
if (cross_type_flag) {
memory->grow(ct, atom->nmax, "rheo:ct");
memory->grow(cgradt, atom->nmax, 3, "rheo:cgradt");
memory->grow(wsame, atom->nmax, "rheo:wsame");
}
nmax_store = atom->nmax; nmax_store = atom->nmax;
} }
@ -224,7 +244,17 @@ void ComputeRHEOVShift::compute_peratom()
} }
} }
if (newton_pair) comm->reverse_comm(this); comm_stage = 0;
if (newton_pair) comm->reverse_comm(this, 3);
// Zero any excluded types
for (i = 0; i < nlocal; i++)
if (!shift_type[type[i]])
for (a = 0; a < dim; a++)
vshift[i][a] = 0.0;
if (cross_type_flag) correct_type_interface();
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -246,7 +276,6 @@ void ComputeRHEOVShift::correct_surfaces()
if (status[i] & PHASECHECK) continue; if (status[i] & PHASECHECK) continue;
//if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
if (status[i] & STATUS_SURFACE) { if (status[i] & STATUS_SURFACE) {
nx = nsurface[i][0]; nx = nsurface[i][0];
ny = nsurface[i][1]; ny = nsurface[i][1];
@ -283,16 +312,278 @@ void ComputeRHEOVShift::correct_surfaces()
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOVShift::correct_type_interface()
{ {
int m, last; int i, j, a, ii, jj, jnum, itype, jtype;
int fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, r, rinv;
double w, wp, dr, w0, prefactor;
double imass, jmass, voli, volj, rhoi, rhoj;
double dx[3];
int dim = domain->dimension;
int *jlist;
int inum, *ilist, *numneigh, **firstneigh;
int *type = atom->type;
int *status = atom->rheo_status;
int *mask = atom->mask;
double **x = atom->x;
double *rho = atom->rho;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
size_t nbytes = nmax_store * sizeof(double);
memset(&ct[0], 0, nbytes);
memset(&wsame[0], 0, nbytes);
memset(&cgradt[0][0], 0, 3 * nbytes);
double ctmp, *dWij, *dWji;
// Calculate color gradient
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
fluidi = !(status[i] & PHASECHECK);
jlist = firstneigh[i];
jnum = numneigh[i];
if (rmass)
imass = rmass[i];
else
imass = mass[itype];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2];
rsq = lensq3(dx);
if (rsq > cutsq) continue;
fluidj = !(status[j] & PHASECHECK);
jtype = type[j];
if (rmass)
jmass = rmass[j];
else
jmass = mass[jtype];
r = sqrt(rsq);
rhoi = rho[i];
rhoj = rho[j];
// Add corrections for walls
if (interface_flag) {
if (fluidi && (!fluidj)) {
rhoj = compute_interface->correct_rho(j);
} else if ((!fluidi) && fluidj) {
rhoi = compute_interface->correct_rho(i);
} else if ((!fluidi) && (!fluidj)) {
rhoi = rho0[itype];
rhoj = rho0[jtype];
}
}
voli = imass / rhoi;
volj = jmass / rhoj;
w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
if (itype != jtype) ctmp = 1;
else ctmp = 0;
ct[i] += ctmp * volj * w;
if (newton_pair || j < nlocal)
ct[j] += ctmp * voli * w;
}
}
comm_stage = 1;
if (newton_pair) comm->reverse_comm(this, 1);
// Calculate color gradient
// Note: in future might want to generalize this so color function can be used
// by other calculations (e.g. surface tension)
// maybe can create custom "calc_grad" method that takes an arbitrary field
// in ComputeRHEOGrad?
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
fluidi = !(status[i] & PHASECHECK);
jlist = firstneigh[i];
jnum = numneigh[i];
imass = mass[itype];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2];
rsq = lensq3(dx);
if (rsq > cutsq) continue;
fluidj = !(status[j] & PHASECHECK);
jtype = type[j];
if (rmass)
jmass = rmass[j];
else
jmass = mass[jtype];
r = sqrt(rsq);
rhoi = rho[i];
rhoj = rho[j];
// Add corrections for walls
if (interface_flag) {
if (fluidi && (!fluidj)) {
rhoj = compute_interface->correct_rho(j);
} else if ((!fluidi) && fluidj) {
rhoi = compute_interface->correct_rho(i);
} else if ((!fluidi) && (!fluidj)) {
rhoi = rho0[itype];
rhoj = rho0[jtype];
}
}
voli = imass / rhoi;
volj = jmass / rhoj;
w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
dWij = compute_kernel->dWij;
dWji = compute_kernel->dWji;
if (itype != jtype) ctmp = 1;
else ctmp = 0;
for (a = 0; a < dim; a++) {
cgradt[i][a] -= ctmp * volj * dWij[a];
if (newton_pair || j < nlocal)
cgradt[j][a] -= ctmp * voli * dWji[a];
}
if (itype == jtype) {
wsame[i] += w * r;
if (newton_pair || j < nlocal)
wsame[j] += w * r;
}
}
}
comm_stage = 2;
if (newton_pair) comm->reverse_comm(this, 4);
comm->forward_comm(this, 1);
// Correct shifting at fluid-fluid interface
// remove normal shifting component for interfacial particles
// Based on Yang, Rakhsha, Hu, & Negrut 2022
double ntmp[3], minv, dot;
for (i = 0; i < nlocal; i++) {
// If isolated, just don't shift
if (wsame[i] < wmin) {
for (a = 0; a < dim; a++)
vshift[i][a] = 0.0;
continue;
}
if (ct[i] < cmin) continue;
minv = 0;
for (a = 0; a < dim; a++)
minv += cgradt[i][a] * cgradt[i][a];
if (minv != 0)
minv = 1 / sqrt(minv);
for (a = 0; a < dim; a++)
ntmp[a] = cgradt[i][a] * minv;
dot = 0.0;
for (a = 0; a < dim; a++)
dot += ntmp[a] * vshift[i][a];
// To allowing shifting into the same phase bulk
// if (dot > 0.0) continue;
for (a = 0; a < dim; a++)
vshift[i][a] -= (1.0 - scale) * ntmp[a] * dot;
}
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOVShift::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = wsame[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0; m = 0;
last = first + n; last = first + n;
for (int i = first; i < last; i++) { for (i = first; i < last; i++)
buf[m++] = vshift[i][0]; wsame[i] = buf[m++];
buf[m++] = vshift[i][1]; }
buf[m++] = vshift[i][2];
/* ---------------------------------------------------------------------- */
int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, a, last;
m = 0;
last = first + n;
if (comm_stage == 0) {
for (i = first; i < last; i++) {
buf[m++] = vshift[i][0];
buf[m++] = vshift[i][1];
buf[m++] = vshift[i][2];
}
} else if (comm_stage == 1) {
for (i = first; i < last; i++)
buf[m++] = ct[i];
} else {
for (i = first; i < last; i++) {
for (a = 0; a < 3; a++)
buf[m++] = cgradt[i][a];
buf[m++] = wsame[i];
}
} }
return m; return m;
} }
@ -301,14 +592,28 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
{ {
int i, j, m; int i, j, a, m;
m = 0; m = 0;
for (i = 0; i < n; i++) { if (comm_stage == 0) {
j = list[i]; for (i = 0; i < n; i++) {
vshift[j][0] += buf[m++]; j = list[i];
vshift[j][1] += buf[m++]; vshift[j][0] += buf[m++];
vshift[j][2] += buf[m++]; vshift[j][1] += buf[m++];
vshift[j][2] += buf[m++];
}
} else if (comm_stage == 1) {
for (i = 0; i < n; i++) {
j = list[i];
ct[j] += buf[m++];
}
} else {
for (i = 0; i < n; i++) {
j = list[i];
for (a = 0; a < 3; a++)
cgradt[j][a] += buf[m++];
wsame[j] += buf[m++];
}
} }
} }
@ -319,5 +624,9 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
double ComputeRHEOVShift::memory_usage() double ComputeRHEOVShift::memory_usage()
{ {
double bytes = 3 * nmax_store * sizeof(double); double bytes = 3 * nmax_store * sizeof(double);
if (cross_type_flag)
bytes += 5 * nmax_store * sizeof(double);
return bytes; return bytes;
} }

10
src/RHEO/compute_rheo_vshift.h
View File

@ -31,19 +31,25 @@ class ComputeRHEOVShift : public Compute {
void init() override; void init() override;
void init_list(int, class NeighList *) override; void init_list(int, class NeighList *) override;
void compute_peratom() override; void compute_peratom() override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override; int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override; void unpack_reverse_comm(int, int *, double *) override;
double memory_usage() override; double memory_usage() override;
void correct_surfaces(); void correct_surfaces();
void correct_type_interface();
double **vshift; double **vshift;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
int nmax_store; int nmax_store, comm_stage;
double dtv, cut, cutsq, cutthird; double dtv, cut, cutsq, cutthird;
int surface_flag, interface_flag; double scale, wmin, cmin;
int surface_flag, interface_flag, cross_type_flag;
double *rho0; double *rho0;
double *wsame, *ct, **cgradt;
int *shift_type;
class NeighList *list; class NeighList *list;
class ComputeRHEOInterface *compute_interface; class ComputeRHEOInterface *compute_interface;

70
src/RHEO/fix_rheo.cpp
View File

@ -38,23 +38,25 @@ using namespace LAMMPS_NS;
using namespace RHEO_NS; using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
#if 0
// publication was removed from documentation
static const char cite_rheo[] = static const char cite_rheo[] =
"@article{PalermoInPrep,\n" "@article{Palermo2024,\n"
" journal = {in prep},\n" " journal = {Physics of Fluids},\n"
" title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n" " title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\n"
" volume = {36},\n"
" number = {11},\n"
" pages = {113337},\n"
" year = {2024},\n" " year = {2024},\n"
" author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n" " issn = {1070-6631},\n"
" doi = {https://doi.org/10.1063/5.0228823},\n"
" author = {Palermo, Eric T. and Wolf, Ki T. and Clemmer, Joel T. and O'Connor, Thomas C.},\n"
"}\n\n"; "}\n\n";
#endif
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), shift_type(nullptr),
compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr), compute_grad(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
compute_rhosum(nullptr), compute_vshift(nullptr) compute_surface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr)
{ {
time_integrate = 1; time_integrate = 1;
@ -68,10 +70,12 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
shift_flag = 0; shift_flag = 0;
interface_flag = 0; interface_flag = 0;
surface_flag = 0; surface_flag = 0;
oxidation_flag = 0; coordination_flag = 0;
self_mass_flag = 0;
int i; rhosum_self_mass_flag = 0;
shift_cross_type_flag = 0;
int i, nlo, nhi;
int n = atom->ntypes; int n = atom->ntypes;
memory->create(rho0, n + 1, "rheo:rho0"); memory->create(rho0, n + 1, "rheo:rho0");
memory->create(csq, n + 1, "rheo:csq"); memory->create(csq, n + 1, "rheo:csq");
@ -111,6 +115,26 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg], "shift") == 0) { if (strcmp(arg[iarg], "shift") == 0) {
shift_flag = 1; shift_flag = 1;
memory->create(shift_type, n + 1, "rheo:shift_type");
for (i = 1; i <= n; i++) shift_type[i] = 1;
while (iarg < narg) { // optional sub-arguments
if (strcmp(arg[iarg], "scale/cross/type") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
shift_cross_type_flag = 1;
shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
shift_cmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
shift_wmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
iarg += 3;
} else if (strcmp(arg[iarg], "exclude/type") == 0) {
if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
for (i = nlo; i <= nhi; i++) shift_type[i] = 0;
iarg += 1;
} else {
break;
}
iarg += 1;
}
} else if (strcmp(arg[iarg], "thermal") == 0) { } else if (strcmp(arg[iarg], "thermal") == 0) {
thermal_flag = 1; thermal_flag = 1;
} else if (strcmp(arg[iarg], "surface/detection") == 0) { } else if (strcmp(arg[iarg], "surface/detection") == 0) {
@ -127,14 +151,19 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
} else { } else {
error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]); error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
} }
iarg += 3; iarg += 3;
} else if (strcmp(arg[iarg], "interface/reconstruct") == 0) { } else if (strcmp(arg[iarg], "interface/reconstruct") == 0) {
interface_flag = 1; interface_flag = 1;
} else if (strcmp(arg[iarg], "rho/sum") == 0) { } else if (strcmp(arg[iarg], "rho/sum") == 0) {
rhosum_flag = 1; rhosum_flag = 1;
} else if (strcmp(arg[iarg], "self/mass") == 0) { while (iarg < narg) { // optional sub-arguments
self_mass_flag = 1; if (strcmp(arg[iarg], "self/mass") == 0) {
rhosum_self_mass_flag = 1;
} else {
break;
}
iarg += 1;
}
} else if (strcmp(arg[iarg], "density") == 0) { } else if (strcmp(arg[iarg], "density") == 0) {
if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error); if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error);
for (i = 1; i <= n; i++) rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp); for (i = 1; i <= n; i++) rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
@ -152,12 +181,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
iarg += 1; iarg += 1;
} }
if (self_mass_flag && (!rhosum_flag))
error->all(FLERR, "Cannot use self/mass setting without rho/sum");
#if 0
if (lmp->citeme) lmp->citeme->add(cite_rheo); if (lmp->citeme) lmp->citeme->add(cite_rheo);
#endif
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -173,6 +197,7 @@ FixRHEO::~FixRHEO()
memory->destroy(csq); memory->destroy(csq);
memory->destroy(rho0); memory->destroy(rho0);
memory->destroy(shift_type);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -192,7 +217,7 @@ void FixRHEO::post_constructor()
if (rhosum_flag) { if (rhosum_flag) {
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>( compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(
modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag))); modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", rhosum_self_mass_flag)));
compute_rhosum->fix_rheo = this; compute_rhosum->fix_rheo = this;
} }
@ -375,7 +400,6 @@ void FixRHEO::initial_integrate(int /*vflag*/)
// Shifting atoms // Shifting atoms
if (shift_flag) { if (shift_flag) {
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (status[i] & STATUS_NO_SHIFT) continue; if (status[i] & STATUS_NO_SHIFT) continue;
if (status[i] & PHASECHECK) continue; if (status[i] & PHASECHECK) continue;

14
src/RHEO/fix_rheo.h
View File

@ -41,19 +41,25 @@ class FixRHEO : public Fix {
// Model parameters // Model parameters
double cut; double cut;
double *rho0, *csq; double *rho0, *csq;
int self_mass_flag;
int zmin_kernel, zmin_surface, zmin_splash; int zmin_kernel, zmin_surface, zmin_splash;
int kernel_style, surface_style; int kernel_style, surface_style;
double divr_surface; double divr_surface;
// Accessory fixes/computes // Settings flags
int thermal_flag; int thermal_flag;
int rhosum_flag; int rhosum_flag;
int shift_flag;
int interface_flag; int interface_flag;
int surface_flag; int surface_flag;
int oxidation_flag; int shift_flag;
int coordination_flag;
// Optional sub-flags/parameters
int rhosum_self_mass_flag;
int *shift_type;
int shift_cross_type_flag;
double shift_scale, shift_wmin, shift_cmin;
// Accessory fixes/computes
int viscosity_fix_defined; int viscosity_fix_defined;
int pressure_fix_defined; int pressure_fix_defined;
int thermal_fix_defined; int thermal_fix_defined;

80
src/RHEO/fix_rheo_pressure.cpp
View File

@ -31,7 +31,7 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL }; enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL , IDEAL };
static constexpr double SEVENTH = 1.0 / 7.0; static constexpr double SEVENTH = 1.0 / 7.0;
@ -39,12 +39,15 @@ static constexpr double SEVENTH = 1.0 / 7.0;
FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr), Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr),
rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr), rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), gamma(nullptr),
fix_rheo(nullptr) pressure_style(nullptr), fix_rheo(nullptr)
{ {
if (narg < 4) error->all(FLERR, "Illegal fix command"); if (narg < 4) error->all(FLERR, "Illegal fix command");
comm_forward = 1; comm_forward = 1;
variable_csq = 0;
invertible_pressure = 1;
background_flag = 0;
// Currently can only have one instance of fix rheo/pressure // Currently can only have one instance of fix rheo/pressure
if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all"); if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all");
@ -55,7 +58,12 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
memory->create(c_cubic, n + 1, "rheo:c_cubic"); memory->create(c_cubic, n + 1, "rheo:c_cubic");
memory->create(tpower, n + 1, "rheo:tpower"); memory->create(tpower, n + 1, "rheo:tpower");
memory->create(pbackground, n + 1, "rheo:pbackground"); memory->create(pbackground, n + 1, "rheo:pbackground");
for (i = 1; i <= n; i++) pressure_style[i] = NONE; memory->create(gamma, n + 1, "rheo:gamma");
for (i = 1; i <= n; i++) {
pressure_style[i] = NONE;
pbackground[i] = 0.0;
}
int iarg = 3; int iarg = 3;
while (iarg < narg) { while (iarg < narg) {
@ -68,16 +76,14 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg + 1], "tait/water") == 0) { } else if (strcmp(arg[iarg + 1], "tait/water") == 0) {
for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER;
} else if (strcmp(arg[iarg + 1], "tait/general") == 0) { } else if (strcmp(arg[iarg + 1], "tait/general") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait", error); if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait/general", error);
double tpower_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); double tpower_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
double pbackground_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg += 1;
iarg += 2;
for (i = nlo; i <= nhi; i++) { for (i = nlo; i <= nhi; i++) {
pressure_style[i] = TAITGENERAL; pressure_style[i] = TAITGENERAL;
tpower[i] = tpower_one; tpower[i] = tpower_one;
pbackground[i] = pbackground_one;
} }
} else if (strcmp(arg[iarg + 1], "cubic") == 0) { } else if (strcmp(arg[iarg + 1], "cubic") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure cubic", error); if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure cubic", error);
@ -89,6 +95,31 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
pressure_style[i] = CUBIC; pressure_style[i] = CUBIC;
c_cubic[i] = c_cubic_one; c_cubic[i] = c_cubic_one;
} }
invertible_pressure = 0;
} else if (strcmp(arg[iarg + 1], "ideal/gas") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure ideal/gas", error);
double c_gamma_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
iarg += 1;
for (i = nlo; i <= nhi; i++) {
pressure_style[i] = IDEAL;
gamma[i] = c_gamma_one;
}
variable_csq = 1;
if (atom->esph_flag != 1)
error->all(FLERR, "fix rheo/pressure ideal gas equation of state requires atom property esph");
} else if (strcmp(arg[iarg + 1], "background") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure background", error);
double pbackground_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
iarg += 1;
for (i = nlo; i <= nhi; i++)
pbackground[i] = pbackground_one;
background_flag = 1;
} else { } else {
error->all(FLERR, "Illegal fix command, {}", arg[iarg]); error->all(FLERR, "Illegal fix command, {}", arg[iarg]);
} }
@ -110,6 +141,7 @@ FixRHEOPressure::~FixRHEOPressure()
memory->destroy(c_cubic); memory->destroy(c_cubic);
memory->destroy(tpower); memory->destroy(tpower);
memory->destroy(pbackground); memory->destroy(pbackground);
memory->destroy(gamma);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -197,10 +229,11 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixRHEOPressure::calc_pressure(double rho, int type) double FixRHEOPressure::calc_pressure(double rho, int i)
{ {
double p = 0.0; double p = 0.0;
double dr, rr3, rho_ratio; double dr, rr3, rho_ratio;
int type = atom->type[i];
if (pressure_style[type] == LINEAR) { if (pressure_style[type] == LINEAR) {
p = csq[type] * (rho - rho0[type]); p = csq[type] * (rho - rho0[type]);
@ -214,16 +247,25 @@ double FixRHEOPressure::calc_pressure(double rho, int type)
} else if (pressure_style[type] == TAITGENERAL) { } else if (pressure_style[type] == TAITGENERAL) {
rho_ratio = rho * rho0inv[type]; rho_ratio = rho * rho0inv[type];
p = csq[type] * rho0[type] * (pow(rho_ratio, tpower[type]) - 1.0) / tpower[type]; p = csq[type] * rho0[type] * (pow(rho_ratio, tpower[type]) - 1.0) / tpower[type];
p += pbackground[type]; } else if (pressure_style[type] == IDEAL) {
p = (gamma[type] - 1.0) * rho * atom->esph[i] / atom->mass[type];
} }
if (background_flag)
p += pbackground[type];
return p; return p;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixRHEOPressure::calc_rho(double p, int type) double FixRHEOPressure::calc_rho(double p, int i)
{ {
double rho = 0.0; double rho = 0.0;
int type = atom->type[i];
if (background_flag)
p -= pbackground[type];
if (pressure_style[type] == LINEAR) { if (pressure_style[type] == LINEAR) {
rho = csqinv[type] * p + rho0[type]; rho = csqinv[type] * p + rho0[type];
@ -235,10 +277,24 @@ double FixRHEOPressure::calc_rho(double p, int type)
rho *= pow(rho0[type], 6.0 * SEVENTH); rho *= pow(rho0[type], 6.0 * SEVENTH);
rho *= pow(csq[type], -SEVENTH); rho *= pow(csq[type], -SEVENTH);
} else if (pressure_style[type] == TAITGENERAL) { } else if (pressure_style[type] == TAITGENERAL) {
p -= pbackground[type];
rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]); rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]);
rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]); rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]);
rho *= pow(csq[type], -1.0 / tpower[type]); rho *= pow(csq[type], -1.0 / tpower[type]);
} else if (pressure_style[type] == IDEAL) {
rho = p * atom->mass[type] / ((gamma[type] - 1.0) * atom->esph[i]);
} }
return rho; return rho;
} }
/* ---------------------------------------------------------------------- */
double FixRHEOPressure::calc_csq(double p, int i)
{
int type = atom->type[i];
double csq2 = csq[type];
if (pressure_style[type] == IDEAL) {
csq2 = (gamma[type] - 1.0) * atom->esph[i] / atom->mass[type];
}
return csq2;
}

7
src/RHEO/fix_rheo_pressure.h
View File

@ -36,10 +36,13 @@ class FixRHEOPressure : public Fix {
void unpack_forward_comm(int, int, double *) override; void unpack_forward_comm(int, int, double *) override;
double calc_pressure(double, int); double calc_pressure(double, int);
double calc_rho(double, int); double calc_rho(double, int);
double calc_csq(double, int);
int variable_csq;
int invertible_pressure;
private: private:
double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground; double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground, *gamma;
int *pressure_style; int *pressure_style, background_flag;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
}; };

26
src/RHEO/pair_rheo.cpp
View File

@ -80,8 +80,8 @@ void PairRHEO::compute(int eflag, int vflag)
int pair_force_flag, pair_rho_flag, pair_avisc_flag; int pair_force_flag, pair_rho_flag, pair_avisc_flag;
int fluidi, fluidj; int fluidi, fluidj;
double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave; double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave;
double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj;
kappa_ave, dT_prefactor; double eta_ave, kappa_ave, dT_prefactor;
double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
double *dWij, *dWji; double *dWij, *dWji;
double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3]; double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3];
@ -185,7 +185,13 @@ void PairRHEO::compute(int eflag, int vflag)
kappaj = conductivity[j]; kappaj = conductivity[j];
} }
cs_ave = 0.5 * (cs[itype] + cs[jtype]); if (!variable_csq) {
cs_ave = 0.5 * (cs[itype] + cs[jtype]);
} else {
csqi = fix_pressure->calc_csq(rhoi, i);
csqj = fix_pressure->calc_csq(rhoj, j);
cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
}
csq_ave = cs_ave * cs_ave; csq_ave = cs_ave * cs_ave;
pair_rho_flag = 0; pair_rho_flag = 0;
@ -221,7 +227,7 @@ void PairRHEO::compute(int eflag, int vflag)
if (fluidi && (!fluidj)) { if (fluidi && (!fluidj)) {
compute_interface->correct_v(vj, vi, j, i); compute_interface->correct_v(vj, vi, j, i);
rhoj = compute_interface->correct_rho(j); rhoj = compute_interface->correct_rho(j);
Pj = fix_pressure->calc_pressure(rhoj, jtype); Pj = fix_pressure->calc_pressure(rhoj, j);
if ((chi[j] > 0.9) && (r < (cutk * 0.5))) if ((chi[j] > 0.9) && (r < (cutk * 0.5)))
fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv; fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv;
@ -229,7 +235,7 @@ void PairRHEO::compute(int eflag, int vflag)
} else if ((!fluidi) && fluidj) { } else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vi, vj, i, j); compute_interface->correct_v(vi, vj, i, j);
rhoi = compute_interface->correct_rho(i); rhoi = compute_interface->correct_rho(i);
Pi = fix_pressure->calc_pressure(rhoi, itype); Pi = fix_pressure->calc_pressure(rhoi, i);
if (chi[i] > 0.9 && r < (cutk * 0.5)) if (chi[i] > 0.9 && r < (cutk * 0.5))
fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv; fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv;
@ -238,6 +244,14 @@ void PairRHEO::compute(int eflag, int vflag)
rhoi = rho0i; rhoi = rho0i;
rhoj = rho0j; rhoj = rho0j;
} }
// recalculate speed of sound, if necessary
if (variable_csq && ((!fluidi) || (!fluidj))) {
csqi = fix_pressure->calc_csq(rhoi, i);
csqj = fix_pressure->calc_csq(rhoj, j);
cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj));
csq_ave = cs_ave * cs_ave;
}
} }
// Repel if close to inner solid particle // Repel if close to inner solid particle
@ -480,6 +494,8 @@ void PairRHEO::setup()
csq = fix_rheo->csq; csq = fix_rheo->csq;
rho0 = fix_rheo->rho0; rho0 = fix_rheo->rho0;
variable_csq = fix_pressure->variable_csq;
if (cutk != fix_rheo->cut) if (cutk != fix_rheo->cut)
error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk,
fix_rheo->cut); fix_rheo->cut);

1
src/RHEO/pair_rheo.h
View File

@ -45,6 +45,7 @@ class PairRHEO : public Pair {
int rho_damp_flag; int rho_damp_flag;
int thermal_flag; int thermal_flag;
int interface_flag; int interface_flag;
int variable_csq;
int harmonic_means_flag; int harmonic_means_flag;

8
src/comm.cpp
View File

@ -228,13 +228,13 @@ void Comm::init()
} }
for (const auto &compute : modify->get_compute_list()) { for (const auto &compute : modify->get_compute_list()) {
maxforward = MAX(maxforward,compute->comm_forward); maxforward = MAX(maxforward, compute->comm_forward);
maxreverse = MAX(maxreverse,compute->comm_reverse); maxreverse = MAX(maxreverse, compute->comm_reverse);
} }
for (const auto &dump: output->get_dump_list()) { for (const auto &dump: output->get_dump_list()) {
maxforward = MAX(maxforward,dump->comm_forward); maxforward = MAX(maxforward, dump->comm_forward);
maxreverse = MAX(maxreverse,dump->comm_reverse); maxreverse = MAX(maxreverse, dump->comm_reverse);
} }
if (force->newton == 0) maxreverse = 0; if (force->newton == 0) maxreverse = 0;

16
src/comm.h
View File

@ -87,17 +87,17 @@ class Comm : protected Pointers {
// forward/reverse comm from a Pair, Bond, Fix, Compute, Dump // forward/reverse comm from a Pair, Bond, Fix, Compute, Dump
virtual void forward_comm(class Pair *) = 0; virtual void forward_comm(class Pair *, int size = 0) = 0;
virtual void reverse_comm(class Pair *) = 0; virtual void reverse_comm(class Pair *, int size = 0) = 0;
virtual void forward_comm(class Bond *) = 0; virtual void forward_comm(class Bond *, int size = 0) = 0;
virtual void reverse_comm(class Bond *) = 0; virtual void reverse_comm(class Bond *, int size = 0) = 0;
virtual void forward_comm(class Fix *, int size = 0) = 0; virtual void forward_comm(class Fix *, int size = 0) = 0;
virtual void reverse_comm(class Fix *, int size = 0) = 0; virtual void reverse_comm(class Fix *, int size = 0) = 0;
virtual void reverse_comm_variable(class Fix *) = 0; virtual void reverse_comm_variable(class Fix *) = 0;
virtual void forward_comm(class Compute *) = 0; virtual void forward_comm(class Compute *, int size = 0) = 0;
virtual void reverse_comm(class Compute *) = 0; virtual void reverse_comm(class Compute *, int size = 0) = 0;
virtual void forward_comm(class Dump *) = 0; virtual void forward_comm(class Dump *, int size = 0) = 0;
virtual void reverse_comm(class Dump *) = 0; virtual void reverse_comm(class Dump *, int size = 0) = 0;
// forward comm of an array // forward comm of an array

108
src/comm_brick.cpp
View File

@ -984,16 +984,21 @@ void CommBrick::borders()
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Pair forward communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Pair
size > 0 -> Pair passes max size per atom
the latter is only useful if Pair does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::forward_comm(Pair *pair) void CommBrick::forward_comm(Pair *pair, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = pair->comm_forward; if (size) nsize = size;
else nsize = pair->comm_forward;
for (iswap = 0; iswap < nswap; iswap++) { for (iswap = 0; iswap < nswap; iswap++) {
@ -1021,16 +1026,21 @@ void CommBrick::forward_comm(Pair *pair)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Pair reverse communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Pair
size > 0 -> Pair passes max size per atom
the latter is only useful if Pair does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::reverse_comm(Pair *pair) void CommBrick::reverse_comm(Pair *pair, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); if (size) nsize = size;
else nsize = MAX(pair->comm_reverse,pair->comm_reverse_off);
for (iswap = nswap-1; iswap >= 0; iswap--) { for (iswap = nswap-1; iswap >= 0; iswap--) {
@ -1058,16 +1068,21 @@ void CommBrick::reverse_comm(Pair *pair)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Bond forward communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::forward_comm(Bond *bond) void CommBrick::forward_comm(Bond *bond, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = bond->comm_forward; if (size) nsize = size;
else nsize = bond->comm_forward;
for (iswap = 0; iswap < nswap; iswap++) { for (iswap = 0; iswap < nswap; iswap++) {
@ -1095,16 +1110,21 @@ void CommBrick::forward_comm(Bond *bond)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Bond reverse communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::reverse_comm(Bond *bond) void CommBrick::reverse_comm(Bond *bond, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); if (size) nsize = size;
else nsize = MAX(bond->comm_reverse,bond->comm_reverse_off);
for (iswap = nswap-1; iswap >= 0; iswap--) { for (iswap = nswap-1; iswap >= 0; iswap--) {
@ -1175,10 +1195,10 @@ void CommBrick::forward_comm(Fix *fix, int size)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Fix reverse communication invoked by a Fix
size/nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix size = 0 (default) -> use comm_reverse from Fix
size > 0 -> Fix passes max size per atom size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes, the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::reverse_comm(Fix *fix, int size) void CommBrick::reverse_comm(Fix *fix, int size)
@ -1259,16 +1279,21 @@ void CommBrick::reverse_comm_variable(Fix *fix)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Compute forward communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::forward_comm(Compute *compute) void CommBrick::forward_comm(Compute *compute, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = compute->comm_forward; if (size) nsize = size;
else nsize = compute->comm_forward;
for (iswap = 0; iswap < nswap; iswap++) { for (iswap = 0; iswap < nswap; iswap++) {
@ -1297,16 +1322,21 @@ void CommBrick::forward_comm(Compute *compute)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Compute reverse communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::reverse_comm(Compute *compute) void CommBrick::reverse_comm(Compute *compute, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = compute->comm_reverse; if (size) nsize = size;
else nsize = compute->comm_reverse;
for (iswap = nswap-1; iswap >= 0; iswap--) { for (iswap = nswap-1; iswap >= 0; iswap--) {
@ -1334,16 +1364,21 @@ void CommBrick::reverse_comm(Compute *compute)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Dump forward communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::forward_comm(Dump *dump) void CommBrick::forward_comm(Dump *dump, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = dump->comm_forward; if (size) nsize = size;
else nsize = dump->comm_forward;
for (iswap = 0; iswap < nswap; iswap++) { for (iswap = 0; iswap < nswap; iswap++) {
@ -1372,16 +1407,21 @@ void CommBrick::forward_comm(Dump *dump)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Dump reverse communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommBrick::reverse_comm(Dump *dump) void CommBrick::reverse_comm(Dump *dump, int size)
{ {
int iswap,n; int iswap,n,nsize;
double *buf; double *buf;
MPI_Request request; MPI_Request request;
int nsize = dump->comm_reverse; if (size) nsize = size;
else nsize = dump->comm_reverse;
for (iswap = nswap-1; iswap >= 0; iswap--) { for (iswap = nswap-1; iswap >= 0; iswap--) {

22
src/comm_brick.h
View File

@ -32,17 +32,17 @@ class CommBrick : public Comm {
void exchange() override; // move atoms to new procs void exchange() override; // move atoms to new procs
void borders() override; // setup list of atoms to comm void borders() override; // setup list of atoms to comm
void forward_comm(class Pair *) override; // forward comm from a Pair void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair
void reverse_comm(class Pair *) override; // reverse comm from a Pair void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair
void forward_comm(class Bond *) override; // forward comm from a Bond void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond
void reverse_comm(class Bond *) override; // reverse comm from a Bond void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond
void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix
void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix
void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix
void forward_comm(class Compute *) override; // forward from a Compute void forward_comm(class Compute *, int size = 0) override; // forward from a Compute
void reverse_comm(class Compute *) override; // reverse from a Compute void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute
void forward_comm(class Dump *) override; // forward comm from a Dump void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump
void reverse_comm(class Dump *) override; // reverse comm from a Dump void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump
void forward_comm_array(int, double **) override; // forward comm of array void forward_comm_array(int, double **) override; // forward comm of array
void *extract(const char *, int &) override; void *extract(const char *, int &) override;

109
src/comm_tiled.cpp
View File

@ -1382,14 +1382,19 @@ void CommTiled::borders()
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Pair forward communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Pair
size > 0 -> Pair passes max size per atom
the latter is only useful if Pair does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::forward_comm(Pair *pair) void CommTiled::forward_comm(Pair *pair, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = pair->comm_forward; if (size) nsize = size;
else nsize = pair->comm_forward;
for (int iswap = 0; iswap < nswap; iswap++) { for (int iswap = 0; iswap < nswap; iswap++) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1426,14 +1431,19 @@ void CommTiled::forward_comm(Pair *pair)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Pair reverse communication invoked by a Pair
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Fix
size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::reverse_comm(Pair *pair) void CommTiled::reverse_comm(Pair *pair, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); if (size) nsize = MAX(pair->comm_reverse, pair->comm_reverse_off);
else nsize = pair->comm_reverse;
for (int iswap = nswap-1; iswap >= 0; iswap--) { for (int iswap = nswap-1; iswap >= 0; iswap--) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1466,14 +1476,18 @@ void CommTiled::reverse_comm(Pair *pair)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Bond forward communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::forward_comm(Bond *bond, int size)
void CommTiled::forward_comm(Bond *bond)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = bond->comm_forward; if (size) nsize = size;
else nsize = bond->comm_forward;
for (int iswap = 0; iswap < nswap; iswap++) { for (int iswap = 0; iswap < nswap; iswap++) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1510,14 +1524,19 @@ void CommTiled::forward_comm(Bond *bond)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Bond reverse communication invoked by a Bond
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Bond
size > 0 -> Bond passes max size per atom
the latter is only useful if Bond does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::reverse_comm(Bond *bond) void CommTiled::reverse_comm(Bond *bond, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); if (size) nsize = size;
else nsize = MAX(bond->comm_reverse,bond->comm_reverse_off);
for (int iswap = nswap-1; iswap >= 0; iswap--) { for (int iswap = nswap-1; iswap >= 0; iswap--) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1598,10 +1617,10 @@ void CommTiled::forward_comm(Fix *fix, int size)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Fix reverse communication invoked by a Fix
size/nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Fix size = 0 (default) -> use comm_reverse from Fix
size > 0 -> Fix passes max size per atom size > 0 -> Fix passes max size per atom
the latter is only useful if Fix does several comm modes, the latter is only useful if Fix does several comm modes,
some are smaller than max stored in its comm_forward some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::reverse_comm(Fix *fix, int size) void CommTiled::reverse_comm(Fix *fix, int size)
@ -1654,14 +1673,19 @@ void CommTiled::reverse_comm_variable(Fix * /*fix*/)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Compute forward communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::forward_comm(Compute *compute) void CommTiled::forward_comm(Compute *compute, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = compute->comm_forward; if (size) nsize = size;
else nsize = compute->comm_forward;
for (int iswap = 0; iswap < nswap; iswap++) { for (int iswap = 0; iswap < nswap; iswap++) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1697,14 +1721,19 @@ void CommTiled::forward_comm(Compute *compute)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Compute reverse communication invoked by a Compute
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Compute
size > 0 -> Compute passes max size per atom
the latter is only useful if Compute does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::reverse_comm(Compute *compute) void CommTiled::reverse_comm(Compute *compute, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = compute->comm_reverse; if (size) nsize = size;
else nsize = compute->comm_reverse;
for (int iswap = nswap-1; iswap >= 0; iswap--) { for (int iswap = nswap-1; iswap >= 0; iswap--) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1738,14 +1767,19 @@ void CommTiled::reverse_comm(Compute *compute)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
forward communication invoked by a Dump forward communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_forward from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_forward
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::forward_comm(Dump *dump) void CommTiled::forward_comm(Dump *dump, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = dump->comm_forward; if (size) nsize = size;
else nsize = dump->comm_forward;
for (int iswap = 0; iswap < nswap; iswap++) { for (int iswap = 0; iswap < nswap; iswap++) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];
@ -1781,14 +1815,19 @@ void CommTiled::forward_comm(Dump *dump)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
reverse communication invoked by a Dump reverse communication invoked by a Dump
nsize used only to set recv buffer limit size/nsize used only to set recv buffer limit
size = 0 (default) -> use comm_reverse from Dump
size > 0 -> Dump passes max size per atom
the latter is only useful if Dump does several comm modes,
some are smaller than max stored in its comm_reverse
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void CommTiled::reverse_comm(Dump *dump) void CommTiled::reverse_comm(Dump *dump, int size)
{ {
int i,irecv,n,nsend,nrecv; int i,irecv,n,nsize,nsend,nrecv;
int nsize = dump->comm_reverse; if (size) nsize = size;
else nsize = dump->comm_reverse;
for (int iswap = nswap-1; iswap >= 0; iswap--) { for (int iswap = nswap-1; iswap >= 0; iswap--) {
nsend = nsendproc[iswap] - sendself[iswap]; nsend = nsendproc[iswap] - sendself[iswap];

22
src/comm_tiled.h
View File

@ -32,17 +32,17 @@ class CommTiled : public Comm {
void exchange() override; // move atoms to new procs void exchange() override; // move atoms to new procs
void borders() override; // setup list of atoms to comm void borders() override; // setup list of atoms to comm
void forward_comm(class Pair *) override; // forward comm from a Pair void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair
void reverse_comm(class Pair *) override; // reverse comm from a Pair void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair
void forward_comm(class Bond *) override; // forward comm from a Bond void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond
void reverse_comm(class Bond *) override; // reverse comm from a Bond void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond
void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix
void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix
void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix
void forward_comm(class Compute *) override; // forward from a Compute void forward_comm(class Compute *, int size = 0) override; // forward from a Compute
void reverse_comm(class Compute *) override; // reverse from a Compute void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute
void forward_comm(class Dump *) override; // forward comm from a Dump void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump
void reverse_comm(class Dump *) override; // reverse comm from a Dump void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump
void forward_comm_array(int, double **) override; // forward comm of array void forward_comm_array(int, double **) override; // forward comm of array

4
src/domain.cpp
View File

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