From 247e997b88e4f25ad0b5b67605dc577389d34ed0 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Fri, 6 Dec 2024 11:39:46 +0100 Subject: [PATCH 1/4] added tolerance (and fixed some typos) --- cpp/memilio/config.h | 10 ++++++++++ cpp/memilio/math/integrator.h | 17 +++++++++-------- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index c888bebeb1..17d1dbf25e 100644 --- a/cpp/memilio/config.h +++ b/cpp/memilio/config.h @@ -27,6 +27,7 @@ #include "memilio/config_internal.h" #include +#include "ad/ad.hpp" using ScalarType = double; @@ -67,6 +68,15 @@ struct Limits { return 1e-12; } }; + +template +struct Limits> { + /// @brief Returns the limit under which an ad::internal::active_type may be rounded down to zero. + static constexpr float 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, From db18780681a90bc0d3ad464e1dc1fe134a00216a Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Fri, 6 Dec 2024 12:33:31 +0100 Subject: [PATCH 2/4] bug fix --- cpp/memilio/config.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index 17d1dbf25e..97676ac1e2 100644 --- a/cpp/memilio/config.h +++ b/cpp/memilio/config.h @@ -72,7 +72,7 @@ struct Limits { template struct Limits> { /// @brief Returns the limit under which an ad::internal::active_type may be rounded down to zero. - static constexpr float zero_tolerance() + static constexpr AD_TAPE_REAL zero_tolerance() { return Limits::zero_tolerance(); } From 343cbe9b3038323643e61c96afb8756341f14a90 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Thu, 12 Dec 2024 15:30:06 +0100 Subject: [PATCH 3/4] Review Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com> --- cpp/memilio/config.h | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index 97676ac1e2..966878d6d8 100644 --- a/cpp/memilio/config.h +++ b/cpp/memilio/config.h @@ -27,10 +27,19 @@ #include "memilio/config_internal.h" #include -#include "ad/ad.hpp" using ScalarType = double; +namespace ad +{ +namespace internal +{ +// Forward declaration used for support of ad types in Limits. +template +struct active_type; +} // namespace internal +} // namespace ad + namespace mio { /** @@ -69,12 +78,12 @@ struct Limits { } }; -template -struct Limits> { - /// @brief Returns the limit under which an ad::internal::active_type may be rounded down to zero. - static constexpr AD_TAPE_REAL zero_tolerance() +template +struct Limits> { + /// @brief Returns the limit under which an AD value may be rounded down to zero. + static constexpr ADTapeReal zero_tolerance() { - return Limits::zero_tolerance(); + return Limits::zero_tolerance(); } }; /** @} */ From d7d2bf70cb7bfef88f63297dc4ad6598867baf46 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Thu, 12 Dec 2024 16:58:16 +0100 Subject: [PATCH 4/4] small corrections --- cpp/memilio/config.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/cpp/memilio/config.h b/cpp/memilio/config.h index 966878d6d8..98c6509f61 100644 --- a/cpp/memilio/config.h +++ b/cpp/memilio/config.h @@ -35,7 +35,8 @@ namespace ad namespace internal { // Forward declaration used for support of ad types in Limits. -template +// These templates are called AD_TAPE_REAL and DATA_HANDLER internally in AD. +template struct active_type; } // namespace internal } // namespace ad @@ -78,12 +79,12 @@ struct Limits { } }; -template -struct Limits> { +template +struct Limits> { /// @brief Returns the limit under which an AD value may be rounded down to zero. - static constexpr ADTapeReal zero_tolerance() + static constexpr FP zero_tolerance() { - return Limits::zero_tolerance(); + return Limits::zero_tolerance(); } }; /** @} */