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