diff --git a/src/compute_state.cc b/src/compute_state.cc index b928618bb57f0ec37161daae3edeb957b2c03932..e4b6d215b79815a462a3a3fb326038bb9e1078ec 100644 --- a/src/compute_state.cc +++ b/src/compute_state.cc @@ -34,31 +34,31 @@ class DecayingExponential { double const tau; }; -double state_update_dieterich_bisection(double tau, double uol, +double state_update_dieterich_bisection(double tau, double VoL, double old_state) { DecayingExponential::VectorType const start(0); DecayingExponential::VectorType const direction(1); DecayingExponential const phi(tau); MyDirectionalConvexFunction<DecayingExponential> const J( - 1.0, old_state - uol, phi, start, direction); + 1.0, old_state - tau * VoL, phi, start, direction); int bisectionsteps = 0; Bisection const bisection(0.0, 1.0, 1e-12, true, 0); // TODO return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO } -double state_update_dieterich(double tau, double uol, double old_state) { - double const ret = state_update_dieterich_bisection(tau, uol, old_state); +double state_update_dieterich(double tau, double VoL, double old_state) { + double const ret = state_update_dieterich_bisection(tau, VoL, old_state); /* We have - ret - old_state + uol = tau*exp(-ret) + ret - old_state + tau * VoL = tau*exp(-ret) or - log((ret - old_state + uol)/tau) = -ret + log((ret - old_state)/tau + VoL) = -ret */ - assert(std::min(std::abs(ret - old_state + uol - tau * std::exp(-ret)), - std::abs(std::log((ret - old_state + uol) / tau) + ret)) < + assert(std::min(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))), + std::abs(std::log((ret - old_state) / tau + VoL) + ret)) < 1e-12); return ret; } diff --git a/src/compute_state.hh b/src/compute_state.hh index 3fedc4fe7953074d006e0a3afef772dfd7672329..49fe3dcc3d7bc112b31d0d8d2193a12811b5121a 100644 --- a/src/compute_state.hh +++ b/src/compute_state.hh @@ -1,5 +1,5 @@ #ifndef COMPUTE_STATE_HH #define COMPUTE_STATE_HH -double state_update_dieterich(double h, double uol, double old_state); +double state_update_dieterich(double h, double VoL, double old_state); #endif diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 596615b16c33d62875475b63eefe44d555cefb72..448f849e654e58594567e384fc8f40744c722a6a 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -366,20 +366,19 @@ int main(int argc, char *argv[]) { // Update the state for (size_t i = 0; i < frictionalNodes.size(); ++i) { if (frictionalNodes[i][0]) { - double const unorm = ud[i].two_norm() * tau; + double const V = ud[i].two_norm(); // // the (logarithmic) steady state corresponding to the // // current velocity - // std::cout << std::log(L/unorm * tau) << std::endl; + // std::cout << std::log(L/V) << std::endl; switch (parset.get<Config::state_model>( "boundary.friction.state.model")) { case Config::Dieterich: - alpha[i] = - state_update_dieterich(tau, unorm / L, alpha_old[i]); + alpha[i] = state_update_dieterich(tau, V / L, alpha_old[i]); break; case Config::Ruina: - alpha[i] = state_update_ruina(tau, unorm / L, alpha_old[i]); + alpha[i] = state_update_ruina(tau, V * tau / L, alpha_old[i]); break; } }