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

Merge pull request #3181 from jibril-b-coulibaly/new_damp

Browse Source 
Damping fixes for finite-size particles simulations
This commit is contained in:
Axel Kohlmeyer
2022-03-24 11:08:02 -04:00
committed by GitHub
parent d8d12ccf1c 07cec37786
commit 81d0165c5c
10 changed files with 849 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_fix.rst
View File

@ -55,6 +55,7 @@ OPT.
* :doc:`cmap <fix_cmap>`
* :doc:`colvars <fix_colvars>`
* :doc:`controller <fix_controller>`
* :doc:`damping/cundall <fix_damping_cundall>`
* :doc:`deform (k) <fix_deform>`
* :doc:`deposit <fix_deposit>`
* :doc:`dpd/energy (k) <fix_dpd_energy>`
@ -243,6 +244,7 @@ OPT.
* :doc:`vector <fix_vector>`
* :doc:`viscosity <fix_viscosity>`
* :doc:`viscous <fix_viscous>`
* :doc:`viscous/sphere <fix_viscous_sphere>`
* :doc:`wall/body/polygon <fix_wall_body_polygon>`
* :doc:`wall/body/polyhedron <fix_wall_body_polyhedron>`
* :doc:`wall/colloid <fix_wall>`

2
doc/src/fix.rst
View File

@ -198,6 +198,7 @@ accelerated styles exist.
* :doc:`cmap <fix_cmap>` - enables CMAP cross-terms of the CHARMM force field
* :doc:`colvars <fix_colvars>` - interface to the collective variables "Colvars" library
* :doc:`controller <fix_controller>` - apply control loop feedback mechanism
* :doc:`damping/cundall <fix_damping_cundall>` - Cundall non-viscous damping for granular simulations
* :doc:`deform <fix_deform>` - change the simulation box size/shape
* :doc:`deposit <fix_deposit>` - add new atoms above a surface
* :doc:`dpd/energy <fix_dpd_energy>` - constant energy dissipative particle dynamics
@ -386,6 +387,7 @@ accelerated styles exist.
* :doc:`vector <fix_vector>` - accumulate a global vector every N timesteps
* :doc:`viscosity <fix_viscosity>` - Muller-Plathe momentum exchange for viscosity calculation
* :doc:`viscous <fix_viscous>` - viscous damping for granular simulations
* :doc:`viscous/sphere <fix_viscous_sphere>` - viscous damping on angular velocity for granular simulations
* :doc:`wall/body/polygon <fix_wall_body_polygon>` -
* :doc:`wall/body/polyhedron <fix_wall_body_polyhedron>` -
* :doc:`wall/colloid <fix_wall>` - Lennard-Jones wall interacting with finite-size particles

149
doc/src/fix_damping_cundall.rst Normal file
View File

@ -0,0 +1,149 @@
.. index:: fix damping/cundall
fix damping/cundall command
===========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID damping/cundall gamma_l gamma_a keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* damping/cundall = style name of this fix command
* gamma_l = linear damping coefficient (dimensionless)
* gamma_a = angular damping coefficient (dimensionless)
* zero or more keyword/value pairs may be appended
.. parsed-literal::
keyword = *scale*
*scale* values = *type ratio* or *v_name*
type = atom type (1-N)
ratio = factor to scale the damping coefficients by
v_name = reference to atom style variable *name*
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all damping/cundall 0.8 0.8
fix 1 all damping/cundall 0.8 0.5 scale 3 2.5
fix a all damping/cundall 0.8 0.5 scale v_radscale
Description
"""""""""""
Add damping force and torque to finite-size spherical particles in the group
following the model of :ref:`Cundall, 1987 <Cundall1987>`, as implemented in other
granular physics code (e.g., :ref:`Yade-DEM <YadeDEM>`, :ref:`PFC <PFC>`).
The damping is constructed to always have negative mechanical power with respect
to the current velocity/angular velocity to ensure dissipation of kinetic energy.
If used without additional thermostatting (to add kinetic energy to the system),
it has the effect of slowly (or rapidly) freezing the system; hence it can also
be used as a simple energy minimization technique.
The magnitude of the damping force/torque :math:`F_d`/:math:`T_d` is a fraction
:math:`\gamma \in [0;1]` of the current force/torque :math:`F`/:math:`T` on the
particle. Damping is applied component-by-component in each direction
:math:`k\in\{x, y, z\}`:
.. math::
{F_d}_k = - \gamma_l \, F_k \, \mathrm{sign}(F_k v_k)
.. math::
{T_d}_k = - \gamma_a \, T_k \, \mathrm{sign}(T_k \omega_k)
The larger the coefficients, the faster the kinetic energy is reduced.
If the optional keyword *scale* is used, :math:`\gamma_l` and :math:`\gamma_a`
can be scaled up or down by the specified factor for atoms. This factor can be
set for different atom types and thus the *scale* keyword used multiple times
followed by the atom type and the associated scale factor. Alternately the
scaling factor can be computed for each atom (e.g. based on its radius) by
using an :doc:`atom-style variable <variable>`.
.. Note::
The damping force/torque is computed based on the force/torque at the moment
this fix is invoked. Any force/torque added after this fix, e.g., by
:doc:`fix addforce <fix_addforce>` or :doc:`fix addtorque <fix_addtorque>`
will not be damped. When performing simulations with gravity, invoking
:doc:`fix gravity <fix_gravity>` after this fix will maintain the specified
gravitational acceleration.
.. Note::
This scheme is dependent on the coordinates system and does not correspond to
realistic physical processes. It is constructed for numerical convenience and
efficacy.
This non-viscous damping presents the following advantages:
1. damping is independent of velocity, equally damping regions with distinct natural frequencies,
2. damping affects acceleration and vanishes for steady uniform motion of the particles,
3. damping parameter :math:`\gamma` is dimensionless and does not require scaling.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix. No global or per-atom quantities are stored by
this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.
The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
integrator the fix is modifying forces/torques. Default is the outermost level.
The forces/torques due to this fix are imposed during an energy minimization,
invoked by the :doc:`minimize <minimize>` command. This fix should only
be used with damped dynamics minimizers that allow for
non-conservative forces. See the :doc:`min_style <min_style>` command
for details.
Restrictions
""""""""""""
This fix is part of the GRANULAR package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
This fix requires that atoms store torque and a radius as defined by the
:doc:`atom_style sphere <atom_style>` command.
Related commands
""""""""""""""""
:doc:`fix viscous <fix_viscous>`, :doc:`fix viscous/sphere <fix_viscous_sphere>`
Default
"""""""
none
References
""""""""""
.. _Cundall1987:
**(Cundall, 1987)** Cundall, P. A. Distinct Element Models of Rock and Soil
Structure, in Analytical and Computational Methods in Engineering Rock
Mechanics, Ch. 4, pp. 129-163. E. T. Brown, ed. London: Allen & Unwin., 1987.
.. _PFC:
**(PFC)** PFC Particle Flow Code 6.0 Documentation. Itasca Consulting Group.
.. _YadeDEM:
**(Yade-DEM)** V. Smilauer et al. (2021), Yade Documentation 3rd ed.
The Yade Project. DOI:10.5281/zenodo.5705394 (https://yade-dem.org/doc/)

7
doc/src/fix_viscous.rst
View File

@ -106,12 +106,15 @@ for details.
Restrictions
""""""""""""
none
none
Related commands
""""""""""""""""
:doc:`fix langevin <fix_langevin>`
:doc:`fix langevin <fix_langevin>`,
:doc:`fix viscous/sphere <fix_viscous_sphere>`,
:doc:`fix damping/cundall <fix_damping_cundall>`
Default
"""""""

111
doc/src/fix_viscous_sphere.rst Normal file
View File

@ -0,0 +1,111 @@
.. index:: fix viscous/sphere
fix viscous/sphere command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID viscous/sphere gamma keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* viscous/sphere = style name of this fix command
* gamma = damping coefficient (torque/angular velocity units)
* zero or more keyword/value pairs may be appended
.. parsed-literal::
keyword = *scale*
*scale* values = *type ratio* or *v_name*
type = atom type (1-N)
ratio = factor to scale the damping coefficients by
v_name = reference to atom style variable *name*
Examples
""""""""
.. code-block:: LAMMPS
fix 1 flow viscous/sphere 0.1
fix 1 damp viscous/sphere 0.5 scale 3 2.5
fix 1 damp viscous/sphere 0.5 scale v_radscale
Description
"""""""""""
Add a viscous damping torque to finite-size spherical particles in the group
that is proportional to the angular velocity of the atom. In granular
simulations this can be useful for draining the rotational kinetic energy from
the system in a controlled fashion. If used without additional thermostatting
(to add kinetic energy to the system), it has the effect of slowly (or rapidly)
freezing the system; hence it can also be used as a simple energy minimization
technique.
The damping torque :math:`T_i` is given by :math:`T_i = - \gamma \omega_i`.
The larger the coefficient, the faster the rotational kinetic energy is reduced.
If the optional keyword *scale* is used, :math:`\gamma` can be scaled up
or down by the specified factor for atoms. This factor can be set for
different atom types and thus the *scale* keyword used multiple times
followed by the atom type and the associated scale factor. Alternately
the scaling factor can be computed for each atom (e.g. based on its
radius) by using an :doc:`atom-style variable <variable>`.
.. note::
You should specify gamma in torque/angular velocity units. This is not
the same as mass/time units, at least for some of the LAMMPS
:doc:`units <units>` options like "real" or "metal" that are not
self-consistent.
In the current implementation, rather than have the user specify a viscosity,
:math:`\gamma` is specified directly in torque/angular velocity units.
If needed, :math:`\gamma` can be adjusted for atoms of different sizes
(i.e. :math:`\sigma`) by using the *scale* keyword.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix. No global or per-atom quantities are stored by
this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.
The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
integrator the fix is modifying torques. Default is the outermost level.
The torques due to this fix are imposed during an energy minimization,
invoked by the :doc:`minimize <minimize>` command. This fix should only
be used with damped dynamics minimizers that allow for
non-conservative forces. See the :doc:`min_style <min_style>` command
for details.
Restrictions
""""""""""""
This fix is part of the EXTRA-FIX package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
This fix requires that atoms store torque and angular velocity (omega)
and a radius as defined by the :doc:`atom_style sphere <atom_style>`
command.
All particles in the group must be finite-size spheres. They cannot
be point particles.
Related commands
""""""""""""""""
:doc:`fix viscous <fix_viscous>`, :doc:`fix damping/cundall <fix_damping_cundall>`
Default
"""""""
none

8
doc/utils/sphinx-config/false_positives.txt
View File

@ -573,6 +573,8 @@ cuFFT
CuH
Cui
Cummins
Cundall
cundall
Curk
Cusentino
customIDs
@ -655,6 +657,7 @@ delocalized
Delong
delr
deltaHf
dem
Dendrimer
dendritic
Denniston
@ -3096,6 +3099,7 @@ smallsmall
smd
SMD
smi
Smilauer
Smirichinski
Smit
smtbq
@ -3481,6 +3485,7 @@ unsplit
unstrained
untar
untilted
Unwin
uparrow
upenn
upto
@ -3701,6 +3706,8 @@ xy
xyz
xz
xzhou
Yade
yade
yaff
YAFF
Yamada
@ -3746,6 +3753,7 @@ zcm
zeeman
Zeeman
Zemer
zenodo
Zepeda
zflag
Zhang

210
src/EXTRA-FIX/fix_viscous_sphere.cpp Normal file
View File

@ -0,0 +1,210 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_viscous_sphere.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "input.h"
#include "memory.h"
#include "modify.h"
#include "respa.h"
#include "update.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
// type of scaling
enum { NONE, TYPE, VARIABLE };
/* ---------------------------------------------------------------------- */
FixViscousSphere::FixViscousSphere(LAMMPS *_lmp, int narg, char **arg) :
Fix(_lmp, narg, arg), scalegamma(nullptr), scaleval(nullptr), scalevarid(nullptr),
scalestyle(NONE)
{
dynamic_group_allow = 1;
if (!atom->sphere_flag) error->all(FLERR, "Fix viscous/sphere requires atom style sphere");
if (narg < 4) error->all(FLERR, "Illegal fix viscous/sphere command");
gamma = utils::numeric(FLERR, arg[3], false, lmp);
// optional args
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg], "scale") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix viscous/sphere command");
if (utils::strmatch(arg[iarg + 1], "^v_")) {
if (scalestyle != NONE) error->all(FLERR, "Must use only one style of scaling");
scalevarid = utils::strdup(arg[iarg + 1] + 2);
int ivariable = input->variable->find(scalevarid);
if (ivariable < 0)
error->all(FLERR, "Variable name {} for fix viscous/sphere does not exist", scalevarid);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR, "Fix viscous/scale variable {} is not atom-style variable", scalevarid);
scalestyle = VARIABLE;
memory->destroy(scaleval);
memory->create(scaleval, atom->nmax, "fix_viscous/sphere:scaleval");
iarg += 2;
} else {
if (scalestyle == VARIABLE) error->all(FLERR, "Must use only one style of scaling");
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix viscous/sphere command");
if (!scalegamma) {
scalegamma = new double[atom->ntypes + 1];
for (int i = 1; i <= atom->ntypes; i++) scalegamma[i] = 1.0;
}
int itype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
double scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if ((itype < 1) || (itype > atom->ntypes))
error->all(FLERR, "Atom type {} out of range for fix viscous/sphere command:", itype);
scalegamma[itype] = scale;
scalestyle = TYPE;
iarg += 3;
}
} else
error->all(FLERR, "Illegal fix viscous/sphere command");
}
respa_level_support = 1;
ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
FixViscousSphere::~FixViscousSphere()
{
memory->destroy(scaleval);
delete[] scalegamma;
delete[] scalevarid;
}
/* ---------------------------------------------------------------------- */
int FixViscousSphere::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::init()
{
int max_respa = 0;
if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = max_respa = ((Respa *) update->integrate)->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level, max_respa);
}
if (scalestyle == VARIABLE) {
int ivariable = input->variable->find(scalevarid);
if (ivariable < 0)
error->all(FLERR, "Variable name {} for fix viscous/sphere does not exist", scalevarid);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR, "Fix viscous/sphere variable {} is not atom-style variable", scalevarid);
scalevar = ivariable;
}
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::setup(int vflag)
{
if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag, ilevel_respa, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::post_force(int /*vflag*/)
{
// apply drag torque to finite-size atoms in group
// direction is opposed to angular velocity vector
// magnitude depends on atom type
double **omega = atom->omega;
double **torque = atom->torque;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
double drag;
if (scalestyle == VARIABLE) {
memory->grow(scaleval, atom->nmax, "fix_viscous/sphere:scaleval");
input->variable->compute_atom(scalevar, igroup, scaleval, 1, 0);
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (scalestyle == TYPE) {
drag = gamma * scalegamma[type[i]];
} else if (scalestyle == VARIABLE) {
drag = gamma * scaleval[i];
} else { // scalestyle == NONE
drag = gamma;
}
torque[i][0] -= drag * omega[i][0];
torque[i][1] -= drag * omega[i][1];
torque[i][2] -= drag * omega[i][2];
}
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixViscousSphere::min_post_force(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
double FixViscousSphere::memory_usage()
{
if (scalestyle == VARIABLE)
return (double) sizeof(double) * atom->nmax;
else if (scalestyle == TYPE)
return (double) sizeof(double) * atom->ntypes;
else
return 0.0;
}

60
src/EXTRA-FIX/fix_viscous_sphere.h Normal file
View File

@ -0,0 +1,60 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(viscous/sphere,FixViscousSphere);
// clang-format on
#else
#ifndef LMP_FIX_VISCOUS_SPHERE_H
#define LMP_FIX_VISCOUS_SPHERE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixViscousSphere : public Fix {
public:
FixViscousSphere(class LAMMPS *, int, char **);
virtual ~FixViscousSphere();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double memory_usage();
protected:
double gamma, *scalegamma, *scaleval;
const char *scalevarid;
int scalestyle, scalevar;
int ilevel_respa;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

242
src/GRANULAR/fix_damping_cundall.cpp Normal file
View File

@ -0,0 +1,242 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_damping_cundall.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "input.h"
#include "memory.h"
#include "modify.h"
#include "respa.h"
#include "update.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
// type of scaling
enum { NONE, TYPE, VARIABLE };
/* ---------------------------------------------------------------------- */
FixDampingCundall::FixDampingCundall(LAMMPS *_lmp, int narg, char **arg) :
Fix(_lmp, narg, arg), scalegamma(nullptr), scaleval(nullptr), scalevarid(nullptr),
scalestyle(NONE)
{
dynamic_group_allow = 1;
if (!atom->sphere_flag) error->all(FLERR, "Fix damping/cundall requires atom style sphere");
if (narg < 5) error->all(FLERR, "Illegal fix damping/cundall command");
gamma_lin = utils::numeric(FLERR, arg[3], false, lmp);
gamma_ang = utils::numeric(FLERR, arg[4], false, lmp);
// optional args
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg], "scale") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix damping/cundall command");
if (utils::strmatch(arg[iarg + 1], "^v_")) {
if (scalestyle != NONE) error->all(FLERR, "Must use only one style of scaling");
scalevarid = utils::strdup(arg[iarg + 1] + 2);
int ivariable = input->variable->find(scalevarid);
if (ivariable < 0)
error->all(FLERR, "Variable name {} for fix damping/cundall does not exist", scalevarid);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR, "Fix viscous/scale variable {} is not atom-style variable", scalevarid);
scalestyle = VARIABLE;
memory->destroy(scaleval);
memory->create(scaleval, atom->nmax, "fix_damping/cundall:scaleval");
iarg += 2;
} else {
if (scalestyle == VARIABLE) error->all(FLERR, "Must use only one style of scaling");
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix damping/cundall command");
if (!scalegamma) {
scalegamma = new double[atom->ntypes + 1];
for (int i = 1; i <= atom->ntypes; i++) scalegamma[i] = 1.0;
}
int itype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
double scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if ((itype < 1) || (itype > atom->ntypes))
error->all(FLERR, "Atom type {} out of range for fix damping/cundall command:", itype);
scalegamma[itype] = scale;
scalestyle = TYPE;
iarg += 3;
}
} else
error->all(FLERR, "Illegal fix damping/cundall command");
}
respa_level_support = 1;
ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
FixDampingCundall::~FixDampingCundall()
{
memory->destroy(scaleval);
delete[] scalegamma;
delete[] scalevarid;
}
/* ---------------------------------------------------------------------- */
int FixDampingCundall::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::init()
{
int max_respa = 0;
if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = max_respa = ((Respa *) update->integrate)->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level, max_respa);
}
bool fflag = false;
for (auto ifix : modify->get_fix_list()) {
if (fflag && (comm->me == 0) && (ifix->setmask() & POST_FORCE))
error->warning(FLERR, "Fix {} alters forces after fix damping/cundall", ifix->id);
if (ifix == this) fflag = true;
}
if (scalestyle == VARIABLE) {
int ivariable = input->variable->find(scalevarid);
if (ivariable < 0)
error->all(FLERR, "Variable name {} for fix damping/cundall does not exist", scalevarid);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR, "Fix damping/cundall variable {} is not atom-style variable", scalevarid);
scalevar = ivariable;
}
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::setup(int vflag)
{
if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag, ilevel_respa, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::post_force(int /*vflag*/)
{
// apply damping force/torque to finite-size atoms in group
// add a fraction of the current force/torque if work is negative
// subtract a fraction of the current force/torque if work is positive
// applied over each component independently (non-objective)
// magnitude depends on atom type
// see, e.g. Yade-DEM implementation of NewtonIntegrator::cundallDamp1st()
// gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/NewtonIntegrator.cpp
double **v = atom->v;
double **omega = atom->omega;
double **f = atom->f;
double **torque = atom->torque;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
double gamma_l, gamma_a;
int signf0, signf1, signf2;
int signt0, signt1, signt2;
if (scalestyle == VARIABLE) {
memory->grow(scaleval, atom->nmax, "fix_damping/cundall:scaleval");
input->variable->compute_atom(scalevar, igroup, scaleval, 1, 0);
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (scalestyle == TYPE) {
gamma_l = gamma_lin * scalegamma[type[i]];
gamma_a = gamma_ang * scalegamma[type[i]];
} else if (scalestyle == VARIABLE) {
gamma_l = gamma_lin * scaleval[i];
gamma_a = gamma_ang * scaleval[i];
} else { // scalestyle NONE
gamma_l = gamma_lin;
gamma_a = gamma_ang;
}
signf0 = (f[i][0] * v[i][0] >= 0.0) ? 1.0 : -1.0;
signf1 = (f[i][1] * v[i][1] >= 0.0) ? 1.0 : -1.0;
signf2 = (f[i][2] * v[i][2] >= 0.0) ? 1.0 : -1.0;
f[i][0] *= 1.0 - gamma_l * signf0;
f[i][1] *= 1.0 - gamma_l * signf1;
f[i][2] *= 1.0 - gamma_l * signf2;
signt0 = (torque[i][0] * omega[i][0] >= 0.0) ? 1.0 : -1.0;
signt1 = (torque[i][1] * omega[i][1] >= 0.0) ? 1.0 : -1.0;
signt2 = (torque[i][2] * omega[i][2] >= 0.0) ? 1.0 : -1.0;
torque[i][0] *= 1.0 - gamma_a * signt0;
torque[i][1] *= 1.0 - gamma_a * signt1;
torque[i][2] *= 1.0 - gamma_a * signt2;
}
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixDampingCundall::min_post_force(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
double FixDampingCundall::memory_usage()
{
if (scalestyle == VARIABLE)
return (double) sizeof(double) * atom->nmax;
else if (scalestyle == TYPE)
return 2.0 * sizeof(double) * atom->ntypes;
else
return 0.0;
}

60
src/GRANULAR/fix_damping_cundall.h Normal file
View File

@ -0,0 +1,60 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(damping/cundall,FixDampingCundall);
// clang-format on
#else
#ifndef LMP_FIX_DAMPING_CUNDALL_H
#define LMP_FIX_DAMPING_CUNDALL_H
#include "fix.h"
namespace LAMMPS_NS {
class FixDampingCundall : public Fix {
public:
FixDampingCundall(class LAMMPS *, int, char **);
virtual ~FixDampingCundall();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double memory_usage();
protected:
double gamma_lin, gamma_ang, *scalegamma, *scaleval;
const char *scalevarid;
int scalestyle, scalevar;
int ilevel_respa;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/