[Algorit] Solve state-ODE analytically

src/Makefile.am

@@ -5,8 +5,6 @@ SOURCES = \
 	assemblers.cc \
 	boundary_writer.cc \
 	friction_writer.cc \
-	state/compute_state_dieterich_euler.cc \
-	state/compute_state_ruina.cc \
 	solverfactory.cc \
 	one-body-sample.cc \
 	timestepping.cc \

src/state/compute_state_dieterich_euler.cc

-#include <cassert>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tnnmg/problem-classes/directionalconvexfunction.hh>
-
-#include "compute_state_dieterich_euler.hh"
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/solvers/common/interval.hh>
-
-/* Dieterich's law:
-
-   \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$
-
-   this is obviously monotone in \f$\alpha_1\f$; we can thus see if as
-   a convex minimisation problem.
-
-   find the antiderivative:
-
-   \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,
-                          Dune::Solvers::Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-std::exp(-u[0]));
-  }
-
-  void directionalDomain(VectorType const &, VectorType const &,
-                         Dune::Solvers::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 logState_o) {
-  DieterichNonlinearity::VectorType const start(0);
-  DieterichNonlinearity::VectorType const direction(1);
-  DieterichNonlinearity phi;
-  DirectionalConvexFunction<DieterichNonlinearity> const J(
-      1.0 / tau, logState_o / tau - VoL, phi, start, direction);
-
-  int bisectionsteps = 0;
-  Bisection const bisection;
-  return bisection.minimize(J, 0.0, 0.0, bisectionsteps);
-}

src/state/compute_state_dieterich_euler.hh

-#ifndef COMPUTE_STATE_DIETERICH_EULER_HH
-#define COMPUTE_STATE_DIETERICH_EULER_HH
-
-double state_update_dieterich_euler(double h, double VoL, double logState_o);
-#endif

src/state/compute_state_ruina.cc

-#include <cassert>
-
-#include <cmath>
-
-#include "compute_state_ruina.hh"
-
-double state_update_ruina(double tau, double uol, double logState_o) {
-  if (uol == 0)
-    return logState_o;
-
-  double const ret = logState_o - uol * std::log(uol / tau);
-  return ret / (1 + uol);
-}

src/state/compute_state_ruina.hh

-#ifndef COMPUTE_STATE_RUINA_HH
-#define COMPUTE_STATE_RUINA_HH
-
-double state_update_ruina(double h, double uol, double logState_o);
-#endif

src/state/dieterichstateupdater.hh

 #ifndef DIETERICH_STATE_UPDATER_HH
 #define DIETERICH_STATE_UPDATER_HH
 
-#include "compute_state_dieterich_euler.hh"
 #include "stateupdater.hh"
 
 template <class ScalarVector, class Vector>
@@ -39,14 +38,33 @@ void DieterichStateUpdater<ScalarVector, Vector>::setup(double _tau) {
   tau = _tau;
 }
 
+/*
+  Approximate (1-exp(-tV))/V, i.e. t*(exp(x) - 1)/x
+  for x = -t*v close to 0
+
+  with t, V >= 0, so that x should be non-positive
+*/
+double approximateMyPowerSeries(double V, double t) {
+  double mVt = -V * t;
+  if (mVt >= 0)
+    return t;
+
+  return -std::expm1(mVt) / V;
+}
+
 template <class ScalarVector, class Vector>
 void DieterichStateUpdater<ScalarVector, Vector>::solve(
     Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i)
-    if (nodes[i][0]) {
-      double const V = velocity_field[i].two_norm();
-      logState[i] = state_update_dieterich_euler(tau, V / L, logState_o[i]);
-    }
+  for (size_t i = 0; i < nodes.size(); ++i) {
+    if (not nodes[i][0])
+      continue;
+
+    double const VoL = velocity_field[i].two_norm() / L;
+    // log( (1-exp(-tV))/V + exp(alpha_old - tV) ) with a
+    // liftable singularity at V=0
+    logState[i] = std::log(approximateMyPowerSeries(VoL, tau) +
+                           std::exp(logState_o[i] - tau * VoL));
+  }
 }
 
 template <class ScalarVector, class Vector>

src/state/ruinastateupdater.hh

 #ifndef RUINA_STATE_UPDATER_HH
 #define RUINA_STATE_UPDATER_HH
 
-#include "compute_state_ruina.hh"
 #include "stateupdater.hh"
 
 template <class ScalarVector, class Vector>
@@ -42,11 +41,16 @@ void RuinaStateUpdater<ScalarVector, Vector>::setup(double _tau) {
 template <class ScalarVector, class Vector>
 void RuinaStateUpdater<ScalarVector, Vector>::solve(
     Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i)
-    if (nodes[i][0]) {
-      double const V = velocity_field[i].two_norm();
-      logState[i] = state_update_ruina(tau, V * tau / L, logState_o[i]);
-    }
+  for (size_t i = 0; i < nodes.size(); ++i) {
+    if (not nodes[i][0])
+      continue;
+
+    double const VoL = velocity_field[i].two_norm() / L;
+    // (exp(-tV) - 1)*log(V) + alpha_old*exp(-tV)
+    logState[i] =
+        (VoL <= 0) ? logState_o[i] : std::expm1(-tau * VoL) * std::log(VoL) +
+                                         logState_o[i] * std::exp(-tau * VoL);
+  }
 }
 
 template <class ScalarVector, class Vector>