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             :  * Third order diagonally implicit Runge Kutta method (Dirk) with three stages.
      16             :  *
      17             :  * The Butcher tableau for this method is:
      18             :  * gamma       | gamma
      19             :  * (1+gamma)/2 | (1-gamma)/2                         gamma
      20             :  * 1           | (1/4)*(-6*gamma**2 + 16*gamma - 1)  (1/4)*(6*gamma**2 - 20*gamma + 5)  gamma
      21             :  * ------------------------------------------------------------------------------------------
      22             :  *             | (1/4)*(-6*gamma**2 + 16*gamma - 1)  (1/4)*(6*gamma**2 - 20*gamma + 5)  gamma
      23             :  *
      24             :  * where gamma = -sqrt(2)*cos(atan(sqrt(2)/4)/3)/2 + sqrt(6)*sin(atan(sqrt(2)/4)/3)/2 + 1  ~
      25             :  * .435866521508459
      26             :  *
      27             :  * The stability function for this method is:
      28             :  * R(z) = (1.90128552647780115*z**2 + 2.46079651620301599*z - 8) /
      29             :  *        (0.662446064957040178*z**3 - 4.55951098972521484*z**2 + 10.460796516203016*z - 8)
      30             :  *
      31             :  * The method is L-stable:
      32             :  * lim R(z), z->oo = 0
      33             :  *
      34             :  * Notes: This method is derived in detail in: R. Alexander,
      35             :  * "Diagonally implicit Runge-Kutta Methods for Stiff ODEs", SIAM
      36             :  * J. Numer. Anal., 14(6), Dec. 1977, pg. 1006-1021.  Unlike BDF3,
      37             :  * this method is L-stable and so may be more suitable for "stiff"
      38             :  * problems.
      39             :  */
      40             : class LStableDirk3 : public TimeIntegrator
      41             : {
      42             : public:
      43             :   static InputParameters validParams();
      44             : 
      45             :   LStableDirk3(const InputParameters & parameters);
      46             : 
      47           0 :   virtual int order() override { return 3; }
      48             :   virtual void computeTimeDerivatives() override;
      49             :   virtual void computeADTimeDerivatives(ADReal & ad_u_dot,
      50             :                                         const dof_id_type & dof,
      51             :                                         ADReal & ad_u_dotdot) const override;
      52             :   virtual void solve() override;
      53             :   virtual void postResidual(NumericVector<Number> & residual) override;
      54        2518 :   virtual bool overridesSolve() const override { return true; }
      55             : 
      56             : protected:
      57             :   /**
      58             :    * Helper function that actually does the math for computing the time derivative
      59             :    */
      60             :   template <typename T, typename T2>
      61             :   void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const;
      62             : 
      63             :   // Indicates the current stage.
      64             :   unsigned int _stage;
      65             : 
      66             :   // Store pointers to the various stage residuals
      67             :   NumericVector<Number> * _stage_residuals[3];
      68             : 
      69             :   // The parameter of the method, set at construction time and cannot be changed.
      70             :   const Real _gamma; // 0.4358665215084589
      71             : 
      72             :   // Butcher tableau "C" parameters derived from _gamma
      73             :   // 0.4358665215084589, 0.7179332607542295, 1.0000000000000000
      74             :   Real _c[3];
      75             : 
      76             :   // Butcher tableau "A" values derived from _gamma.  We only use the
      77             :   // lower triangle of this.
      78             :   // 0.4358665215084589
      79             :   // 0.2820667392457705,  0.4358665215084589
      80             :   // 1.2084966491760099, -0.6443631706844688, 0.4358665215084589
      81             :   Real _a[3][3];
      82             : };
      83             : 
      84             : template <typename T, typename T2>
      85             : void
      86       65992 : LStableDirk3::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const
      87             : {
      88       65992 :   u_dot -= u_old;
      89       65992 :   u_dot *= 1. / _dt;
      90       65992 : }