Adding custom svector gransubmod quantities

doc/src/fix_wall_gran.rst

@@ -222,10 +222,10 @@ restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
 
 If the :code:`contacts` option is used, this fix generates a per-atom array
-with 8 columns as output, containing the contact information for owned
+with at least 8 columns as output, containing the contact information for owned
 particles (nlocal on each processor). All columns in this per-atom array will
-be zero if no contact has occurred.  The values of these columns are listed in
-the following table:
+be zero if no contact has occurred.  The first 8 values of these columns are
+listed in the following table.
 
 +-------+----------------------------------------------------+----------------+
 | Index | Value                                              | Units          |
@@ -248,6 +248,14 @@ the following table:
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
+If a granular submodel calculates additional contact information (e.g. the
+heat submodels calculate the amount of heat exchanged), these quantities
+are appended to the end of this array. First, any extra values from the
+normal submodel are appended followed by the damping, tangential, rolling,
+twisting, then heat models. See the descriptions of granular submodels in
+the :doc:`pair granular <pair_granular>` page for information on any extra
+quantities.
+
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
 No parameter of this fix can be used with the *start/stop* keywords of the
 :doc:`run <run>` command. This fix is not invoked during :doc:`energy

doc/src/fix_wall_gran_region.rst

@@ -243,10 +243,10 @@ uninterrupted fashion.
    with a different region ID.
 
 If the :code:`contacts` option is used, this fix generates a per-atom array
-with 8 columns as output, containing the contact information for owned
+with at least 8 columns as output, containing the contact information for owned
 particles (nlocal on each processor). All columns in this per-atom array will
-be zero if no contact has occurred. The values of these columns are listed in
-the following table:
+be zero if no contact has occurred.  The first 8 values of these columns are
+listed in the following table.
 
 +-------+----------------------------------------------------+----------------+
 | Index | Value                                              | Units          |
@@ -269,6 +269,14 @@ the following table:
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
+If a granular submodel calculates additional contact information (e.g. the
+heat submodels calculate the amount of heat exchanged), these quantities
+are appended to the end of this array. First, any extra values from the
+normal submodel are appended followed by the damping, tangential, rolling,
+twisting, then heat models. See the descriptions of granular submodels in
+the :doc:`pair granular <pair_granular>` page for information on any extra
+quantities.
+
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
 No parameter of this fix can be used with the *start/stop* keywords of the
 :doc:`run <run>` command. This fix is not invoked during :doc:`energy

doc/src/pair_granular.rst

@@ -166,15 +166,15 @@ initially will not experience force until they come into contact
 experience a tensile force up to :math:`3\pi\gamma R`, at which point they
 lose contact.
 
-The *mdr* model is a mechanically-derived contact model designed to capture the 
+The *mdr* model is a mechanically-derived contact model designed to capture the
 contact response between adhesive elastic-plastic particles into large deformation.
 The theoretical foundations of the *mdr* model are detailed in the
-two-part series :ref:`Zunker and Kamrin Part I <Zunker2024I>` and 
-:ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development 
+two-part series :ref:`Zunker and Kamrin Part I <Zunker2024I>` and
+:ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development
 and demonstrations of its application to industrially relevant powder
 compaction processes are presented in :ref:`Zunker et al. <Zunker2025>`.
 
-The model requires the following inputs: 
+The model requires the following inputs:
 
    1. *Young's modulus* :math:`E > 0` : The Young's modulus is commonly reported
    for various powders.
@@ -197,12 +197,12 @@ The model requires the following inputs:
    response.
 
    6. *Coefficient of restitution* :math:`0 \le e \le 1` : The coefficient of
-   restitution is a tunable parameter that controls damping in the normal direction. 
+   restitution is a tunable parameter that controls damping in the normal direction.
 
 .. note::
 
    The values for :math:`E`, :math:`\nu`, :math:`Y`, and :math:`\Delta\gamma` (i.e.,
-   :math:`K_{Ic}`) should be selected for zero porosity to reflect the intrinsic 
+   :math:`K_{Ic}`) should be selected for zero porosity to reflect the intrinsic
    material property rather than the bulk powder property.
 
 The *mdr* model produces a nonlinear force-displacement response, therefore the
@@ -218,7 +218,7 @@ Here, :math:`\kappa = E/(3(1-2\nu))` is the bulk modulus and
 .. note::
 
    The *mdr* model requires some specific settings to function properly,
-   please read the following text carefully to ensure all requirements are 
+   please read the following text carefully to ensure all requirements are
    followed.
 
 The *atom_style* must be set to *sphere 1* to enable dynamic particle
@@ -230,7 +230,7 @@ deformation is accumulated, a new enlarged apparent radius is defined to
 ensure that that volume change due to plastic deformation is not lost.
 This apparent radius is stored as the *atom radius* meaning it is used
 for subsequent neighbor list builds and contact detection checks. The
-advantage of this is that multi-neighbor dependent effects such as 
+advantage of this is that multi-neighbor dependent effects such as
 formation of secondary contacts caused by radial expansion are captured
 by the *mdr* model. Setting *atom_style sphere 1* ensures that updates to
 the particle radii are properly reflected throughout the simulation.
@@ -243,14 +243,14 @@ Newton's third law must be set to *off*. This ensures that the neighbor lists
 are constructed properly for the topological penalty algorithm used to screen
 for non-physical contacts occurring through obstructing particles, an issue
 prevalent under large deformation conditions. For more information on this
-algorithm see :ref:`Zunker et al. <Zunker2025>`. 
+algorithm see :ref:`Zunker et al. <Zunker2025>`.
 
 .. code-block:: LAMMPS
 
    newton off
 
-The damping model must be set to *none*. The *mdr* model already has a built 
-in damping model. 
+The damping model must be set to *none*. The *mdr* model already has a built
+in damping model.
 
 .. code-block:: LAMMPS
 
@@ -268,7 +268,7 @@ Additionally, the following *mdr* inputs must match between the
 *pair_style* and *fix wall/gran/region* definitions: :math:`E`,
 :math:`\nu`, :math:`Y`, :math:`\psi_b`, and :math:`e`. The exception
 is :math:`\Delta\gamma`, which may vary, permitting different
-adhesive behaviors between particle-particle and particle-wall interactions.  
+adhesive behaviors between particle-particle and particle-wall interactions.
 
 .. note::
 
@@ -297,7 +297,7 @@ contact. In the input script, these quantities are initialized by calling
    *examples/granular* directory. The first is a die compaction
    simulation involving 200 particles named *in.tableting.200*.
    The second is a triaxial compaction simulation involving 12
-   particles named *in.triaxial.compaction.12*. 
+   particles named *in.triaxial.compaction.12*.
 ----------
 
 In addition, the normal force is augmented by a damping term of the
@@ -806,7 +806,10 @@ supported are:
 2. *radius* : :math:`k_{s}`
 3. *area* : :math:`h_{s}`
 
-If the *heat* keyword is not specified, the model defaults to *none*.
+If the *heat* keyword is not specified, the model defaults to *none*. All
+heat models calculate an additional pairwise quantity accessible by the
+single() function (described below) which is the heat conducted between the
+two particles.
 
 For *heat* *radius*, the heat
 :math:`Q` conducted between two particles is given by
@@ -921,7 +924,7 @@ The single() function of these pair styles returns 0.0 for the energy
 of a pairwise interaction, since energy is not conserved in these
 dissipative potentials.  It also returns only the normal component of
 the pairwise interaction force.  However, the single() function also
-calculates 13 extra pairwise quantities.  The first 3 are the
+calculates at least 13 extra pairwise quantities.  The first 3 are the
 components of the tangential force between particles I and J, acting
 on particle I.  The fourth is the magnitude of this tangential force.
 The next 3 (5-7) are the components of the rolling torque acting on
@@ -929,9 +932,13 @@ particle I. The next entry (8) is the magnitude of the rolling torque.
 The next entry (9) is the magnitude of the twisting torque acting
 about the vector connecting the two particle centers.
 The next 3 (10-12) are the components of the vector connecting
-the centers of the two particles (x_I - x_J). The last quantity (13)
-is the heat flow between the two particles, set to 0 if no heat model
-is active.
+the centers of the two particles (x_I - x_J). If a granular submodel
+calculates additional contact information (e.g. the heat submodels
+calculate the amount of heat exchanged), these quantities are appended
+to the end of this list. First, any extra values from the normal submodel
+are appended followed by the damping, tangential, rolling, twisting, then
+heat models. See the descriptions of specific granular submodels above
+for information on any extra quantities.
 
 These extra quantities can be accessed by the :doc:`compute pair/local <compute_pair_local>` command, as *p1*, *p2*, ...,
 *p12*\ .

src/GRANULAR/fix_wall_gran.cpp

@@ -199,7 +199,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
       iarg += 3;
     } else if (strcmp(arg[iarg],"contacts") == 0) {
       peratom_flag = 1;
-      size_peratom_cols = 8;
+      size_peratom_cols = 8 + model->nsvector;
       peratom_freq = 1;
       iarg += 1;
     } else if (strcmp(arg[iarg],"temperature") == 0) {
@@ -363,7 +363,7 @@ void FixWallGran::setup(int vflag)
 
 void FixWallGran::post_force(int /*vflag*/)
 {
-  int i,j;
+  int i,j,n;
   double dx,dy,dz,del1,del2,delxy,delr,rwall,meff;
   double *forces, *torquesi;
   double vwall[3];
@@ -523,7 +523,9 @@ void FixWallGran::post_force(int /*vflag*/)
     if (use_history) model->history = history_one[i];
     if (heat_flag) model->Ti = temperature[i];
 
+    if (peratom_flag) model->calculate_svector = 1;
     model->calculate_forces();
+    if (peratom_flag) model->calculate_svector = 0;
 
     forces = model->forces;
     torquesi = model->torquesi;
@@ -544,6 +546,9 @@ void FixWallGran::post_force(int /*vflag*/)
       array_atom[i][5] = x[i][1] - dy;
       array_atom[i][6] = x[i][2] - dz;
       array_atom[i][7] = radius[i];
+
+      for (n = 0; n < model->nsvector; n++)
+        array_atom[i][8 + n] = model->svector[n];
     }
   }
 }

src/GRANULAR/fix_wall_gran.h

@@ -84,7 +84,7 @@ class FixWallGran : public Fix {
 
   // store particle interactions
 
-  int store;
+  int nsvector;
 
   void clear_stored_contacts();
 };

src/GRANULAR/fix_wall_gran_region.cpp

@@ -118,7 +118,7 @@ void FixWallGranRegion::init()
 
 void FixWallGranRegion::post_force(int /*vflag*/)
 {
-  int i, m, nc, iwall;
+  int i, n, m, nc, iwall;
   double *forces, *torquesi;
   double meff, vwall[3], w0[3] = {0.0, 0.0, 0.0};
   bool touchflag = false;
@@ -262,7 +262,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
       if (use_history) model->history = history_many[i][c2r[ic]];
       if (heat_flag) model->Ti = temperature[i];
 
+      if (peratom_flag) model->calculate_svector = 1;
       model->calculate_forces();
+      if (peratom_flag) model->calculate_svector = 0;
 
       forces = model->forces;
       torquesi = model->torquesi;
@@ -283,6 +285,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
         array_atom[i][5] = x[i][1] - model->dx[1];
         array_atom[i][6] = x[i][2] - model->dx[2];
         array_atom[i][7] = radius[i];
+
+        for (n = 0; n < model->nsvector; n++)
+          array_atom[i][8 + n] = model->svector[n];
       }
     }
   }

src/GRANULAR/gran_sub_mod.cpp

@@ -42,6 +42,7 @@ GranSubMod::GranSubMod(class GranularModel *gm, LAMMPS *lmp) : Pointers(lmp)
   beyond_contact = 0;
   num_coeffs = 0;
   contact_radius_flag = 0;
+  nsvector = 0;
 
   nondefault_history_transfer = 0;
   transfer_history_factor = nullptr;

src/GRANULAR/gran_sub_mod.h

@@ -46,6 +46,8 @@ namespace Granular_NS {
 
     GranularModel *gm;
 
+    int nsvector, index_svector;
+
    protected:
     int allocated;
 

src/GRANULAR/gran_sub_mod_heat.cpp

@@ -50,6 +50,7 @@ GranSubModHeatRadius::GranSubModHeatRadius(GranularModel *gm, LAMMPS *lmp) : Gra
   num_coeffs = 1;
   contact_radius_flag = 1;
   conductivity = 0.0;
+  nsvector = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -65,7 +66,9 @@ void GranSubModHeatRadius::coeffs_to_local()
 
 double GranSubModHeatRadius::calculate_heat()
 {
-  return 2 * conductivity * gm->contact_radius * (gm->Tj - gm->Ti);
+  double heat = 2 * conductivity * gm->contact_radius * (gm->Tj - gm->Ti);
+  if (gm->calculate_svector) gm->svector[index_svector] = heat;
+  return heat;
 }
 
 
@@ -78,6 +81,7 @@ GranSubModHeatArea::GranSubModHeatArea(GranularModel *gm, LAMMPS *lmp) : GranSub
   num_coeffs = 1;
   contact_radius_flag = 1;
   heat_transfer_coeff = 0.0;
+  nsvector = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -93,5 +97,7 @@ void GranSubModHeatArea::coeffs_to_local()
 
 double GranSubModHeatArea::calculate_heat()
 {
-  return heat_transfer_coeff * MY_PI * gm->contact_radius * gm->contact_radius * (gm->Tj - gm->Ti);
+  double heat = heat_transfer_coeff * MY_PI * gm->contact_radius * gm->contact_radius * (gm->Tj - gm->Ti);
+  if (gm->calculate_svector) gm->svector[index_svector] = heat;
+  return heat;
 }

src/GRANULAR/granular_model.cpp

@@ -27,6 +27,7 @@
 #include "force.h"
 #include "gran_sub_mod.h"
 #include "math_extra.h"
+#include "memory.h"
 
 #include "style_gran_sub_mod.h"    // IWYU pragma: keep
 
@@ -64,6 +65,10 @@ GranularModel::GranularModel(LAMMPS *lmp) : Pointers(lmp)
   twisting_model = nullptr;
   heat_model = nullptr;
 
+  calculate_svector = 0;
+  nsvector = 0;
+  svector = nullptr;
+
   for (int i = 0; i < NSUBMODELS; i++) sub_models[i] = nullptr;
   transfer_history_factor = nullptr;
 
@@ -100,6 +105,7 @@ GranularModel::~GranularModel()
   delete[] gran_sub_mod_class;
   delete[] gran_sub_mod_names;
   delete[] gran_sub_mod_types;
+  delete[] svector;
 
   for (int i = 0; i < NSUBMODELS; i++) delete sub_models[i];
 }
@@ -298,6 +304,18 @@ void GranularModel::init()
   }
 
   for (int i = 0; i < NSUBMODELS; i++) sub_models[i]->init();
+
+  int index_svector = 0;
+  for (int i = 0; i < NSUBMODELS; i++) {
+    if (sub_models[i]->nsvector != 0) {
+      sub_models[i]->index_svector = index_svector;
+      nsvector += sub_models[i]->nsvector;
+      index_svector += sub_models[i]->nsvector;
+    }
+  }
+
+  if (nsvector != 0)
+    svector = new double[nsvector];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -498,9 +516,8 @@ void GranularModel::calculate_forces()
     if (contact_type == PAIR) sub3(torquesj, tortwist, torquesj);
   }
 
-  if (heat_defined) {
+  if (heat_defined)
     dq = heat_model->calculate_heat();
-  }
 }
 
 /* ----------------------------------------------------------------------

src/GRANULAR/granular_model.h

@@ -96,6 +96,10 @@ class GranularModel : protected Pointers {
   double magtwist;
   bool touch;
 
+  // Extra output
+  int calculate_svector, nsvector;
+  double *svector;
+
  protected:
   int rolling_defined, twisting_defined, heat_defined; // Flag optional sub models
   int classic_model;                                   // Flag original pair/gran calculations

src/GRANULAR/pair_granular.cpp

@@ -417,8 +417,10 @@ void PairGranular::init_style()
       error->all(FLERR,"Heat conduction in pair granular requires atom style with heatflow property");
   }
 
-  // allocate history and initialize models
+  // allocate history and aggregate model information
   class GranularModel* model;
+  double nsvector_total;
+  extra_svector = 0;
   int size_max[NSUBMODELS] = {0};
   for (int n = 0; n < nmodels; n++) {
     model = models_list[n];
@@ -429,13 +431,23 @@ void PairGranular::init_style()
     }
     if (model->size_history != 0) use_history = 1;
 
-    for (int i = 0; i < NSUBMODELS; i++)
+    nsvector_total = 0;
+    for (int i = 0; i < NSUBMODELS; i++) {
+      nsvector_total += model->sub_models[i]->nsvector;
       if (model->sub_models[i]->size_history > size_max[i])
         size_max[i] = model->sub_models[i]->size_history;
+    }
+    extra_svector = MAX(extra_svector, nsvector_total);
 
     if (model->nondefault_history_transfer) nondefault_history_transfer = 1;
   }
 
+  if (extra_svector != 0) {
+    single_extra = 12 + extra_svector;
+    delete[] svector;
+    svector = new double[single_extra];
+  }
+
   size_history = 0;
   if (use_history) {
     for (int i = 0; i < NSUBMODELS; i++) size_history += size_max[i];
@@ -770,7 +782,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   model->omegaj = omega[j];
   model->history = history;
 
+  model->calculate_svector = 1;
   model->calculate_forces();
+  model->calculate_svector = 0;
 
   // apply forces & torques
   // Calculate normal component, normalized by r
@@ -790,6 +804,14 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   svector[10] = model->dx[1];
   svector[11] = model->dx[2];
 
+  // add submodel-specific quantities
+  for (int n = 0; n < model->nsvector; n++)
+    svector[12 + n] = model->svector[n];
+
+  // zero any values unused by this specific model
+  for (int n = 12 + model->nsvector; n < single_extra; n++)
+    svector[n] = 0.0;
+
   return 0.0;
 }
 

src/GRANULAR/pair_granular.h

@@ -82,6 +82,8 @@ class PairGranular : public Pair {
   // optional user-specified global cutoff, per-type user-specified cutoffs
   double **cutoff_type;
   double cutoff_global;
+
+  int extra_svector;
 };
 
 }    // namespace LAMMPS_NS