RungeKutta3.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 RungeKutta3.
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 TimeId;
26 namespace TimeSteppers {
27 template <typename Vars, typename DerivVars>
28 class History;
29 } // namespace TimeSteppers
30 /// \endcond
31 
32 namespace TimeSteppers {
33 
34 /// \ingroup TimeSteppersGroup
35 ///
36 /// A "strong stability-preserving" 3rd-order Runge-Kutta time-stepper.
37 /// Major reference: J. Hesthaven & T. Warburton, Nodal Discontinuous
38 /// Galerkin Methods. section 5.7
39 class RungeKutta3 : public TimeStepper::Inherit {
40  public:
41  using options = tmpl::list<>;
42  static constexpr OptionString help = {
43  "A third-order strong stability-preserving Runge-Kutta time-stepper."};
44 
45  RungeKutta3() = default;
46  RungeKutta3(const RungeKutta3&) noexcept = default;
47  RungeKutta3& operator=(const RungeKutta3&) noexcept = default;
48  RungeKutta3(RungeKutta3&&) noexcept = default;
49  RungeKutta3& operator=(RungeKutta3&&) noexcept = default;
50  ~RungeKutta3() noexcept override = default;
51 
52  template <typename Vars, typename DerivVars>
53  void update_u(gsl::not_null<Vars*> u,
55  const TimeDelta& time_step) const noexcept;
56 
57  template <typename Vars, typename DerivVars>
58  void dense_update_u(gsl::not_null<Vars*> u,
59  const History<Vars, DerivVars>& history,
60  double time) const noexcept;
61 
62  uint64_t number_of_substeps() const noexcept override;
63 
64  size_t number_of_past_steps() const noexcept override;
65 
66  double stable_step() const noexcept override;
67 
68  TimeId next_time_id(const TimeId& current_id,
69  const TimeDelta& time_step) const noexcept override;
70 
71  WRAPPED_PUPable_decl_template(RungeKutta3); // NOLINT
72 
73  explicit RungeKutta3(CkMigrateMessage* /*unused*/) noexcept {}
74 
75  // clang-tidy: do not pass by non-const reference
76  void pup(PUP::er& p) noexcept override { // NOLINT
77  TimeStepper::Inherit::pup(p);
78  }
79 };
80 
81 inline bool constexpr operator==(const RungeKutta3& /*lhs*/,
82  const RungeKutta3& /*rhs*/) noexcept {
83  return true;
84 }
85 
86 inline bool constexpr operator!=(const RungeKutta3& /*lhs*/,
87  const RungeKutta3& /*rhs*/) noexcept {
88  return false;
89 }
90 
91 template <typename Vars, typename DerivVars>
92 void RungeKutta3::update_u(
93  const gsl::not_null<Vars*> u,
95  const TimeDelta& time_step) const noexcept {
96  const size_t substep = history->size() - 1;
97  const auto& vars = (history->end() - 1).value();
98  const auto& dt_vars = (history->end() - 1).derivative();
99  const auto& U0 = history->begin().value();
100 
101  switch (substep) {
102  case 0: {
103  // from (5.32) of Hesthaven
104  // v^(1) = u^n + dt*RHS(u^n,t^n)
105  // On entry V = u^n, U0 = u^n, rhs0 = RHS(u^n,t^n),
106  // time = t^n
107  *u += time_step.value() * dt_vars;
108  // On exit v = v^(1), time = t^n + dt
109  break;
110  }
111  case 1: {
112  // from (5.32) of Hesthaven
113  // v^(2) = (1/4)*( 3*u^n + v^(1) + dt*RHS(v^(1),t^n + dt) )
114  // On entry V = v^(1), U0 = u^n, rhs0 = RHS(v^(1),t^n + dt),
115  // time = t^n + dt
116  *u += 0.25 * (3.0 * (U0 - vars) + time_step.value() * dt_vars);
117  // On exit v = v^(2), time = t^n + (1/2)*dt
118  break;
119  }
120  case 2: {
121  // from (5.32) of Hesthaven
122  // u^(n+1) = (1/3)*( u^n + 2*v^(2) + 2*dt*RHS(v^(2),t^n + (1/2)*dt) )
123  // On entry V = v^(2), U0 = u^n, rhs0 = RHS(v^(2),t^n + (1/2)*dt),
124  // time = t^n + (1/2)*dt
125  *u += (1.0 / 3.0) * (U0 - vars + 2.0 * time_step.value() * dt_vars);
126  // On exit v = u^(n+1), time = t^n + dt
127  break;
128  }
129  default:
130  ERROR("Bad substep value in RK3: " << substep);
131  }
132 
133  // Clean up old history
134  if (history->size() == number_of_substeps()) {
135  history->mark_unneeded(history->end());
136  }
137 }
138 
139 template <typename Vars, typename DerivVars>
140 void RungeKutta3::dense_update_u(gsl::not_null<Vars*> u,
141  const History<Vars, DerivVars>& history,
142  const double time) const noexcept {
143  ASSERT(history.size() == 3, "Can only dense output on last substep");
144  const double step_start = history[0].value();
145  const double step_end = history[1].value();
146  const double time_step = step_end - step_start;
147  const double output_fraction = (time - step_start) / time_step;
148  ASSERT(output_fraction >= 0, "Attempting dense output at time " << time
149  << ", but already progressed past " << step_start);
150  ASSERT(output_fraction <= 1,
151  "Requested time (" << time << " not within step [" << step_start
152  << ", " << step_end << "]");
153 
154  // arXiv:1605.02429
155  *u += (1.0 - output_fraction * (1.0 - output_fraction / 3.0)) *
156  history.begin().value() +
157  output_fraction * (1.0 - output_fraction) *
158  (history.begin() + 1).value() +
159  (2.0 / 3.0 * square(output_fraction) - 1.0) *
160  (history.begin() + 2).value() +
161  2.0 / 3.0 * square(output_fraction) * time_step *
162  (history.begin() + 2).derivative();
163 }
164 } // namespace TimeSteppers
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
A unique identifier for the temporal state of an integrated system.
Definition: TimeId.hpp:25
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:49
History data used by a TimeStepper.
Definition: History.hpp:25
Defines Time and TimeDelta.
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
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:61
Represents an interval of time within a single slab.
Definition: Time.hpp:108
Define simple functions for constant expressions.
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
A "strong stability-preserving" 3rd-order Runge-Kutta time-stepper. Major reference: J...
Definition: RungeKutta3.hpp:39
Defines functions and classes from the GSL.
Holds classes that take time steps.
Definition: BoundaryHistory.hpp:23
Defines macro ERROR.
double value() const noexcept
Approximate numerical length of the interval.
Definition: Time.hpp:131
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12