diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index c888bebeb1..98c6509f61 100644 --- a/cpp/memilio/config.h +++ b/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 +struct active_type; +} // namespace internal +} // namespace ad + namespace mio { /** @@ -67,6 +78,15 @@ struct Limits { return 1e-12; } }; + +template +struct Limits> { + /// @brief Returns the limit under which an AD value may be rounded down to zero. + static constexpr FP zero_tolerance() + { + return Limits::zero_tolerance(); + } +}; /** @} */ } // namespace mio diff --git a/cpp/memilio/math/integrator.h b/cpp/memilio/math/integrator.h index d5de016007..94214f529f 100644 --- a/cpp/memilio/math/integrator.h +++ b/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 using DerivFunction = std::function> y, FP t, Eigen::Ref> 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 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(dt, dt_copy, Limits::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,