1038 odeintegrator can run indefinitely

cpp/memilio/math/adapt_rk.h

@@ -208,19 +208,22 @@ class RKIntegratorCore : public IntegratorCore<FP>
      * @param[in,out] t current time
      * @param[in,out] dt current time step size h=dt
      * @param[out] ytp1 approximated value y(t+1)
+     * @param[in] force_step_size Forces a step with exactly size dt.
      */
     bool step(const DerivFunction<FP>& f, Eigen::Ref<Eigen::VectorXd const> yt, double& t, double& dt,
-              Eigen::Ref<Eigen::VectorXd> ytp1) const override
+              Eigen::Ref<Eigen::VectorXd> ytp1, const bool force_step_size = false) const override
     {
         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);
-        }
+        if (!force_step_size) {
+            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);
+            dt = std::min(dt, m_dt_max);
+        }
 
         double t_eval; // shifted time for evaluating yt
         double dt_new; // updated dt
@@ -236,7 +239,7 @@ class RKIntegratorCore : public IntegratorCore<FP>
         m_yt_eval = yt;
 
         while (!converged && !dt_is_invalid) {
-            if (dt < m_dt_min) {
+            if (dt < m_dt_min && !force_step_size) {
                 dt_is_invalid = true;
                 dt            = m_dt_min;
             }
@@ -264,7 +267,7 @@ class RKIntegratorCore : public IntegratorCore<FP>
             // calculate mixed tolerance
             m_eps = m_abs_tol + ytp1.array().abs() * m_rel_tol;
 
-            converged = (m_error_estimate <= m_eps).all(); // convergence criterion
+            converged = (m_error_estimate <= m_eps).all() || force_step_size; // convergence criterion
 
             if (converged || dt_is_invalid) {
                 // if sufficiently exact, return ytp1, which currently contains the lower order approximation

cpp/memilio/math/euler.h

@@ -22,6 +22,7 @@
 
 #include "memilio/config.h"
 #include "memilio/math/integrator.h"
+#include "memilio/utils/compiler_diagnostics.h"
 
 namespace mio
 {
@@ -41,10 +42,12 @@ class EulerIntegratorCore : public IntegratorCore<FP>
      * @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)
+     * @param[in] force_step_size This argument is unused, as the step width is always fixed.
      */
-    bool step(const DerivFunction<FP>& f, Eigen::Ref<const Vector<FP>> yt, FP& t, FP& dt,
-              Eigen::Ref<Vector<FP>> ytp1) const override
+    bool step(const DerivFunction<FP>& f, Eigen::Ref<const Vector<FP>> yt, FP& t, FP& dt, Eigen::Ref<Vector<FP>> ytp1,
+              const bool force_step_size = false) const override
     {
+        mio::unused(force_step_size);
         // we are misusing the next step y as temporary space to store the derivative
         f(yt, t, ytp1);
         ytp1 = yt + dt * ytp1;

cpp/memilio/math/integrator.h

@@ -60,11 +60,14 @@ class IntegratorCore
      * @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.
+     * @param[in] force_step_size This value may be ignored by non adaptive integrators. False by default.
+     *     (If adaptive, setting this to true allows making integration steps with any dt, ignoring the restriction to
+     *      [dt_min, dt_max]. This is meant to enable making a final integration step that is smaller than dt_min.)
+     * @return Always true for non adaptive methods.
      *     (If adaptive, returns whether the adaptive step sizing was successful.)
      */
     virtual bool step(const DerivFunction<FP>& f, Eigen::Ref<const Vector<FP>> yt, FP& t, FP& dt,
-                      Eigen::Ref<Vector<FP>> ytp1) const = 0;
+                      Eigen::Ref<Vector<FP>> ytp1, const bool force_step_size = false) const = 0;
 };
 
 /**
@@ -109,7 +112,7 @@ class OdeIntegrator
 
         results.reserve(results.get_num_time_points() + num_steps);
 
-        bool step_okay = true;
+        bool is_step_sizing_okay = true;
 
         FP dt_copy; // used to check whether step sizing is adaptive
         FP dt_restore = 0; // used to restore dt if dt was decreased to reach tmax
@@ -127,19 +130,30 @@ class OdeIntegrator
             dt_copy = dt;
 
             results.add_time_point();
-            step_okay &= m_core->step(f, results[i], t, dt, results[i + 1]);
-            results.get_last_time() = t;
+            bool step_okay = m_core->step(f, results[i], t, dt, results[i + 1]);
 
             // if dt has been changed (even slighly) by step, register the current m_core as adaptive
             m_is_adaptive |= !floating_point_equal(dt, dt_copy);
+
+            // catch case where dt_min of an adaptive stepper is larger than tmax - t
+            if (t > tmax) {
+                log_info("Forcing last step size.");
+                // reset time and step size
+                dt = dt_copy;
+                t  = tmax - dt_copy;
+                // recalculate the erroneous last step and overwrite its results
+                step_okay = m_core->step(f, results[i], t, dt, results[i + 1], true);
+            }
+            is_step_sizing_okay &= step_okay;
+            results.get_last_time() = t;
         }
         // if dt was decreased to reach tmax in the last time iteration,
         // we restore it as it is now probably smaller than required for tolerances
         dt = max(dt, dt_restore);
 
         if (m_is_adaptive) {
-            if (!step_okay) {
-                log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
+            if (!is_step_sizing_okay) {
+                log_warning("Adaptive step sizing failed. Forced at least one integration step of size dt_min.");
             }
             else if (fabs((tmax - t) / (tmax - t0)) > 1e-14) {
                 log_warning("Last time step too small. Could not reach tmax exactly.");

cpp/memilio/math/stepper_wrapper.h

@@ -27,6 +27,8 @@
 #include "boost/numeric/odeint/stepper/controlled_runge_kutta.hpp"
 #include "boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp"
 #include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp"
+
+#include <limits>
 // #include "boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp" // TODO: reenable once boost bug is fixed
 
 namespace mio
@@ -70,24 +72,27 @@ class ControlledStepperWrapper : public mio::IntegratorCore<FP>
     }
 
     /**
-    * @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}).
-    */
+     * @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}).
+     * @param[in] force_step_size Forces a step with exactly size dt.
+     */
     bool step(const mio::DerivFunction<FP>& f, Eigen::Ref<Vector<FP> const> yt, FP& t, FP& dt,
-              Eigen::Ref<Vector<FP>> ytp1) const override
+              Eigen::Ref<Vector<FP>> ytp1, const bool force_step_size = false) const override
     {
         using boost::numeric::odeint::fail;
         using std::max;
         assert(0 <= m_dt_min);
         assert(m_dt_min <= m_dt_max);
 
-        if (dt < m_dt_min || dt > m_dt_max) {
+        const FP dt_forced = dt;
+        if ((dt < m_dt_min || dt > m_dt_max) && !force_step_size) {
             mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min,
                              m_dt_max);
+            dt = std::min(dt, m_dt_max);
         }
         // set initial values for exit conditions
         auto step_result = fail;
@@ -98,14 +103,14 @@ class ControlledStepperWrapper : public mio::IntegratorCore<FP>
         // 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) {
+        do {
+            if (dt < m_dt_min && !force_step_size) {
                 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.
-
+            // note that the step_adjuster_type never allows or returns dt >= dt_max
             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
@@ -115,18 +120,18 @@ class ControlledStepperWrapper : public mio::IntegratorCore<FP>
                     },
                     m_yt, t, m_ytp1, dt);
             }
-        }
+        } while (step_result == fail && is_dt_valid && !force_step_size);
         // 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 = max(dt, m_dt_min);
+        // clamp dt (again)
+        // the last adaptive step (successful or not) may have calculated a new step size smaller than m_dt_min,
+        // and when forcing step size, dt may have been larger than m_dt_max
+        dt = std::clamp(dt, m_dt_min, m_dt_max);
         // 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;
+            t += force_step_size ? dt_forced : dt;
+            return force_step_size;
         }
         else {
             // successfully made an integration step
@@ -157,8 +162,7 @@ class ControlledStepperWrapper : public mio::IntegratorCore<FP>
     /// @param dt_max sets the maximum step size
     void set_dt_max(FP dt_max)
     {
-        m_dt_max  = dt_max;
-        m_stepper = create_stepper();
+        m_dt_max = dt_max;
     }
 
 private:
@@ -167,7 +171,7 @@ class ControlledStepperWrapper : public mio::IntegratorCore<FP>
     {
         // for more options see: boost/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
         return Stepper(typename Stepper::error_checker_type(m_abs_tol, m_rel_tol),
-                       typename Stepper::step_adjuster_type(m_dt_max));
+                       typename Stepper::step_adjuster_type(std::numeric_limits<FP>::max())); // handle dt_max manually
     }
 
     FP m_abs_tol, m_rel_tol; ///< Absolute and relative tolerances for integration.