Rescale functional; Stricter assertions

src/compute_state.cc

@@ -15,9 +15,11 @@ class DecayingExponential {
   typedef Dune::FieldVector<double, 1> VectorType;
   typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
 
+  DecayingExponential(double h) : h(h) {}
+
   void directionalSubDiff(VectorType const &u, VectorType const &v,
                           Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-std::exp(-u[0]));
+    D[0] = D[1] = v[0] * (-h * std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
@@ -25,12 +27,15 @@ class DecayingExponential {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
+
+private:
+  double const h;
 };
 
 double compute_state_update_bisection(double h, double unorm, double L,
                                       double old_state) {
   MyDirectionalConvexFunction<DecayingExponential> const J(
-      1.0 / h, (old_state - unorm / L) / h, DecayingExponential(), 0, 1);
+      1.0, old_state - unorm / L, DecayingExponential(h), 0, 1);
   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
@@ -39,7 +44,7 @@ double compute_state_update_bisection(double h, double unorm, double L,
 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-10);
+  assert(std::abs(ret1 - old_state + unorm / L - h * std::exp(-ret1)) < 1e-12);
 
   return ret1;
 }