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