Prepare for Ruina nonlinearity; u/L -> uol

src/compute_state.cc

@@ -10,7 +10,8 @@
 
 #include "compute_state.hh"
 
-// The nonlinearity exp(-x)
+// The nonlinearity exp(-x) [ Note: the Dieterich law has an additional linear
+// term! ]
 class DecayingExponential {
 public:
   typedef Dune::FieldVector<double, 1> VectorType;
@@ -33,22 +34,23 @@ class DecayingExponential {
   double const h;
 };
 
-double compute_state_update_bisection(double h, double unorm, double L,
-                                      double old_state) {
+double state_update_dieterich_bisection(double h, double uol,
+                                        double old_state) {
   DecayingExponential::VectorType const start(0);
   DecayingExponential::VectorType const direction(1);
   DecayingExponential const phi(h);
   MyDirectionalConvexFunction<DecayingExponential> const J(
-      1.0, old_state - unorm / L, phi, start, direction);
+      1.0, old_state - uol, 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 compute_state_update(double h, double unorm, double L,
-                            double old_state) {
-  double ret1 = compute_state_update_bisection(h, unorm, L, old_state);
-  assert(std::abs(ret1 - old_state + unorm / L - h * std::exp(-ret1)) < 1e-12);
+double state_update_dieterich(double h, double uol, // unorm / L
+                              double old_state) {
+  double const ret1 = state_update_dieterich_bisection(h, uol, old_state);
+  assert(std::abs(ret1 - old_state + uol - h * std::exp(-ret1)) < 1e-12);
 
   return ret1;
 }

src/compute_state.hh

 #ifndef COMPUTE_STATE_HH
 #define COMPUTE_STATE_HH
 
-double compute_state_update(double h, double unorm, double L, double old_state);
+double state_update_dieterich(double h, double uol, double old_state);
 
 #endif

src/one-body-sample.cc

@@ -277,7 +277,7 @@ int main(int argc, char *argv[]) {
               // velocity
               // std::cout << std::log(L/unorm * h) << std::endl;
 
-              (*s4_new)[i] = compute_state_update(h, unorm, L, s4_old[i]);
+              (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]);
             }
           }
           if (parset.get<bool>("printProgress")) {