diff --git a/dune/tectonic/frictiondata.hh b/dune/tectonic/frictiondata.hh index 8cc724698a9c163e7a5c157dbef4a2179e1f3817..656683f90283355a77e9855677270fe58c95ef25 100644 --- a/dune/tectonic/frictiondata.hh +++ b/dune/tectonic/frictiondata.hh @@ -9,6 +9,7 @@ struct FrictionData { a(parset.get<double>("a")), b(parset.get<double>("b")), mu0(parset.get<double>("mu0")), + mumin(parset.get<double>("mumin")), normalStress(normalStress) {} double const L; @@ -16,6 +17,7 @@ struct FrictionData { double const a; double const b; double const mu0; + double const mumin; double const normalStress; }; #endif diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh index bc2b22739bc28281e3caaee9e209be0faca79c48..73d17dbaffedc515245d625142008988e81cb161 100644 --- a/dune/tectonic/frictionpotential.hh +++ b/dune/tectonic/frictionpotential.hh @@ -34,23 +34,15 @@ class FrictionPotential : public FrictionPotentialWrapper { FrictionPotential(double coefficient, FrictionData const &fd) : fd(fd), weightTimesNormalStress(coefficient * fd.normalStress) {} - double evaluate(double V) const { - assert(V >= 0); - if (V <= V_m) - return 0; - - // V log(V/V_m) - V + V_m - return weightTimesNormalStress * fd.a * (V * std::log(V / V_m) - V + V_m); - } - // log(V/V_m) if V >= V_0 // 0 otherwise double differential(double V) const { - assert(V >= 0); - if (V <= V_m) - return 0; + assert(V >= 0.0); + if (V <= V_cutoff) { + return fd.mumin; + } - return weightTimesNormalStress * fd.a * std::log(V / V_m); + return weightTimesNormalStress * fd.a * (std::log(V) - logV_m); } // 1/V if V > V_0 @@ -58,7 +50,7 @@ class FrictionPotential : public FrictionPotentialWrapper { // 0 if V < V_0 double second_deriv(double V) const { assert(V >= 0); - if (V <= V_m) + if (V <= V_cutoff) return 0; return weightTimesNormalStress * fd.a / V; @@ -67,7 +59,7 @@ class FrictionPotential : public FrictionPotentialWrapper { double regularity(double V) const { assert(V >= 0); // TODO: Make this controllable - if (std::abs(V - V_m) < 1e-14) + if (std::abs(V - V_cutoff) < 1e-14) return std::numeric_limits<double>::infinity(); return std::abs(second_deriv(V)); @@ -77,14 +69,16 @@ class FrictionPotential : public FrictionPotentialWrapper { // with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ] void updateState(double state) { // state is assumed to be logarithmic - V_m = fd.V0 * - std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a); + logV_m = std::log(fd.V0) + + (-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a); + V_cutoff = std::exp(logV_m + fd.mumin / fd.a); } private: FrictionData const &fd; double const weightTimesNormalStress; - double V_m; + double logV_m; + double V_cutoff; // velocity at which mu falls short of mumin }; class TrivialFunction : public FrictionPotentialWrapper { diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset index 1f819b2408ba7f4b41c404f6addd5a4bf47e0941..d3eecdd8b256d56b35cfcac88cc1b212738c9f9a 100644 --- a/src/one-body-sample.parset +++ b/src/one-body-sample.parset @@ -51,6 +51,7 @@ requiredResidual = 1e-12 [boundary.friction] mu0 = 0.6 +mumin = 0.0 a = 0.010 b = 0.015 V0 = 1e-6