From 7ea0dc3996db3e9315d06fae83162ec2990ad114 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 29 Mar 2024 19:00:14 -0600 Subject: [PATCH] Adding more oxidation features + doc pages --- doc/src/Howto_rheo.rst | 99 +++++++++++++++++++++- doc/src/Packages_details.rst | 36 ++++++++ doc/src/Packages_list.rst | 5 ++ doc/src/fix_rheo.rst | 108 ++++++++++++++++++++---- doc/src/pair_rheo.rst | 26 +++++- src/RHEO/atom_vec_rheo.cpp | 2 +- src/RHEO/atom_vec_rheo_thermal.cpp | 2 +- src/RHEO/bond_rheo_shell.cpp | 80 ++++++++++++++++++ src/RHEO/bond_rheo_shell.h | 6 ++ src/RHEO/compute_rheo_grad.cpp | 2 +- src/RHEO/compute_rheo_kernel.cpp | 2 +- src/RHEO/compute_rheo_property_atom.cpp | 7 +- src/RHEO/compute_rheo_property_atom.h | 3 +- src/RHEO/compute_rheo_rho_sum.cpp | 5 ++ src/RHEO/compute_rheo_vshift.cpp | 2 + src/RHEO/fix_rheo.cpp | 2 +- src/RHEO/fix_rheo.h | 1 + src/RHEO/fix_rheo_oxidation.cpp | 29 ++++++- src/RHEO/fix_rheo_oxidation.h | 3 + src/RHEO/fix_rheo_pressure.cpp | 2 +- src/RHEO/fix_rheo_thermal.cpp | 35 +++++--- src/RHEO/fix_rheo_viscosity.cpp | 2 +- src/RHEO/pair_rheo.cpp | 2 +- src/RHEO/pair_rheo_solid.cpp | 5 ++ 24 files changed, 423 insertions(+), 43 deletions(-) diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index c55631455b..146716ba18 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -1,4 +1,99 @@ Reproducing hydrodynamics and elastic objects (RHEO) -====================== +==================================================== -Text +The RHEO package is built around an implementation of smoothed particle +hydrodynamics (SPH) coupled to the :doc:`BPM package ` to model +solid elements of a system. The SPH solver supports many advanced options +including reproducing kernels, particle shifting, free surface identification, +and solid surface reconstruction. To model fluid-solid systems, the status of +particles can dynamically change between a fluid and solid state, e.g. during +melting/solidification, which determines how they interact and their physical +behavior. The package is designed with modularity in mind, so one can easily +turn various features on/off, adjust physical details of the system, or +develop new capabilities. Additional numerical details can be found in +:ref:`(Palermo) ` and :ref:`(Clemmer) `. + +---------- + +At the core of the package is :doc:`fix rheo ` which integrates +particle trajectories and controls many optional features (e.g. the use +of reproducing kernels). In conjunction to fix rheo, one must specify an +instance of :doc:`fix rheo/pressure ` and +:doc:`fix rheo/viscosity ` to define a pressure equation +of state and viscosity model, respectively. Optionally, one can model +a heat equation with :doc:`fix rheo/thermal`, which also allows the user +to specify equations for a particle's thermal conductivity, specific heat, +latent heat, and melting temperature. Fix rheo must be defined prior to all +other RHEO fixes. + +Typically, RHEO requires atom style rheo. In addition to typical atom +properties like positions and forces, particles store a local density, +viscosity, pressure, and status. If thermal evolution is modeled, one must +use atom style rheo/thermal which also include a local temperature and +conductivity. The status variable uses bitmasking to track various +properties of a particle such as its current phase (fluid or solid) and its +location relative to a surface. Many of these properties (and others) can +be easily accessed using +:doc:`compute rheo/property/atom `. + +Fluid interactions, including pressure forces, viscous forces, and heat exchange, +are calculated using :doc:`pair rheo `. Unlike typical pair styles, +pair rheo ignores the :doc:`special bond ` settings. Instead, +it determines whether to calculate forces based on the status of particles: +hydrodynamic forces are only calculated if a fluid particle is involved. + +---------- + +To model elastic objects, there are current two mechanisms in RHEO, one designed +for bulk solid bodies and the other for thin shells. Both mechanisms rely on +overlaying bonds and therefore require a hybrid of atom style bond and rheo +(or rheo/thermal). + +To create an elastic solid body, one has to (a) change the status of constituent +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, +looks at a particles fluid/solid status 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. + +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. Using the *react* option, one can specify a maximum +bond length and a bond type. Then, when solidifying, particles will 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 then delete +any bonds of the specified type when reverting to a fluid state. Special bonds are +updated as bonds are created/broken. + +The other option for elastic objects is an elastic shell that is nominally much +thinner than a particle diameter, e.g. a oxide skin which gradually forms over time +on the surface of a fluid. Currently, this is implemented using +:doc:`fix rheo/oxidaton ` and bond style +:doc:`rheo/shell `. Essentially, fix rheo/oxidaton creates candidate +bonds of a specified type between surface fluid particles within a specified distance. +a newly created rheo/shell bond will then start a timer. While the timer is counting +down, the bond will delete itself if particles move too far apart or move away from the +surface. However, if the timer reaches a user-defined threshold, then the bond will +activate and apply additional forces to the fluid particles. Bond style rheo/shell +then operates very similarly to a BPM bond style, storing a reference length and +breaking if stretched too far. Unlike the above method, this option does not remove +the underlying fluid interactions (although particle shifting is turned off) and does +not modify special bond settings of particles. + +While these two options are not expected to be appropriate for every multiphase system, +either framework can be modified to create more suitable models (e.g. by changing the +criteria for creating/deleting a bond or altering force calculations). + +---------- + +.. _howto-howto_rheo_palermo: + +**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation. + +.. _howto-howto_rheo_clemmer: + +**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index a3d65d9d65..39a9deef63 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -101,6 +101,7 @@ page gives those details. * :ref:`QEQ ` * :ref:`QMMM ` * :ref:`QTB ` + * :ref:`RHEO ` * :ref:`REACTION ` * :ref:`REAXFF ` * :ref:`REPLICA ` @@ -2571,6 +2572,41 @@ another set. ---------- +.. _PKG-RHEO: + +RHEO package +------------ + +**Contents:** + +Pair styles, bond styles, fixes, and computes for reproducing +hydrodynamics and elastic objects. See the +:doc:`Howto rheo ` page for an overview. + +**Authors:** Joel T. Clemmer (Sandia National Labs), +Thomas C. O'Connor (Carnegie Mellon University) + +.. versionadded:: TBD + +**Supporting info:** + +* src/RHEO filenames -> commands +* :doc:`Howto_rheo ` +* :doc:`atom_style rheo ` +* :doc:`atom_style rheo/thermal ` +* :doc:`bond_style rheo/shell ` +* :doc:`compute rheo/property/atom ` +* :doc:`fix rheo ` +* :doc:`fix rheo/oxidation ` +* :doc:`fix rheo/pressure ` +* :doc:`fix rheo/thermal ` +* :doc:`fix rheo/viscosity ` +* :doc:`pair_style rheo ` +* :doc:`pair_style rheo/solid ` +* examples/rheo + +---------- + .. _PKG-RIGID: RIGID package diff --git a/doc/src/Packages_list.rst b/doc/src/Packages_list.rst index c0a1164513..4e1d32385f 100644 --- a/doc/src/Packages_list.rst +++ b/doc/src/Packages_list.rst @@ -403,6 +403,11 @@ whether an extra library is needed to build and use the package: - :doc:`fix qtb ` :doc:`fix qbmsst ` - qtb - no + * - :ref:`RHEO ` + - reproducing hydrodynamics and elastic objects + - :doc:`Howto rheo ` + - rheo + - no * - :ref:`REACTION ` - chemical reactions in classical MD - :doc:`fix bond/react ` diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index c61d1939db..76d4ae3972 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -8,40 +8,107 @@ Syntax .. parsed-literal:: - fix ID group-ID rheo cut kstyle keyword values... + fix ID group-ID rheo cut kstyle zmin keyword values... * ID, group-ID are documented in :doc:`fix ` command * rheo = style name of this fix command -* cut = *quintic* or *RK0* or *RK1* or *RK2* +* cut = cutoff for the kernel (distance) +* 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 = *shift* or *thermal* or *surface/detection* or *interface/reconstruction* or - *rho/sum* or *density* or *sound/squared* +* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or + *shift* or *rho/sum* or *density* or *speed/sound* .. parsed-literal:: - *shift* values = none, turns on velocity shifting *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 (unitless) - *limit/splash* = threshold for splash particles (unitless) - *interface/reconstruct* values = none, reconstructs interfaces with solid particles + *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 - *density* values = *rho0* (density) - *sound/squared* values = *csq* (velocity\^2) + *density* values = *rho01*, ... *rho0N* (density) + *speed/sound* values = *cs0*, ... *csN* (velocity) Examples """""""" .. code-block:: LAMMPS - fix 1 all rheo 1.0 quintic thermal density 0.1 sound/squared 10.0 + fix 1 all rheo 1.0 quintic thermal density 0.1 speed/sound 10.0 fix 1 all rheo 1.0 RK1 shift surface/detection coordination 40 Description """"""""""" -Fix description... +Perform time integration for RHEO particles, updating positions, velocities, +and densities. For a detailed breakdown of the integration timestep and +numerical details, see :ref:`(Palermo) `. For an +overview 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 +spline function commonly used in SPH. The other options, *RK0*, *RK1*, and +*RK2*, are zeroth, first, and second order reproducing. To generate a reproducing kernel, a particle must have sufficient neighbors inside the +kernel cutoff distance (a coordination number) to accurately calculate +moments. This threshold is set by *zmin*. If reproducing kernels are +requested but a particle has fewer neighbors, then it will revert to a +non-reproducing quintic kernel until it gains more neighbors. + +To model temperature evolution, one must specify the *thermal* keyword, +define a separate instance of :doc:`fix rheo/thermal `, +and use atom style rheo/thermal. + +By default, the density of solid RHEO particles does not evolve and forces +with fluid particles are calculated using the current velocity of the solid +particle. If the *interface/reconstruct* keyword is used, then the density +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. 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. 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 free surface, or isolated +from other particles in a splash or droplet. Shifting is then disabled in +the direction away from the free surface to prevent it from diffusing +particles away from the bulk fluid. Surface detection can also be used +to control surface-nucleated effects like oxidation when used in combination +with :doc:`fix rheo/oxidation `. + +The *surface/detection* keyword takes three arguments: *sdstyle*, *limit*, +and *limi/splash*. The first, *sdstyle*, specifies whether surface particles +are identified using a coordination number (*coordination*) or the divergence +of the local particle positions (*divergence*). The threshold value for a +surface particle for either of these criteria is set by the numerical value +of *limit*. Additionally, if a particle's coordination number is too low, +i.e. if it has separated off from the bulk in a droplet, it is not possible +to define surfaces and a particle is classified as a splash. The coordination +threshold for this classification is set by the numerical value of +*limit/splash*. + +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 *density* is used to specify the equilbrium density of each of the N +particle types. It must be followed by N numerical values specifying each +type's equilibrium density *rho0*. + +The *density* is used to specify the speed of sound of each of the N particle +types. It must be followed by N numerical values specifying each type's speed +of sound *cs*. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -55,13 +122,14 @@ the :doc:`run ` command. This fix is not invoked during :doc:`energy minim Restrictions """""""""""" -This fix must be used with an atom style that includes density -such as atom_style rheo or rheo/thermal. This fix must be used in -conjuction with :doc:`fix rheo/pressure `. and +This fix must be used with atom style rheo or rheo/thermal. +This fix must be used in conjuction with +:doc:`fix rheo/pressure `. and :doc:`fix rheo/viscosity `, If the *thermal* setting is used, there must also be an instance of :doc:`fix rheo/thermal `. The fix group must be -set to all. +set to all. Only one instance of fix rheo may be defined and it +must be defined prior to all other RHEO fixes. This fix is part of the RHEO package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. @@ -78,4 +146,10 @@ Related commands Default """"""" -*rho0* and *csq* are set to 1.0. +*rho0* and *cs* are set to 1.0 for all atom types. + +---------- + +.. _howto-howto_rheo_palermo: + +**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation. diff --git a/doc/src/pair_rheo.rst b/doc/src/pair_rheo.rst index 6f706a77ac..f6c3d9e3ba 100644 --- a/doc/src/pair_rheo.rst +++ b/doc/src/pair_rheo.rst @@ -16,8 +16,8 @@ Syntax .. parsed-literal:: - *rho/damp* args = density damping prefactor :math:`\xi` (units?) - *artificial/visc* args = artificial viscosity prefactor :math:`\zeta` (units?) + *rho/damp* args = density damping prefactor :math:`\xi` + *artificial/visc* args = artificial viscosity prefactor :math:`\zeta` *harmonic/means* args = none Examples @@ -31,7 +31,27 @@ Examples Description """"""""""" -pair style... +Pair style *rheo* computes pressure and viscous forces between particles +in the :doc:`rheo package `. If thermal evolution is turned +on in :doc:`fix rheo `, then the pair style also calculates +heat exchanged between particles. + +The *artificial/viscosity* keyword is used to specify the magnitude +:math:`\zeta` of an optional artificial viscosity contribution to forces. +This factor can help stabilize simulations by smoothing out small length +scale variations in velocity fields. + +The *rho/damp* keyword is used to specify the magnitude :math:`\xi` of +an optional pairwise damping term between the density of particles. This +factor can help stabilize simulations by smoothing out small length +scale variations in density fields. + +If particles have different viscosities or conductivities, the +*harmonic/means* keyword changes how they are averaged before calculating +pairwise forces or heat exchanges. By default, an arithmetic averaged is +used, however, a harmonic mean may improve stability in multiphase systems +with large disparities in viscosities. This keyword has no effect on +results if viscosities and conductivities are constant. No coefficients are defined for each pair of atoms types via the :doc:`pair_coeff ` command as in the examples diff --git a/src/RHEO/atom_vec_rheo.cpp b/src/RHEO/atom_vec_rheo.cpp index ec44a230ec..843269a717 100644 --- a/src/RHEO/atom_vec_rheo.cpp +++ b/src/RHEO/atom_vec_rheo.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL) ----------------------------------------------------------------------- */ #include "atom_vec_rheo.h" diff --git a/src/RHEO/atom_vec_rheo_thermal.cpp b/src/RHEO/atom_vec_rheo_thermal.cpp index 26394c9175..7174d4cb66 100644 --- a/src/RHEO/atom_vec_rheo_thermal.cpp +++ b/src/RHEO/atom_vec_rheo_thermal.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL) ----------------------------------------------------------------------- */ #include "atom_vec_rheo_thermal.h" diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index 6a9f99d8b1..6a71136d9d 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL) +----------------------------------------------------------------------- */ + #include "bond_rheo_shell.h" #include "atom.h" @@ -40,17 +45,38 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr) { partial_flag = 1; + comm_reverse = 1; tform = rmax = -1; single_extra = 1; svector = new double[1]; + + // For nbond, create an instance of fix property atom + // Need restarts + exchanging with neighbors since it needs to persist + // between timesteps (fix property atom will handle callbacks) + + int tmp1, tmp2; + index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); + if (index_nb == -1) { + id_fix = utils::strdup("bond_rheo_shell_fix_property_atom"); + modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix)); + index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); + } + nbond = atom->ivector[index_nb]; + + //Store non-persistent per atom quantities, intermediate + + nmax_store = atom->nmax; + memory->create(dbond, nmax_store, "rheo/react:dbond"); } /* ---------------------------------------------------------------------- */ BondRHEOShell::~BondRHEOShell() { + if (modify->nfix) modify->delete_fix(id_fix); + delete[] id_fix; delete[] svector; if (allocated) { @@ -59,6 +85,8 @@ BondRHEOShell::~BondRHEOShell() memory->destroy(ecrit); memory->destroy(gamma); } + + memory->destroy(dbond); } /* ---------------------------------------------------------------------- @@ -151,6 +179,15 @@ void BondRHEOShell::compute(int eflag, int vflag) double **bondstore = fix_bond_history->bondstore; + if (atom->nmax > nmax_store){ + nmax_store = atom->nmax; + memory->destroy(dbond); + memory->create(dbond, nmax_store, "rheo/shell:dbond"); + } + + size_t nbytes = nmax_store * sizeof(int); + memset(&dbond[0], 0, nbytes); + for (n = 0; n < nbondlist; n++) { // skip bond if already broken @@ -196,6 +233,8 @@ void BondRHEOShell::compute(int eflag, int vflag) if (bondstore[n][1] >= tform) { bondstore[n][0] = r; r0 = r; + if (newton_bond || i1 < nlocal) dbond[i1] ++; + if (newton_bond || i2 < nlocal) dbond[i2] ++; } else { continue; } @@ -205,6 +244,8 @@ void BondRHEOShell::compute(int eflag, int vflag) if (fabs(e) > ecrit[type]) { bondlist[n][2] = 0; process_broken(i1, i2); + if (newton_bond || i1 < nlocal) dbond[i1] --; + if (newton_bond || i2 < nlocal) dbond[i2] --; continue; } @@ -233,6 +274,17 @@ void BondRHEOShell::compute(int eflag, int vflag) if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz); } + + + // Communicate changes in nbond + if (newton_bond) comm->reverse_comm(this); + + for(i = 0; i < nlocal; i++) { + nbond[i] += dbond[i]; + + // If it has bonds, no shifting + if (nbond[i] != 0) status[i] |= STATUS_NO_SHIFT; + } } /* ---------------------------------------------------------------------- */ @@ -419,6 +471,34 @@ void BondRHEOShell::read_restart_settings(FILE *fp) MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, world); } + +/* ---------------------------------------------------------------------- */ + +int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last; + m = 0; + last = first + n; + + for (i = first; i < last; i++) { + buf[m++] = dbond[i]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void BondRHEOShell::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + dbond[j] += buf[m++]; + } +} + /* ---------------------------------------------------------------------- */ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce) diff --git a/src/RHEO/bond_rheo_shell.h b/src/RHEO/bond_rheo_shell.h index 562be6d9a6..513d481eeb 100644 --- a/src/RHEO/bond_rheo_shell.h +++ b/src/RHEO/bond_rheo_shell.h @@ -37,12 +37,18 @@ class BondRHEOShell : public BondBPM { void read_restart(FILE *) override; void write_restart_settings(FILE *) override; void read_restart_settings(FILE *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; double single(int, double, int, int, double &) override; protected: double *k, *ecrit, *gamma; double tform, rmax; + int *dbond, *nbond; + int index_nb, nmax_store; + char *id_fix; + void process_ineligibility(int, int); void allocate(); void store_data(); diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index acfc01d793..5ed43a421a 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL), Thomas O'Connor (CMU) ----------------------------------------------------------------------- */ #include "compute_rheo_grad.h" diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 6f58d79243..bd16a89b6a 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL), Thomas O'Connor (CMU) ----------------------------------------------------------------------- */ #include "compute_rheo_kernel.h" diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 380ff398d8..2d6ff1e55a 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -56,7 +56,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (nvalues == 1) size_peratom_cols = 0; else size_peratom_cols = nvalues; - pressure_flag = thermal_flag = interface_flag = surface_flag = shift_flag = 0; + pressure_flag = thermal_flag = interface_flag = surface_flag = shift_flag = shell_flag = 0; // parse input values // customize a new keyword by adding to if statement @@ -112,6 +112,9 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a arg[iarg], atom->get_style()); pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; thermal_flag = 1; + } else if (strcmp(arg[iarg],"nbond/shell") == 0) { + shell_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_nbond_shell; } else { avec_index[i] = atom->avec->property_atom(arg[iarg]); if (avec_index[i] < 0) @@ -156,6 +159,8 @@ void ComputeRHEOPropertyAtom::init() error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo"); if (thermal_flag && !(fix_rheo->thermal_flag)) error->all(FLERR, "Cannot request thermal property without fix rheo/thermal"); + if (shell_flag && !(fix_rheo->oxidation_flag)) + error->all(FLERR, "Cannot request number of shell bonds without fix rheo/oxidation"); compute_interface = fix_rheo->compute_interface; compute_kernel = fix_rheo->compute_kernel; diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index f3596fbbf9..fd73b5883f 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -34,7 +34,7 @@ class ComputeRHEOPropertyAtom : public Compute { private: int nvalues, nmax; - int pressure_flag, thermal_flag, interface_flag, surface_flag, shift_flag; + int pressure_flag, thermal_flag, interface_flag, surface_flag, shift_flag, shell_flag; int *avec_index; int *col_index; double *buf; @@ -54,6 +54,7 @@ class ComputeRHEOPropertyAtom : public Compute { void pack_shift_v(int); void pack_gradv(int); void pack_pressure(int); + void pack_nbond_shell(int); void pack_atom_style(int); int get_vector_index(char*); diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index 82d3aa4bc6..8322ad4ad4 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL) +----------------------------------------------------------------------- */ + #include "compute_rheo_rho_sum.h" #include "atom.h" diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 569c8569f7..533287911a 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -147,6 +147,8 @@ void ComputeRHEOVShift::compute_peratom() fluidj = !(status[j] & PHASECHECK); if ((!fluidi) && (!fluidj)) continue; + + // Will skip shifting in FixRHEO initial integrate, but also skip here to save time if ((status[i] & STATUS_NO_SHIFT) && (status[j] & STATUS_NO_SHIFT)) continue; dx[0] = xtmp - x[j][0]; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index beba940174..b0cb1513fc 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -221,7 +221,7 @@ void FixRHEO::setup_pre_force(int /*vflag*/) { // Check to confirm accessory fixes do not preceed FixRHEO // Note: fixes set this flag in setup_pre_force() - if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined) + if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || oxidation_fix_defined) error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes"); // Calculate surfaces diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 8ec28c7d0e..33fd0084db 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -55,6 +55,7 @@ class FixRHEO : public Fix { int viscosity_fix_defined; int pressure_fix_defined; int thermal_fix_defined; + int oxidation_fix_defined; class ComputeRHEOGrad *compute_grad; class ComputeRHEOKernel *compute_kernel; diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index bd7babedbb..1628cef13f 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -60,6 +60,7 @@ FixRHEOOxidation::~FixRHEOOxidation() int FixRHEOOxidation::setmask() { int mask = 0; + mask |= PRE_FORCE; mask |= POST_INTEGRATE; return mask; } @@ -70,7 +71,7 @@ void FixRHEOOxidation::init() { auto fixes = modify->get_fix_by_style("^rheo$"); if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/oxidation"); - class FixRHEO *fix_rheo = dynamic_cast(fixes[0]); + fix_rheo = dynamic_cast(fixes[0]); double cut_kernel = fix_rheo->h; if (cut > cut_kernel) @@ -102,6 +103,26 @@ void FixRHEOOxidation::init_list(int /*id*/, NeighList *ptr) /* ---------------------------------------------------------------------- */ +void FixRHEOOxidation::setup_pre_force(int /*vflag*/) +{ + // Not strictly required that this fix be after FixRHEO, + // but enforce to be consistent with other RHEO fixes + fix_rheo->oxidation_fix_defined = 1; + + if (!fix_rheo->surface_flag) error->all(FLERR, + "fix rheo/oxidation requires surface calculation in fix rheo"); + + pre_force(0); +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOThermal::pre_force(int /*vflag*/) +{ +} + +/* ---------------------------------------------------------------------- */ + void FixRHEOOxidation::post_integrate() { int i, j, n, ii, jj, inum, jnum, bflag; @@ -187,4 +208,10 @@ void FixRHEOOxidation::post_integrate() } } } + + //todo: + // allow option to create near-surface bonds + // extract # of bonds in property/atom + // check bond style shell used, get index to bonds, share with compute property atom + // add type option to compute nbond/atom } diff --git a/src/RHEO/fix_rheo_oxidation.h b/src/RHEO/fix_rheo_oxidation.h index 4c81605611..ca36a6fdf9 100644 --- a/src/RHEO/fix_rheo_oxidation.h +++ b/src/RHEO/fix_rheo_oxidation.h @@ -33,6 +33,8 @@ class FixRHEOOxidation : public Fix { int setmask() override; void init() override; void init_list(int, class NeighList *) override; + void setup_pre_force(int) override; + void pre_force(int) override; void post_integrate() override; private: @@ -40,6 +42,7 @@ class FixRHEOOxidation : public Fix { double cut, cutsq; class NeighList *list; + class FixRHEO *fix_rheo; //class FixBondHistory *fix_bond_history; }; diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 8c523b2b35..d6dea8aa1a 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL), Thomas O'Connor (CMU) ----------------------------------------------------------------------- */ #include "fix_rheo_pressure.h" diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index c6d10bcc79..7d9aff5424 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -498,21 +498,36 @@ void FixRHEOThermal::break_bonds() int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - // Rapidly delete all bonds for local atoms that melt (no shifting) + // Rapidly delete all bonds for local atoms that melt of a given type for (int i = 0; i < nlocal; i++) { if (!(status[i] & STATUS_MELTING)) continue; - for (m = 0; m < num_bond[i]; m++) { - j = atom->map(bond_atom[i][m]); - bond_type[i][m] = 0; + for (m = (num_bond[i] - 1); m >= 0; m--) { + if (bond_type[i][m] != btype) continue; - if (n_histories > 0) - for (auto &ihistory: histories) - dynamic_cast(ihistory)->delete_history(i, m); + j = atom->map(bond_atom[i][m]); + + nmax = num_bond[i] - 1; + if (m == nmax) { + if (n_histories > 0) + for (auto &ihistory: histories) + dynamic_cast(ihistory)->delete_history(i, m); + } else { + bond_type[i][m] = bond_type[i][nmax]; + bond_atom[i][m] = bond_atom[i][nmax]; + if (n_histories > 0) { + for (auto &ihistory: histories) { + auto fix_bond_history = dynamic_cast (ihistory); + fix_bond_history->shift_history(i, m, nmax); + fix_bond_history->delete_history(i, nmax); + } + } + } + bond_type[i][nmax] = 0; + num_bond[i]--; if (fix_update_special_bonds) fix_update_special_bonds->add_broken_bond(i, j); } - num_bond[i] = 0; } // Update bond list and break solid-melted bonds @@ -530,7 +545,7 @@ void FixRHEOThermal::break_bonds() // Delete bonds for non-melted local atoms (shifting) if (i < nlocal) { for (m = 0; m < num_bond[i]; m++) { - if (bond_atom[i][m] == tag[j]) { + if (bond_atom[i][m] == tag[j] && bond_type[i][m] == btype) { nmax = num_bond[i] - 1; bond_type[i][m] = bond_type[i][nmax]; bond_atom[i][m] = bond_atom[i][nmax]; @@ -549,7 +564,7 @@ void FixRHEOThermal::break_bonds() if (j < nlocal) { for (m = 0; m < num_bond[j]; m++) { - if (bond_atom[j][m] == tag[i]) { + if (bond_atom[j][m] == tag[i] && bond_type[j][m] == btype) { nmax = num_bond[j] - 1; bond_type[j][m] = bond_type[j][nmax]; bond_atom[j][m] = bond_atom[j][nmax]; diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index 2fffa8b29c..abaa55ca70 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL), Thomas O'Connor (CMU) ----------------------------------------------------------------------- */ #include "fix_rheo_viscosity.h" diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index b07e914af1..fe5a1d6dec 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -13,7 +13,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) + Joel Clemmer (SNL), Thomas O'Connor (CMU) ----------------------------------------------------------------------- */ #include "pair_rheo.h" diff --git a/src/RHEO/pair_rheo_solid.cpp b/src/RHEO/pair_rheo_solid.cpp index d0a68d5230..f6dcd95879 100644 --- a/src/RHEO/pair_rheo_solid.cpp +++ b/src/RHEO/pair_rheo_solid.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL) +----------------------------------------------------------------------- */ + #include "pair_rheo_solid.h" #include "atom.h"