[Extend] Add regularised rate/state law

dune/tectonic/frictionpotential.hh

@@ -69,6 +69,38 @@ class TruncatedRateState : public FrictionPotential {
   double Vmin;
 };
 
+class RegularisedRateState : public FrictionPotential {
+public:
+  RegularisedRateState(double coefficient, double _normalStress,
+                       FrictionData _fd)
+      : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
+
+  double coefficientOfFriction(double V) const {
+    return fd.a * std::asinh(0.5 * V / Vmin);
+  }
+
+  double differential(double V) const {
+    return weight * (fd.C - normalStress * coefficientOfFriction(V));
+  }
+
+  double second_deriv(double V) const {
+    return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V);
+  }
+
+  double regularity(double V) const { return std::abs(second_deriv(V)); }
+
+  void updateAlpha(double alpha) {
+    double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
+    Vmin = fd.V0 / std::exp(logrest);
+  }
+
+private:
+  FrictionData const fd;
+  double const weight;
+  double const normalStress;
+  double Vmin;
+};
+
 class TrivialFunction : public FrictionPotential {
 public:
   double evaluate(double) const { return 0; }

src/assemblers.cc

@@ -122,6 +122,7 @@ void MyAssembler<GridView, dimension>::assembleNormalStress(
 
 template <class GridView, int dimension>
 auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+    Config::FrictionModel frictionModel,
     BoundaryPatch<GridView> const &frictionalBoundary,
     GlobalFrictionData<dimension> const &frictionInfo,
     ScalarVector const &normalStress)
@@ -136,9 +137,18 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
-  return std::make_shared<
-      GlobalRateStateFriction<Matrix, Vector, TruncatedRateState, GridView>>(
-      frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+  switch (frictionModel) {
+    case Config::Truncated:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, TruncatedRateState, GridView>>(
+          frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+    case Config::Regularised:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, RegularisedRateState, GridView>>(
+          frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+    default:
+      assert(false);
+  }
 }
 
 template <class GridView, int dimension>

src/assemblers.hh

@@ -16,6 +16,8 @@
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
 
+#include "enums.hh"
+
 template <class GridView, int dimension> class MyAssembler {
 public:
   using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
@@ -76,6 +78,7 @@ template <class GridView, int dimension> class MyAssembler {
                             double poissonRatio, Vector const &displacement);
 
   std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
+      Config::FrictionModel frictionModel,
       BoundaryPatch<GridView> const &frictionalBoundary,
       GlobalFrictionData<dimension> const &frictionInfo,
       ScalarVector const &normalStress);

src/enum_friction_model.cc

+#include <dune/common/exceptions.hh>
+
+#include "enums.hh"
+
+template <> struct StringToEnum<Config::FrictionModel> {
+  static Config::FrictionModel convert(std::string const &s) {
+    if (s == "Truncated")
+      return Config::Truncated;
+
+    if (s == "Regularised")
+      return Config::Regularised;
+
+    DUNE_THROW(Dune::Exception, "failed to parse enum");
+  }
+};

src/enums.hh

@@ -2,6 +2,10 @@
 #define SRC_ENUMS_HH
 
 struct Config {
+  enum FrictionModel {
+    Truncated,
+    Regularised
+  };
   enum stateModel {
     AgeingLaw,
     SlipLaw

src/sand-wedge-data/parset.cfg

@@ -27,6 +27,7 @@ V0              = 5e-5  # [m/s]
 L               = 2e-5  # [m]      ?
 initialAlpha    = 0.916290731874155 # [ ]      ?
 stateModel      = AgeingLaw
+frictionModel   = Truncated
 [boundary.friction.weakening]
 a               = 0.015 # [1]      ?
 b               = 0.030 # [1]      ?

src/sand-wedge.cc

@@ -68,6 +68,7 @@
 #include "assemblers.hh"
 #include "tobool.hh"
 #include "enum_parser.cc"
+#include "enum_friction_model.cc"
 #include "enum_scheme.cc"
 #include "enum_state_model.cc"
 #include "enum_verbosity.cc"
@@ -248,6 +249,7 @@ int main(int argc, char *argv[]) {
                                      body.getPoissonRatio(), u_initial);
     MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
     auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
+        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
         frictionalBoundary, frictionInfo, normalStress);
     myGlobalFriction->updateAlpha(alpha_initial);