[Cleanup] Rework Dieterich state update

* Merge three files into one
* Merge two functions into one (*_bisection is gone)
* Move linear term out of the nonlinearity
* Use DirectionalConvexFunction from dune-tnnmg

src/Makefile.am

@@ -8,7 +8,6 @@ bin_PROGRAMS = \
 SOURCES = \
 	assemblers.cc \
 	state/compute_state_dieterich_euler.cc \
-	state/compute_state_dieterich_common.cc \
 	state/compute_state_ruina.cc \
 	solverfactory.cc \
 	one-body-sample.cc \

src/state/compute_state_dieterich_common.cc

-#include "compute_state_dieterich_common.hh"
-
-DieterichNonlinearity::DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
-
-void DieterichNonlinearity::directionalSubDiff(VectorType const &u,
-                                               VectorType const &v,
-                                               Interval<double> &D) const {
-  D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
-}
-
-void DieterichNonlinearity::directionalDomain(VectorType const &,
-                                              VectorType const &,
-                                              Interval<double> &dom) const {
-  dom[0] = -std::numeric_limits<double>::max();
-  dom[1] = std::numeric_limits<double>::max();
-}

src/state/compute_state_dieterich_common.hh

-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
-
-class DieterichNonlinearity {
-public:
-  using VectorType = Dune::FieldVector<double, 1>;
-  using MatrixType = Dune::FieldMatrix<double, 1, 1>;
-
-  DieterichNonlinearity(double _VoL);
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const;
-  void directionalDomain(VectorType const &, VectorType const &,
-                         Interval<double> &dom) const;
-
-private:
-  double const VoL;
-};

src/state/compute_state_dieterich_euler.cc

@@ -2,36 +2,72 @@
 
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include <dune/tnnmg/problem-classes/directionalconvexfunction.hh>
 
-#include "compute_state_dieterich_common.hh"
 #include "compute_state_dieterich_euler.hh"
 
-double state_update_dieterich_euler_bisection(double tau, double VoL,
-                                              double old_state) {
-  DieterichNonlinearity::VectorType const start(0);
-  DieterichNonlinearity::VectorType const direction(1);
-  DieterichNonlinearity const phi(VoL);
-  MyDirectionalConvexFunction<DieterichNonlinearity> const J(
-      1.0 / tau, old_state / tau, phi, start, direction);
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/fufem/interval.hh>
 
-  int bisectionsteps = 0;
-  Bisection const bisection;
-  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
-}
+/* Dieterich's law:
 
-double state_update_dieterich_euler(double tau, double VoL, double old_state) {
-  double const ret =
-      state_update_dieterich_euler_bisection(tau, VoL, old_state);
-  /* We have
+   \f$ \dot \theta = 1 - V \theta/L \f$
+
+   Transformed using \alpha = \log \theta:
+
+   \f$ \dot \alpha = \exp(-\alpha) - V/L \f$
+
+   discretised in time (using backward euler):
+
+   \f$\delta \alpha_1 = \exp(-\alpha_1) - V/L\f$
+
+   or
+
+   \f$ 0 = ( \alpha_1 - \alpha_0 )/\tau - ( \exp(-\alpha_1) - V/L ) \f$
 
-     ret - old_state + tau * VoL = tau*exp(-ret)
+   this is obviously monotone in \f$\alpha_1\f$; we can thus see if as
+   a convex minimisation problem.
 
-     or
+   find the antiderivative:
 
-     log((ret - old_state)/tau + VoL) = -ret
-  */
-  assert(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))) < 1e-12 ||
-         std::abs(std::log((ret - old_state) / tau + VoL) + ret) < 1e-12);
-  return ret;
+   \f$ J(\alpha_1)
+   = 0.5/\tau * \alpha_1^2 - ( \alpha_0/\tau - V/L ) \alpha_1 +
+   \exp(-\alpha_1)\f$
+
+   Consequence: linear part is \f$ \alpha_0/\tau - V/L \f$; nonlinear
+   part is \f$\exp(-\alpha_1)\f$ which has derivative \f$
+   -\exp(-\alpha_1) \f$.
+*/
+
+//! The nonlinearity \f$f(x) = exp(-x)\f$
+class DieterichNonlinearity {
+public:
+  using VectorType = Dune::FieldVector<double, 1>;
+  using MatrixType = Dune::FieldMatrix<double, 1, 1>;
+
+  DieterichNonlinearity() {}
+
+  void directionalSubDiff(VectorType const &u, VectorType const &v,
+                          Interval<double> &D) const {
+    D[0] = D[1] = v[0] * (-std::exp(-u[0]));
+  }
+
+  void directionalDomain(VectorType const &, VectorType const &,
+                         Interval<double> &dom) const {
+    dom[0] = -std::numeric_limits<double>::max();
+    dom[1] = std::numeric_limits<double>::max();
+  }
+};
+
+double state_update_dieterich_euler(double tau, double VoL, double old_state) {
+  DieterichNonlinearity::VectorType const start(0);
+  DieterichNonlinearity::VectorType const direction(1);
+  DieterichNonlinearity phi;
+  DirectionalConvexFunction<DieterichNonlinearity> const J(
+      1.0 / tau, old_state / tau - VoL, phi, start, direction);
+
+  int bisectionsteps = 0;
+  Bisection const bisection;
+  return bisection.minimize(J, 0.0, 0.0, bisectionsteps);
 }