Proper time discretisation

dune/tectonic/globalruinanonlinearity.hh

@@ -26,13 +26,14 @@ class GlobalRuinaNonlinearity
       shared_ptr<BlockVector<FieldVector<double, 1>> const> mu,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> eta,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress,
-      shared_ptr<BlockVector<FieldVector<double, 1>> const> state)
+      shared_ptr<BlockVector<FieldVector<double, 1>> const> state, double h)
       : nodalIntegrals(nodalIntegrals),
         a(a),
         mu(mu),
         eta(eta),
         normalStress(normalStress),
         state(state),
+        h(h),
         trivialNonlinearity(
             new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
         restrictions(nodalIntegrals->size()) // TODO: can we get the size from
@@ -54,7 +55,7 @@ class GlobalRuinaNonlinearity
 
     auto const func = make_shared<RuinaFunction const>(
         (*nodalIntegrals)[i][0], (*a)[i][0], (*mu)[i][0], (*eta)[i][0],
-        (*normalStress)[i][0], (*state)[i][0]);
+        (*normalStress)[i][0], (*state)[i][0], h);
     restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
     return restrictions[i];
   }
@@ -68,6 +69,7 @@ class GlobalRuinaNonlinearity
   shared_ptr<BlockVector<FieldVector<double, 1>> const> eta;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> state;
+  double const h;
 
   std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
 };

dune/tectonic/nicefunction.hh

@@ -31,14 +31,15 @@ class NiceFunction : public VirtualFunction<double, double> {
 class RuinaFunction : public NiceFunction {
 public:
   RuinaFunction(double coefficient, double a, double mu, double eta,
-                double normalStress, double state)
+                double normalStress, double state, double h)
       : coefficient(coefficient),
         a(a),
         mu(mu),
         rho(exp(-mu / a)),
         eta(eta),
         normalStress(normalStress),
-        state(state) {}
+        state(state),
+        h(h) {}
 
   /*
     If mu and sigma_n denote the coefficient of friction and the
@@ -48,17 +49,18 @@ class RuinaFunction : public NiceFunction {
 
     with the constants a and b from Ruina's model and
 
-    h(beta) = beta (a (log beta - 1) + mu)
+    r(beta) = beta (a (log beta - 1) + mu)
 
     as well as
 
     rho = exp(-mu/a)
   */
   void virtual evaluate(double const &x, double &y) const {
-    double const arg = std::max(eta * x, rho);
-    double const h = arg * (a * (std::log(arg) - 1) + mu + state);
-    y = 1 / eta * normalStress * h + a * rho + mu + state;
+    double const arg = std::max(eta * x / h, rho);
+    double const r = arg * (a * (std::log(arg) - 1) + mu + state);
+    y = 1 / eta * normalStress * r + a * rho + mu + state;
     y *= coefficient;
+    y *= h;
   }
 
   /*
@@ -69,18 +71,20 @@ class RuinaFunction : public NiceFunction {
     = a * log(eta x) + mu
   */
   double virtual leftDifferential(double s) const {
-    if (eta * s <= rho)
+    double const arg = eta * s / h;
+    if (arg <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
+    return coefficient * normalStress * (a * std::log(arg) + mu + state);
   }
 
   /* see above */
   double virtual rightDifferential(double s) const {
-    if (eta * s <= rho)
+    double const arg = eta * s / h;
+    if (arg <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
+    return coefficient * normalStress * (a * std::log(arg) + mu + state);
   }
 
   /*
@@ -90,10 +94,10 @@ class RuinaFunction : public NiceFunction {
   */
   double virtual second_deriv(double s) const {
     // includes the case eta * s = rho for which there is no second derivative
-    if (eta * s <= rho)
+    if (eta * s / h <= rho)
       return 0;
     else
-      return coefficient * normalStress * (a / s);
+      return coefficient * normalStress * a / s;
   }
 
   double virtual regularity(double s) const {
@@ -110,6 +114,7 @@ class RuinaFunction : public NiceFunction {
   double const eta;
   double const normalStress;
   double const state;
+  double const h;
 
   double const rho;
 };

src/one-body-sample.cc

@@ -103,11 +103,12 @@ template <class GridType, class GridView, class LocalVectorType, class FEBasis>
 void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                       Dune::BitSetVector<1> const &neumannNodes,
                       Dune::BlockVector<LocalVectorType> &f,
-                      int run) { // constant sample function on neumann boundary
+                      double time) { // constant sample function on neumann
+                                     // boundary
   BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
   LocalVectorType SampleVector(0);
   // FIXME: random values (time-dependent)
-  SampleVector[0] = run / 10.0;
+  SampleVector[0] = time;
   SampleVector[1] = 0;
   ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
   NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
@@ -145,7 +146,8 @@ assemble_nonlinearity(
     int size, Dune::ParameterTree const &parset,
     Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
         nodalIntegrals,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state) {
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+    double h) {
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
 
   // {{{ Assemble terms for the nonlinearity
@@ -166,7 +168,7 @@ assemble_nonlinearity(
 
     return Dune::make_shared<
         Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
-        nodalIntegrals, a, mu, eta, normalStress, state);
+        nodalIntegrals, a, mu, eta, normalStress, state, h);
   } else if (friction_model == std::string("Laursen")) {
     return Dune::make_shared<Dune::GlobalLaursenNonlinearity<
         Dune::LinearFunction, VectorType, MatrixType> const>(mu, normalStress,
@@ -312,6 +314,7 @@ int main(int argc, char *argv[]) {
 
     timer.reset();
     auto const timesteps = parset.get<size_t>("timesteps");
+    double h = 1.0 / timesteps;
     for (size_t run = 1; run <= timesteps; ++run) {
       if (parset.get<bool>("printProgress")) {
         std::cout << ".";
@@ -320,7 +323,7 @@ int main(int argc, char *argv[]) {
 
       // b = neumann
       assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
-          leafView, p1Basis, neumannNodes, b1, run);
+          leafView, p1Basis, neumannNodes, b1, h * run);
       b2 = b1;
       b3 = b1;
       b4 = b1;
@@ -339,7 +342,7 @@ int main(int argc, char *argv[]) {
         auto myGlobalNonlinearity =
             assemble_nonlinearity<VectorType, OperatorType>(
                 grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
-                state);
+                state, h);
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                             *myGlobalNonlinearity, b1);
         MyBlockProblemType myBlockProblem(parset, myConvexProblem);
@@ -362,7 +365,7 @@ int main(int argc, char *argv[]) {
         auto myGlobalNonlinearity =
             assemble_nonlinearity<VectorType, OperatorType>(
                 grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
-                state);
+                state, h);
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                             *myGlobalNonlinearity, b4);