1154 integrator is handled as adaptive if using small fixed step sizes (#1155)

Add a tolerance in OdeIntegrator when detecting whether the current IntegratorCore is adaptive. 

Co-authored-by: reneSchm <[email protected]>

Author: lenaploetzke

cpp/memilio/config.h

@@ -30,6 +30,17 @@
 
 using ScalarType = double;
 
+namespace ad
+{
+namespace internal
+{
+// Forward declaration used for support of ad types in Limits.
+// These templates are called AD_TAPE_REAL and DATA_HANDLER internally in AD.
+template <class FP, class DataHandler>
+struct active_type;
+} // namespace internal
+} // namespace ad
+
 namespace mio
 {
 /**
@@ -67,6 +78,15 @@ struct Limits<double> {
         return 1e-12;
     }
 };
+
+template <class FP, class DataHandler>
+struct Limits<ad::internal::active_type<FP, DataHandler>> {
+    /// @brief Returns the limit under which an AD value may be rounded down to zero.
+    static constexpr FP zero_tolerance()
+    {
+        return Limits<FP>::zero_tolerance();
+    }
+};
 /** @} */
 
 } // namespace mio

cpp/memilio/math/integrator.h

@@ -20,6 +20,7 @@
 #ifndef INTEGRATOR_H
 #define INTEGRATOR_H
 
+#include "memilio/config.h"
 #include "memilio/math/floating_point.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/utils/logging.h"
@@ -30,7 +31,7 @@ namespace mio
 {
 
 /**
- * Function template to be integrated
+ * Function template to be integrated.
  */
 template <typename FP = double>
 using DerivFunction = std::function<void(Eigen::Ref<const mio::Vector<FP>> y, FP t, Eigen::Ref<mio::Vector<FP>> dydt)>;
@@ -118,8 +119,8 @@ class IntegratorCore
 };
 
 /**
- * Integrate initial value problems (IVP) of ordinary differential equations (ODE) of the form y' = f(y, t), y(t0) = y0.
- * tparam FP a floating point type accepted by Eigen
+ * @brief Integrate initial value problems (IVP) of ordinary differential equations (ODE) of the form y' = f(y, t), y(t0) = y0.
+ * @tparam FP a floating point type accepted by Eigen
  */
 template <typename FP = double>
 class OdeIntegrator
@@ -141,7 +142,7 @@ class OdeIntegrator
      * @param[in] tmax Time end point. Must be greater than results.get_last_time().
      * @param[in, out] dt Initial integration step size. May be changed by the IntegratorCore.
      * @param[in, out] results List of results. Must contain at least one time point. The last entry is used as
-     * intitial time and value. A new entry is added for each integration step.
+     * initial time and value. A new entry is added for each integration step.
      * @return A reference to the last value in the results time series.
      */
 
@@ -168,9 +169,9 @@ class OdeIntegrator
         FP t              = t0;
 
         for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > 1e-10; ++i) {
-            //we don't make timesteps too small as the error estimator of an adaptive integrator
+            // We don't make time steps too small as the error estimator of an adaptive integrator
             //may not be able to handle it. this is very conservative and maybe unnecessary,
-            //but also unlikely to happen. may need to be reevaluated
+            //but also unlikely to happen. may need to be reevaluated.
 
             if (dt > tmax - t) {
                 dt_restore = dt;
@@ -187,8 +188,8 @@ class OdeIntegrator
             step_okay &= m_core->step(f, results[i], t, dt, results[i + 1]);
             results.get_last_time() = t;
 
-            // 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);
+            // if dt has been changed by step, register the current m_core as adaptive.
+            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, Limits<FP>::zero_tolerance());
         }
         m_core->get_dt_min() = dt_min_restore; // restore dt_min
         // if dt was decreased to reach tmax in the last time iteration,