Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <cstdint>
8 : #include <optional>
9 : #include <vector>
10 :
11 : #include "Time/StepperErrorEstimate.hpp"
12 : #include "Time/TimeStepId.hpp"
13 : #include "Time/TimeSteppers/TimeStepper.hpp"
14 : #include "Utilities/Gsl.hpp"
15 :
16 : /// \cond
17 : struct StepperErrorTolerances;
18 : class TimeDelta;
19 : namespace TimeSteppers {
20 : template <typename T>
21 : class ConstUntypedHistory;
22 : template <typename T>
23 : class MutableUntypedHistory;
24 : } // namespace TimeSteppers
25 : /// \endcond
26 :
27 : namespace TimeSteppers {
28 : /*!
29 : * \ingroup TimeSteppersGroup
30 : * Intermediate base class implementing a generic Runge-Kutta scheme.
31 : *
32 : * Implements most of the virtual methods of TimeStepper for a generic
33 : * Runge-Kutta method. From the TimeStepper interface, derived
34 : * classes need only implement `order` and `stable_step`, and should
35 : * not include the `TIME_STEPPER_*` macros. All other methods are
36 : * implemented in terms of a Butcher tableau returned by the
37 : * `butcher_tableau` function.
38 : */
39 1 : class RungeKutta : public virtual TimeStepper {
40 : public:
41 0 : struct ButcherTableau {
42 : /*!
43 : * The times of the substeps, excluding the initial time step.
44 : * Often called \f$c\f$ in the literature.
45 : */
46 1 : std::vector<double> substep_times;
47 : /*!
48 : * The coefficient matrix of the substeps. Do not include the
49 : * initial empty row or the coefficients for the full step. Often
50 : * called \f$A\f$ in the literature.
51 : */
52 1 : std::vector<std::vector<double>> substep_coefficients;
53 : /*!
54 : * The coefficients for the final result. Often called \f$b\f$ in
55 : * the literature.
56 : *
57 : * If the number of coefficients is smaller than the number of
58 : * substeps defined in `substep_coefficients`, the extra substeps
59 : * will only be performed when an error estimate is requested.
60 : */
61 1 : std::vector<double> result_coefficients;
62 : /*!
63 : * The coefficients for an error estimate. Often called
64 : * \f$b^*\f$ or \f$\hat{b}\f$ in the literature. This is assumed
65 : * to be one order less accurate than the main method.
66 : */
67 1 : std::vector<double> error_coefficients;
68 : /*!
69 : * Coefficient polynomials for dense output. Each entry is the
70 : * coefficients of a polynomial that will be evaluated with the
71 : * fraction of the way through the step at which output is
72 : * desired:
73 : *
74 : * \f{equation}{
75 : * y = y_0 + \sum_i y'_i \sum_j p_{ij} x^j \qquad 0 \le x \le 1
76 : * \f}
77 : *
78 : * The derivative at the start of the next step is available as an
79 : * additional substep after the final real substep. Trailing zero
80 : * polynomials can be omitted.
81 : */
82 1 : std::vector<std::vector<double>> dense_coefficients;
83 : };
84 :
85 1 : uint64_t number_of_substeps() const override;
86 :
87 1 : uint64_t number_of_substeps_for_error() const override;
88 :
89 1 : size_t number_of_past_steps() const override;
90 :
91 1 : bool monotonic() const override;
92 :
93 1 : bool dense_output_uses_fsal() const override;
94 :
95 1 : TimeStepId next_time_id(const TimeStepId& current_id,
96 : const TimeDelta& time_step) const override;
97 :
98 1 : TimeStepId next_time_id_for_error(const TimeStepId& current_id,
99 : const TimeDelta& time_step) const override;
100 :
101 0 : virtual const ButcherTableau& butcher_tableau() const = 0;
102 :
103 : private:
104 : template <typename T>
105 0 : void update_u_impl(gsl::not_null<T*> u, const ConstUntypedHistory<T>& history,
106 : const TimeDelta& time_step) const;
107 :
108 : template <typename T>
109 0 : std::optional<StepperErrorEstimate> update_u_impl(
110 : gsl::not_null<T*> u, const ConstUntypedHistory<T>& history,
111 : const TimeDelta& time_step,
112 : const StepperErrorTolerances& tolerances) const;
113 :
114 : template <typename T>
115 0 : void clean_history_impl(const MutableUntypedHistory<T>& history) const;
116 :
117 : template <typename T>
118 0 : bool dense_update_u_impl(gsl::not_null<T*> u,
119 : const ConstUntypedHistory<T>& history,
120 : double time) const;
121 :
122 : template <typename T>
123 0 : bool can_change_step_size_impl(const TimeStepId& time_id,
124 : const ConstUntypedHistory<T>& history) const;
125 :
126 : TIME_STEPPER_DECLARE_OVERLOADS
127 : };
128 : } // namespace TimeSteppers
|
