diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index 6d2c0154db7226136e34cb8bb5ffbb5609d33c05..59e7ec53aecc5fd23b96ed22ec46a2a598721cc4 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -24,7 +24,7 @@ class FrictionPotentialWrapper {
DUNE_THROW(Dune::NotImplemented, "evaluation not implemented");
}
- void virtual updateLogState(double) = 0;
+ void virtual updateAlpha(double) = 0;
};
class FrictionPotential : public FrictionPotentialWrapper {
@@ -60,8 +60,8 @@ class FrictionPotential : public FrictionPotentialWrapper {
return std::abs(second_deriv(V));
}
- void updateLogState(double logState) {
- logrest = (fd.mu0 + fd.b * (logState + std::log(fd.V0 / fd.L))) / fd.a;
+ void updateAlpha(double alpha) {
+ logrest = (fd.mu0 + fd.b * alpha) / fd.a;
Vmin = fd.V0 / std::exp(logrest);
}
@@ -85,6 +85,6 @@ class TrivialFunction : public FrictionPotentialWrapper {
double regularity(double) const { return 0; }
- void updateLogState(double) {}
+ void updateAlpha(double) {}
};
#endif
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 2f00da1b4091585baf39065dadd988e21d28aa93..60dde9b4bb1ed1501d7c4fba27069bf11656e6b2 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -82,6 +82,6 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
return res->regularity(x);
}
- void virtual updateLogState(ScalarVector const &logState) = 0;
+ void virtual updateAlpha(ScalarVector const &alpha) = 0;
};
#endif
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 4c82a517e057cc22d779293c3d301ff844d5935b..69c751b46ae3e4ee150a553c97233c0ef061ed09 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -51,9 +51,9 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
}
}
- void updateLogState(ScalarVector const &logState) override {
+ void updateAlpha(ScalarVector const &alpha) override {
for (size_t i = 0; i < restrictions.size(); ++i)
- restrictions[i]->updateLogState(logState[i]);
+ restrictions[i]->updateAlpha(alpha[i]);
}
void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index eaea43c4f0f5850e23c96acee0a7d840ace61ed8..ed96fdd5000aa8b7fd1a49065854ac1201678644 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -24,7 +24,7 @@ template <size_t dimension> class LocalFriction {
return func->evaluate(x.two_norm());
}
- void updateLogState(double logState) { func->updateLogState(logState); }
+ void updateAlpha(double alpha) { func->updateAlpha(alpha); }
double regularity(VectorType const &x) const {
double const xnorm = x.two_norm();
diff --git a/src/friction_writer.cc b/src/friction_writer.cc
index 175f8e6a94bfbd8c6ccede17bbe0d702b07ef51a..cc24c2e4ea739d00ba2862b1aa3ee530b82924fc 100644
--- a/src/friction_writer.cc
+++ b/src/friction_writer.cc
@@ -12,7 +12,7 @@ FrictionWriter<ScalarVector, Vector>::FrictionWriter(
typename BW::Projector projector)
: BW(vertexCoordinates, _boundaryNodes, prefix, projector),
coefficientWriter(prefix + "Coefficients", std::fstream::out),
- stateWriter(prefix + "LogStates", std::fstream::out) {}
+ stateWriter(prefix + "Alpha", std::fstream::out) {}
template <class ScalarVector, class Vector>
FrictionWriter<ScalarVector, Vector>::~FrictionWriter() {
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 0becf2b53bd71ae9765a8b5df39cbdb7be77ba96..e1adc85ac6eeb66f318fd83e09e12631d0f4470f 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/src/sand-wedge-data/parset.cfg
@@ -25,7 +25,7 @@ C = 10 # [Pa]
mu0 = 0.7 # [1]
V0 = 5e-5 # [m/s]
L = 2e-5 # [m] ?
-initialLogState = 0 # [ ] ?
+initialAlpha = 0.916290731874155 # [ ] ?
stateModel = Dieterich
[boundary.friction.weakening]
a = 0.015 # [1] ?
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 65f2ec000da34f3f63e5b0a3959827f42b2f9b65..ec95af912e27fa9a4a55b5b90e29ca44260e9ab2 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -243,7 +243,7 @@ int main(int argc, char *argv[]) {
ur_initial = 0.0;
ScalarVector alpha_initial(fineVertexCount);
- alpha_initial = parset.get<double>("boundary.friction.initialLogState");
+ alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
ScalarVector normalStress(fineVertexCount);
myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
body.getYoungModulus(),
@@ -251,7 +251,7 @@ int main(int argc, char *argv[]) {
MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
frictionalBoundary, frictionInfo, normalStress);
- myGlobalNonlinearity->updateLogState(alpha_initial);
+ myGlobalNonlinearity->updateAlpha(alpha_initial);
Vector v_initial(fineVertexCount);
v_initial = 0.0;
@@ -350,7 +350,8 @@ int main(int argc, char *argv[]) {
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes,
- parset.get<double>("boundary.friction.L"));
+ parset.get<double>("boundary.friction.L"),
+ parset.get<double>("boundary.friction.V0"));
Vector v = v_initial;
Vector v_m(fineVertexCount);
@@ -395,7 +396,7 @@ int main(int argc, char *argv[]) {
size_t iterationCounter;
auto solveVelocityProblem = [&](Vector &_velocityIterate,
ScalarVector const &_alpha) {
- myGlobalNonlinearity->updateLogState(_alpha);
+ myGlobalNonlinearity->updateAlpha(_alpha);
// NIT: Do we really need to pass u here?
typename NonlinearFactory::ConvexProblem convexProblem(
@@ -420,7 +421,7 @@ int main(int argc, char *argv[]) {
Arithmetic::addProduct(v_m, lambda, v);
stateUpdater->solve(v_m);
- stateUpdater->extractLogState(alpha);
+ stateUpdater->extractAlpha(alpha);
if (stateFPI == 1)
relaxationWriter << "N ";
diff --git a/src/state.hh b/src/state.hh
index 9d752b651a2dc9eebec20b85eea41fbbc048a886..57c24f622095712d3de33c71fe979390592bf177 100644
--- a/src/state.hh
+++ b/src/state.hh
@@ -11,14 +11,14 @@
template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
Config::stateModel model, ScalarVector const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, double L) {
+ Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) {
switch (model) {
case Config::Dieterich:
return std::make_shared<DieterichStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, L);
+ alpha_initial, frictionalNodes, L, V0);
case Config::Ruina:
return std::make_shared<RuinaStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, L);
+ alpha_initial, frictionalNodes, L, V0);
default:
assert(false);
}
diff --git a/src/state/dieterichstateupdater.hh b/src/state/dieterichstateupdater.hh
index 30a9c6dc57a2ede08a836e5857fe715356231061..f8414bd4468fc9fbdc2db009b6ed0c6ef120ba9c 100644
--- a/src/state/dieterichstateupdater.hh
+++ b/src/state/dieterichstateupdater.hh
@@ -7,31 +7,33 @@
template <class ScalarVector, class Vector>
class DieterichStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- DieterichStateUpdater(ScalarVector _logState_initial,
- Dune::BitSetVector<1> const &_nodes, double _L);
+ DieterichStateUpdater(ScalarVector _alpha_initial,
+ Dune::BitSetVector<1> const &_nodes, double _L,
+ double _V0);
void virtual nextTimeStep();
void virtual setup(double _tau);
void virtual solve(Vector const &velocity_field);
- void virtual extractLogState(ScalarVector &);
+ void virtual extractAlpha(ScalarVector &);
private:
- ScalarVector logState_o;
- ScalarVector logState;
+ ScalarVector alpha_o;
+ ScalarVector alpha;
Dune::BitSetVector<1> const &nodes;
double const L;
+ double const V0;
double tau;
};
template <class ScalarVector, class Vector>
DieterichStateUpdater<ScalarVector, Vector>::DieterichStateUpdater(
- ScalarVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
- double _L)
- : logState(_logState_initial), nodes(_nodes), L(_L) {}
+ ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L,
+ double _V0)
+ : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector>
void DieterichStateUpdater<ScalarVector, Vector>::nextTimeStep() {
- logState_o = logState;
+ alpha_o = alpha;
}
template <class ScalarVector, class Vector>
@@ -40,17 +42,14 @@ void DieterichStateUpdater<ScalarVector, Vector>::setup(double _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
+ Compute [ 1-\exp(c*x) ] / x under the assumption that x is
+ non-negative
*/
-double approximateMyPowerSeries(double V, double t) {
- double mVt = -V * t;
- if (mVt >= 0)
- return t;
-
- return -std::expm1(mVt) / V;
+double liftSingularity(double c, double x) {
+ if (x <= 0)
+ return -c;
+ else
+ return -std::expm1(c * x) / x;
}
template <class ScalarVector, class Vector>
@@ -60,18 +59,17 @@ void DieterichStateUpdater<ScalarVector, Vector>::solve(
if (not toBool(nodes[i]))
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));
+ double const V = velocity_field[i].two_norm();
+ double const mtoL = -tau / L;
+ alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) +
+ V0 * liftSingularity(mtoL, V));
}
}
template <class ScalarVector, class Vector>
-void DieterichStateUpdater<ScalarVector, Vector>::extractLogState(
- ScalarVector &_logState) {
- _logState = logState;
+void DieterichStateUpdater<ScalarVector, Vector>::extractAlpha(
+ ScalarVector &_alpha) {
+ _alpha = alpha;
}
#endif
diff --git a/src/state/explicit.tex b/src/state/explicit.tex
new file mode 100644
index 0000000000000000000000000000000000000000..4e60ce3470930938f74d4e2c62fdc25049565fe0
--- /dev/null
+++ b/src/state/explicit.tex
@@ -0,0 +1,16 @@
+\documentclass[preview]{standalone}
+\providecommand\Vref{V_\ast}
+\usepackage{amsmath}
+\usepackage{mathtools}
+\DeclarePairedDelimiter\paren{\lparen}{\rparen}
+\begin{document}
+\begin{equation*}
+ \alpha(t) =
+ \begin{cases*}
+ \log\paren*{e^{\alpha_0-tV/L} + \Vref \frac{1 - e^{-tV/L}}V}
+ &for the ageing law,\\
+ (e^{-\frac {tV}L} - 1) \log \frac V\Vref + \alpha_0 e^{-\frac {tV}L}
+ &for the slip law.
+ \end{cases*}
+\end{equation*}
+\end{document}
diff --git a/src/state/ruinastateupdater.hh b/src/state/ruinastateupdater.hh
index 578f1e9fecd517dc3ce02182c4a417705a723517..616a837b620eaca1e8d7452e22ead8cc35229acd 100644
--- a/src/state/ruinastateupdater.hh
+++ b/src/state/ruinastateupdater.hh
@@ -7,31 +7,32 @@
template <class ScalarVector, class Vector>
class RuinaStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- RuinaStateUpdater(ScalarVector _logState_initial,
- Dune::BitSetVector<1> const &_nodes, double _L);
+ RuinaStateUpdater(ScalarVector _alpha_initial,
+ Dune::BitSetVector<1> const &_nodes, double _L, double _V0);
void virtual nextTimeStep();
void virtual setup(double _tau);
void virtual solve(Vector const &velocity_field);
- void virtual extractLogState(ScalarVector &);
+ void virtual extractAlpha(ScalarVector &);
private:
- ScalarVector logState_o;
- ScalarVector logState;
+ ScalarVector alpha_o;
+ ScalarVector alpha;
Dune::BitSetVector<1> const &nodes;
double const L;
+ double const V0;
double tau;
};
template <class ScalarVector, class Vector>
RuinaStateUpdater<ScalarVector, Vector>::RuinaStateUpdater(
- ScalarVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
- double _L)
- : logState(_logState_initial), nodes(_nodes), L(_L) {}
+ ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L,
+ double _V0)
+ : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector>
void RuinaStateUpdater<ScalarVector, Vector>::nextTimeStep() {
- logState_o = logState;
+ alpha_o = alpha;
}
template <class ScalarVector, class Vector>
@@ -46,18 +47,17 @@ void RuinaStateUpdater<ScalarVector, Vector>::solve(
if (not toBool(nodes[i]))
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);
+ double const V = velocity_field[i].two_norm();
+ double const mtVoL = -tau * V / L;
+ alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) +
+ alpha_o[i] * std::exp(mtVoL);
}
}
template <class ScalarVector, class Vector>
-void RuinaStateUpdater<ScalarVector, Vector>::extractLogState(
- ScalarVector &_logState) {
- _logState = logState;
+void RuinaStateUpdater<ScalarVector, Vector>::extractAlpha(
+ ScalarVector &_alpha) {
+ _alpha = alpha;
}
#endif
diff --git a/src/state/stateupdater.hh b/src/state/stateupdater.hh
index 9ff937564e9d63c4213adb021c71a48f51725c8b..0ebf4774c11bf7825954ec4006a0b870645e7241 100644
--- a/src/state/stateupdater.hh
+++ b/src/state/stateupdater.hh
@@ -6,7 +6,7 @@ template <class ScalarVector, class Vector> class StateUpdater {
void virtual nextTimeStep() = 0;
void virtual setup(double _tau) = 0;
void virtual solve(Vector const &velocity_field) = 0;
- void virtual extractLogState(ScalarVector &logState) = 0;
+ void virtual extractAlpha(ScalarVector &alpha) = 0;
};
#endif
diff --git a/src/vtk.cc b/src/vtk.cc
index 0155384e070b215d12aa28de23e842c2715cf53c..6755b8b7b122f1312ab1372f10788e7ec9d2cc9d 100644
--- a/src/vtk.cc
+++ b/src/vtk.cc
@@ -32,10 +32,10 @@ void MyVTKWriter<VertexBasis, CellBasis>::write(
vertexBasis, v, "velocity");
writer.addVertexData(velocityPointer);
- auto const logStatePointer =
+ auto const AlphaPointer =
std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
- vertexBasis, alpha, "logState");
- writer.addVertexData(logStatePointer);
+ vertexBasis, alpha, "Alpha");
+ writer.addVertexData(AlphaPointer);
auto const stressPointer =
std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(