DormandPrince5.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines class DormandPrince5.
6 
7 #pragma once
8 
9 #include <cstddef>
10 #include <cstdint>
11 #include <ostream>
12 #include <pup.h>
13 
14 #include "Options/Options.hpp"
16 #include "Time/Time.hpp"
17 #include "Time/TimeSteppers/TimeStepper.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/TMPL.hpp"
23 
24 /// \cond
25 struct TimeStepId;
26 namespace TimeSteppers {
27 template <typename Vars, typename DerivVars>
28 class History;
29 } // namespace TimeSteppers
30 /// \endcond
31 
32 namespace TimeSteppers {
33 
34 /*!
35  * \ingroup TimeSteppersGroup
36  *
37  * The standard 5th-order Dormand-Prince time stepping method, given e.g. in
38  * Sec. 7.2 of \cite NumericalRecipes.
39  *
40  * \f{eqnarray}{
41  * \frac{du}{dt} & = & \mathcal{L}(t,u).
42  * \f}
43  * Given a solution \f$u(t^n)=u^n\f$, this stepper computes
44  * \f$u(t^{n+1})=u^{n+1}\f$ using the following equations:
45  *
46  * \f{align}{
47  * k^{(1)} & = dt \mathcal{L}(t^n, u^n),\\
48  * k^{(i)} & = dt \mathcal{L}(t^n + c_i dt,
49  * u^n + \sum_{j=1}^{i-1} a_{ij} k^{(j)}),
50  * \mbox{ } 2 \leq i \leq 6,\\
51  * u^{n+1} & = u^n + \sum_{i=1}^{6} b_i k^{(i)}.
52  * \f}
53  *
54  * Here the coefficients \f$a_{ij}\f$, \f$b_i\f$, and \f$c_i\f$ are given
55  * in e.g. Sec. 7.2 of \cite NumericalRecipes. Note that \f$c_1 = 0\f$.
56  *
57  */
58 class DormandPrince5 : public TimeStepper::Inherit {
59  public:
60  using options = tmpl::list<>;
61  static constexpr Options::String help = {
62  "The standard Dormand-Prince 5th-order time stepper."};
63 
64  DormandPrince5() = default;
65  DormandPrince5(const DormandPrince5&) noexcept = default;
66  DormandPrince5& operator=(const DormandPrince5&) noexcept = default;
67  DormandPrince5(DormandPrince5&&) noexcept = default;
68  DormandPrince5& operator=(DormandPrince5&&) noexcept = default;
69  ~DormandPrince5() noexcept override = default;
70 
71  template <typename Vars, typename DerivVars>
72  void update_u(gsl::not_null<Vars*> u,
74  const TimeDelta& time_step) const noexcept;
75 
76  template <typename Vars, typename ErrVars, typename DerivVars>
77  bool update_u(gsl::not_null<Vars*> u, gsl::not_null<ErrVars*> u_error,
79  const TimeDelta& time_step) const noexcept;
80 
81  template <typename Vars, typename DerivVars>
82  void dense_update_u(gsl::not_null<Vars*> u,
83  const History<Vars, DerivVars>& history,
84  double time) const noexcept;
85 
86  size_t order() const noexcept override;
87 
88  size_t error_estimate_order() const noexcept override;
89 
90  uint64_t number_of_substeps() const noexcept override;
91 
92  uint64_t number_of_substeps_for_error() const noexcept override;
93 
94  size_t number_of_past_steps() const noexcept override;
95 
96  double stable_step() const noexcept override;
97 
98  TimeStepId next_time_id(const TimeStepId& current_id,
99  const TimeDelta& time_step) const noexcept override;
100 
101  TimeStepId next_time_id_for_error(
102  const TimeStepId& current_id,
103  const TimeDelta& time_step) const noexcept override;
104 
105  template <typename Vars, typename DerivVars>
106  bool can_change_step_size(
107  const TimeStepId& time_id,
108  const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
109  noexcept {
110  return time_id.substep() == 0;
111  }
112 
113  WRAPPED_PUPable_decl_template(DormandPrince5); // NOLINT
114 
115  explicit DormandPrince5(CkMigrateMessage* /*unused*/) noexcept {}
116 
117  // clang-tidy: do not pass by non-const reference
118  void pup(PUP::er& p) noexcept override { // NOLINT
119  TimeStepper::Inherit::pup(p);
120  }
121 
122  private:
123  // Coefficients from the Dormand-Prince 5 Butcher tableau (e.g. Sec. 7.2
124  // of \cite NumericalRecipes).
125  static constexpr double a2_ = 0.2;
126  static constexpr std::array<double, 2> a3_{{3.0 / 40.0, 9.0 / 40.0}};
127  static constexpr std::array<double, 3> a4_{
128  {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0}};
129  static constexpr std::array<double, 4> a5_{
130  {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0}};
131  static constexpr std::array<double, 5> a6_{{9017.0 / 3168.0, -355.0 / 33.0,
132  46732.0 / 5247.0, 49.0 / 176.0,
133  -5103.0 / 18656.0}};
134  static constexpr std::array<double, 6> b_{{35.0 / 384.0, 0.0, 500.0 / 1113.0,
135  125.0 / 192.0, -2187.0 / 6784.0,
136  11.0 / 84.0}};
137  // Coefficients for the embedded method, for generating an error measure
138  // during adaptive stepping (e.g. Sec. 7.2 of \cite NumericalRecipes).
139  static constexpr std::array<double, 7> b_alt_{
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}};
142 
143  static const std::array<Time::rational_t, 6> c_;
144 
145  // Coefficients for dense output, taken from Sec. 7.2 of
146  // \cite NumericalRecipes
147  static constexpr std::array<double, 6> d_{
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}};
151 };
152 
153 inline bool constexpr operator==(const DormandPrince5& /*lhs*/,
154  const DormandPrince5& /*rhs*/) noexcept {
155  return true;
156 }
157 
158 inline bool constexpr operator!=(const DormandPrince5& /*lhs*/,
159  const DormandPrince5& /*rhs*/) noexcept {
160  return false;
161 }
162 
163 template <typename Vars, typename DerivVars>
164 void DormandPrince5::update_u(
165  const gsl::not_null<Vars*> u,
166  const gsl::not_null<History<Vars, DerivVars>*> history,
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();
174 
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();
179  }
180  };
181 
182  if (substep == 0) {
183  *u = (history->end() - 1).value() +
184  (a2_ * dt) * history->begin().derivative();
185  } else if (substep < 6) {
186  *u = u0;
187  if (substep == 1) {
188  increment_u(a3_);
189  } else if (substep == 2) {
190  increment_u(a4_);
191  } else if (substep == 3) {
192  increment_u(a5_);
193  } else if (substep == 4) {
194  increment_u(a6_);
195  } else {
196  increment_u(b_);
197  }
198  } else {
199  ERROR("Substep in DP5 should be one of 0,1,2,3,4,5, not " << substep);
200  }
201 
202  // Clean up old history
203  if (history->size() == number_of_substeps()) {
204  history->mark_unneeded(history->end());
205  }
206 }
207 
208 template <typename Vars, typename ErrVars, typename DerivVars>
209 bool DormandPrince5::update_u(gsl::not_null<Vars*> u,
210  gsl::not_null<ErrVars*> u_error,
211  gsl::not_null<History<Vars, DerivVars>*> history,
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();
219 
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();
225  }
226  };
227  if (substep == 0) {
228  *u = (history->end() - 1).value() +
229  (a2_ * dt) * history->begin().derivative();
230  } else if (substep < 7) {
231  *u = u0;
232  if (substep == 1) {
233  increment_u(a3_, u);
234  } else if (substep == 2) {
235  increment_u(a4_, u);
236  } else if (substep == 3) {
237  increment_u(a5_, u);
238  } else if (substep == 4) {
239  increment_u(a6_, u);
240  } else if (substep == 5) {
241  increment_u(b_, u);
242  } else {
243  increment_u(b_, u);
244  *u_error = u0;
245  increment_u(b_alt_, u_error);
246  *u_error = *u - *u_error;
247  }
248  } else {
249  ERROR("Substep in adaptive DP5 should be one of 0,1,2,3,4,5,6, not "
250  << substep);
251  }
252 
253  // Clean up old history
254  if (history->size() == number_of_substeps_for_error()) {
255  history->mark_unneeded(history->end());
256  }
257  return substep == 6;
258 }
259 
260 template <typename Vars, typename DerivVars>
261 void DormandPrince5::dense_update_u(const gsl::not_null<Vars*> u,
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();
270  // The history does not contain the final step; specifically,
271  // step_end = t0 + c[4] * dt, so dt = (t_end - t0) / c[4]. But since
272  // c[4] = 1.0 for DP5, we don't need to divide by c[4] here.
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 "
277  << t0);
278  ASSERT(output_fraction <= 1.0, "Requested time ("
279  << time << ") not within step [" << t0
280  << ", " << t0 + dt << "]");
281 
282  // Get the solution u1 at the final time
283  const auto& u0 = history.begin().value();
284  Vars u1 = u0;
285  for (size_t i = 0; i < 6; ++i) {
286  u1 += (gsl::at(b_, i) * dt) *
287  (history.begin() + static_cast<int>(i)).derivative();
288  }
289 
290  // We need the following: k1, k3, k4, k5, k6.
291  // Here k1 = dt * l1, k3 = dt * l3, etc.
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();
297 
298  // Compute the updating coefficents, called rcontN in Numerical recipes,
299  // that will be reused, so I don't have to compute them more than once.
300  const Vars rcont2 = u1 - u0;
301  const Vars rcont3 = dt * l1 - rcont2;
302 
303  // The formula for dense output is given in Numerical Recipes Sec. 17.2.3.
304  // Note: L(t+dt, u1) after the step is unavailable in the history; so here,
305  // approximate L(t+dt, u1) by l6.
306  *u = u0 + output_fraction *
307  (rcont2 +
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))));
314 }
315 } // namespace TimeSteppers
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
CharmPupable.hpp
Options.hpp
Error.hpp
TimeSteppers::DormandPrince5
Definition: DormandPrince5.hpp:58
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
cstddef
Assert.hpp
std::array< double, 2 >
TimeDelta
Definition: Time.hpp:88
cstdint
ConstantExpressions.hpp
Time.hpp
TimeStepId
Definition: TimeStepId.hpp:25
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
TimeSteppers::History
Definition: History.hpp:33
TimeSteppers
Definition: BoundaryHistory.hpp:23
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
std::time
T time(T... args)
ostream
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13