Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <vector>
       7             : 
       8             : #include "Time/TimeSteppers/ImexTimeStepper.hpp"
       9             : #include "Time/TimeSteppers/RungeKutta.hpp"
      10             : #include "Utilities/Gsl.hpp"
      11             : 
      12             : /// \cond
      13             : class TimeDelta;
      14             : namespace TimeSteppers {
      15             : template <typename T>
      16             : class ConstUntypedHistory;
      17             : template <typename T>
      18             : class MutableUntypedHistory;
      19             : }  // namespace TimeSteppers
      20             : /// \endcond
      21             : 
      22             : namespace TimeSteppers {
      23             : /*!
      24             :  * \ingroup TimeSteppersGroup
      25             :  * Intermediate base class implementing a generic IMEX Runge-Kutta
      26             :  * scheme.
      27             :  *
      28             :  * Implements most of the virtual methods of ImexTimeStepper for a
      29             :  * generic Runge-Kutta method.  Derived classes must implement the
      30             :  * requirements of the `RungeKutta` base class, as well as
      31             :  * `imex_order()`, `implicit_stage_order()`, and
      32             :  * `implicit_butcher_tableau()`.
      33             :  */
      34           1 : class ImexRungeKutta : public virtual RungeKutta,
      35             :                        public virtual ImexTimeStepper {
      36             :  public:
      37             :   /// Implicit part of the Butcher tableau.  Most parts of the tableau
      38             :   /// must be the same as the explicit part, and so are omitted.
      39           1 :   struct ImplicitButcherTableau {
      40             :     /*!
      41             :      * The coefficient matrix of the substeps.  We can only reasonably
      42             :      * support EDIRK methods (including special cases such as ESDIRK
      43             :      * and QESDIRK, see \cite Kennedy2016), so the tableau must be
      44             :      * lower-triangular with an empty first row.  As with the explicit
      45             :      * tableau, the initial blank row should be omitted.  For a
      46             :      * stiffly-accurate method, the final row must be the same as the
      47             :      * result coefficients in the explicit tableau.
      48             :      *
      49             :      * More general DIRK methods can be implemented inefficiently by
      50             :      * adding an unused initial substep to convert them to EDIRK form.
      51             :      */
      52           1 :     std::vector<std::vector<double>> substep_coefficients;
      53             :   };
      54             : 
      55             :   /*!
      56             :    * Smallest order of the intermediate result at any substep.  For
      57             :    * the methods supported by this class, this cannot exceed 2.
      58             :    */
      59           1 :   virtual size_t implicit_stage_order() const = 0;
      60             : 
      61           0 :   virtual const ImplicitButcherTableau& implicit_butcher_tableau() const = 0;
      62             : 
      63             :  private:
      64             :   template <typename T>
      65           0 :   void add_inhomogeneous_implicit_terms_impl(
      66             :       gsl::not_null<T*> u, const MutableUntypedHistory<T>& implicit_history,
      67             :       const TimeDelta& time_step) const;
      68             : 
      69             :   template <typename T>
      70           0 :   double implicit_weight_impl(const MutableUntypedHistory<T>& implicit_history,
      71             :                               const TimeDelta& time_step) const;
      72             : 
      73             :   IMEX_TIME_STEPPER_DECLARE_OVERLOADS
      74             : };
      75             : }  // namespace TimeSteppers