From 5dcbbc084ba0c42f8ed2be3f0f07ac6b9f575cb6 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 1 Nov 2022 13:30:24 -0600 Subject: [PATCH] Adding developers documentation page --- doc/src/Howto_granular.rst | 11 +++ doc/src/Modify_gsm.rst | 183 +++++++++++++++++++++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 doc/src/Modify_gsm.rst diff --git a/doc/src/Howto_granular.rst b/doc/src/Howto_granular.rst index 9492e5755e..af52904204 100644 --- a/doc/src/Howto_granular.rst +++ b/doc/src/Howto_granular.rst @@ -43,6 +43,14 @@ The fix style *freeze* zeroes both the force and torque of frozen atoms, and should be used for granular system instead of the fix style *setforce*\ . +To model heat conduction, one must use the atom style: +* :doc:`atom_style sphere/temp ` +a temperature integration fix +* :doc:`fix temp/integrate ` +and a heat conduction option defined in both +* :doc:`pair_style granular ` +* :doc:`fix wall/gran ` + For computational efficiency, you can eliminate needless pairwise computations between frozen atoms by using this command: @@ -55,3 +63,6 @@ computations between frozen atoms by using this command: will be the same as in 3d. If you wish to model granular particles in 2d as 2d discs, see the note on this topic on the :doc:`Howto 2d ` doc page, where 2d simulations are discussed. + +To add custom granular contact models, see the +:doc:`modifying granular submodels page `. diff --git a/doc/src/Modify_gsm.rst b/doc/src/Modify_gsm.rst new file mode 100644 index 0000000000..06172b4dec --- /dev/null +++ b/doc/src/Modify_gsm.rst @@ -0,0 +1,183 @@ +Granular Submodel (GSM) styles +=========== + +In granular models, particles are spheres with a finite radius and rotational +degrees of freedom as further described in the +:doc:`Howto granular page `. Interactions between pair of +particles or particles and walls may therefore depend on many different modes +of motion as described in :doc:`pair granular ` and +:doc:`fix wall/gran `. In both cases, the exchange of forces, +torques, and heat flow between two types of bodies is defined using a +GranularModel class. The GranularModel class organizes the details of an +interaction using a series of granular submodels each of which describe a +particular interaction mode (e.g. normal forces or rolling friction). From a +parent GSM class, several types of submodel classes are derived: + +* GSMNormal: normal force submodel +* GSMDamping: normal damping submodel +* GSMTangential: tangential forces and sliding friction submodel +* GSMRolling: rolling friction submodel +* GSMTwisting: twisting friction submodel +* GSMHeat: heat conduction submodel + +For each type of submodel, more classes are further derived, each describing +a specific implementation. For instance, from the GSMNormal class the +GSMNormalHooke, GSMNormalHertz, and GSMNormalJKR classes are derived which +calculate Hookean, Hertzian, or JKR normal forces, respectively. This modular +structure simplifies the addition of new granualar contact models as as one only +needs to create a new GSM class without having to modify the more complex +PairGranular, FixGranWall, and GranularModel classes. Most GSM methods are also +already defined by the parent classes so new contact models typically only require +edits to a few relevant methods (e.g. methods that define coefficients and +calculate forces). + +Each GSM class has a pointer to both the LAMMPS class and the GranularModel +class which owns it, ``lmp`` and ``gm``, respectively. The GranularModel class +includes several public variables that describe the geometry/dynamics of the +contact such as + +.. list-table:: + + * - ``xi`` and ``xj`` + - Positions of the two contacting bodies + * - ``vi`` and ``vj`` + - Velocities of the two contacting bodies + * - ``omegai`` and ``omegaj`` + - Angular velocities of the two contacting bodies + * - ``dx`` and ``nx`` + - The displacement and normalized displacement vectors + * - ``r``, ``rsq``, and ``rinv`` + - The distance, distance squared, and inverse distance + * - ``radsum`` + - The sum of particle radii + * - ``vr``, ``vn``, and ``vt`` + - The relative velocity vector and its normal and tangential components + * - ``wr`` + - The relative rotational velocity + +These quantities, among others, are calculated in the ``GranularModel->check_contact()`` +and ``GranularModel->calculate_forces()`` methods which can be referred to for more +details. + +To create a new GSM class, it is recommended that one first looks at similar GSM +classes. All GSM classes share several general methods which may need to be defined + +.. list-table:: + + * - ``GSM->mix_coeff()`` + - Optional method to define how coefficients are mixed for different atom types. By default, coefficients are mixed using a geometric mean. + * - ``GSM->coeffs_to_local()`` + - Parses coefficients to define local variables, run once at model construction. + * - ``GSM->init()`` + - Optional method to define local variables after other GSM types were created. For instance, this method may be used by a tangential model that derives parameters from the normal stiffness. + +There are also several type-specific methods + +.. list-table:: + + * - ``GSMNormal->touch()`` + - Optional method to test when particles are in contact. By default, this is when particles overlap. + * - ``GSMNormal->pulloff_distance()`` + - Optional method to return the distance at which particles stop interacting. By default, this is when particles no longer overlap. + * - ``GSMNormal->calculate_area()`` + - Optional method to return the surface area of the contact. By default, this returns the geometric cross section. + * - ``GSMNormal->set_fncrit()`` + - Optional method that defines the critical force to break the contact used by some tangential, rolling, and twisting submodels. By default, this is the current total normal force, including damping. + * - ``GSMNormal->calculate_forces()`` + - Required method that returns the normal contact force + * - ``GSMDamping->calculate_forces()`` + - Required method that returns the normal damping force + * - ``GSMTangential->calculate_forces()`` + - Required method that calculates tangential forces/torques + * - ``GSMTwisting->calculate_forces()`` + - Required method that calculates twisting friction forces/torques + * - ``GSMRolling->calculate_forces()`` + - Required method that calculates rolling friction forces/torques + * - ``GSMHeat->calculate_heat()`` + - Required method that returns the rate of heat flow + +As an example, say one wanted to create a new normal force option that consisted +of a Hookean force with a piecewise stiffness. This could be done by adding a new +set of files ``gsm_custom.h``: + +.. code-block:: c++ + + #ifdef GSM_CLASS + // clang-format off + GSMStyle(hooke/piecewise, + GSMNormalHookePiecewise, + NORMAL); + // clang-format on + #else + + #ifndef GSM_CUSTOM_H_ + #define GSM_CUSTOM_H_ + + #include "gsm.h" + #include "gsm_normal.h" + + namespace LAMMPS_NS { + namespace Granular_NS { + + class GSMNormalHookePiecewise : public GSMNormal { + public: + GSMNormalHookePiecewise(class GranularModel *, class LAMMPS *); + void coeffs_to_local() override; + void set_knfac(); + double calculate_forces(); + protected: + double k1, k2, delta_switch; + }; + + } // namespace Granular_NS + } // namespace LAMMPS_NS + + #endif /*GSM_CUSTOM_H_ */ + #endif /*GSM_CLASS_H_ */ + + +and ``gsm_custom.cpp`` + +.. code-block:: c++ + + #include "gsm_custom.h" + #include "gsm_normal.h" + #include "granular_model.h" + + using namespace LAMMPS_NS; + using namespace Granular_NS; + + GSMNormalHookePiecewise::GSMNormalHookePiecewise(GranularModel *gm, LAMMPS *lmp) : GSMNormal(gm, lmp) + { + num_coeffs = 4; + } + + /* ---------------------------------------------------------------------- */ + + void GSMNormalHookePiecewise::coeffs_to_local() + { + k1 = coeffs[0]; + k2 = coeffs[1]; + damp = coeffs[2]; + delta_switch = coeffs[3]; + } + + /* ---------------------------------------------------------------------- */ + + double GSMNormalHookePiecewise::calculate_forces() + { + if (gm->delta >= delta_switch { + Fne = k1 * delta_switch + k2 * (gm->delta - delta_switch); + } else { + Fne = k1 * gm->delta; + } + return Fne; + } + + /* ---------------------------------------------------------------------- */ + + void GSMNormalHookePiecewise::set_knfac() + { + if (gm->delta < delta_switch) knfac = k1; + else knfac = k2; + }