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             : class LStableDirk4;
      15             : 
      16             : /**
      17             :  * Fourth-order diagonally implicit Runge Kutta method (Dirk) with
      18             :  * three stages plus an update.
      19             :  *
      20             :  * The Butcher tableau for this method is:
      21             :  * gamma   | gamma                   0                          0
      22             :  * 1/2     | 1/2-gamma               gamma                      0
      23             :  * 1-gamma | 2*gamma                 1-4*gamma                  gamma
      24             :  * --------|-------------------------------------------------------------------------
      25             :  *         | 1/(24*(1/2-gamma)**2)   1 - 1/(12*(1/2-gamma)**2)  1/(24*(1/2-gamma)**2)
      26             :  *
      27             :  * where gamma = 1/2 + sqrt(3)/3 * cos(pi/18) ~ 1.06857902130162881
      28             :  *
      29             :  * The stability function for this method is:
      30             :  * R(z) = (-0.76921266689461*z**3 - 0.719846311954034*z**2 + 2.20573706179581*z - 1.0)/
      31             :  *        (1.22016884572716*z**3 - 3.42558337748497*z**2 + 3.20573706551682*z - 1.0)
      32             :  *
      33             :  * The method is *not* L-stable, it is only A-stable:
      34             :  # lim R(z), z->oo = -0.630414937726258
      35             :  *
      36             :  * Notes:
      37             :  * 1.) Method is originally due to:
      38             :  *     M. Crouzeix, "Sur l'approximation des equations differentielles
      39             :  *     operationelles lineaires par des methodes de Runge Kutta",
      40             :  *     Ph.D. thesis, Universite Paris VI, Paris, 1975.
      41             :  * 2.) Since gamma is slightly larger than 1, the first stage involves
      42             :  *     evaluation of the non-time residuals for t > t_n+dt, while the
      43             :  *     third stage involves evaluation of the non-time residual for t
      44             :  *     < t_n, which may present an issue for the first timestep (if
      45             :  *     e.g. material properties or forcing functions are not defined
      46             :  *     for t<0. We could handle this by using an alternate (more
      47             :  *     expensive) method in the first timestep, or by using a
      48             :  *     lower-order method for the first timestep and then switching to
      49             :  *     this method for subsequent timesteps.
      50             :  */
      51             : class AStableDirk4 : public TimeIntegrator
      52             : {
      53             : public:
      54             :   static InputParameters validParams();
      55             : 
      56             :   AStableDirk4(const InputParameters & parameters);
      57             : 
      58           0 :   virtual int order() override { return 4; }
      59             :   virtual void computeTimeDerivatives() override;
      60             :   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        2804 :   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             :   // Store pointers to the various stage residuals
      78             :   NumericVector<Number> * _stage_residuals[3];
      79             : 
      80             :   // The parameter of the method, set at construction time and cannot be changed.
      81             :   const Real _gamma; // 1.06857902130162881
      82             : 
      83             :   // Butcher tableau "C" parameters derived from _gamma
      84             :   // 1.06857902130162881, 0.5, -.06857902130162881
      85             :   Real _c[3];
      86             : 
      87             :   // Butcher tableau "A" values derived from _gamma.  We only use the
      88             :   // lower triangle of this.
      89             :   // 1.06857902130162881
      90             :   // -.56857902130162881, 1.06857902130162881
      91             :   // 2.13715804260325762, -3.27431608520651524, 1.06857902130162881
      92             :   Real _a[3][3];
      93             : 
      94             :   // The Butcher tableau "b" parameters derived from _gamma;
      95             :   // 1.2888640051572051e-01, 7.4222719896855893e-01, 1.2888640051572051e-01
      96             :   Real _b[3];
      97             : 
      98             :   // If true, we use a more expensive method (LStableDirk4) to
      99             :   // "bootstrap" the first timestep of this method and avoid
     100             :   // evaluating residuals before the initial time.
     101             :   bool _safe_start;
     102             : 
     103             :   // A pointer to the "bootstrapping" method to use if _safe_start==true.
     104             :   std::shared_ptr<LStableDirk4> _bootstrap_method;
     105             : };
     106             : 
     107             : template <typename T, typename T2>
     108             : void
     109       89874 : AStableDirk4::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const
     110             : {
     111       89874 :   u_dot -= u_old;
     112       89874 :   u_dot *= 1. / _dt;
     113       89874 : }