Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "TimeIntegrator.h"
      13             : 
      14             : /**
      15             :  * Fourth-order diagonally implicit Runge Kutta method (Dirk) with five stages.
      16             :  *
      17             :  * The Butcher tableau for this method is:
      18             :  * 1/4 |  1/4
      19             :  * 0   | -1/4 1/4
      20             :  * 1/2 |  1/8 1/8 1/4
      21             :  * 1   | -3/2 3/4 3/2   1/4
      22             :  * 1   |    0 1/6 2/3 -1/12 1/4
      23             :  * ----------------------------
      24             :  * 1   |    0 1/6 2/3 -1/12 1/4
      25             :  *
      26             :  *
      27             :  * The stability function for this method is:
      28             :  * R(z) = -(28*z**4 + 32*z**3 - 384*z**2 - 768*z + 3072)/
      29             :  *         (3*z**5 - 60*z**4 + 480*z**3 - 1920*z**2 + 3840*z - 3072)
      30             :  *
      31             :  * The method is L-stable:
      32             :  * lim R(z), z->oo = 0
      33             :  *
      34             :  * Notes:
      35             :  * I found the method here:
      36             :  *
      37             :  * L. M. Skvortsov, "Diagonally implicit Runge-Kutta Methods
      38             :  * for Stiff Problems", Computational Mathematics and Mathematical
      39             :  * Physics vol 46, no 12, pp. 2110-2123, 2006.
      40             :  *
      41             :  * but it may not be the original source.  There is also a 4th-order
      42             :  * rule with 5 stages on page 107 of:
      43             :  *
      44             :  * E. Hairer and G. Wanner, Solving Ordinary Differential Equations,
      45             :  * Vol. 2: Stiff and Differential-Algebraic Problems (Springer, Berlin,
      46             :  * 1987-1991; Mir, Moscow, 1999).
      47             :  *
      48             :  * but its coefficients have less favorable "amplification factors"
      49             :  * than the present rule.
      50             :  */
      51             : class LStableDirk4 : public TimeIntegrator
      52             : {
      53             : public:
      54             :   static InputParameters validParams();
      55             : 
      56             :   LStableDirk4(const InputParameters & parameters);
      57             : 
      58           0 :   virtual int order() override { return 4; }
      59             :   virtual void computeTimeDerivatives() override;
      60             :   virtual void computeADTimeDerivatives(ADReal & ad_u_dot,
      61             :                                         const dof_id_type & dof,
      62             :                                         ADReal & ad_u_dotdot) const override;
      63             :   virtual void solve() override;
      64             :   virtual void postResidual(NumericVector<Number> & residual) override;
      65         362 :   virtual bool overridesSolve() const override { return true; }
      66             : 
      67             : protected:
      68             :   /**
      69             :    * Helper function that actually does the math for computing the time derivative
      70             :    */
      71             :   template <typename T, typename T2>
      72             :   void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const;
      73             : 
      74             :   // Indicates the current stage.
      75             :   unsigned int _stage;
      76             : 
      77             :   // The number of stages in the method.  According to S9.4.2/4 of the
      78             :   // standard, we can specify a constant initializer like this for
      79             :   // integral types, it does not have to appear outside the class
      80             :   // definition.
      81             :   static const unsigned int _n_stages = 5;
      82             : 
      83             :   // Store pointers to the various stage residuals
      84             :   NumericVector<Number> * _stage_residuals[_n_stages];
      85             : 
      86             :   // Butcher tableau "C" parameters derived from _gamma
      87             :   static const Real _c[_n_stages];
      88             : 
      89             :   // Butcher tableau "A" values derived from _gamma.  We only use the
      90             :   // lower triangle of this.
      91             :   static const Real _a[_n_stages][_n_stages];
      92             : };
      93             : 
      94             : template <typename T, typename T2>
      95             : void
      96       18558 : LStableDirk4::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const
      97             : {
      98       18558 :   u_dot -= u_old;
      99       18558 :   u_dot *= 1. / _dt;
     100       18558 : }