949 Correct handling of dt_min in integrators

cpp/benchmarks/graph_simulation.cpp

@@ -169,8 +169,8 @@ BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::EulerIntegratorCore)->Name("Graph Si
 BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::RKIntegratorCore)->Name("Graph Simulation - adapt_rk");
 BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
     ->Name("Graph Simulation - rk_ck54 (boost)");
-BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
-    ->Name("Graph Simulation - rk_dopri5 (boost)");
+// BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
+// ->Name("Graph Simulation - rk_dopri5 (boost)"); // TODO: reenable once boost bug is fixed
 BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
     ->Name("Graph Simulation - rkf78 (boost)");
 // run all benchmarks

cpp/benchmarks/integrator_step.config

@@ -11,7 +11,9 @@
                 6377.873644,
                 35.24915600,
                 30.02961100,
+                0.000000000,
                 182.1458650,
+                0.000000000,
                 66.15305900,
                 79.53062100,
                 3069.383604,

cpp/benchmarks/integrator_step.cpp

@@ -58,8 +58,8 @@ BENCHMARK_TEMPLATE(integrator_step, mio::RKIntegratorCore)->Name("Dummy 3/3");
 BENCHMARK_TEMPLATE(integrator_step, mio::RKIntegratorCore)->Name("simulate SecirModel adapt_rk");
 BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
     ->Name("simulate SecirModel boost rk_ck54");
-BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
-    ->Name("simulate SecirModel boost rk_dopri5");
+// BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
+// ->Name("simulate SecirModel boost rk_dopri5"); // TODO: reenable once boost bug is fixed
 BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
     ->Name("simulate SecirModel boost rkf78");
 // run all benchmarks

cpp/benchmarks/simulation.cpp

@@ -49,8 +49,8 @@ BENCHMARK_TEMPLATE(simulation, mio::RKIntegratorCore)->Name("Dummy 3/3");
 BENCHMARK_TEMPLATE(simulation, mio::RKIntegratorCore)->Name("simulate SecirModel adapt_rk");
 BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
     ->Name("simulate SecirModel boost rk_ck54");
-BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
-    ->Name("simulate SecirModel boost rk_dopri5");
+// BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
+//     ->Name("simulate SecirModel boost rk_dopri5"); // TODO: reenable once boost bug is fixed
 BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
     ->Name("simulate SecirModel boost rkf78");
 // run all benchmarks

cpp/memilio/math/adapt_rk.cpp

@@ -18,6 +18,7 @@
 * limitations under the License.
 */
 #include "memilio/math/adapt_rk.h"
+#include "memilio/utils/logging.h"
 
 namespace mio
 {
@@ -73,11 +74,20 @@ Tableau::Tableau()
 bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::VectorXd> yt, double& t, double& dt,
                             Eigen::Ref<Eigen::VectorXd> ytp1) const
 {
+    assert(0 <= m_dt_min);
+    assert(m_dt_min <= m_dt_max);
+
+    if (dt < m_dt_min || dt > m_dt_max) {
+        mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min, m_dt_max);
+    }
+
+    dt = std::min(dt, m_dt_max);
+
     double t_eval; // shifted time for evaluating yt
     double dt_new; // updated dt
 
-    bool converged              = false; // carry for convergence criterion
-    bool failed_step_size_adapt = false;
+    bool converged     = false; // carry for convergence criterion
+    bool dt_is_invalid = false;
 
     if (m_yt_eval.size() != yt.size()) {
         m_yt_eval.resize(yt.size());
@@ -86,7 +96,11 @@ bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::Vect
 
     m_yt_eval = yt;
 
-    while (!converged && !failed_step_size_adapt) {
+    while (!converged && !dt_is_invalid) {
+        if (dt < m_dt_min) {
+            dt_is_invalid = true;
+            dt            = m_dt_min;
+        }
         // compute first column of kt, i.e. kt_0 for each y in yt_eval
         f(m_yt_eval, t, m_kt_values.col(0));
 
@@ -113,7 +127,7 @@ bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::Vect
 
         converged = (m_error_estimate <= m_eps).all(); // convergence criterion
 
-        if (converged) {
+        if (converged || dt_is_invalid) {
             // if sufficiently exact, return ytp1, which currently contains the lower order approximation
             // (higher order is not always higher accuracy)
             t += dt; // this is the t where ytp1 belongs to
@@ -128,14 +142,12 @@ bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::Vect
         // and to avoid dt_new -> dt for step decreases when |error_estimate - eps| -> 0
         dt_new *= 0.9;
         // check if updated dt stays within desired bounds and update dt for next step
-        if (m_dt_min < dt_new) {
-            dt = std::min(dt_new, m_dt_max);
-        }
-        else {
-            failed_step_size_adapt = true;
-        }
+        dt = std::min(dt_new, m_dt_max);
     }
-    return !failed_step_size_adapt;
+    dt = std::max(dt, m_dt_min);
+    // return 'converged' in favor of '!dt_is_invalid', as these values only differ if step sizing failed,
+    // but the step with size dt_min was accepted.
+    return converged;
 }
 
 } // namespace mio

cpp/memilio/math/integrator.cpp

@@ -58,7 +58,7 @@ Eigen::Ref<Eigen::VectorXd> OdeIntegrator::advance(const DerivFunction& f, const
     }
 
     if (!step_okay) {
-        log_warning("Adaptive step sizing failed.");
+        log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
     }
     else if (std::abs((tmax - t) / (tmax - t0)) > 1e-14) {
         log_warning("Last time step too small. Could not reach tmax exactly.");

cpp/memilio/math/integrator.h

@@ -22,11 +22,8 @@
 
 #include "memilio/utils/time_series.h"
 
-#include "memilio/math/eigen.h"
 #include <memory>
-#include <vector>
 #include <functional>
-#include <algorithm>
 
 namespace mio
 {
@@ -43,13 +40,27 @@ class IntegratorCore
     virtual ~IntegratorCore(){};
 
     /**
-     * @brief Step of the integration with possible adaptive with
+     * @brief Make a single integration step.
      *
-     * @param[in] f Right hand side of ODE
-     * @param[in] yt value of y at t, y(t)
-     * @param[in,out] t current time step h=dt
-     * @param[in,out] dt current time step h=dt
-     * @param[out] ytp1 approximated value y(t+1)
+     * The behaviour of this method changes when the integration scheme has adaptive step sizing. 
+     * These changes are noted in the parentheses (...) below.
+     * Adaptive integrators must have bounds dt_min and dt_max for dt.
+     * The adaptive step sizing is considered to be successful, if a step of at least size dt_min sufficed tolerances.
+     * Tolerances are defined in each implementation, usually using a criterion with absolute and relative tolerances.
+     * Even if the step sizing failed, the integrator will make a step of at least size dt_min.
+     *
+     * @param[in] f Right hand side of the ODE. May be called multiple times with different arguments.
+     * @param[in] yt The known value of y at time t.
+     * @param[in,out] t The current time. It will be increased by dt.
+     *     (If adaptive, the increment is instead within [dt_min, dt].)
+     * @param[in,out] dt The current step size h=dt. Will not be changed.
+     *     (If adaptive, the given dt is used as the maximum step size, and must be within [dt_min, dt_max].
+     *      During integration, dt is adjusted in [dt_min, dt_max] to have an optimal size for the next step.)
+     * @param[out] ytp1 Set to the approximated value of y at time t + dt.
+     *     (If adaptive, this time may be smaller, but it is at least t + dt_min, at most t + dt_max.
+     *      Note that the increment on t may be different from the returned value of dt.)
+     * @return Always true for nonadaptive methods.
+     *     (If adaptive, returns whether the adaptive step sizing was successful.)
      */
     virtual bool step(const DerivFunction& f, Eigen::Ref<const Eigen::VectorXd> yt, double& t, double& dt,
                       Eigen::Ref<Eigen::VectorXd> ytp1) const = 0;

cpp/memilio/math/stepper_wrapper.h

@@ -22,6 +22,7 @@
 
 #include "memilio/utils/compiler_diagnostics.h"
 #include "memilio/math/integrator.h"
+#include "memilio/utils/logging.h"
 
 GCC_CLANG_DIAGNOSTIC(push)
 GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow")
@@ -32,7 +33,7 @@ MSVC_WARNING_DISABLE_PUSH(4127)
 #include "boost/numeric/odeint/stepper/runge_kutta4.hpp"
 #include "boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp"
 #include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp"
-#include "boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp"
+// #include "boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp" // TODO: reenable once boost bug is fixed
 MSVC_WARNING_POP()
 GCC_CLANG_DIAGNOSTIC(pop)
 
@@ -50,6 +51,12 @@ template <template <class State = Eigen::VectorXd, class Value = double, class D
           class ControlledStepper>
 class ControlledStepperWrapper : public mio::IntegratorCore
 {
+    using Stepper                         = boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>>;
+    static constexpr bool is_fsal_stepper = std::is_same_v<typename Stepper::stepper_type::stepper_category,
+                                                           boost::numeric::odeint::explicit_error_stepper_fsal_tag>;
+    static_assert(!is_fsal_stepper,
+                  "FSAL steppers cannot be used until https://github.com/boostorg/odeint/issues/72 is resolved.");
+
 public:
     /**
      * @brief Set up the integrator
@@ -70,32 +77,67 @@ class ControlledStepperWrapper : public mio::IntegratorCore
     }
 
     /**
-    * @brief Make a single integration step of a system of ODEs and adapt step width
-    * @param[in] yt value of y at t, y(t)
-    * @param[in,out] t current time step h=dt
-    * @param[in,out] dt current time step h=dt
-    * @param[out] ytp1 approximated value y(t+1)
+    * @brief Make a single integration step on a system of ODEs and adapt the step size dt.
+
+    * @param[in] yt Value of y at t_{k}, y(t_{k}).
+    * @param[in,out] t Current time step t_{k} for some k. Will be set to t_{k+1} in [t_{k} + dt_min, t_{k} + dt].
+    * @param[in,out] dt Current time step size h=dt. Overwritten by an estimated optimal step size for the next step.
+    * @param[out] ytp1 The approximated value of y(t_{k+1}).
     */
     bool step(const mio::DerivFunction& f, Eigen::Ref<Eigen::VectorXd const> yt, double& t, double& dt,
               Eigen::Ref<Eigen::VectorXd> ytp1) const override
     {
-        // copy y(t) to dydt, to retrieve the VectorXd from the Ref
-        dydt               = yt;
-        const double t_old = t; // t is updated by try_step on a successfull step
-        do {
-            // we use the scheme try_step(sys, inout, t, dt) with sys=f, inout=y(t) for
-            // in-place computation. This is similiar to do_step, but it can update t and dt
-            m_stepper.try_step(
-                // reorder arguments of the DerivFunction f for the stepper
-                [&](const Eigen::VectorXd& x, Eigen::VectorXd& dxds, double s) {
-                    dxds.resizeLike(x); // try_step calls sys with a vector of size 0 for some reason
-                    f(x, s, dxds);
-                },
-                dydt, t, dt);
-            // stop on a successfull step or a failed step size adaption (w.r.t. the minimal step size)
-        } while (t == t_old && dt > m_dt_min);
-        ytp1 = dydt; // output new y(t)
-        return dt > m_dt_min;
+        using boost::numeric::odeint::fail;
+        assert(0 <= m_dt_min);
+        assert(m_dt_min <= m_dt_max);
+
+        if (dt < m_dt_min || dt > m_dt_max) {
+            mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min,
+                             m_dt_max);
+        }
+        // set initial values for exit conditions
+        auto step_result = fail;
+        bool is_dt_valid = true;
+        // copy vectors from the references, since the stepper cannot (trivially) handle Eigen::Ref
+        m_ytp1 = ytp1;
+        m_yt   = yt;
+        // make a integration step, adapting dt to a possibly larger value on success,
+        // or a strictly smaller value on fail.
+        // stop only on a successful step or a failed step size adaption (w.r.t. the minimal step size m_dt_min)
+        while (step_result == fail && is_dt_valid) {
+            if (dt < m_dt_min) {
+                is_dt_valid = false;
+                dt          = m_dt_min;
+            }
+            // we use the scheme try_step(sys, in, t, out, dt) with sys=f, in=y(t_{k}), out=y(t_{k+1}).
+            // this is similiar to do_step, but it can adapt the step size dt. If successful, it also updates t.
+
+            if constexpr (!is_fsal_stepper) { // prevent compile time errors with fsal steppers
+                step_result = m_stepper.try_step(
+                    // reorder arguments of the DerivFunction f for the stepper
+                    [&](const Eigen::VectorXd& x, Eigen::VectorXd& dxds, double s) {
+                        dxds.resizeLike(x); // boost resizers cannot resize Eigen::Vector, hence we need to do that here
+                        f(x, s, dxds);
+                    },
+                    m_yt, t, m_ytp1, dt);
+            }
+        }
+        // output the new value by copying it back to the reference
+        ytp1 = m_ytp1;
+        // bound dt from below
+        // the last adaptive step (successful or not) may have calculated a new step size smaller than m_dt_min
+
+        dt = std::max(dt, m_dt_min);
+        // check whether the last step failed (which means that m_dt_min was still too large to suffice tolerances)
+        if (step_result == fail) {
+            // adaptive stepping failed, but we still return the result of the last attempt
+            t += m_dt_min;
+            return false;
+        }
+        else {
+            // successfully made an integration step
+            return true;
+        }
     }
 
     /// @param tol the required absolute tolerance for comparison of the iterative approximation
@@ -126,21 +168,18 @@ class ControlledStepperWrapper : public mio::IntegratorCore
     }
 
 private:
-    boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>> create_stepper()
+    /// @brief (Re)initialize the internal stepper.
+    Stepper create_stepper()
     {
         // for more options see: boost/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
-        return boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>>(
-            boost::numeric::odeint::default_error_checker<typename ControlledStepper<>::value_type,
-                                                          typename ControlledStepper<>::algebra_type,
-                                                          typename ControlledStepper<>::operations_type>(m_abs_tol,
-                                                                                                         m_rel_tol),
-            boost::numeric::odeint::default_step_adjuster<typename ControlledStepper<>::value_type,
-                                                          typename ControlledStepper<>::time_type>(m_dt_max));
+        return Stepper(typename Stepper::error_checker_type(m_abs_tol, m_rel_tol),
+                       typename Stepper::step_adjuster_type(m_dt_max));
     }
 
-    double m_abs_tol, m_rel_tol, m_dt_min, m_dt_max; // integrator parameters
-    mutable Eigen::VectorXd dydt;
-    mutable boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>> m_stepper;
+    double m_abs_tol, m_rel_tol; ///< Absolute and relative tolerances for integration.
+    double m_dt_min, m_dt_max; ///< Lower and upper bound to the step size dt.
+    mutable Eigen::VectorXd m_ytp1, m_yt; ///< Temporary storage to avoid allocations in step function.
+    mutable Stepper m_stepper; ///< A stepper instance used for integration.
 };
 
 } // namespace mio