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