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 "LStableDirk3.h"
      11             : #include "NonlinearSystem.h"
      12             : #include "FEProblem.h"
      13             : #include "PetscSupport.h"
      14             : using namespace libMesh;
      15             : 
      16             : registerMooseObject("MooseApp", LStableDirk3);
      17             : 
      18             : InputParameters
      19       14436 : LStableDirk3::validParams()
      20             : {
      21       14436 :   InputParameters params = TimeIntegrator::validParams();
      22       14436 :   params.addClassDescription(
      23             :       "Third order diagonally implicit Runge Kutta method (Dirk) with three stages.");
      24       14436 :   return params;
      25           0 : }
      26             : 
      27         114 : LStableDirk3::LStableDirk3(const InputParameters & parameters)
      28             :   : TimeIntegrator(parameters),
      29         114 :     _stage(1),
      30         114 :     _gamma(-std::sqrt(2.) * std::cos(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. +
      31         114 :            std::sqrt(6.) * std::sin(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. + 1.)
      32             : {
      33         114 :   mooseInfo("LStableDirk3 and other multistage TimeIntegrators are known not to work with "
      34             :             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
      35             : 
      36             :   // Name the stage residuals "residual_stage1", "residual_stage2", etc.
      37         456 :   for (unsigned int stage = 0; stage < 3; ++stage)
      38             :   {
      39         342 :     std::ostringstream oss;
      40         342 :     oss << "residual_stage" << stage + 1;
      41         342 :     _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
      42         342 :   }
      43             : 
      44             :   // Initialize parameters
      45         114 :   _c[0] = _gamma;
      46         114 :   _c[1] = .5 * (1 + _gamma);
      47         114 :   _c[2] = 1.0;
      48             : 
      49         114 :   _a[0][0] = _gamma;
      50         114 :   _a[1][0] = .5 * (1 - _gamma); /**/
      51         114 :   _a[1][1] = _gamma;
      52         114 :   _a[2][0] = .25 * (-6 * _gamma * _gamma + 16 * _gamma - 1); /**/
      53         114 :   _a[2][1] = .25 * (6 * _gamma * _gamma - 20 * _gamma + 5);  /**/
      54         114 :   _a[2][2] = _gamma;
      55         114 : }
      56             : 
      57             : void
      58       65992 : LStableDirk3::computeTimeDerivatives()
      59             : {
      60             :   // We are multiplying by the method coefficients in postResidual(), so
      61             :   // the time derivatives are of the same form at every stage although
      62             :   // the current solution varies depending on the stage.
      63       65992 :   if (!_sys.solutionUDot())
      64           0 :     mooseError("LStableDirk3: Time derivative of solution (`u_dot`) is not stored. Please set "
      65             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      66             : 
      67       65992 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      68       65992 :   u_dot = *_solution;
      69       65992 :   computeTimeDerivativeHelper(u_dot, _solution_old);
      70       65992 :   u_dot.close();
      71       65992 :   computeDuDotDu();
      72       65992 : }
      73             : 
      74             : void
      75           0 : LStableDirk3::computeADTimeDerivatives(ADReal & ad_u_dot,
      76             :                                        const dof_id_type & dof,
      77             :                                        ADReal & /*ad_u_dotdot*/) const
      78             : {
      79           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
      80           0 : }
      81             : 
      82             : void
      83        1259 : LStableDirk3::solve()
      84             : {
      85             :   // Time at end of step
      86        1259 :   Real time_old = _fe_problem.timeOld();
      87             : 
      88             :   // Reset iteration counts
      89        1259 :   _n_nonlinear_iterations = 0;
      90        1259 :   _n_linear_iterations = 0;
      91             : 
      92             :   // A for-loop would increment _stage too far, so we use an extra
      93             :   // loop counter.
      94        5036 :   for (unsigned int current_stage = 1; current_stage < 4; ++current_stage)
      95             :   {
      96             :     // Set the current stage value
      97        3777 :     _stage = current_stage;
      98             : 
      99             :     // This ensures that all the Output objects in the OutputWarehouse
     100             :     // have had solveSetup() called, and sets the default solver
     101             :     // parameters for PETSc.
     102        3777 :     _fe_problem.initPetscOutputAndSomeSolverSettings();
     103             : 
     104        3777 :     _console << "Stage " << _stage << std::endl;
     105             : 
     106             :     // Set the time for this stage
     107        3777 :     _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
     108             : 
     109             :     // If we previously used coloring, destroy the old object so it doesn't leak when we allocate a
     110             :     // new object in the following lines
     111        3777 :     _nl->destroyColoring();
     112             : 
     113             :     // Potentially setup finite differencing contexts for the solve
     114        3777 :     _nl->potentiallySetupFiniteDifferencing();
     115             : 
     116             :     // Do the solve
     117        3777 :     _nl->system().solve();
     118             : 
     119             :     // Update the iteration counts
     120        3777 :     _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     121        3777 :     _n_linear_iterations += getNumLinearIterationsLastSolve();
     122             : 
     123             :     // Abort time step immediately on stage failure - see TimeIntegrator doc page
     124        3777 :     if (!_fe_problem.converged(_nl->number()))
     125           0 :       return;
     126             :   }
     127             : }
     128             : 
     129             : void
     130       65948 : LStableDirk3::postResidual(NumericVector<Number> & residual)
     131             : {
     132             :   // Error if _stage got messed up somehow.
     133       65948 :   if (_stage > 3)
     134           0 :     mooseError(
     135             :         "LStableDirk3::postResidual(): Member variable _stage can only have values 1, 2, or 3.");
     136             : 
     137             :   // In the standard RK notation, the residual of stage 1 of s is given by:
     138             :   //
     139             :   // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
     140             :   //
     141             :   // where:
     142             :   // .) M is the mass matrix
     143             :   // .) Y_i is the stage solution
     144             :   // .) dt is the timestep, and is accounted for in the _Re_time residual.
     145             :   // .) f are the "non-time" residuals evaluated for a given stage solution.
     146             :   // .) The minus signs are already "baked in" to the residuals and so do not appear below.
     147             : 
     148             :   // Store this stage's non-time residual.  We are calling operator=
     149             :   // here, and that calls close().
     150       65948 :   *_stage_residuals[_stage - 1] = *_Re_non_time;
     151             : 
     152             :   // Build up the residual for this stage.
     153       65948 :   residual.add(1., *_Re_time);
     154      197878 :   for (unsigned int j = 0; j < _stage; ++j)
     155      131930 :     residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
     156       65948 :   residual.close();
     157       65948 : }