Skip to content
Snippets Groups Projects
Commit 50b8e604 authored by Elias Pipping's avatar Elias Pipping
Browse files

[Cleanup] Call state logState everywhere

parent 2bb07e1f
Branches
No related tags found
No related merge requests found
Showing with 41 additions and 42 deletions
......@@ -24,7 +24,7 @@ class FrictionPotentialWrapper {
DUNE_THROW(NotImplemented, "evaluation not implemented");
}
void virtual updateState(double) = 0;
void virtual updateLogState(double) = 0;
};
// phi(V) = V log(V/V_m) - V + V_m if V >= V_m
......@@ -67,10 +67,9 @@ class FrictionPotential : public FrictionPotentialWrapper {
// V_m = V_0 exp(-K/a),
// 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
void updateLogState(double logState) {
logV_m = std::log(fd.V0) +
(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
(-(fd.mu0 + fd.b * (logState + std::log(fd.V0 / fd.L))) / fd.a);
V_cutoff = std::exp(logV_m + fd.mumin / fd.a);
}
......@@ -91,7 +90,7 @@ class TrivialFunction : public FrictionPotentialWrapper {
double regularity(double) const { return 0; }
void updateState(double) {}
void updateLogState(double) {}
};
}
#endif
......@@ -82,7 +82,7 @@ class GlobalNonlinearity {
return res->regularity(x);
}
virtual void updateState(SingletonVectorType const &state) = 0;
virtual void updateLogState(SingletonVectorType const &logState) = 0;
};
}
#endif
......@@ -42,9 +42,9 @@ class GlobalRuinaNonlinearity
}
}
void updateState(SingletonVectorType const &state) override {
void updateLogState(SingletonVectorType const &logState) override {
for (size_t i = 0; i < restrictions.size(); ++i)
restrictions[i]->updateState(state[i]);
restrictions[i]->updateLogState(logState[i]);
}
/*
......
......@@ -43,7 +43,7 @@ template <size_t dimension> class LocalFriction {
return func->evaluate(x.two_norm());
}
void updateState(double state) { func->updateState(state); }
void updateLogState(double logState) { func->updateLogState(logState); }
double regularity(VectorType const &x) const {
double const xnorm = x.two_norm();
......
......@@ -4,9 +4,9 @@
#include "friction_writer.hh"
double computeCOF(FrictionData const &fd, double V, double log_state) {
double computeCOF(FrictionData const &fd, double V, double logState) {
double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
fd.b * (log_state + std::log(fd.V0 / fd.L));
fd.b * (logState + std::log(fd.V0 / fd.L));
return (mu > 0) ? mu : 0;
}
......
......@@ -411,7 +411,7 @@ int main(int argc, char *argv[]) {
accelerationRHS = 0.0;
Arithmetic::addProduct(accelerationRHS, A, u_initial);
Arithmetic::addProduct(accelerationRHS, C, v_initial);
myGlobalNonlinearity->updateState(alpha_initial);
myGlobalNonlinearity->updateLogState(alpha_initial);
// NOTE: We assume differentiability of Psi at v0 here!
myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
accelerationRHS -= ell;
......@@ -511,7 +511,7 @@ int main(int argc, char *argv[]) {
size_t iterationCounter;
auto solveVelocityProblem = [&](VectorType &_velocityIterate,
SingletonVectorType const &_alpha) {
myGlobalNonlinearity->updateState(_alpha);
myGlobalNonlinearity->updateLogState(_alpha);
// NIT: Do we really need to pass u here?
typename NonlinearFactoryType::ConvexProblemType const convexProblem(
......@@ -535,7 +535,7 @@ int main(int argc, char *argv[]) {
double lastStateCorrection;
for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
stateUpdater->solve(v);
stateUpdater->extractState(alpha);
stateUpdater->extractLogState(alpha);
if (stateFPI == 1)
relaxationWriter << "N ";
......
......@@ -60,12 +60,12 @@ class DieterichNonlinearity {
}
};
double state_update_dieterich_euler(double tau, double VoL, double old_state) {
double state_update_dieterich_euler(double tau, double VoL, double logState_o) {
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);
1.0 / tau, logState_o / tau - VoL, phi, start, direction);
int bisectionsteps = 0;
Bisection const bisection;
......
#ifndef COMPUTE_STATE_DIETERICH_EULER_HH
#define COMPUTE_STATE_DIETERICH_EULER_HH
double state_update_dieterich_euler(double h, double VoL, double old_state);
double state_update_dieterich_euler(double h, double VoL, double logState_o);
#endif
......@@ -4,10 +4,10 @@
#include "compute_state_ruina.hh"
double state_update_ruina(double tau, double uol, double old_state) {
double state_update_ruina(double tau, double uol, double logState_o) {
if (uol == 0)
return old_state;
return logState_o;
double const ret = old_state - uol * std::log(uol / tau);
double const ret = logState_o - uol * std::log(uol / tau);
return ret / (1 + uol);
}
#ifndef COMPUTE_STATE_RUINA_HH
#define COMPUTE_STATE_RUINA_HH
double state_update_ruina(double h, double uol, double old_state);
double state_update_ruina(double h, double uol, double logState_o);
#endif
......@@ -14,10 +14,10 @@ class DieterichStateUpdater
virtual void nextTimeStep();
virtual void setup(double _tau);
virtual void solve(VectorType const &velocity_field);
virtual void extractState(SingletonVectorType &state);
virtual void extractLogState(SingletonVectorType &);
private:
SingletonVectorType alpha_o;
SingletonVectorType logState_o;
SingletonVectorType alpha;
Dune::BitSetVector<1> const &nodes;
double const L;
......@@ -32,7 +32,7 @@ DieterichStateUpdater<SingletonVectorType, VectorType>::DieterichStateUpdater(
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
alpha_o = alpha;
logState_o = alpha;
}
template <class SingletonVectorType, class VectorType>
......@@ -47,14 +47,14 @@ void DieterichStateUpdater<SingletonVectorType, VectorType>::solve(
for (size_t i = 0; i < nodes.size(); ++i)
if (nodes[i][0]) {
double const V = velocity_field[i].two_norm();
alpha[i] = state_update_dieterich_euler(tau, V / L, alpha_o[i]);
alpha[i] = state_update_dieterich_euler(tau, V / L, logState_o[i]);
}
}
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::extractState(
SingletonVectorType &state) {
state = alpha;
void DieterichStateUpdater<SingletonVectorType, VectorType>::extractLogState(
SingletonVectorType &logState) {
logState = alpha;
}
#endif
......@@ -13,10 +13,10 @@ class RuinaStateUpdater : public StateUpdater<SingletonVectorType, VectorType> {
virtual void nextTimeStep();
virtual void setup(double _tau);
virtual void solve(VectorType const &velocity_field);
virtual void extractState(SingletonVectorType &state);
virtual void extractLogState(SingletonVectorType &);
private:
SingletonVectorType alpha_o;
SingletonVectorType logState_o;
SingletonVectorType alpha;
Dune::BitSetVector<1> const &nodes;
double const L;
......@@ -31,7 +31,7 @@ RuinaStateUpdater<SingletonVectorType, VectorType>::RuinaStateUpdater(
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
alpha_o = alpha;
logState_o = alpha;
}
template <class SingletonVectorType, class VectorType>
......@@ -45,14 +45,14 @@ void RuinaStateUpdater<SingletonVectorType, VectorType>::solve(
for (size_t i = 0; i < nodes.size(); ++i)
if (nodes[i][0]) {
double const V = velocity_field[i].two_norm();
alpha[i] = state_update_ruina(tau, V * tau / L, alpha_o[i]);
alpha[i] = state_update_ruina(tau, V * tau / L, logState_o[i]);
}
}
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::extractState(
SingletonVectorType &state) {
state = alpha;
void RuinaStateUpdater<SingletonVectorType, VectorType>::extractLogState(
SingletonVectorType &logState) {
logState = alpha;
}
#endif
......@@ -6,7 +6,7 @@ template <class SingletonVectorType, class VectorType> class StateUpdater {
virtual void nextTimeStep() = 0;
virtual void setup(double _tau) = 0;
virtual void solve(VectorType const &velocity_field) = 0;
virtual void extractState(SingletonVectorType &state) = 0;
virtual void extractLogState(SingletonVectorType &logState) = 0;
};
#endif
......@@ -7,7 +7,7 @@
template <class VertexBasis, class CellBasis, class VectorType,
class SingletonVectorType>
void writeVtk(VertexBasis const &vertexBasis, VectorType const &displacement,
SingletonVectorType const &state, CellBasis const &cellBasis,
SingletonVectorType const &logState, CellBasis const &cellBasis,
SingletonVectorType const &stress, std::string const &filename) {
Dune::VTKWriter<typename VertexBasis::GridView> writer(
vertexBasis.getGridView());
......@@ -17,10 +17,10 @@ void writeVtk(VertexBasis const &vertexBasis, VectorType const &displacement,
vertexBasis, displacement, "displacement");
writer.addVertexData(displacementPointer);
auto const statePointer = Dune::make_shared<
auto const logStatePointer = Dune::make_shared<
VTKBasisGridFunction<VertexBasis, SingletonVectorType> const>(
vertexBasis, state, "state");
writer.addVertexData(statePointer);
vertexBasis, logState, "logState");
writer.addVertexData(logStatePointer);
auto const vonmisesPointer = Dune::make_shared<
VTKBasisGridFunction<CellBasis, SingletonVectorType> const>(
......
......@@ -9,7 +9,7 @@
template <class VertexBasis, class CellBasis, class VectorType,
class SingletonVectorType>
void writeVtk(VertexBasis const &vertexBasis, VectorType const &displacement,
SingletonVectorType const &state, CellBasis const &cellBasis,
SingletonVectorType const &logState, CellBasis const &cellBasis,
SingletonVectorType const &stress, std::string const &filename);
#endif
......@@ -28,5 +28,5 @@ using MyP0Basis = P0Basis<GridView, double>;
template void writeVtk<P1Basis, MyP0Basis, VectorType, SingletonVectorType>(
P1Basis const &vertexBasis, VectorType const &displacement,
SingletonVectorType const &state, MyP0Basis const &cellBasis,
SingletonVectorType const &logState, MyP0Basis const &cellBasis,
SingletonVectorType const &stress, std::string const &filename);
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment