diff --git a/src/Makefile.am b/src/Makefile.am index 4fa02ef8b5e8d1d63aee750ca1220075a72c4b75..89b6e40d91526597e9a4c465dc04fe0ef693cc66 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 449bcbe356cfda2636cd21b24d7e0ffaa0470ad0..0000000000000000000000000000000000000000 --- 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 3528e8206cb43d04908addd60beb9b9660f2308a..0000000000000000000000000000000000000000 --- 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 1a551016c8f2836b0f37b8c10587d6a8bad522fd..6a1b421480b36e9da407f89c45cf6ae4240c6a95 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); }