added MDR damping method

src/GRANULAR/fix_granular_mdr.h

@@ -51,7 +51,8 @@ namespace Granular_MDR_NS {
     PENALTY,           // contact penalty
     DELTA_MAX,
     DELTAP_0,
-    DELTAP_1
+    DELTAP_1,
+    DAMP_SCALE
   };
 
 }    // namespace Granular_MDR_NS
@@ -97,7 +98,7 @@ class FixGranularMDR : public Fix {
   int index_sigmayy;            // yy-component of the stress tensor, not necessary forforce calculation
   int index_sigmazz;            // zz-component of the stress tensor, not necessary forforce calculation
   int index_history_setup_flag; // flag to check if history variables have beeninitialized
-  int index_dRavg;           // total contacts on particle
+  int index_dRavg;              // average radius update increment
 };
 
 }    // namespace LAMMPS_NS

src/GRANULAR/gran_sub_mod_damping.cpp

@@ -14,6 +14,7 @@
 #include "gran_sub_mod_damping.h"
 
 #include "gran_sub_mod_normal.h"
+#include "fix_granular_mdr.h"
 #include "granular_model.h"
 #include "math_special.h"
 #include "math_const.h"
@@ -176,3 +177,28 @@ void GranSubModDampingCoeffRestitution::init()
     damp /= sqrt(MY_PI * MY_PI + logcor * logcor);
   }
 }
+
+/* ----------------------------------------------------------------------
+   MDR damping
+------------------------------------------------------------------------- */
+
+GranSubModDampingMDR::GranSubModDampingMDR(GranularModel *gm, LAMMPS *lmp) :
+    GranSubModDamping(gm, lmp)
+{
+  contact_radius_flag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double GranSubModDampingMDR::calculate_forces()
+{
+  using namespace Granular_MDR_NS;
+  double *history = & gm->history[gm->normal_model->history_index]; // load in all history variables
+  //printf("%p %d damping\n", (void*)& history[DAMP_SCALE], gm->normal_model->history_index);
+  damp_prefactor = damp * history[DAMP_SCALE];
+  return -damp_prefactor * gm->vnnr;
+}
+
+  
+
+  //printf("DAMP_SCALE = %d, damp_scale_history = %e, damp = %e, damp_prefactor = %e, F_DAMP = %e  \n", DAMP_SCALE, history[DAMP_SCALE], damp, damp_prefactor, -damp_prefactor * gm->vnnr);

src/GRANULAR/gran_sub_mod_damping.h

@@ -19,6 +19,7 @@ GranSubModStyle(mass_velocity,GranSubModDampingMassVelocity,DAMPING);
 GranSubModStyle(viscoelastic,GranSubModDampingViscoelastic,DAMPING);
 GranSubModStyle(tsuji,GranSubModDampingTsuji,DAMPING);
 GranSubModStyle(coeff_restitution,GranSubModDampingCoeffRestitution,DAMPING);
+GranSubModStyle(mdr,GranSubModDampingMDR,DAMPING);
 // clang-format on
 #else
 
@@ -95,6 +96,14 @@ namespace Granular_NS {
 
   /* ---------------------------------------------------------------------- */
 
+  class GranSubModDampingMDR : public GranSubModDamping {
+   public:
+    GranSubModDampingMDR(class GranularModel *, class LAMMPS *);
+    double calculate_forces() override;
+  };
+
+  /* ---------------------------------------------------------------------- */
+
 }    // namespace Granular_NS
 }    // namespace LAMMPS_NS
 

src/GRANULAR/gran_sub_mod_normal.cpp

@@ -445,7 +445,7 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
 
   num_coeffs = 6;
   contact_radius_flag = 1;
-  size_history = 26;
+  size_history = 27;
   nsvector = 1;
   fix_mdr_flag = 0;
   id_fix = nullptr;
@@ -484,7 +484,7 @@ void GranSubModNormalMDR::coeffs_to_local()
   if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0");
   if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0");
   //if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0");
-  if (damp < 0.0) error->all(FLERR, "Illegal MDR normal model, damping coefficent must be greater than 1");
+  if (damp < 0.0) error->all(FLERR, "Illegal MDR normal model, damping coefficent must be greater than or equal to 0");
 
   G = E / (2.0 * (1.0 + nu));            // shear modulus
   kappa = E / (3.0 * (1.0 - 2.0 * nu));  // bulk modulus
@@ -1000,19 +1000,27 @@ double GranSubModNormalMDR::calculate_forces()
   const double wij = MAX(1.0 - pij, 0.0);
 
   // assign final force
-  double damp_prefactor;
+  double damp_scale;
   if (gm->contact_type != PAIR) {
     a_damp = a_damp/2.0;
-    damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff2particle * a_damp);
+    damp_scale = 0.5*sqrt(gm->meff * 2.0 * Eeff2particle * a_damp);
     double *deltao_offset = &history[DELTAO_0];
     const double wfm = std::exp(10.7*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0; // wall force magnifier
     //const double wfm = 1.0;
     F = wij * F0 * wfm;
   } else {
-    damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff * a_damp);
+    damp_scale = 0.5*sqrt(gm->meff * 2.0 * Eeff * a_damp);
     F = wij * (F0 + F1) * 0.5;
   }
 
+  if (history_update) {
+    double *damp_scale_offset = & history[DAMP_SCALE];
+    //printf("damp_scale_history = %e  \n", history[DAMP_SCALE]);
+    (a_damp < 0.0) ? *damp_scale_offset = 0.0 : *damp_scale_offset = damp_scale;
+    //printf("DAMP_SCALE = %d, damp_scale_history = %e, damp = %e, a_damp = %e, damp_scale = %e  \n", DAMP_SCALE, history[DAMP_SCALE], damp, a_damp, damp_scale);
+    //printf("%p %d normal\n", (void*)& history[DAMP_SCALE], history_index);
+  }
+
   // calculate damping force
   //if (F > 0.0) {
   //  double Eeff2;
@@ -1033,7 +1041,7 @@ double GranSubModNormalMDR::calculate_forces()
   //}
 
   // calculate damping force
-  if (F > 0.0) F += -wij * damp * gm->vnnr;
+  // if (F > 0.0) F += -wij * damp * gm->vnnr;
 
   //F += -wij * damp_prefactor * gm->vnnr;