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             :  * Second-order Runge-Kutta (implicit midpoint) time integration.
      16             :  *
      17             :  * The Butcher tableau for this method is:
      18             :  * 1/2 | 1/2
      19             :  * ----------
      20             :  *     | 1
      21             :  *
      22             :  * The stability function for this method is:
      23             :  * R(z) = -(z + 2)/(z - 2)
      24             :  * which is the same as Crank-Nicolson.
      25             :  *
      26             :  * The method is A-stable, but not L-stable,
      27             :  * lim R(z), z->oo = -1
      28             :  *
      29             :  * Although strictly speaking this is a "one stage" method, we treat
      30             :  * the "update" step as a second stage, since in finite element
      31             :  * analysis the update step requires a mass matrix solve.  The
      32             :  * implicit midpoint method can also be written in a single stage as:
      33             :  *
      34             :  * M*y_{n+1} = M*y_n + h*f(t_n + h/2, (y_{n+1} + y_n)/2)
      35             :  *
      36             :  * However this is less convenient to implement in MOOSE, since it
      37             :  * requires evaluating the time and non-time residuals on different
      38             :  * solution vectors, and the Jacobian contributions have an extra
      39             :  * factor of 1/2.
      40             :  */
      41             : class ImplicitMidpoint : public TimeIntegrator
      42             : {
      43             : public:
      44             :   static InputParameters validParams();
      45             : 
      46             :   ImplicitMidpoint(const InputParameters & parameters);
      47             : 
      48           0 :   virtual int order() override { return 2; }
      49             : 
      50             :   virtual void computeTimeDerivatives() override;
      51             :   void computeADTimeDerivatives(ADReal & ad_u_dot,
      52             :                                 const dof_id_type & dof,
      53             :                                 ADReal & ad_u_dotdot) const override;
      54             :   virtual void solve() override;
      55             :   virtual void postResidual(NumericVector<Number> & residual) override;
      56        2496 :   virtual bool overridesSolve() const override { return true; }
      57             : 
      58             : protected:
      59             :   /**
      60             :    * Helper function that actually does the math for computing the time derivative
      61             :    */
      62             :   template <typename T, typename T2>
      63             :   void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const;
      64             : 
      65             :   unsigned int _stage;
      66             : 
      67             :   /// Buffer to store non-time residual from the first stage.
      68             :   NumericVector<Number> * _residual_stage1;
      69             : };
      70             : 
      71             : template <typename T, typename T2>
      72             : void
      73       56339 : ImplicitMidpoint::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const
      74             : {
      75       56339 :   u_dot -= u_old;
      76       56339 :   u_dot *= 1. / _dt;
      77       56339 : }