diff --git a/src/Makefile.am b/src/Makefile.am
index 4cca1056223006b91a0d55c175242bca7d0a0bb5..b648151d932a2f0dff6b7b3d15e6bd29030ebb45 100644
--- a/src/Makefile.am
+++ b/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 \
diff --git a/src/state/compute_state_dieterich_euler.cc b/src/state/compute_state_dieterich_euler.cc
deleted file mode 100644
index 103c56f2c3d514c086ad91cfdf58bf71a4191220..0000000000000000000000000000000000000000
--- a/src/state/compute_state_dieterich_euler.cc
+++ /dev/null
@@ -1,73 +0,0 @@
-#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);
-}
diff --git a/src/state/compute_state_dieterich_euler.hh b/src/state/compute_state_dieterich_euler.hh
deleted file mode 100644
index 95690601f2d6caaa602afcb2cba4c2ae22cbd279..0000000000000000000000000000000000000000
--- a/src/state/compute_state_dieterich_euler.hh
+++ /dev/null
@@ -1,5 +0,0 @@
-#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
diff --git a/src/state/compute_state_ruina.cc b/src/state/compute_state_ruina.cc
deleted file mode 100644
index 4ad016936a1ec4c46d228a5d07b7c47c85145f04..0000000000000000000000000000000000000000
--- a/src/state/compute_state_ruina.cc
+++ /dev/null
@@ -1,13 +0,0 @@
-#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);
-}
diff --git a/src/state/compute_state_ruina.hh b/src/state/compute_state_ruina.hh
deleted file mode 100644
index 07dc4d7549fe91a0e8e5294c7bec037ace6b1df9..0000000000000000000000000000000000000000
--- a/src/state/compute_state_ruina.hh
+++ /dev/null
@@ -1,5 +0,0 @@
-#ifndef COMPUTE_STATE_RUINA_HH
-#define COMPUTE_STATE_RUINA_HH
-
-double state_update_ruina(double h, double uol, double logState_o);
-#endif
diff --git a/src/state/dieterichstateupdater.hh b/src/state/dieterichstateupdater.hh
index cd28d97810b046077904047281df23ee5a788083..ffb417a08d39a3e3b5ebe1d040869664dbfbd75e 100644
--- a/src/state/dieterichstateupdater.hh
+++ b/src/state/dieterichstateupdater.hh
@@ -1,7 +1,6 @@
#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>
diff --git a/src/state/ruinastateupdater.hh b/src/state/ruinastateupdater.hh
index 18604f991a105f8e5f871fc02a7f5c056e4f9be9..d87870bb0716508d32176a46e08265562770e333 100644
--- a/src/state/ruinastateupdater.hh
+++ b/src/state/ruinastateupdater.hh
@@ -1,7 +1,6 @@
#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>