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 "ErrorHandling/Assert.hpp"
15 #include "ErrorHandling/Error.hpp"
16 #include "Options/Options.hpp"
17 #include "Parallel/CharmPupable.hpp"
18 #include "Time/Time.hpp"
19 #include "Time/TimeSteppers/TimeStepper.hpp"
20 #include "Utilities/ConstantExpressions.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,
73  gsl::not_null<History<Vars, DerivVars>*> history,
74  const TimeDelta& time_step) const noexcept;
75 
76  template <typename Vars, typename DerivVars>
77  void dense_update_u(gsl::not_null<Vars*> u,
78  const History<Vars, DerivVars>& history,
79  double time) const noexcept;
80 
81  uint64_t number_of_substeps() const noexcept override;
82 
83  size_t number_of_past_steps() const noexcept override;
84 
85  double stable_step() const noexcept override;
86 
87  TimeStepId next_time_id(const TimeStepId& current_id,
88  const TimeDelta& time_step) const noexcept override;
89 
90  template <typename Vars, typename DerivVars>
91  bool can_change_step_size(
92  const TimeStepId& time_id,
93  const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
94  noexcept {
95  return time_id.substep() == 0;
96  }
97 
98  WRAPPED_PUPable_decl_template(DormandPrince5); // NOLINT
99 
100  explicit DormandPrince5(CkMigrateMessage* /*unused*/) noexcept {}
101 
102  // clang-tidy: do not pass by non-const reference
103  void pup(PUP::er& p) noexcept override { // NOLINT
104  TimeStepper::Inherit::pup(p);
105  }
106 
107  private:
108  // Coefficients from the Dormand-Prince 5 Butcher tableau (e.g. Sec. 7.2
109  // of \cite NumericalRecipes).
110  static constexpr double a2_ = 0.2;
111  static constexpr std::array<double, 2> a3_{{3.0 / 40.0, 9.0 / 40.0}};
112  static constexpr std::array<double, 3> a4_{
113  {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0}};
114  static constexpr std::array<double, 4> a5_{
115  {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0}};
116  static constexpr std::array<double, 5> a6_{{9017.0 / 3168.0, -355.0 / 33.0,
117  46732.0 / 5247.0, 49.0 / 176.0,
118  -5103.0 / 18656.0}};
119  static constexpr std::array<double, 6> b_{{35.0 / 384.0, 0.0, 500.0 / 1113.0,
120  125.0 / 192.0, -2187.0 / 6784.0,
121  11.0 / 84.0}};
122  static const std::array<Time::rational_t, 5> c_;
123 
124  // Coefficients for dense output, taken from Sec. 7.2 of
125  // \cite NumericalRecipes
126  static constexpr std::array<double, 6> d_{
127  {-12715105075.0 / 11282082432.0, 87487479700.0 / 32700410799.0,
128  -10690763975.0 / 1880347072.0, 701980252875.0 / 199316789632.0,
129  -1453857185.0 / 822651844.0, 69997945.0 / 29380423.0}};
130 };
131 
132 inline bool constexpr operator==(const DormandPrince5& /*lhs*/,
133  const DormandPrince5& /*rhs*/) noexcept {
134  return true;
135 }
136 
137 inline bool constexpr operator!=(const DormandPrince5& /*lhs*/,
138  const DormandPrince5& /*rhs*/) noexcept {
139  return false;
140 }
141 
142 template <typename Vars, typename DerivVars>
143 void DormandPrince5::update_u(
144  const gsl::not_null<Vars*> u,
145  const gsl::not_null<History<Vars, DerivVars>*> history,
146  const TimeDelta& time_step) const noexcept {
147  const size_t substep = history->size() - 1;
148  const auto& u0 = history->begin().value();
149  const double dt = time_step.value();
150 
151  const auto increment_u = [&u, &history, &dt](const auto& coeffs) noexcept {
152  for (size_t i = 0; i < coeffs.size(); ++i) {
153  *u += (gsl::at(coeffs, i) * dt) *
154  (history->begin() + static_cast<int>(i)).derivative();
155  }
156  };
157 
158  if (substep == 0) {
159  *u += (a2_ * dt) * history->begin().derivative();
160  } else if (substep < 6) {
161  *u = u0;
162  if (substep == 1) {
163  increment_u(a3_);
164  } else if (substep == 2) {
165  increment_u(a4_);
166  } else if (substep == 3) {
167  increment_u(a5_);
168  } else if (substep == 4) {
169  increment_u(a6_);
170  } else {
171  increment_u(b_);
172  }
173  } else {
174  ERROR("Substep in DP5 should be one of 0,1,2,3,4,5, not " << substep);
175  }
176 
177  // Clean up old history
178  if (history->size() == number_of_substeps()) {
179  history->mark_unneeded(history->end());
180  }
181 }
182 
183 template <typename Vars, typename DerivVars>
184 void DormandPrince5::dense_update_u(const gsl::not_null<Vars*> u,
185  const History<Vars, DerivVars>& history,
186  const double time) const noexcept {
187  ASSERT(history.size() == number_of_substeps(),
188  "DP5 can only dense output on last substep ("
189  << number_of_substeps() - 1 << "), not substep "
190  << history.size() - 1);
191  const double t0 = history[0].value();
192  const double t_end = history[history.size() - 1].value();
193  // The history does not contain the final step; specifically,
194  // step_end = t0 + c[4] * dt, so dt = (t_end - t0) / c[4]. But since
195  // c[4] = 1.0 for DP5, we don't need to divide by c[4] here.
196  const double dt = t_end - t0;
197  const double output_fraction = (time - t0) / dt;
198  ASSERT(output_fraction >= 0.0, "Attempting dense output at time "
199  << time << ", but already progressed past "
200  << t0);
201  ASSERT(output_fraction <= 1.0, "Requested time ("
202  << time << ") not within step [" << t0
203  << ", " << t0 + dt << "]");
204 
205  // Get the solution u1 at the final time
206  const auto& u0 = history.begin().value();
207  Vars u1 = u0;
208  for (size_t i = 0; i < 6; ++i) {
209  u1 += (gsl::at(b_, i) * dt) *
210  (history.begin() + static_cast<int>(i)).derivative();
211  }
212 
213  // We need the following: k1, k3, k4, k5, k6.
214  // Here k1 = dt * l1, k3 = dt * l3, etc.
215  const auto& l1 = history.begin().derivative();
216  const auto& l3 = (history.begin() + 2).derivative();
217  const auto& l4 = (history.begin() + 3).derivative();
218  const auto& l5 = (history.begin() + 4).derivative();
219  const auto& l6 = (history.begin() + 5).derivative();
220 
221  // Compute the updating coefficents, called rcontN in Numerical recipes,
222  // that will be reused, so I don't have to compute them more than once.
223  const Vars rcont2 = u1 - u0;
224  const Vars rcont3 = dt * l1 - rcont2;
225 
226  // The formula for dense output is given in Numerical Recipes Sec. 17.2.3.
227  // Note: L(t+dt, u1) after the step is unavailable in the history; so here,
228  // approximate L(t+dt, u1) by l6.
229  *u = u0 + output_fraction *
230  (rcont2 +
231  (1.0 - output_fraction) *
232  (rcont3 + output_fraction *
233  ((rcont2 - dt * l6 - rcont3) +
234  ((1.0 - output_fraction) * dt) *
235  (d_[0] * l1 + d_[1] * l3 + d_[2] * l4 +
236  d_[3] * l5 + (d_[4] + d_[5]) * l6))));
237 }
238 } // 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:51
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: Gsl.hpp:183
© Copyright 2017 - 2020 SXS Collaboration, MIT License