From 4e4fca3701b9991f0f6782decd769b6bb547cede Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 18 Sep 2012 12:56:05 +0200
Subject: [PATCH] Do not split the Dieterich nonlinearity

It makes the code harder to read and the computational benefit is
questionable
---
 src/compute_state_dieterich.cc | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/src/compute_state_dieterich.cc b/src/compute_state_dieterich.cc
index b69a121b..3b8441cf 100644
--- a/src/compute_state_dieterich.cc
+++ b/src/compute_state_dieterich.cc
@@ -10,18 +10,17 @@
 
 #include "compute_state_dieterich.hh"
 
-// The nonlinearity exp(-x)
-// NOTE: the Dieterich law has an additional linear term!
-class DecayingExponential {
+// psi(beta) = V/L beta + e^(-beta)
+class DieterichNonlinearity {
 public:
   typedef Dune::FieldVector<double, 1> VectorType;
   typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
 
-  DecayingExponential(double tau) : tau(tau) {}
+  DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
 
   void directionalSubDiff(VectorType const &u, VectorType const &v,
                           Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-tau * std::exp(-u[0]));
+    D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
@@ -31,16 +30,16 @@ class DecayingExponential {
   }
 
 private:
-  double const tau;
+  double const VoL;
 };
 
 double state_update_dieterich_bisection(double tau, double VoL,
                                         double old_state) {
-  DecayingExponential::VectorType const start(0);
-  DecayingExponential::VectorType const direction(1);
-  DecayingExponential const phi(tau);
-  MyDirectionalConvexFunction<DecayingExponential> const J(
-      1.0, old_state - tau * VoL, phi, start, direction);
+  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);
 
   int bisectionsteps = 0;
   Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
-- 
GitLab