unorm/tau -> VoL

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;
 }

src/compute_state.hh

 #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

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;
               }
             }