From 41b575c5aa4d787d739e95a9c07fbd1c474c1a01 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 21 Jul 2013 20:30:38 +0200
Subject: [PATCH] [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 | 1 -
src/state/compute_state_dieterich_common.cc | 16 ----
src/state/compute_state_dieterich_common.hh | 18 -----
src/state/compute_state_dieterich_euler.cc | 84 +++++++++++++++------
4 files changed, 60 insertions(+), 59 deletions(-)
delete mode 100644 src/state/compute_state_dieterich_common.cc
delete mode 100644 src/state/compute_state_dieterich_common.hh
diff --git a/src/Makefile.am b/src/Makefile.am
index 4fa02ef8..89b6e40d 100644
--- a/src/Makefile.am
+++ b/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 \
diff --git a/src/state/compute_state_dieterich_common.cc b/src/state/compute_state_dieterich_common.cc
deleted file mode 100644
index 449bcbe3..00000000
--- a/src/state/compute_state_dieterich_common.cc
+++ /dev/null
@@ -1,16 +0,0 @@
-#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();
-}
diff --git a/src/state/compute_state_dieterich_common.hh b/src/state/compute_state_dieterich_common.hh
deleted file mode 100644
index 3528e820..00000000
--- a/src/state/compute_state_dieterich_common.hh
+++ /dev/null
@@ -1,18 +0,0 @@
-#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;
-};
diff --git a/src/state/compute_state_dieterich_euler.cc b/src/state/compute_state_dieterich_euler.cc
index 1a551016..6a1b4214 100644
--- a/src/state/compute_state_dieterich_euler.cc
+++ b/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);
}
--
GitLab