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             : #include "LStableDirk2.h"
      11             : #include "NonlinearSystem.h"
      12             : #include "FEProblem.h"
      13             : #include "PetscSupport.h"
      14             : 
      15             : using namespace libMesh;
      16             : 
      17             : registerMooseObject("MooseApp", LStableDirk2);
      18             : 
      19             : InputParameters
      20       14484 : LStableDirk2::validParams()
      21             : {
      22       14484 :   InputParameters params = TimeIntegrator::validParams();
      23       14484 :   params.addClassDescription(
      24             :       "Second order diagonally implicit Runge Kutta method (Dirk) with two stages.");
      25       14484 :   return params;
      26           0 : }
      27             : 
      28         146 : LStableDirk2::LStableDirk2(const InputParameters & parameters)
      29             :   : TimeIntegrator(parameters),
      30         146 :     _stage(1),
      31         146 :     _residual_stage1(addVector("residual_stage1", false, GHOSTED)),
      32         146 :     _residual_stage2(addVector("residual_stage2", false, GHOSTED)),
      33         292 :     _alpha(1. - 0.5 * std::sqrt(2))
      34             : {
      35         146 :   mooseInfo("LStableDirk2 and other multistage TimeIntegrators are known not to work with "
      36             :             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
      37         146 : }
      38             : 
      39             : void
      40        7327 : LStableDirk2::computeTimeDerivatives()
      41             : {
      42             :   // We are multiplying by the method coefficients in postResidual(), so
      43             :   // the time derivatives are of the same form at every stage although
      44             :   // the current solution varies depending on the stage.
      45        7327 :   if (!_sys.solutionUDot())
      46           0 :     mooseError("LStableDirk2: Time derivative of solution (`u_dot`) is not stored. Please set "
      47             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      48             : 
      49        7327 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      50        7327 :   u_dot = *_solution;
      51        7327 :   computeTimeDerivativeHelper(u_dot, _solution_old);
      52        7327 :   u_dot.close();
      53        7327 :   computeDuDotDu();
      54        7327 : }
      55             : 
      56             : void
      57           0 : LStableDirk2::computeADTimeDerivatives(ADReal & ad_u_dot,
      58             :                                        const dof_id_type & dof,
      59             :                                        ADReal & /*ad_u_dotdot*/) const
      60             : {
      61           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
      62           0 : }
      63             : 
      64             : void
      65         231 : LStableDirk2::solve()
      66             : {
      67             :   // Time at end of step
      68         231 :   Real time_new = _fe_problem.time();
      69             : 
      70             :   // Time at beginning of step
      71         231 :   Real time_old = _fe_problem.timeOld();
      72             : 
      73             :   // Time at stage 1
      74         231 :   Real time_stage1 = time_old + _alpha * _dt;
      75             : 
      76             :   // Reset iteration counts
      77         231 :   _n_nonlinear_iterations = 0;
      78         231 :   _n_linear_iterations = 0;
      79             : 
      80             :   // Compute first stage
      81         231 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
      82         231 :   _console << "1st stage" << std::endl;
      83         231 :   _stage = 1;
      84         231 :   _fe_problem.time() = time_stage1;
      85         231 :   _nl->system().solve();
      86         231 :   _nl->destroyColoring();
      87         231 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
      88         231 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
      89             : 
      90             :   // Abort time step immediately on stage failure - see TimeIntegrator doc page
      91         231 :   if (!_fe_problem.converged(_nl->number()))
      92          10 :     return;
      93             : 
      94             :   // Compute second stage
      95         221 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
      96         221 :   _console << "2nd stage" << std::endl;
      97         221 :   _stage = 2;
      98         221 :   _fe_problem.timeOld() = time_stage1;
      99         221 :   _fe_problem.time() = time_new;
     100         221 :   _nl->potentiallySetupFiniteDifferencing();
     101         221 :   _nl->system().solve();
     102         221 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     103         221 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
     104             : 
     105             :   // Reset time at beginning of step to its original value
     106         221 :   _fe_problem.timeOld() = time_old;
     107             : }
     108             : 
     109             : void
     110        7305 : LStableDirk2::postResidual(NumericVector<Number> & residual)
     111             : {
     112        7305 :   if (_stage == 1)
     113             :   {
     114             :     // In the standard RK notation, the stage 1 residual is given by:
     115             :     //
     116             :     // R := (Y_1 - y_n)/dt - alpha*f(t_n + alpha*dt, Y_1) = 0
     117             :     //
     118             :     // where:
     119             :     // .) f(t_n + alpha*dt, Y_1) corresponds to the residuals of the
     120             :     //    non-time kernels.  We'll save this as "_residual_stage1" to use
     121             :     //    later.
     122             :     // .) (Y_1 - y_n)/dt corresponds to the residual of the time kernels.
     123             :     // .) The minus sign in front of alpha is already "baked in" to
     124             :     //    the non-time residuals, so it does not appear here.
     125        3690 :     *_residual_stage1 = *_Re_non_time;
     126        3690 :     _residual_stage1->close();
     127             : 
     128        3690 :     residual.add(1., *_Re_time);
     129        3690 :     residual.add(_alpha, *_residual_stage1);
     130        3690 :     residual.close();
     131             :   }
     132        3615 :   else if (_stage == 2)
     133             :   {
     134             :     // In the standard RK notation, the stage 2 residual is given by:
     135             :     //
     136             :     // R := (Y_2 - y_n)/dt - (1-alpha)*f(t_n + alpha*dt, Y_1) - alpha*f(t_n + dt, Y_2) = 0
     137             :     //
     138             :     // where:
     139             :     // .) f(t_n + alpha*dt, Y_1) has already been saved as _residual_stage1.
     140             :     // .) f(t_n + dt, Y_2) will now be saved as "_residual_stage2".
     141             :     // .) (Y_2 - y_n)/dt corresponds to the residual of the time kernels.
     142             :     // .) The minus signs are once again "baked in" to the non-time
     143             :     //    residuals, so they do not appear here.
     144             :     //
     145             :     // The solution at the end of stage 2, i.e. Y_2, is also the final solution.
     146        3615 :     *_residual_stage2 = *_Re_non_time;
     147        3615 :     _residual_stage2->close();
     148             : 
     149        3615 :     residual.add(1., *_Re_time);
     150        3615 :     residual.add(1. - _alpha, *_residual_stage1);
     151        3615 :     residual.add(_alpha, *_residual_stage2);
     152        3615 :     residual.close();
     153             :   }
     154             :   else
     155           0 :     mooseError("LStableDirk2::postResidual(): Member variable _stage can only have values 1 or 2.");
     156        7305 : }