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/EvolutionOrdering.hpp"
17 #include "Time/Time.hpp"
18 #include "Time/TimeSteppers/TimeStepper.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/TMPL.hpp"
24 
25 /// \cond
26 struct TimeStepId;
27 namespace TimeSteppers {
28 template <typename Vars, typename DerivVars>
29 class History;
30 } // namespace TimeSteppers
31 /// \endcond
32 
33 namespace TimeSteppers {
34 
35 /*!
36  * \ingroup TimeSteppersGroup
37  *
38  * The standard 5th-order Dormand-Prince time stepping method, given e.g. in
39  * Sec. 7.2 of \cite NumericalRecipes.
40  *
41  * \f{eqnarray}{
42  * \frac{du}{dt} & = & \mathcal{L}(t,u).
43  * \f}
44  * Given a solution \f$u(t^n)=u^n\f$, this stepper computes
45  * \f$u(t^{n+1})=u^{n+1}\f$ using the following equations:
46  *
47  * \f{align}{
48  * k^{(1)} & = dt \mathcal{L}(t^n, u^n),\\
49  * k^{(i)} & = dt \mathcal{L}(t^n + c_i dt,
50  * u^n + \sum_{j=1}^{i-1} a_{ij} k^{(j)}),
51  * \mbox{ } 2 \leq i \leq 6,\\
52  * u^{n+1} & = u^n + \sum_{i=1}^{6} b_i k^{(i)}.
53  * \f}
54  *
55  * Here the coefficients \f$a_{ij}\f$, \f$b_i\f$, and \f$c_i\f$ are given
56  * in e.g. Sec. 7.2 of \cite NumericalRecipes. Note that \f$c_1 = 0\f$.
57  *
58  */
59 class DormandPrince5 : public TimeStepper::Inherit {
60  public:
61  using options = tmpl::list<>;
62  static constexpr Options::String help = {
63  "The standard Dormand-Prince 5th-order time stepper."};
64 
65  DormandPrince5() = default;
66  DormandPrince5(const DormandPrince5&) noexcept = default;
67  DormandPrince5& operator=(const DormandPrince5&) noexcept = default;
68  DormandPrince5(DormandPrince5&&) noexcept = default;
69  DormandPrince5& operator=(DormandPrince5&&) noexcept = default;
70  ~DormandPrince5() noexcept override = default;
71 
72  template <typename Vars, typename DerivVars>
73  void update_u(gsl::not_null<Vars*> u,
75  const TimeDelta& time_step) const noexcept;
76 
77  template <typename Vars, typename ErrVars, typename DerivVars>
78  bool update_u(gsl::not_null<Vars*> u, gsl::not_null<ErrVars*> u_error,
80  const TimeDelta& time_step) const noexcept;
81 
82  template <typename Vars, typename DerivVars>
83  bool dense_update_u(gsl::not_null<Vars*> u,
84  const History<Vars, DerivVars>& history,
85  double time) const noexcept;
86 
87  size_t order() const noexcept override;
88 
89  size_t error_estimate_order() const noexcept override;
90 
91  uint64_t number_of_substeps() const noexcept override;
92 
93  uint64_t number_of_substeps_for_error() const noexcept override;
94 
95  size_t number_of_past_steps() const noexcept override;
96 
97  double stable_step() const noexcept override;
98 
99  TimeStepId next_time_id(const TimeStepId& current_id,
100  const TimeDelta& time_step) const noexcept override;
101 
102  TimeStepId next_time_id_for_error(
103  const TimeStepId& current_id,
104  const TimeDelta& time_step) const noexcept override;
105 
106  template <typename Vars, typename DerivVars>
107  bool can_change_step_size(
108  const TimeStepId& time_id,
109  const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
110  noexcept {
111  return time_id.substep() == 0;
112  }
113 
114  WRAPPED_PUPable_decl_template(DormandPrince5); // NOLINT
115 
116  explicit DormandPrince5(CkMigrateMessage* /*unused*/) noexcept {}
117 
118  // clang-tidy: do not pass by non-const reference
119  void pup(PUP::er& p) noexcept override { // NOLINT
120  TimeStepper::Inherit::pup(p);
121  }
122 
123  private:
124  // Coefficients from the Dormand-Prince 5 Butcher tableau (e.g. Sec. 7.2
125  // of \cite NumericalRecipes).
126  static constexpr double a2_ = 0.2;
127  static constexpr std::array<double, 2> a3_{{3.0 / 40.0, 9.0 / 40.0}};
128  static constexpr std::array<double, 3> a4_{
129  {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0}};
130  static constexpr std::array<double, 4> a5_{
131  {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0}};
132  static constexpr std::array<double, 5> a6_{{9017.0 / 3168.0, -355.0 / 33.0,
133  46732.0 / 5247.0, 49.0 / 176.0,
134  -5103.0 / 18656.0}};
135  static constexpr std::array<double, 6> b_{{35.0 / 384.0, 0.0, 500.0 / 1113.0,
136  125.0 / 192.0, -2187.0 / 6784.0,
137  11.0 / 84.0}};
138  // Coefficients for the embedded method, for generating an error measure
139  // during adaptive stepping (e.g. Sec. 7.2 of \cite NumericalRecipes).
140  static constexpr std::array<double, 7> b_alt_{
141  {5179.0 / 57600.0, 0.0, 7571.0 / 16695.0, 393.0 / 640.0,
142  -92097.0 / 339200.0, 187.0 / 2100.0, 1.0 / 40.0}};
143 
144  static const std::array<Time::rational_t, 6> c_;
145 
146  // Coefficients for dense output, taken from Sec. 7.2 of
147  // \cite NumericalRecipes
148  static constexpr std::array<double, 6> d_{
149  {-12715105075.0 / 11282082432.0, 87487479700.0 / 32700410799.0,
150  -10690763975.0 / 1880347072.0, 701980252875.0 / 199316789632.0,
151  -1453857185.0 / 822651844.0, 69997945.0 / 29380423.0}};
152 };
153 
154 inline bool constexpr operator==(const DormandPrince5& /*lhs*/,
155  const DormandPrince5& /*rhs*/) noexcept {
156  return true;
157 }
158 
159 inline bool constexpr operator!=(const DormandPrince5& /*lhs*/,
160  const DormandPrince5& /*rhs*/) noexcept {
161  return false;
162 }
163 
164 template <typename Vars, typename DerivVars>
165 void DormandPrince5::update_u(
166  const gsl::not_null<Vars*> u,
167  const gsl::not_null<History<Vars, DerivVars>*> history,
168  const TimeDelta& time_step) const noexcept {
169  ASSERT(history->integration_order() == 5,
170  "Fixed-order stepper cannot run at order "
171  << history->integration_order());
172  const size_t substep = (history->end() - 1).time_step_id().substep();
173 
174  // Clean up old history
175  if (substep == 0) {
176  history->mark_unneeded(history->end() - 1);
177  }
178 
179  const auto& u0 = history->begin().value();
180  const double dt = time_step.value();
181 
182  const auto increment_u = [&u, &history, &dt](const auto& coeffs) noexcept {
183  for (size_t i = 0; i < coeffs.size(); ++i) {
184  *u += (gsl::at(coeffs, i) * dt) *
185  (history->begin() + static_cast<int>(i)).derivative();
186  }
187  };
188 
189  if (substep == 0) {
190  *u = (history->end() - 1).value() +
191  (a2_ * dt) * history->begin().derivative();
192  } else if (substep < 6) {
193  *u = u0;
194  if (substep == 1) {
195  increment_u(a3_);
196  } else if (substep == 2) {
197  increment_u(a4_);
198  } else if (substep == 3) {
199  increment_u(a5_);
200  } else if (substep == 4) {
201  increment_u(a6_);
202  } else {
203  increment_u(b_);
204  }
205  } else {
206  ERROR("Substep in DP5 should be one of 0,1,2,3,4,5, not " << substep);
207  }
208 }
209 
210 template <typename Vars, typename ErrVars, typename DerivVars>
211 bool DormandPrince5::update_u(
212  const gsl::not_null<Vars*> u, const gsl::not_null<ErrVars*> u_error,
213  const gsl::not_null<History<Vars, DerivVars>*> history,
214  const TimeDelta& time_step) const noexcept {
215  ASSERT(history->integration_order() == 5,
216  "Fixed-order stepper cannot run at order "
217  << history->integration_order());
218  const size_t substep = (history->end() - 1).time_step_id().substep();
219 
220  // Clean up old history
221  if (substep == 0) {
222  history->mark_unneeded(history->end() - 1);
223  }
224 
225  const auto& u0 = history->begin().value();
226  const double dt = time_step.value();
227 
228  const auto increment_u = [&history, &dt](const auto& coeffs,
229  auto local_u) noexcept {
230  for (size_t i = 0; i < coeffs.size(); ++i) {
231  *local_u += (gsl::at(coeffs, i) * dt) *
232  (history->begin() + static_cast<int>(i)).derivative();
233  }
234  };
235  if (substep == 0) {
236  *u = (history->end() - 1).value() +
237  (a2_ * dt) * history->begin().derivative();
238  } else if (substep < 7) {
239  *u = u0;
240  if (substep == 1) {
241  increment_u(a3_, u);
242  } else if (substep == 2) {
243  increment_u(a4_, u);
244  } else if (substep == 3) {
245  increment_u(a5_, u);
246  } else if (substep == 4) {
247  increment_u(a6_, u);
248  } else if (substep == 5) {
249  increment_u(b_, u);
250  } else {
251  increment_u(b_, u);
252  *u_error = u0;
253  increment_u(b_alt_, u_error);
254  *u_error = *u - *u_error;
255  }
256  } else {
257  ERROR("Substep in adaptive DP5 should be one of 0,1,2,3,4,5,6, not "
258  << substep);
259  }
260  return substep == 6;
261 }
262 
263 template <typename Vars, typename DerivVars>
264 bool DormandPrince5::dense_update_u(const gsl::not_null<Vars*> u,
265  const History<Vars, DerivVars>& history,
266  const double time) const noexcept {
267  if ((history.end() - 1).time_step_id().substep() != 0) {
268  return false;
269  }
270  const double t0 = history.front().value();
271  const double t_end = history.back().value();
272  if (time == t_end) {
273  // Special case necessary for dense output at the initial time,
274  // before taking a step.
275  *u = (history.end() - 1).value();
276  return true;
277  }
278  const evolution_less<double> before{t_end > t0};
279  if (history.size() == 1 or before(t_end, time)) {
280  return false;
281  }
282  const double dt = t_end - t0;
283  const double output_fraction = (time - t0) / dt;
284  ASSERT(output_fraction >= 0.0, "Attempting dense output at time "
285  << time << ", but already progressed past "
286  << t0);
287  ASSERT(output_fraction <= 1.0, "Requested time ("
288  << time << ") not within step [" << t0
289  << ", " << t0 + dt << "]");
290 
291  const auto& u0 = history.begin().value();
292  const auto& u1 = (history.end() - 1).value();
293 
294  // We need the following: k1, k3, k4, k5, k6.
295  // Here k1 = dt * l1, k3 = dt * l3, etc.
296  const auto& l1 = history.begin().derivative();
297  const auto& l3 = (history.begin() + 2).derivative();
298  const auto& l4 = (history.begin() + 3).derivative();
299  const auto& l5 = (history.begin() + 4).derivative();
300  const auto& l6 = (history.begin() + 5).derivative();
301  const auto& l7 = (history.begin() + 6).derivative();
302 
303  // Compute the updating coefficents, called rcontN in Numerical recipes,
304  // that will be reused, so I don't have to compute them more than once.
305  const Vars rcont2 = u1 - u0;
306  const Vars rcont3 = dt * l1 - rcont2;
307 
308  // The formula for dense output is given in Numerical Recipes Sec. 17.2.3.
309  *u = u0 + output_fraction *
310  (rcont2 + (1.0 - output_fraction) *
311  (rcont3 + output_fraction *
312  ((rcont2 - dt * l7 - rcont3) +
313  ((1.0 - output_fraction) * dt) *
314  (d_[0] * l1 + d_[1] * l3 +
315  d_[2] * l4 + d_[3] * l5 +
316  d_[4] * l6 + d_[5] * l7))));
317  return true;
318 }
319 } // 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
evolution_comparator
Definition: EvolutionOrdering.hpp:20
TimeSteppers::DormandPrince5
Definition: DormandPrince5.hpp:59
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
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