diff --git a/doc/src/Developer_comm_ops.rst b/doc/src/Developer_comm_ops.rst index 116b4d6ad8..66947de74f 100644 --- a/doc/src/Developer_comm_ops.rst +++ b/doc/src/Developer_comm_ops.rst @@ -79,19 +79,19 @@ containing ``double`` values. To correctly store integers that may be 64-bit (bigint, tagint, imageint) in the buffer, you need to use the :ref:`ubuf union ` construct. -The *Fix*, *Compute*, and *Dump* classes can also invoke the same kind -of forward and reverse communication operations using the same *Comm* -class methods. Likewise, the same pack/unpack methods and +The *Fix*, *Bond*, *Compute*, and *Dump* classes can also invoke the +same kind of forward and reverse communication operations using the +same *Comm* class methods. Likewise, the same pack/unpack methods and comm_forward/comm_reverse variables must be defined by the calling -*Fix*, *Compute*, or *Dump* class. +*Fix*, *Bond*, *Compute*, or *Dump* class. -For *Fix* classes, there is an optional second argument to the +For all of these classes, there is an optional second argument to the *forward_comm()* and *reverse_comm()* call which can be used when the -fix performs multiple modes of communication, with different numbers -of values per atom. The fix should set the *comm_forward* and +class performs multiple modes of communication, with different numbers +of values per atom. The class should set the *comm_forward* and *comm_reverse* variables to the maximum value, but can invoke the communication for a particular mode with a smaller value. For this -to work, the *pack_forward_comm()*, etc methods typically use a class +to work, the *pack_forward_comm()*, etc. methods typically use a class member variable to choose which values to pack/unpack into/from the buffer. diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index 90f9ff82d0..2407fd484d 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -15,8 +15,9 @@ details of the system, or develop new capabilities. For instance, the numerics associated with calculating gradients, reproducing kernels, etc. are separated into distinct classes to simplify the development of new integration schemes which can call these calculations. Additional numerical details can be found in -:ref:`(Clemmer) `. Example movies illustrating some of these -capabilities are found at https://www.lammps.org/movies.html#rheopackage. +:ref:`(Palermo) ` and :ref:`(Clemmer) `. +Example movies illustrating some of these capabilities are found at +https://www.lammps.org/movies.html#rheopackage. Note, if you simply want to run a traditional SPH simulation, the :ref:`SPH package ` 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 ` command), (b) create bpm bonds between the particles (see the :doc:`bpm howto ` page for more details), and (c) use :doc:`pair rheo/solid ` to apply repulsive contact forces between distinct solid bodies. Akin to pair rheo, -pair rheo/solid considers a particles fluid/solid phase to determine whether to +pair rheo/solid considers a particle's fluid/solid phase to determine whether to apply forces. However, unlike pair rheo, pair rheo/solid does obey special bond settings such that contact forces do not have to be calculated between two bonded solid particles in the same elastic body. @@ -79,10 +80,10 @@ In systems with thermal evolution, fix rheo/thermal can optionally set a melting/solidification temperature allowing particles to dynamically swap their state between fluid and solid when the temperature exceeds or drops below the critical temperature, respectively. Using the *react* option, one can specify a maximum -bond length and a bond type. Then, when solidifying, particles will search their +bond length and a bond type. Then, when solidifying, particles search their local neighbors and automatically create bonds with any neighboring solid particles -in range. For BPM bond styles, bonds will then use the immediate position of the two -particles to calculate a reference state. When melting, particles will delete any +in range. For BPM bond styles, bonds then use the immediate position of the two +particles to calculate a reference state. When melting, particles delete any bonds of the specified type when reverting to a fluid state. Special bonds are updated as bonds are created/broken. @@ -107,6 +108,10 @@ criteria for creating/deleting a bond or altering force calculations). ---------- +.. _howto_rheo_palermo: + +**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024). + .. _howto_rheo_clemmer: **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index eb88ef0536..0214289dee 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -16,21 +16,36 @@ Syntax * kstyle = *quintic* or *RK0* or *RK1* or *RK2* * zmin = minimal number of neighbors for reproducing kernels * zero or more keyword/value pairs may be appended to args -* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *self/mass* or *speed/sound* +* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *speed/sound* .. parsed-literal:: - *thermal* values = none, turns on thermal evolution - *interface/reconstruct* values = none, reconstructs interfaces with solid particles - *surface/detection* values = *sdstyle* *limit* *limit/splash* - *sdstyle* = *coordination* or *divergence* - *limit* = threshold for surface particles - *limit/splash* = threshold for splash particles - *shift* values = none, turns on velocity shifting - *rho/sum* values = none, uses the kernel to compute the density of particles - *self/mass* values = none, a particle uses its own mass in a rho summation - *density* values = *rho01*, ... *rho0N* (density) - *speed/sound* values = *cs0*, ... *csN* (velocity) + *thermal* turns on thermal evolution + values = none + *interface/reconstruct* reconstructs interfaces with solid particles + values = none + *surface/detection* detects free-surfaces with an absence of particles + values = *sdstyle* *limit* *limit/splash* + *sdstyle* = *coordination* or *divergence* + *limit* = threshold for surface particles + *limit/splash* = threshold for splash particles (unitless) + *shift* turns on velocity shifting + values = none + optional args = *exclude/type* or *scale/cross/type* + *exclude/type* values = *types* + *types* = list of types + *scale/cross/type* values = *shiftscale* *cmin* *wmin* + *shiftscale* = fraction of shifting in normal direction to preserve (unitless) + *cmin* = minimum color function value required for scaling (unitless) + *wmin* = minimum local same-type support required for any shifting (unitless) + *rho/sum* density evolution performed by a kernel summation + values = none + optional args = *self/mass* + *self/mass* values = none, a particle uses its own mass in summation + *density* specify equilibrium densities for each atom type + values = *rho01*, ... *rho0N* (density) + *speed/sound* specify speeds of sound for each atom type + values = *cs0*, ... *csN* (velocity) Examples """""""" @@ -39,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 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 """"""""""" @@ -46,8 +63,10 @@ Description .. versionadded:: 29Aug2024 Perform time integration for RHEO particles, updating positions, velocities, -and densities. For an overview of other features available in the RHEO package, -see :doc:`the RHEO howto `. +and densities. For a detailed breakdown of the integration timestep and +numerical details, see :ref:`(Palermo) `. For an overview +and list of other features available in the RHEO package, see +:doc:`the RHEO howto `. The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four kernels are currently available. The *quintic* kernel is a standard quintic @@ -70,16 +89,51 @@ and velocity of solid particles are alternatively reconstructed for every fluid-solid interaction to ensure no-slip and pressure-balanced boundaries. This is done by estimating the location of the fluid-solid interface and extrapolating fluid particle properties across the interface to calculate a -temporary apparent density and velocity for a solid particle. +temporary apparent density and velocity for a solid particle. The numerical +details are the same as those described in +:ref:`(Palermo) ` except there is an additional +restriction that the reconstructed solid density cannot be less than the +equilibrium density. This prevents fluid particles from sticking to solid +surfaces. A modified form of Fickian particle shifting can be enabled with the *shift* keyword. This effectively shifts particle positions to generate a -more uniform spatial distribution. Shifting currently does not consider the +more uniform spatial distribution. By default, shifting does not consider the type of a particle and therefore may be inappropriate in systems consisting -of multiple fluid phases. +of multiple atom types representing multiple fluid phases. However, two +optional subarguments can follow the *shift* keyword, *exclude/type* and +*scale/cross/type* to adjust shifting at fluid interfaces. -In systems with free surfaces, the *surface/detection* keyword can be used -to classify the location of particles as being within the bulk fluid, on a +The *exclude/type* option lets the user specify a list of atom types which +are not shifted, *types*. A wild-card asterisk can be used in place +of or in conjunction with the *types* argument to toggle shifting for +multiple atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with +no numeric values means all types from 1 to :math:`N`. A leading asterisk +means all types from 1 to n (inclusive). A trailing asterisk means all types +from m to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). + +The *scale/cross/type* option is designed to handle interfaces between fluids +made up of different atom types. Similar to the method by +:ref:`(Yang) `, 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. Shifting is then disabled in the normal direction away from the free surface to prevent particles from diffusing away. Surface detection can also be used @@ -101,10 +155,9 @@ threshold for this classification is set by the numerical value of By default, RHEO integrates particles' densities using a mass diffusion equation. Alternatively, one can update densities every timestep by performing a kernel summation of the masses of neighboring particles by specifying the *rho/sum* -keyword. - -The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*. -Typically, the density :math:`\rho` of a particle is calculated as the sum over neighbors +keyword. Following this keyword, one may include the optional *self/mass* subargument +which modifies the behavior of the density summation. Typically, the density +:math:`\rho` of a particle is calculated as the sum over neighbors .. math:: \rho_i = \sum_{j} W_{ij} M_j @@ -120,7 +173,9 @@ equilibrium density *rho0*. The *speed/sound* keyword is used to specify the speed of sound of each of the N particle types. It must be followed by N numerical values specifying each type's -speed of sound *cs*. +speed of sound *cs*. These values may be ignored if the pressure equation of +state has a non-constant speed of sound, as discussed further in +:doc:`fix rheo/pressure `. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -163,6 +218,14 @@ Default ---------- +.. _rheo_palermo: + +**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024). + +.. _rheo_yang: + +**(Yang)** Yang, Rakhsha, Hu, Negrut, J. Comp. Physics, 458, 111079 (2022). + .. _fix_rheo_hu: -**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006). +**(Hu)** Hu, and Adams, J. Comp. Physics, 213, 844-861 (2006). diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index 2a714b298b..7d5d08ddd1 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -14,13 +14,16 @@ Syntax * rheo/pressure = style name of this fix command * one or more types and pressure styles must be appended * types = lists of types (see below) -* pstyle = *linear* or *taitwater* or *cubic* +* pstyle = *linear* or *tait/water* or *tait/general* or *cubic* or *ideal/gas* or *background* .. parsed-literal:: *linear* args = none - *taitwater* args = none + *tait/water* args = none + *tait/general* args = exponent :math:`gamma` (unitless) *cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2) + *ideal/gas* args = heat capacity ratio :math:`gamma` (unitless) + *background* args = background pressure :math:`P[b]` (pressure) Examples """""""" @@ -29,6 +32,7 @@ Examples fix 1 all rheo/pressure * linear fix 1 all rheo/pressure 1 linear 2 cubic 10.0 + fix 1 all rheo/pressure * linear * background 0.1 Description """"""""""" @@ -40,13 +44,12 @@ define different equations of state for different atom types. An equation must be specified for every atom type. One first defines the atom *types*. A wild-card asterisk can be used in place -of or in conjunction with the *types* argument to set the coefficients for -multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" -or "m\*n". If :math:`N` is the number of atom types, then an asterisk with -no numeric values means all types from 1 to :math:`N`. A leading asterisk -means all types from 1 to n (inclusive). A trailing asterisk means all types -from m to :math:`N` (inclusive). A middle asterisk means all types from m to n -(inclusive). +of or in conjunction with the *types* argument to set values for multiple atom +types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is +the number of atom types, then an asterisk with no numeric values means all types +from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive). +A trailing asterisk means all types from m to :math:`N` (inclusive). A middle +asterisk means all types from m to n (inclusive). The *types* definition is followed by the pressure style, *pstyle*. Current options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear @@ -54,7 +57,7 @@ equation of state with a particle pressure :math:`P` calculated as .. math:: - P = c (\rho - \rho_0) + P = c^2 (\rho - \rho_0) where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density, and :math:`\rho` is the current density of a particle. The numerical values of @@ -63,14 +66,39 @@ is a cubic equation of state which has an extra argument :math:`A_3`, .. math:: - P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . + P = c^2 ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . -Style *taitwater* is Tait's equation of state: +Style *tait/water* is Tait's equation of state: .. math:: P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr]. +Style *tait/general* generalizes this equation of state + +.. math:: + + P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr]. + +where :math:`\gamma` is an exponent. + +Style *ideal/gas* is the ideal gas equation of state + +.. math:: + + P = (\gamma - 1) \rho e + +where :math:`\gamma` is the heat capacity ratio and :math:`e` is the internal energy of +a particle per unit mass. This style is only compatible with atom style rheo/thermal. +Note that when using this style, the speed of sound is no longer constant such that the +value of :math:`c` specified in :doc:`fix rheo ` 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 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/fix_rheo_thermal.rst b/doc/src/fix_rheo_thermal.rst index 214bc1db86..bc0e8dbfd0 100644 --- a/doc/src/fix_rheo_thermal.rst +++ b/doc/src/fix_rheo_thermal.rst @@ -70,13 +70,13 @@ of the energy is used to shift energies. This may be inappropriate in systems with multiple atom types with different specific heats. For each property, one must first define a list of atom types. A wild-card -asterisk can be used in place of or in conjunction with the *types* argument -to set the coefficients for multiple pairs of atom types. This takes the -form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom -types, then an asterisk with no numeric values means all types from 1 to -:math:`N`. A leading asterisk means all types from 1 to n (inclusive). -A trailing asterisk means all types from m to :math:`N` (inclusive). A -middle asterisk means all types from m to n (inclusive). +asterisk can be used in place of or in conjunction with the *types* argument to +set values for multiple atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with no +numeric values means all types from 1 to :math:`N`. A leading asterisk means +all types from 1 to n (inclusive). A trailing asterisk means all types from m +to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). The *types* definition for each property is followed by the style. Currently, the only option is *constant*. Style *constant* simply applies a constant value diff --git a/doc/src/fix_rheo_viscosity.rst b/doc/src/fix_rheo_viscosity.rst index 804059e6f8..4d72d9b209 100644 --- a/doc/src/fix_rheo_viscosity.rst +++ b/doc/src/fix_rheo_viscosity.rst @@ -45,13 +45,12 @@ viscosities for different atom types, but a viscosity must be specified for every atom type. One first defines the atom *types*. A wild-card asterisk can be used in place -of or in conjunction with the *types* argument to set the coefficients for -multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" -or "m\*n". If :math:`N` is the number of atom types, then an asterisk with -no numeric values means all types from 1 to :math:`N`. A leading asterisk -means all types from 1 to n (inclusive). A trailing asterisk means all types -from m to :math:`N` (inclusive). A middle asterisk means all types from m to n -(inclusive). +of or in conjunction with the *types* argument to set values for multiple atom +types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is +the number of atom types, then an asterisk with no numeric values means all types +from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive). +A trailing asterisk means all types from m to :math:`N` (inclusive). A middle +asterisk means all types from m to n (inclusive). The *types* definition is followed by the viscosity style, *vstyle*. Two options are available, *constant* and *power*. Style *constant* simply diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp index 4cebf34eb2..5bf7464ac0 100644 --- a/src/KOKKOS/comm_kokkos.cpp +++ b/src/KOKKOS/comm_kokkos.cpp @@ -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) { k_sendlist.sync(); - CommBrick::forward_comm(fix,size); + CommBrick::forward_comm(fix, size); } else { k_sendlist.sync(); - forward_comm_device(fix,size); + forward_comm_device(fix, size); } } @@ -455,10 +455,10 @@ void CommKokkos::forward_comm_device(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix 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 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) @@ -482,72 +482,94 @@ void CommKokkos::reverse_comm_variable(Fix *fix) /* ---------------------------------------------------------------------- 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(); - CommBrick::forward_comm(compute); + CommBrick::forward_comm(compute, size); } /* ---------------------------------------------------------------------- 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 - 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 - 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(); - CommBrick::reverse_comm(compute); + CommBrick::reverse_comm(compute, size); } /* ---------------------------------------------------------------------- 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) { k_sendlist.sync(); - CommBrick::forward_comm(pair); + CommBrick::forward_comm(pair, size); } else { k_sendlist.sync(); - forward_comm_device(pair); + forward_comm_device(pair, size); } } /* ---------------------------------------------------------------------- */ template -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; 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(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) { k_sendlist.sync(); - CommBrick::reverse_comm(pair); + CommBrick::reverse_comm(pair, size); } else { k_sendlist.sync(); - reverse_comm_device(pair); + reverse_comm_device(pair, size); } } /* ---------------------------------------------------------------------- */ template -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; DAT::tdual_xfloat_1d k_buf_tmp; KokkosBase* pairKKBase = dynamic_cast(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; 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(); - 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(); - CommBrick::reverse_comm(dump); + CommBrick::reverse_comm(dump, size); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/comm_kokkos.h b/src/KOKKOS/comm_kokkos.h index 4fb4dfbe29..42941ff517 100644 --- a/src/KOKKOS/comm_kokkos.h +++ b/src/KOKKOS/comm_kokkos.h @@ -45,24 +45,24 @@ class CommKokkos : public CommBrick { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - 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_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward 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 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 forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward 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 template void forward_comm_device(); template void reverse_comm_device(); - template void forward_comm_device(Pair *pair); - template void reverse_comm_device(Pair *pair); + template void forward_comm_device(Pair *pair, int size=0); + template void reverse_comm_device(Pair *pair, int size=0); template void forward_comm_device(Fix *fix, int size=0); template void exchange_device(); template void borders_device(); diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp index 2e4ca30bed..afddc079f4 100644 --- a/src/KOKKOS/comm_tiled_kokkos.cpp +++ b/src/KOKKOS/comm_tiled_kokkos.cpp @@ -417,42 +417,58 @@ void CommTiledKokkos::borders() /* ---------------------------------------------------------------------- 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 - 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 - 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 - 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) { - CommTiled::forward_comm(fix,size); + CommTiled::forward_comm(fix, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Fix 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 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) { - 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 - 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 - 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 - 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 - 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); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/comm_tiled_kokkos.h b/src/KOKKOS/comm_tiled_kokkos.h index 9033714796..ef226489c8 100644 --- a/src/KOKKOS/comm_tiled_kokkos.h +++ b/src/KOKKOS/comm_tiled_kokkos.h @@ -46,17 +46,17 @@ class CommTiledKokkos : public CommTiled { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - 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_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward 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 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 forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward 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 diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index 63c77e3f3b..589f85e195 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -290,8 +290,7 @@ void ComputeRHEOGrad::compute_peratom() void ComputeRHEOGrad::forward_gradients() { comm_stage = COMMGRAD; - comm_forward = ncomm_grad; - comm->forward_comm(this); + comm->forward_comm(this, ncomm_grad); } /* ---------------------------------------------------------------------- */ @@ -299,8 +298,7 @@ void ComputeRHEOGrad::forward_gradients() void ComputeRHEOGrad::forward_fields() { comm_stage = COMMFIELD; - comm_forward = ncomm_field; - comm->forward_comm(this); + comm->forward_comm(this, ncomm_field); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index eec5659051..bbc00bd8d8 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -97,6 +97,9 @@ void ComputeRHEOInterface::init() auto fixes = modify->get_fix_by_style("rheo/pressure"); fix_pressure = dynamic_cast(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); } @@ -178,7 +181,7 @@ void ComputeRHEOInterface::compute_peratom() dot = 0; 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; } } @@ -192,7 +195,7 @@ void ComputeRHEOInterface::compute_peratom() dot = 0; 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; } } @@ -210,7 +213,7 @@ void ComputeRHEOInterface::compute_peratom() if (status[i] & PHASECHECK) { if (normwf[i] != 0.0) { // 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 { rho[i] = rho0[itype]; } @@ -218,8 +221,7 @@ void ComputeRHEOInterface::compute_peratom() } comm_stage = 1; - comm_forward = 2; - comm->forward_comm(this); + comm->forward_comm(this, 2); } /* ---------------------------------------------------------------------- */ @@ -369,9 +371,8 @@ void ComputeRHEOInterface::store_forces() } // Forward comm forces - comm_forward = 3; comm_stage = 0; - comm->forward_comm(this); + comm->forward_comm(this, 3); } /* ---------------------------------------------------------------------- diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 31b1f65fc2..3a6629bfef 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -90,7 +90,6 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : comm_forward = ncor * Mdim; } - comm_forward_save = comm_forward; corrections_calculated = 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; memory->create(coordination, nmax_store, "rheo:coordination"); 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) { + if (kernel_style == WENDLANDC4) + return calc_w_wendlandc4(r); + if (kernel_style == QUINTIC) + return calc_w_quintic(r); + double w = 0.0; - int corrections_i, corrections_j, corrections; - - if (kernel_style == WENDLANDC4) return calc_w_wendlandc4(r); - - if (kernel_style != QUINTIC) { - corrections_i = check_corrections(i); - corrections_j = check_corrections(j); - corrections = corrections_i & corrections_j; - } else { - corrections = 0; - } + int corrections_i = check_corrections(i); + int corrections_j = check_corrections(j); + int corrections = corrections_i && corrections_j; if (!corrections) - w = calc_w_quintic(r); - else if (kernel_style == RK0) + return calc_w_quintic(r); + + double dx[3] = {delx, dely, delz}; + if (kernel_style == RK0) w = calc_w_rk0(i, j, r); 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) - w = calc_w_rk2(i, j, delx, dely, delz, r); + w = calc_w_rk2(i, j, dx, r); 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) { + 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; - 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); - - 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); + wp = calc_dw_scalar_quintic(delx, dely, delz, r); // Overwrite if there are corrections + double dxij[3] = {delx, dely, delz}; + double dxji[3] = {-delx, -dely, -delz}; + if (kernel_style == RK1) { - if (corrections_i) calc_dw_rk1(i, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk1(j, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk1(i, dxij, r, dWij); + if (corrections_j) calc_dw_rk1(j, dxji, r, dWji); } else if (kernel_style == RK2) { - if (corrections_i) calc_dw_rk2(i, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk2(j, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk2(i, dxij, r, dWij); + if (corrections_j) calc_dw_rk2(j, dxji, r, dWji); } 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 *dW1, double *dW2) +double ComputeRHEOKernel::calc_dw_scalar_quintic(double delx, double dely, double delz, double r) { - double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv; + double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s; s = r * 3.0 * cutinv; @@ -300,14 +301,23 @@ double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, } 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[1] = dely * wprinv; dW1[2] = delz * wprinv; - dW2[0] = -delx * wprinv; - dW2[1] = -dely * wprinv; - dW2[2] = -delz * wprinv; + scale3(-1.0, dW1, dW2); return wp; } @@ -361,9 +371,7 @@ double ComputeRHEOKernel::calc_dw_wendlandc4(double delx, double dely, double de dW1[1] = dely * wprinv; dW1[2] = delz * wprinv; - dW2[0] = -delx * wprinv; - dW2[1] = -dely * wprinv; - dW2[2] = -delz * wprinv; + scale3(-1.0, dW1, dW2); 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; - 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); 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; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; + w = calc_w_quintic(r); 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, - double *dW) +void ComputeRHEOKernel::calc_dw_rk1(int i, double *dx, double r, double *dW) { int a, b; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; 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) //compute derivative operators + zero3(dW); for (a = 0; a < dim; a++) { - dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels 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, - double *dW) +void ComputeRHEOKernel::calc_dw_rk2(int i, double *dx, double r, double *dW) { int a, b; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; 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) //compute derivative operators + zero3(dW); for (a = 0; a < dim; a++) { - dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels 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_tags.clear(); - if (kernel_style == QUINTIC) return; + if (correction_order == -1) return; corrections_calculated = 1; int i, j, ii, jj, inum, jnum, a, b, lapack_error; @@ -666,6 +661,7 @@ void ComputeRHEOKernel::compute_peratom() w = calc_w_quintic(r); rhoj = rho[j]; + if (interface_flag) if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j); @@ -781,7 +777,6 @@ void ComputeRHEOKernel::compute_peratom() // communicate calculated quantities comm_stage = 1; - comm_forward = comm_forward_save; comm->forward_comm(this); } @@ -829,8 +824,7 @@ void ComputeRHEOKernel::compute_coordination() // communicate calculated quantities comm_stage = 0; - comm_forward = 1; - comm->forward_comm(this); + comm->forward_comm(this, 1); } /* ---------------------------------------------------------------------- */ @@ -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 * /*pbc*/) { + int a, b; int m = 0; + for (int i = 0; i < n; i++) { int j = list[i]; 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) { buf[m++] = C0[j]; } else { - for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) buf[m++] = C[j][a][b]; + for (a = 0; a < ncor; a++) + 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) { + int a, b; int m = 0; int last = first + n; + for (int i = first; i < last; i++) { if (comm_stage == 0) { coordination[i] = buf[m++]; @@ -883,8 +881,8 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) if (kernel_style == RK0) { C0[i] = buf[m++]; } else { - for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; + for (a = 0; a < ncor; a++) + for (b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; } } } diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index 8b70509e6a..037e9e3683 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -40,6 +40,7 @@ class ComputeRHEOKernel : public Compute { double calc_w(int, int, double, double, double, double); double calc_dw(int, int, double, double, double, 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_w_wendlandc4(double); double calc_dw_wendlandc4(double, double, double, double, double *, double *); @@ -51,7 +52,7 @@ class ComputeRHEOKernel : public Compute { class FixRHEO *fix_rheo; private: - int comm_stage, comm_forward_save; + int comm_stage; int interface_flag; int lapack_error_flag; std::unordered_set lapack_error_tags; @@ -69,10 +70,10 @@ class ComputeRHEOKernel : public Compute { int check_corrections(int); double calc_w_rk0(int, int, double); - double calc_w_rk1(int, int, double, double, double, double); - double calc_w_rk2(int, int, double, double, double, double); - void calc_dw_rk1(int, double, double, double, double, double *); - void calc_dw_rk2(int, double, double, double, double, double *); + double calc_w_rk1(int, int, double*, double); + double calc_w_rk2(int, int, double*, double); + void calc_dw_rk1(int, double*, double, double*); + void calc_dw_rk2(int, double*, double, double*); }; } // namespace LAMMPS_NS #endif diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index fb5848ac75..4805c21903 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -79,7 +79,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a size_peratom_cols = nvalues; 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 // customize a new keyword by adding to if statement @@ -109,6 +109,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a surface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; } else if (strcmp(arg[iarg], "coordination") == 0) { + coordination_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; } else if (strcmp(arg[iarg], "pressure") == 0) { pressure_flag = 1; @@ -223,6 +224,11 @@ void ComputeRHEOPropertyAtom::compute_peratom() { 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 if (atom->nmax > nmax) { @@ -434,7 +440,7 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); + buf[n] = fix_pressure->calc_pressure(rho[i], i); else buf[n] = 0.0; n += nvalues; @@ -484,7 +490,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (index == index_transpose) - p = fix_pressure->calc_pressure(rho[i], type[i]); + p = fix_pressure->calc_pressure(rho[i], i); else p = 0.0; buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p; diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index 14f28075a3..745423b4ae 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -36,7 +36,7 @@ class ComputeRHEOPropertyAtom : public Compute { private: int nvalues, nmax; 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 *col_index, *col_t_index; double *buf; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 8af60afeda..74193f28f5 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -81,6 +81,8 @@ void ComputeRHEOSurface::init() threshold_splash = fix_rheo->zmin_splash; interface_flag = fix_rheo->interface_flag; + fix_rheo->coordination_flag = 1; + cutsq = cut * cut; // need an occasional half neighbor list @@ -197,10 +199,8 @@ void ComputeRHEOSurface::compute_peratom() // reverse gradC and divr, forward divr comm_stage = 0; - comm_reverse = dim * dim + 1; - comm_forward = 1; - if (newton) comm->reverse_comm(this); - comm->forward_comm(this); + if (newton) comm->reverse_comm(this, dim * dim + 1); + comm->forward_comm(this, 1); // calculate nsurface for local atoms // Note, this isn't forwarded to ghosts @@ -297,10 +297,8 @@ void ComputeRHEOSurface::compute_peratom() // forward/reverse status and rsurface comm_stage = 1; - comm_reverse = 2; - comm_forward = 2; - if (newton) comm->reverse_comm(this); - comm->forward_comm(this); + if (newton) comm->reverse_comm(this, 2); + comm->forward_comm(this, 2); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 8581d9ea8f..248e857fb3 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -27,6 +27,7 @@ #include "error.h" #include "fix_rheo.h" #include "force.h" +#include "math_extra.h" #include "memory.h" #include "neigh_list.h" #include "neigh_request.h" @@ -36,20 +37,22 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; +using namespace MathExtra; /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr), - compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) + Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(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"); + comm_forward = 0; comm_reverse = 3; surface_flag = 0; - nmax_store = atom->nmax; - memory->create(vshift, nmax_store, 3, "rheo:vshift"); + nmax_store = 0; } /* ---------------------------------------------------------------------- */ @@ -73,9 +76,19 @@ void ComputeRHEOVShift::init() compute_surface = fix_rheo->compute_surface; rho0 = fix_rheo->rho0; + shift_type = fix_rheo->shift_type; cut = fix_rheo->cut; cutsq = cut * cut; 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) { 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; } @@ -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] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) { if (status[i] & STATUS_SURFACE) { nx = nsurface[i][0]; 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; last = first + n; - for (int i = first; i < last; i++) { - buf[m++] = vshift[i][0]; - buf[m++] = vshift[i][1]; - buf[m++] = vshift[i][2]; + for (i = first; i < last; i++) + wsame[i] = buf[m++]; +} + +/* ---------------------------------------------------------------------- */ + +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; } @@ -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) { - int i, j, m; + int i, j, a, m; m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - vshift[j][0] += buf[m++]; - vshift[j][1] += buf[m++]; - vshift[j][2] += buf[m++]; + if (comm_stage == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + vshift[j][0] += 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 bytes = 3 * nmax_store * sizeof(double); + + if (cross_type_flag) + bytes += 5 * nmax_store * sizeof(double); + return bytes; } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 485c6525f3..b76d48d5e3 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -31,19 +31,25 @@ class ComputeRHEOVShift : public Compute { void init() override; void init_list(int, class NeighList *) 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; void unpack_reverse_comm(int, int *, double *) override; double memory_usage() override; void correct_surfaces(); + void correct_type_interface(); double **vshift; class FixRHEO *fix_rheo; private: - int nmax_store; + int nmax_store, comm_stage; 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 *wsame, *ct, **cgradt; + int *shift_type; class NeighList *list; class ComputeRHEOInterface *compute_interface; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 4a25a46501..56738c7a0c 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -38,23 +38,25 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; -#if 0 -// publication was removed from documentation static const char cite_rheo[] = - "@article{PalermoInPrep,\n" - " journal = {in prep},\n" - " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n" + "@article{Palermo2024,\n" + " journal = {Physics of Fluids},\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" - " 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"; -#endif /* ---------------------------------------------------------------------- */ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), - compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr), - compute_rhosum(nullptr), compute_vshift(nullptr) + Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), shift_type(nullptr), + compute_grad(nullptr), compute_kernel(nullptr), compute_interface(nullptr), + compute_surface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr) { time_integrate = 1; @@ -68,10 +70,12 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : shift_flag = 0; interface_flag = 0; surface_flag = 0; - oxidation_flag = 0; - self_mass_flag = 0; + coordination_flag = 0; - int i; + rhosum_self_mass_flag = 0; + shift_cross_type_flag = 0; + + int i, nlo, nhi; int n = atom->ntypes; memory->create(rho0, n + 1, "rheo:rho0"); memory->create(csq, n + 1, "rheo:csq"); @@ -111,6 +115,26 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg], "shift") == 0) { 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) { thermal_flag = 1; } else if (strcmp(arg[iarg], "surface/detection") == 0) { @@ -127,14 +151,19 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : } else { error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]); } - iarg += 3; } else if (strcmp(arg[iarg], "interface/reconstruct") == 0) { interface_flag = 1; } else if (strcmp(arg[iarg], "rho/sum") == 0) { rhosum_flag = 1; - } else if (strcmp(arg[iarg], "self/mass") == 0) { - self_mass_flag = 1; + while (iarg < narg) { // optional sub-arguments + if (strcmp(arg[iarg], "self/mass") == 0) { + rhosum_self_mass_flag = 1; + } else { + break; + } + iarg += 1; + } } else if (strcmp(arg[iarg], "density") == 0) { 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); @@ -152,12 +181,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : iarg += 1; } - if (self_mass_flag && (!rhosum_flag)) - error->all(FLERR, "Cannot use self/mass setting without rho/sum"); - -#if 0 if (lmp->citeme) lmp->citeme->add(cite_rheo); -#endif } /* ---------------------------------------------------------------------- */ @@ -173,6 +197,7 @@ FixRHEO::~FixRHEO() memory->destroy(csq); memory->destroy(rho0); + memory->destroy(shift_type); } /* ---------------------------------------------------------------------- @@ -192,7 +217,7 @@ void FixRHEO::post_constructor() if (rhosum_flag) { compute_rhosum = dynamic_cast( - 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; } @@ -375,7 +400,6 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Shifting atoms if (shift_flag) { for (i = 0; i < nlocal; i++) { - if (status[i] & STATUS_NO_SHIFT) continue; if (status[i] & PHASECHECK) continue; diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 24e86b21d6..aa80bc938e 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -41,19 +41,25 @@ class FixRHEO : public Fix { // Model parameters double cut; double *rho0, *csq; - int self_mass_flag; int zmin_kernel, zmin_surface, zmin_splash; int kernel_style, surface_style; double divr_surface; - // Accessory fixes/computes + // Settings flags int thermal_flag; int rhosum_flag; - int shift_flag; int interface_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 pressure_fix_defined; int thermal_fix_defined; diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 10bf6ff2c9..d972f19ffa 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL }; +enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL , IDEAL }; 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) : Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr), - rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr), - fix_rheo(nullptr) + rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), gamma(nullptr), + pressure_style(nullptr), fix_rheo(nullptr) { if (narg < 4) error->all(FLERR, "Illegal fix command"); comm_forward = 1; + variable_csq = 0; + invertible_pressure = 1; + background_flag = 0; // Currently can only have one instance of fix rheo/pressure 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(tpower, n + 1, "rheo:tpower"); 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; while (iarg < narg) { @@ -68,16 +76,14 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg + 1], "tait/water") == 0) { for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; } 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 pbackground_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 2; + iarg += 1; for (i = nlo; i <= nhi; i++) { pressure_style[i] = TAITGENERAL; tpower[i] = tpower_one; - pbackground[i] = pbackground_one; } } else if (strcmp(arg[iarg + 1], "cubic") == 0) { 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; 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 { error->all(FLERR, "Illegal fix command, {}", arg[iarg]); } @@ -110,6 +141,7 @@ FixRHEOPressure::~FixRHEOPressure() memory->destroy(c_cubic); memory->destroy(tpower); 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 dr, rr3, rho_ratio; + int type = atom->type[i]; if (pressure_style[type] == LINEAR) { p = csq[type] * (rho - rho0[type]); @@ -214,16 +247,25 @@ double FixRHEOPressure::calc_pressure(double rho, int type) } else if (pressure_style[type] == TAITGENERAL) { rho_ratio = rho * rho0inv[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; } /* ---------------------------------------------------------------------- */ -double FixRHEOPressure::calc_rho(double p, int type) +double FixRHEOPressure::calc_rho(double p, int i) { double rho = 0.0; + int type = atom->type[i]; + + if (background_flag) + p -= pbackground[type]; if (pressure_style[type] == LINEAR) { 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(csq[type], -SEVENTH); } else if (pressure_style[type] == TAITGENERAL) { - p -= pbackground[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(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; } + +/* ---------------------------------------------------------------------- */ + +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; +} diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index ca165b1ed5..c15e9b4ddc 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -36,10 +36,13 @@ class FixRHEOPressure : public Fix { void unpack_forward_comm(int, int, double *) override; double calc_pressure(double, int); double calc_rho(double, int); + double calc_csq(double, int); + int variable_csq; + int invertible_pressure; private: - double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground; - int *pressure_style; + double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground, *gamma; + int *pressure_style, background_flag; class FixRHEO *fix_rheo; }; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index aa0c685aaa..5e06c45a3c 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -80,8 +80,8 @@ void PairRHEO::compute(int eflag, int vflag) int pair_force_flag, pair_rho_flag, pair_avisc_flag; int fluidi, fluidj; 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, - kappa_ave, dT_prefactor; + double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj; + double eta_ave, kappa_ave, dT_prefactor; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; double *dWij, *dWji; 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]; } - 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; pair_rho_flag = 0; @@ -221,7 +227,7 @@ void PairRHEO::compute(int eflag, int vflag) if (fluidi && (!fluidj)) { compute_interface->correct_v(vj, vi, j, i); 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))) 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) { compute_interface->correct_v(vi, vj, i, j); 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)) 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; 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 @@ -480,6 +494,8 @@ void PairRHEO::setup() csq = fix_rheo->csq; rho0 = fix_rheo->rho0; + variable_csq = fix_pressure->variable_csq; + if (cutk != fix_rheo->cut) error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, fix_rheo->cut); diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index 444fcc2cb4..f1fcd10bf8 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -45,6 +45,7 @@ class PairRHEO : public Pair { int rho_damp_flag; int thermal_flag; int interface_flag; + int variable_csq; int harmonic_means_flag; diff --git a/src/comm.cpp b/src/comm.cpp index 1b1546f893..98f1ebf28a 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -228,13 +228,13 @@ void Comm::init() } for (const auto &compute : modify->get_compute_list()) { - maxforward = MAX(maxforward,compute->comm_forward); - maxreverse = MAX(maxreverse,compute->comm_reverse); + maxforward = MAX(maxforward, compute->comm_forward); + maxreverse = MAX(maxreverse, compute->comm_reverse); } for (const auto &dump: output->get_dump_list()) { - maxforward = MAX(maxforward,dump->comm_forward); - maxreverse = MAX(maxreverse,dump->comm_reverse); + maxforward = MAX(maxforward, dump->comm_forward); + maxreverse = MAX(maxreverse, dump->comm_reverse); } if (force->newton == 0) maxreverse = 0; diff --git a/src/comm.h b/src/comm.h index 8515866dff..13e44eda57 100644 --- a/src/comm.h +++ b/src/comm.h @@ -87,17 +87,17 @@ class Comm : protected Pointers { // forward/reverse comm from a Pair, Bond, Fix, Compute, Dump - virtual void forward_comm(class Pair *) = 0; - virtual void reverse_comm(class Pair *) = 0; - virtual void forward_comm(class Bond *) = 0; - virtual void reverse_comm(class Bond *) = 0; + virtual void forward_comm(class Pair *, int size = 0) = 0; + virtual void reverse_comm(class Pair *, int size = 0) = 0; + virtual void forward_comm(class Bond *, int size = 0) = 0; + virtual void reverse_comm(class Bond *, 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_variable(class Fix *) = 0; - virtual void forward_comm(class Compute *) = 0; - virtual void reverse_comm(class Compute *) = 0; - virtual void forward_comm(class Dump *) = 0; - virtual void reverse_comm(class Dump *) = 0; + virtual void forward_comm(class Compute *, int size = 0) = 0; + virtual void reverse_comm(class Compute *, int size = 0) = 0; + virtual void forward_comm(class Dump *, int size = 0) = 0; + virtual void reverse_comm(class Dump *, int size = 0) = 0; // forward comm of an array diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index c75b183584..1e9bfe6aca 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -984,16 +984,21 @@ void CommBrick::borders() /* ---------------------------------------------------------------------- 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; MPI_Request request; - int nsize = pair->comm_forward; + if (size) nsize = size; + else nsize = pair->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1021,16 +1026,21 @@ void CommBrick::forward_comm(Pair *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; 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--) { @@ -1058,16 +1068,21 @@ void CommBrick::reverse_comm(Pair *pair) /* ---------------------------------------------------------------------- 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; MPI_Request request; - int nsize = bond->comm_forward; + if (size) nsize = size; + else nsize = bond->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1095,16 +1110,21 @@ void CommBrick::forward_comm(Bond *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; 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--) { @@ -1175,10 +1195,10 @@ void CommBrick::forward_comm(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix 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 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) @@ -1259,16 +1279,21 @@ void CommBrick::reverse_comm_variable(Fix *fix) /* ---------------------------------------------------------------------- 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; MPI_Request request; - int nsize = compute->comm_forward; + if (size) nsize = size; + else nsize = compute->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1297,16 +1322,21 @@ void CommBrick::forward_comm(Compute *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; MPI_Request request; - int nsize = compute->comm_reverse; + if (size) nsize = size; + else nsize = compute->comm_reverse; for (iswap = nswap-1; iswap >= 0; iswap--) { @@ -1334,16 +1364,21 @@ void CommBrick::reverse_comm(Compute *compute) /* ---------------------------------------------------------------------- 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; MPI_Request request; - int nsize = dump->comm_forward; + if (size) nsize = size; + else nsize = dump->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1372,16 +1407,21 @@ void CommBrick::forward_comm(Dump *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; MPI_Request request; - int nsize = dump->comm_reverse; + if (size) nsize = size; + else nsize = dump->comm_reverse; for (iswap = nswap-1; iswap >= 0; iswap--) { diff --git a/src/comm_brick.h b/src/comm_brick.h index 809b3856e4..8878b44f4f 100644 --- a/src/comm_brick.h +++ b/src/comm_brick.h @@ -32,17 +32,17 @@ class CommBrick : public Comm { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - 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_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward 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 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 forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward 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 *extract(const char *, int &) override; diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index f4657c4790..e8b5d19fa5 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -1382,14 +1382,19 @@ void CommTiled::borders() /* ---------------------------------------------------------------------- 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++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1426,14 +1431,19 @@ void CommTiled::forward_comm(Pair *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--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1466,14 +1476,18 @@ void CommTiled::reverse_comm(Pair *pair) /* ---------------------------------------------------------------------- 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) +void CommTiled::forward_comm(Bond *bond, int size) { - 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++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1510,14 +1524,19 @@ void CommTiled::forward_comm(Bond *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--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1598,10 +1617,10 @@ void CommTiled::forward_comm(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix 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 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) @@ -1654,14 +1673,19 @@ void CommTiled::reverse_comm_variable(Fix * /*fix*/) /* ---------------------------------------------------------------------- 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++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1697,14 +1721,19 @@ void CommTiled::forward_comm(Compute *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--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1738,14 +1767,19 @@ void CommTiled::reverse_comm(Compute *compute) /* ---------------------------------------------------------------------- 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++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1781,14 +1815,19 @@ void CommTiled::forward_comm(Dump *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--) { nsend = nsendproc[iswap] - sendself[iswap]; diff --git a/src/comm_tiled.h b/src/comm_tiled.h index e150c710fc..64b80d8d18 100644 --- a/src/comm_tiled.h +++ b/src/comm_tiled.h @@ -32,17 +32,17 @@ class CommTiled : public Comm { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - 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_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward 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 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 forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward 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