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 : TimeStepId next_time_id(const TimeStepId& current_id, 94 : const TimeDelta& time_step) const override; 95 : 96 1 : TimeStepId next_time_id_for_error(const TimeStepId& current_id, 97 : const TimeDelta& time_step) const override; 98 : 99 0 : virtual const ButcherTableau& butcher_tableau() const = 0; 100 : 101 : private: 102 : template <typename T> 103 0 : void update_u_impl(gsl::not_null<T*> u, const ConstUntypedHistory<T>& history, 104 : const TimeDelta& time_step) const; 105 : 106 : template <typename T> 107 0 : std::optional<StepperErrorEstimate> update_u_impl( 108 : gsl::not_null<T*> u, const ConstUntypedHistory<T>& history, 109 : const TimeDelta& time_step, 110 : const std::optional<StepperErrorTolerances>& tolerances) const; 111 : 112 : template <typename T> 113 0 : void clean_history_impl(const MutableUntypedHistory<T>& history) const; 114 : 115 : template <typename T> 116 0 : bool dense_update_u_impl(gsl::not_null<T*> u, 117 : const ConstUntypedHistory<T>& history, 118 : double time) const; 119 : 120 : template <typename T> 121 0 : bool can_change_step_size_impl(const TimeStepId& time_id, 122 : const ConstUntypedHistory<T>& history) const; 123 : 124 : TIME_STEPPER_DECLARE_OVERLOADS 125 : }; 126 : } // namespace TimeSteppers