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