RungeKutta4.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 RungeKutta4.
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"
18 #include "Time/Time.hpp"
19 #include "Time/TimeSteppers/TimeStepper.hpp" // IWYU pragma: keep
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 4th-order Runge-Kutta method, given e.g. in
38  * https://en.wikipedia.org/wiki/Runge-Kutta_methods
39  * that solves equations of the form
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{eqnarray}{
48  * v^{(1)} & = & u^n + dt\cdot \mathcal{L}(t^n, u^n)/2,\\
49  * v^{(2)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt/2, v^{(1)})/2,\\
50  * v^{(3)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt/2, v^{(2)}),\\
51  * v^{(4)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt, v^{(3)}),\\
52  * u^{n+1} & = & (2v^{(1)} + 4v^{(2)} + 2v^{(3)} + v^{(4)} - 3 u^n)/6.
53  * \f}
54  *
55  * Note that in the implementation, the expression for \f$u^{n+1}\f$ is
56  * computed simultaneously with \f$v^{(4)}\f$, so that there are
57  * actually only four substeps per step.
58  */
59 class RungeKutta4 : public TimeStepper::Inherit {
60  public:
61  using options = tmpl::list<>;
62  static constexpr OptionString help = {
63  "The standard fourth-order Runge-Kutta time-stepper."};
64 
65  RungeKutta4() = default;
66  RungeKutta4(const RungeKutta4&) noexcept = default;
67  RungeKutta4& operator=(const RungeKutta4&) noexcept = default;
68  RungeKutta4(RungeKutta4&&) noexcept = default;
69  RungeKutta4& operator=(RungeKutta4&&) noexcept = default;
70  ~RungeKutta4() 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 DerivVars>
78  void dense_update_u(gsl::not_null<Vars*> u,
79  const History<Vars, DerivVars>& history,
80  double time) const noexcept;
81 
82  uint64_t number_of_substeps() const noexcept override;
83 
84  size_t number_of_past_steps() const noexcept override;
85 
86  double stable_step() const noexcept override;
87 
88  TimeStepId next_time_id(const TimeStepId& current_id,
89  const TimeDelta& time_step) const noexcept override;
90 
91  template <typename Vars, typename DerivVars>
92  bool can_change_step_size(
93  const TimeStepId& time_id,
94  const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
95  noexcept {
96  return time_id.substep() == 0;
97  }
98 
99  WRAPPED_PUPable_decl_template(RungeKutta4); // NOLINT
100 
101  explicit RungeKutta4(CkMigrateMessage* /*unused*/) noexcept {}
102 
103  // clang-tidy: do not pass by non-const reference
104  void pup(PUP::er& p) noexcept override { // NOLINT
105  TimeStepper::Inherit::pup(p);
106  }
107 };
108 
109 inline bool constexpr operator==(const RungeKutta4& /*lhs*/,
110  const RungeKutta4& /*rhs*/) noexcept {
111  return true;
112 }
113 
114 inline bool constexpr operator!=(const RungeKutta4& /*lhs*/,
115  const RungeKutta4& /*rhs*/) noexcept {
116  return false;
117 }
118 
119 template <typename Vars, typename DerivVars>
120 void RungeKutta4::update_u(
121  const gsl::not_null<Vars*> u,
122  const gsl::not_null<History<Vars, DerivVars>*> history,
123  const TimeDelta& time_step) const noexcept {
124  const size_t substep = history->size() - 1;
125  const auto& dt_vars = (history->end() - 1).derivative();
126  const auto& u0 = history->begin().value();
127 
128  switch (substep) {
129  case 0: {
130  // from (17.1.3) of Numerical Recipes 3rd Edition
131  // v^(1) = u^n + dt * \mathcal{L}(u^n,t^n)/2
132  // On entry V = u^n, u0 = u^n, rhs0 = \mathcal{L}(u^n, t^n),
133  // time = t^n
134  *u += 0.5 * time_step.value() * dt_vars;
135  // On exit v = v^(1), time = t^n + (1/2)*dt
136  break;
137  }
138  case 1: {
139  // from (17.1.3) of Numerical Recipes 3rd Edition
140  // v^(2) = u^n + dt * \mathcal{L}(v^(1), t^n + (1/2)*dt)/2
141  // On entry V = v^(1), u0 = u^n, rhs0 = \mathcal{L}(v^(1), t^n + dt/2),
142  // time = t^n + dt
143  *u = u0 + 0.5 * time_step.value() * dt_vars;
144  // On exit v = v^(2), time = t^n + (1/2)*dt
145  break;
146  }
147  case 2: {
148  // from (17.1.3) of Numerical Recipes 3rd Edition
149  // v^(3) = u^n + dt * \mathcal{L}(v^(2), t^n + (1/2)*dt))
150  // On entry V = v^(2), u0 = u^n,
151  // rhs0 = \mathcal{L}(v^(2), t^n + (1/2)*dt), time = t^n + (1/2)*dt
152  *u = u0 + time_step.value() * dt_vars;
153  // On exit v = v^(3), time = t^n + dt
154  break;
155  }
156  case 3: {
157  // from (17.1.3) of Numerical Recipes 3rd Edition
158  // u^(n+1) = (2v^(1) + 4*v^(2) + 2*v^(3) + v^(4) - 3*u0)/6
159  // On entry V = v^(3), u0 = u^n, rhs0 = \mathcal{L}(v^(3), t^n + dt),
160  // time = t^n + dt
161  // Note: v^(4) = u0 + dt * \mathcal{L}(t+dt, v^(3)); inserting this gives
162  // u^(n+1) = (2v^(1) + 4*v^(2) + 2*v^(3)
163  // + dt*\mathcal{L}(t+dt,v^(3)) - 2*u0)/6
164  constexpr double one_sixth = 1.0 / 6.0;
165  *u = (2.0 * (history->begin() + 1).value() +
166  4.0 * (history->begin() + 2).value() +
167  2.0 * (history->begin() + 3).value() +
168  (time_step.value() * dt_vars - 2.0 * u0)) *
169  one_sixth;
170  // On exit v = u^(n+1), time = t^n + dt
171  break;
172  }
173  default:
174  ERROR("Substep in RK4 should be one of 0,1,2,3, 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 RungeKutta4::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  "RK4 can only dense output on last substep ("
189  << number_of_substeps() - 1 << "), not substep "
190  << history.size() - 1);
191  const double step_start = history[0].value();
192  const double step_end = history[history.size() - 1].value();
193  const double time_step = step_end - step_start;
194  const double output_fraction = (time - step_start) / time_step;
195  ASSERT(output_fraction >= 0.0, "Attempting dense output at time "
196  << time << ", but already progressed past "
197  << step_start);
198  ASSERT(output_fraction <= 1.0, "Requested time ("
199  << time << ") not within step ["
200  << step_start << ", " << step_end << "]");
201 
202  // Numerical Recipes Eq. (17.2.15). This implements cubic interpolation
203  // throughout the step. Because the history only is available through the
204  // penultimate step, i) the value after the step is computed algebraically
205  // from previous substeps, and ii) the derivative at the final step is
206  // approximated as the derivative at the penultimate substep.
207  constexpr double one_sixth = 1.0 / 6.0;
208  const auto& u0 = history.begin().value();
209  const auto& dt_vars = (history.end() - 1).derivative();
210  const Vars u1 =
211  (2.0 * (history.begin() + 1).value() +
212  4.0 * (history.begin() + 2).value() +
213  2.0 * (history.begin() + 3).value() + (time_step * dt_vars - 2.0 * u0)) *
214  one_sixth;
215  *u = (1.0 - output_fraction) * u0 + output_fraction * u1 +
216  (output_fraction) * (output_fraction - 1.0) *
217  ((1.0 - 2.0 * output_fraction) * (u1 - u0) +
218  (output_fraction - 1.0) * time_step * history.begin().derivative() +
219  output_fraction * time_step * dt_vars);
220 }
221 } // namespace TimeSteppers
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
Defines classes and functions for making classes creatable from input files.
Defines macros to allow serialization of abstract template base classes.
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
History data used by a TimeStepper.
Definition: History.hpp:34
Defines Time and TimeDelta.
A unique identifier for the temporal state of an integrated system.
Definition: TimeStepId.hpp:25
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
const_iterator begin() const noexcept
These iterators directly return the Time of the past values. The other data can be accessed through t...
Definition: History.hpp:68
Represents an interval of time within a single slab.
Definition: Time.hpp:88
Define simple functions for constant expressions.
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
Defines functions and classes from the GSL.
const_iterator end() const noexcept
These iterators directly return the Time of the past values. The other data can be accessed through t...
Definition: History.hpp:73
Holds classes that take time steps.
Definition: BoundaryHistory.hpp:24
Defines macro ERROR.
double value() const noexcept
Approximate numerical length of the interval.
Definition: Time.cpp:111
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182
Definition: RungeKutta4.hpp:59