17 #include "Time/TimeSteppers/TimeStepper.hpp"
27 template <
typename Vars,
typename DerivVars>
60 using options = tmpl::list<>;
62 "The standard Dormand-Prince 5th-order time stepper."};
71 template <
typename Vars,
typename DerivVars>
74 const TimeDelta& time_step)
const noexcept;
76 template <
typename Vars,
typename ErrVars,
typename DerivVars>
79 const TimeDelta& time_step)
const noexcept;
81 template <
typename Vars,
typename DerivVars>
84 double time)
const noexcept;
86 size_t order()
const noexcept
override;
88 size_t error_estimate_order()
const noexcept
override;
90 uint64_t number_of_substeps()
const noexcept
override;
92 uint64_t number_of_substeps_for_error()
const noexcept
override;
94 size_t number_of_past_steps()
const noexcept
override;
96 double stable_step()
const noexcept
override;
99 const TimeDelta& time_step)
const noexcept
override;
103 const TimeDelta& time_step)
const noexcept
override;
105 template <
typename Vars,
typename DerivVars>
106 bool can_change_step_size(
110 return time_id.substep() == 0;
118 void pup(PUP::er& p) noexcept
override {
119 TimeStepper::Inherit::pup(p);
125 static constexpr
double a2_ = 0.2;
128 {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0}};
130 {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0}};
132 46732.0 / 5247.0, 49.0 / 176.0,
135 125.0 / 192.0, -2187.0 / 6784.0,
140 {5179.0 / 57600.0, 0.0, 7571.0 / 16695.0, 393.0 / 640.0,
141 -92097.0 / 339200.0, 187.0 / 2100.0, 1.0 / 40.0}};
148 {-12715105075.0 / 11282082432.0, 87487479700.0 / 32700410799.0,
149 -10690763975.0 / 1880347072.0, 701980252875.0 / 199316789632.0,
150 -1453857185.0 / 822651844.0, 69997945.0 / 29380423.0}};
158 inline bool constexpr operator!=(
const DormandPrince5& ,
159 const DormandPrince5& ) noexcept {
163 template <
typename Vars,
typename DerivVars>
164 void DormandPrince5::update_u(
167 const TimeDelta& time_step)
const noexcept {
168 ASSERT(history->integration_order() == 5,
169 "Fixed-order stepper cannot run at order "
170 << history->integration_order());
171 const size_t substep = history->size() - 1;
172 const auto& u0 = history->begin().value();
173 const double dt = time_step.value();
175 const auto increment_u = [&u, &history, &dt](
const auto& coeffs) noexcept {
176 for (
size_t i = 0; i < coeffs.size(); ++i) {
177 *u += (
gsl::at(coeffs, i) * dt) *
178 (history->begin() +
static_cast<int>(i)).derivative();
183 *u = (history->end() - 1).value() +
184 (a2_ * dt) * history->begin().derivative();
185 }
else if (substep < 6) {
189 }
else if (substep == 2) {
191 }
else if (substep == 3) {
193 }
else if (substep == 4) {
199 ERROR(
"Substep in DP5 should be one of 0,1,2,3,4,5, not " << substep);
203 if (history->size() == number_of_substeps()) {
204 history->mark_unneeded(history->end());
208 template <
typename Vars,
typename ErrVars,
typename DerivVars>
212 const TimeDelta& time_step)
const noexcept {
213 ASSERT(history->integration_order() == 5,
214 "Fixed-order stepper cannot run at order "
215 << history->integration_order());
216 const size_t substep = history->size() - 1;
217 const auto& u0 = history->begin().value();
218 const double dt = time_step.value();
220 const auto increment_u = [&history, &dt](
const auto& coeffs,
221 auto local_u) noexcept {
222 for (
size_t i = 0; i < coeffs.size(); ++i) {
223 *local_u += (
gsl::at(coeffs, i) * dt) *
224 (history->begin() +
static_cast<int>(i)).derivative();
228 *u = (history->end() - 1).value() +
229 (a2_ * dt) * history->begin().derivative();
230 }
else if (substep < 7) {
234 }
else if (substep == 2) {
236 }
else if (substep == 3) {
238 }
else if (substep == 4) {
240 }
else if (substep == 5) {
245 increment_u(b_alt_, u_error);
246 *u_error = *u - *u_error;
249 ERROR(
"Substep in adaptive DP5 should be one of 0,1,2,3,4,5,6, not "
254 if (history->size() == number_of_substeps_for_error()) {
255 history->mark_unneeded(history->end());
260 template <
typename Vars,
typename DerivVars>
262 const History<Vars, DerivVars>& history,
263 const double time)
const noexcept {
264 ASSERT(history.size() == number_of_substeps(),
265 "DP5 can only dense output on last substep ("
266 << number_of_substeps() - 1 <<
"), not substep "
267 << history.size() - 1);
268 const double t0 = history[0].value();
269 const double t_end = history[history.size() - 1].value();
273 const double dt = t_end - t0;
274 const double output_fraction = (
time - t0) / dt;
275 ASSERT(output_fraction >= 0.0,
"Attempting dense output at time "
276 << time <<
", but already progressed past "
278 ASSERT(output_fraction <= 1.0,
"Requested time ("
279 << time <<
") not within step [" << t0
280 <<
", " << t0 + dt <<
"]");
283 const auto& u0 = history.begin().value();
285 for (
size_t i = 0; i < 6; ++i) {
287 (history.begin() +
static_cast<int>(i)).derivative();
292 const auto& l1 = history.begin().derivative();
293 const auto& l3 = (history.begin() + 2).derivative();
294 const auto& l4 = (history.begin() + 3).derivative();
295 const auto& l5 = (history.begin() + 4).derivative();
296 const auto& l6 = (history.begin() + 5).derivative();
300 const Vars rcont2 = u1 - u0;
301 const Vars rcont3 = dt * l1 - rcont2;
306 *u = u0 + output_fraction *
308 (1.0 - output_fraction) *
309 (rcont3 + output_fraction *
310 ((rcont2 - dt * l6 - rcont3) +
311 ((1.0 - output_fraction) * dt) *
312 (d_[0] * l1 + d_[1] * l3 + d_[2] * l4 +
313 d_[3] * l5 + (d_[4] + d_[5]) * l6))));