Skip to content

Commit 69a6b68

Browse files
reneSchmmknaranjaHenrZujubicker
authored
949 Correct handling of dt_min in integrators (#960)
- Enforce time step of (adaptive) ODE integrators to be in [dtmin, dtmax]. - Temporarily deactivate FSAL stepper due to a bug in boost. - Document expected behavior of IntegratorCore::step implementations Co-authored-by: Martin J. Kühn <[email protected]> Co-authored-by: Henrik Zunker <[email protected]> Co-authored-by: Julia Bicker <[email protected]>
1 parent afd4982 commit 69a6b68

12 files changed

Lines changed: 318 additions & 173 deletions

‎cpp/benchmarks/graph_simulation.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -169,8 +169,8 @@ BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::EulerIntegratorCore)->Name("Graph Si
169169
BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::RKIntegratorCore)->Name("Graph Simulation - adapt_rk");
170170
BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
171171
->Name("Graph Simulation - rk_ck54 (boost)");
172-
BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
173-
->Name("Graph Simulation - rk_dopri5 (boost)");
172+
// BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
173+
// ->Name("Graph Simulation - rk_dopri5 (boost)"); // TODO: reenable once boost bug is fixed
174174
BENCHMARK_TEMPLATE(graph_sim_secirvvs, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
175175
->Name("Graph Simulation - rkf78 (boost)");
176176
// run all benchmarks

‎cpp/benchmarks/integrator_step.config‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,9 @@
1111
6377.873644,
1212
35.24915600,
1313
30.02961100,
14+
0.000000000,
1415
182.1458650,
16+
0.000000000,
1517
66.15305900,
1618
79.53062100,
1719
3069.383604,

‎cpp/benchmarks/integrator_step.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,8 @@ BENCHMARK_TEMPLATE(integrator_step, mio::RKIntegratorCore)->Name("Dummy 3/3");
5858
BENCHMARK_TEMPLATE(integrator_step, mio::RKIntegratorCore)->Name("simulate SecirModel adapt_rk");
5959
BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
6060
->Name("simulate SecirModel boost rk_ck54");
61-
BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
62-
->Name("simulate SecirModel boost rk_dopri5");
61+
// BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
62+
// ->Name("simulate SecirModel boost rk_dopri5"); // TODO: reenable once boost bug is fixed
6363
BENCHMARK_TEMPLATE(integrator_step, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
6464
->Name("simulate SecirModel boost rkf78");
6565
// run all benchmarks

‎cpp/benchmarks/simulation.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,8 @@ BENCHMARK_TEMPLATE(simulation, mio::RKIntegratorCore)->Name("Dummy 3/3");
4949
BENCHMARK_TEMPLATE(simulation, mio::RKIntegratorCore)->Name("simulate SecirModel adapt_rk");
5050
BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_cash_karp54>)
5151
->Name("simulate SecirModel boost rk_ck54");
52-
BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
53-
->Name("simulate SecirModel boost rk_dopri5");
52+
// BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_dopri5>)
53+
// ->Name("simulate SecirModel boost rk_dopri5"); // TODO: reenable once boost bug is fixed
5454
BENCHMARK_TEMPLATE(simulation, mio::ControlledStepperWrapper<boost::numeric::odeint::runge_kutta_fehlberg78>)
5555
->Name("simulate SecirModel boost rkf78");
5656
// run all benchmarks

‎cpp/memilio/math/adapt_rk.cpp‎

Lines changed: 23 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
* limitations under the License.
1919
*/
2020
#include "memilio/math/adapt_rk.h"
21+
#include "memilio/utils/logging.h"
2122

2223
namespace mio
2324
{
@@ -73,11 +74,20 @@ Tableau::Tableau()
7374
bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::VectorXd> yt, double& t, double& dt,
7475
Eigen::Ref<Eigen::VectorXd> ytp1) const
7576
{
77+
assert(0 <= m_dt_min);
78+
assert(m_dt_min <= m_dt_max);
79+
80+
if (dt < m_dt_min || dt > m_dt_max) {
81+
mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min, m_dt_max);
82+
}
83+
84+
dt = std::min(dt, m_dt_max);
85+
7686
double t_eval; // shifted time for evaluating yt
7787
double dt_new; // updated dt
7888

79-
bool converged = false; // carry for convergence criterion
80-
bool failed_step_size_adapt = false;
89+
bool converged = false; // carry for convergence criterion
90+
bool dt_is_invalid = false;
8191

8292
if (m_yt_eval.size() != yt.size()) {
8393
m_yt_eval.resize(yt.size());
@@ -86,7 +96,11 @@ bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::Vect
8696

8797
m_yt_eval = yt;
8898

89-
while (!converged && !failed_step_size_adapt) {
99+
while (!converged && !dt_is_invalid) {
100+
if (dt < m_dt_min) {
101+
dt_is_invalid = true;
102+
dt = m_dt_min;
103+
}
90104
// compute first column of kt, i.e. kt_0 for each y in yt_eval
91105
f(m_yt_eval, t, m_kt_values.col(0));
92106

@@ -113,7 +127,7 @@ bool RKIntegratorCore::step(const DerivFunction& f, Eigen::Ref<const Eigen::Vect
113127

114128
converged = (m_error_estimate <= m_eps).all(); // convergence criterion
115129

116-
if (converged) {
130+
if (converged || dt_is_invalid) {
117131
// if sufficiently exact, return ytp1, which currently contains the lower order approximation
118132
// (higher order is not always higher accuracy)
119133
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
128142
// and to avoid dt_new -> dt for step decreases when |error_estimate - eps| -> 0
129143
dt_new *= 0.9;
130144
// check if updated dt stays within desired bounds and update dt for next step
131-
if (m_dt_min < dt_new) {
132-
dt = std::min(dt_new, m_dt_max);
133-
}
134-
else {
135-
failed_step_size_adapt = true;
136-
}
145+
dt = std::min(dt_new, m_dt_max);
137146
}
138-
return !failed_step_size_adapt;
147+
dt = std::max(dt, m_dt_min);
148+
// return 'converged' in favor of '!dt_is_invalid', as these values only differ if step sizing failed,
149+
// but the step with size dt_min was accepted.
150+
return converged;
139151
}
140152

141153
} // namespace mio

‎cpp/memilio/math/integrator.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ Eigen::Ref<Eigen::VectorXd> OdeIntegrator::advance(const DerivFunction& f, const
5858
}
5959

6060
if (!step_okay) {
61-
log_warning("Adaptive step sizing failed.");
61+
log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
6262
}
6363
else if (std::abs((tmax - t) / (tmax - t0)) > 1e-14) {
6464
log_warning("Last time step too small. Could not reach tmax exactly.");

‎cpp/memilio/math/integrator.h‎

Lines changed: 20 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -22,11 +22,8 @@
2222

2323
#include "memilio/utils/time_series.h"
2424

25-
#include "memilio/math/eigen.h"
2625
#include <memory>
27-
#include <vector>
2826
#include <functional>
29-
#include <algorithm>
3027

3128
namespace mio
3229
{
@@ -43,13 +40,27 @@ class IntegratorCore
4340
virtual ~IntegratorCore(){};
4441

4542
/**
46-
* @brief Step of the integration with possible adaptive with
43+
* @brief Make a single integration step.
4744
*
48-
* @param[in] f Right hand side of ODE
49-
* @param[in] yt value of y at t, y(t)
50-
* @param[in,out] t current time step h=dt
51-
* @param[in,out] dt current time step h=dt
52-
* @param[out] ytp1 approximated value y(t+1)
45+
* The behaviour of this method changes when the integration scheme has adaptive step sizing.
46+
* These changes are noted in the parentheses (...) below.
47+
* Adaptive integrators must have bounds dt_min and dt_max for dt.
48+
* The adaptive step sizing is considered to be successful, if a step of at least size dt_min sufficed tolerances.
49+
* Tolerances are defined in each implementation, usually using a criterion with absolute and relative tolerances.
50+
* Even if the step sizing failed, the integrator will make a step of at least size dt_min.
51+
*
52+
* @param[in] f Right hand side of the ODE. May be called multiple times with different arguments.
53+
* @param[in] yt The known value of y at time t.
54+
* @param[in,out] t The current time. It will be increased by dt.
55+
* (If adaptive, the increment is instead within [dt_min, dt].)
56+
* @param[in,out] dt The current step size h=dt. Will not be changed.
57+
* (If adaptive, the given dt is used as the maximum step size, and must be within [dt_min, dt_max].
58+
* During integration, dt is adjusted in [dt_min, dt_max] to have an optimal size for the next step.)
59+
* @param[out] ytp1 Set to the approximated value of y at time t + dt.
60+
* (If adaptive, this time may be smaller, but it is at least t + dt_min, at most t + dt_max.
61+
* Note that the increment on t may be different from the returned value of dt.)
62+
* @return Always true for nonadaptive methods.
63+
* (If adaptive, returns whether the adaptive step sizing was successful.)
5364
*/
5465
virtual bool step(const DerivFunction& f, Eigen::Ref<const Eigen::VectorXd> yt, double& t, double& dt,
5566
Eigen::Ref<Eigen::VectorXd> ytp1) const = 0;

‎cpp/memilio/math/stepper_wrapper.h‎

Lines changed: 73 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222

2323
#include "memilio/utils/compiler_diagnostics.h"
2424
#include "memilio/math/integrator.h"
25+
#include "memilio/utils/logging.h"
2526

2627
GCC_CLANG_DIAGNOSTIC(push)
2728
GCC_CLANG_DIAGNOSTIC(ignored "-Wshadow")
@@ -32,7 +33,7 @@ MSVC_WARNING_DISABLE_PUSH(4127)
3233
#include "boost/numeric/odeint/stepper/runge_kutta4.hpp"
3334
#include "boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp"
3435
#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp"
35-
#include "boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp"
36+
// #include "boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp" // TODO: reenable once boost bug is fixed
3637
MSVC_WARNING_POP()
3738
GCC_CLANG_DIAGNOSTIC(pop)
3839

@@ -50,6 +51,12 @@ template <template <class State = Eigen::VectorXd, class Value = double, class D
5051
class ControlledStepper>
5152
class ControlledStepperWrapper : public mio::IntegratorCore
5253
{
54+
using Stepper = boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>>;
55+
static constexpr bool is_fsal_stepper = std::is_same_v<typename Stepper::stepper_type::stepper_category,
56+
boost::numeric::odeint::explicit_error_stepper_fsal_tag>;
57+
static_assert(!is_fsal_stepper,
58+
"FSAL steppers cannot be used until https://github.com/boostorg/odeint/issues/72 is resolved.");
59+
5360
public:
5461
/**
5562
* @brief Set up the integrator
@@ -70,32 +77,67 @@ class ControlledStepperWrapper : public mio::IntegratorCore
7077
}
7178

7279
/**
73-
* @brief Make a single integration step of a system of ODEs and adapt step width
74-
* @param[in] yt value of y at t, y(t)
75-
* @param[in,out] t current time step h=dt
76-
* @param[in,out] dt current time step h=dt
77-
* @param[out] ytp1 approximated value y(t+1)
80+
* @brief Make a single integration step on a system of ODEs and adapt the step size dt.
81+
82+
* @param[in] yt Value of y at t_{k}, y(t_{k}).
83+
* @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].
84+
* @param[in,out] dt Current time step size h=dt. Overwritten by an estimated optimal step size for the next step.
85+
* @param[out] ytp1 The approximated value of y(t_{k+1}).
7886
*/
7987
bool step(const mio::DerivFunction& f, Eigen::Ref<Eigen::VectorXd const> yt, double& t, double& dt,
8088
Eigen::Ref<Eigen::VectorXd> ytp1) const override
8189
{
82-
// copy y(t) to dydt, to retrieve the VectorXd from the Ref
83-
dydt = yt;
84-
const double t_old = t; // t is updated by try_step on a successfull step
85-
do {
86-
// we use the scheme try_step(sys, inout, t, dt) with sys=f, inout=y(t) for
87-
// in-place computation. This is similiar to do_step, but it can update t and dt
88-
m_stepper.try_step(
89-
// reorder arguments of the DerivFunction f for the stepper
90-
[&](const Eigen::VectorXd& x, Eigen::VectorXd& dxds, double s) {
91-
dxds.resizeLike(x); // try_step calls sys with a vector of size 0 for some reason
92-
f(x, s, dxds);
93-
},
94-
dydt, t, dt);
95-
// stop on a successfull step or a failed step size adaption (w.r.t. the minimal step size)
96-
} while (t == t_old && dt > m_dt_min);
97-
ytp1 = dydt; // output new y(t)
98-
return dt > m_dt_min;
90+
using boost::numeric::odeint::fail;
91+
assert(0 <= m_dt_min);
92+
assert(m_dt_min <= m_dt_max);
93+
94+
if (dt < m_dt_min || dt > m_dt_max) {
95+
mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min,
96+
m_dt_max);
97+
}
98+
// set initial values for exit conditions
99+
auto step_result = fail;
100+
bool is_dt_valid = true;
101+
// copy vectors from the references, since the stepper cannot (trivially) handle Eigen::Ref
102+
m_ytp1 = ytp1;
103+
m_yt = yt;
104+
// make a integration step, adapting dt to a possibly larger value on success,
105+
// or a strictly smaller value on fail.
106+
// stop only on a successful step or a failed step size adaption (w.r.t. the minimal step size m_dt_min)
107+
while (step_result == fail && is_dt_valid) {
108+
if (dt < m_dt_min) {
109+
is_dt_valid = false;
110+
dt = m_dt_min;
111+
}
112+
// we use the scheme try_step(sys, in, t, out, dt) with sys=f, in=y(t_{k}), out=y(t_{k+1}).
113+
// this is similiar to do_step, but it can adapt the step size dt. If successful, it also updates t.
114+
115+
if constexpr (!is_fsal_stepper) { // prevent compile time errors with fsal steppers
116+
step_result = m_stepper.try_step(
117+
// reorder arguments of the DerivFunction f for the stepper
118+
[&](const Eigen::VectorXd& x, Eigen::VectorXd& dxds, double s) {
119+
dxds.resizeLike(x); // boost resizers cannot resize Eigen::Vector, hence we need to do that here
120+
f(x, s, dxds);
121+
},
122+
m_yt, t, m_ytp1, dt);
123+
}
124+
}
125+
// output the new value by copying it back to the reference
126+
ytp1 = m_ytp1;
127+
// bound dt from below
128+
// the last adaptive step (successful or not) may have calculated a new step size smaller than m_dt_min
129+
130+
dt = std::max(dt, m_dt_min);
131+
// check whether the last step failed (which means that m_dt_min was still too large to suffice tolerances)
132+
if (step_result == fail) {
133+
// adaptive stepping failed, but we still return the result of the last attempt
134+
t += m_dt_min;
135+
return false;
136+
}
137+
else {
138+
// successfully made an integration step
139+
return true;
140+
}
99141
}
100142

101143
/// @param tol the required absolute tolerance for comparison of the iterative approximation
@@ -126,21 +168,18 @@ class ControlledStepperWrapper : public mio::IntegratorCore
126168
}
127169

128170
private:
129-
boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>> create_stepper()
171+
/// @brief (Re)initialize the internal stepper.
172+
Stepper create_stepper()
130173
{
131174
// for more options see: boost/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
132-
return boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>>(
133-
boost::numeric::odeint::default_error_checker<typename ControlledStepper<>::value_type,
134-
typename ControlledStepper<>::algebra_type,
135-
typename ControlledStepper<>::operations_type>(m_abs_tol,
136-
m_rel_tol),
137-
boost::numeric::odeint::default_step_adjuster<typename ControlledStepper<>::value_type,
138-
typename ControlledStepper<>::time_type>(m_dt_max));
175+
return Stepper(typename Stepper::error_checker_type(m_abs_tol, m_rel_tol),
176+
typename Stepper::step_adjuster_type(m_dt_max));
139177
}
140178

141-
double m_abs_tol, m_rel_tol, m_dt_min, m_dt_max; // integrator parameters
142-
mutable Eigen::VectorXd dydt;
143-
mutable boost::numeric::odeint::controlled_runge_kutta<ControlledStepper<>> m_stepper;
179+
double m_abs_tol, m_rel_tol; ///< Absolute and relative tolerances for integration.
180+
double m_dt_min, m_dt_max; ///< Lower and upper bound to the step size dt.
181+
mutable Eigen::VectorXd m_ytp1, m_yt; ///< Temporary storage to avoid allocations in step function.
182+
mutable Stepper m_stepper; ///< A stepper instance used for integration.
144183
};
145184

146185
} // namespace mio

0 commit comments

Comments
 (0)