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             : // MOOSE includes
      11             : #include "AStableDirk4.h"
      12             : #include "NonlinearSystem.h"
      13             : #include "FEProblem.h"
      14             : #include "PetscSupport.h"
      15             : #include "LStableDirk4.h"
      16             : 
      17             : using namespace libMesh;
      18             : 
      19             : registerMooseObject("MooseApp", AStableDirk4);
      20             : 
      21             : InputParameters
      22       14508 : AStableDirk4::validParams()
      23             : {
      24       14508 :   InputParameters params = TimeIntegrator::validParams();
      25       14508 :   params.addClassDescription("Fourth-order diagonally implicit Runge Kutta method (Dirk) with "
      26             :                              "three stages plus an update.");
      27       14508 :   params.addParam<bool>("safe_start", true, "If true, use LStableDirk4 to bootstrap this method.");
      28       14508 :   return params;
      29           0 : }
      30             : 
      31         162 : AStableDirk4::AStableDirk4(const InputParameters & parameters)
      32             :   : TimeIntegrator(parameters),
      33         162 :     _stage(1),
      34         324 :     _gamma(0.5 + std::sqrt(3) / 3. * std::cos(libMesh::pi / 18.)),
      35         162 :     _safe_start(getParam<bool>("safe_start"))
      36             : {
      37         162 :   mooseInfo("AStableDirk4 and other multistage TimeIntegrators are known not to work with "
      38             :             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
      39             : 
      40             :   // Name the stage residuals "residual_stage1", "residual_stage2", etc.
      41         648 :   for (unsigned int stage = 0; stage < 3; ++stage)
      42             :   {
      43         486 :     std::ostringstream oss;
      44         486 :     oss << "residual_stage" << stage + 1;
      45         486 :     _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
      46         486 :   }
      47             : 
      48             :   // Initialize parameters
      49         162 :   _c[0] = _gamma;
      50         162 :   _c[1] = .5;
      51         162 :   _c[2] = 1.0 - _gamma;
      52             : 
      53         162 :   _a[0][0] = _gamma;
      54         162 :   _a[1][0] = .5 - _gamma; /**/
      55         162 :   _a[1][1] = _gamma;
      56         162 :   _a[2][0] = 2. * _gamma;     /**/
      57         162 :   _a[2][1] = 1 - 4. * _gamma; /**/
      58         162 :   _a[2][2] = _gamma;
      59             : 
      60         162 :   _b[0] = 1. / (24. * (.5 - _gamma) * (.5 - _gamma));
      61         162 :   _b[1] = 1. - 1. / (12. * (.5 - _gamma) * (.5 - _gamma));
      62         162 :   _b[2] = _b[0];
      63             : 
      64             :   // If doing a _safe_start, construct the bootstrapping
      65             :   // TimeIntegrator.  Note that this method will also add
      66             :   // residual_stage vectors to the system, but since they have the
      67             :   // same name as the vectors we added, they won't be duplicated.
      68         162 :   if (_safe_start)
      69             :   {
      70          90 :     Factory & factory = _app.getFactory();
      71          90 :     InputParameters params = factory.getValidParams("LStableDirk4");
      72             : 
      73             :     // We need to set some parameters that are normally set in
      74             :     // FEProblemBase::addTimeIntegrator() to ensure that the
      75             :     // getCheckedPointerParam() sanity checking is happy.  This is why
      76             :     // constructing MOOSE objects "manually" is generally frowned upon.
      77          90 :     params.set<SystemBase *>("_sys") = &_sys;
      78             : 
      79          90 :     _bootstrap_method = factory.create<LStableDirk4>("LStableDirk4", name() + "_bootstrap", params);
      80          90 :   }
      81         162 : }
      82             : 
      83             : void
      84       89874 : AStableDirk4::computeTimeDerivatives()
      85             : {
      86             :   // We are multiplying by the method coefficients in postResidual(), so
      87             :   // the time derivatives are of the same form at every stage although
      88             :   // the current solution varies depending on the stage.
      89       89874 :   if (!_sys.solutionUDot())
      90           0 :     mooseError("AStableDirk4: Time derivative of solution (`u_dot`) is not stored. Please set "
      91             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      92             : 
      93       89874 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      94       89874 :   u_dot = *_solution;
      95       89874 :   computeTimeDerivativeHelper(u_dot, _solution_old);
      96       89874 :   u_dot.close();
      97       89874 :   computeDuDotDu();
      98       89874 : }
      99             : 
     100             : void
     101           0 : AStableDirk4::computeADTimeDerivatives(ADReal & ad_u_dot,
     102             :                                        const dof_id_type & dof,
     103             :                                        ADReal & /*ad_u_dotdot*/) const
     104             : {
     105           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
     106           0 : }
     107             : 
     108             : void
     109        1402 : AStableDirk4::solve()
     110             : {
     111             :   // Reset iteration counts
     112        1402 :   _n_nonlinear_iterations = 0;
     113        1402 :   _n_linear_iterations = 0;
     114             : 
     115        1402 :   if (_t_step == 1 && _safe_start)
     116             :   {
     117          41 :     _bootstrap_method->solve();
     118          41 :     _n_nonlinear_iterations = _bootstrap_method->getNumNonlinearIterations();
     119          41 :     _n_linear_iterations = _bootstrap_method->getNumLinearIterations();
     120             :   }
     121             :   else
     122             :   {
     123             :     // Time at end of step
     124        1361 :     Real time_old = _fe_problem.timeOld();
     125             : 
     126             :     // A for-loop would increment _stage too far, so we use an extra
     127             :     // loop counter.  The method has three stages and an update stage,
     128             :     // which we treat as just one more (special) stage in the implementation.
     129        6805 :     for (unsigned int current_stage = 1; current_stage < 5; ++current_stage)
     130             :     {
     131             :       // Set the current stage value
     132        5444 :       _stage = current_stage;
     133             : 
     134             :       // This ensures that all the Output objects in the OutputWarehouse
     135             :       // have had solveSetup() called, and sets the default solver
     136             :       // parameters for PETSc.
     137        5444 :       _fe_problem.initPetscOutputAndSomeSolverSettings();
     138             : 
     139        5444 :       if (current_stage < 4)
     140             :       {
     141        4083 :         _console << "Stage " << _stage << std::endl;
     142        4083 :         _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
     143             :       }
     144             :       else
     145             :       {
     146        1361 :         _console << "Update Stage." << std::endl;
     147        1361 :         _fe_problem.time() = time_old + _dt;
     148             :       }
     149             : 
     150             :       // Do the solve
     151        5444 :       _nl->system().solve();
     152             : 
     153             :       // Update the iteration counts
     154        5444 :       _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     155        5444 :       _n_linear_iterations += getNumLinearIterationsLastSolve();
     156             : 
     157             :       // Abort time step immediately on stage failure - see TimeIntegrator doc page
     158        5444 :       if (!_fe_problem.converged(_nl->number()))
     159           0 :         return;
     160             :     }
     161             :   }
     162             : }
     163             : 
     164             : void
     165       89874 : AStableDirk4::postResidual(NumericVector<Number> & residual)
     166             : {
     167       89874 :   if (_t_step == 1 && _safe_start)
     168        3406 :     _bootstrap_method->postResidual(residual);
     169             : 
     170             :   else
     171             :   {
     172             :     // Error if _stage got messed up somehow.
     173       86468 :     if (_stage > 4)
     174           0 :       mooseError("AStableDirk4::postResidual(): Member variable _stage can only have values 1-4.");
     175             : 
     176       86468 :     if (_stage < 4)
     177             :     {
     178             :       // In the standard RK notation, the residual of stage 1 of s is given by:
     179             :       //
     180             :       // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
     181             :       //
     182             :       // where:
     183             :       // .) M is the mass matrix
     184             :       // .) Y_i is the stage solution
     185             :       // .) dt is the timestep, and is accounted for in the _Re_time residual.
     186             :       // .) f are the "non-time" residuals evaluated for a given stage solution.
     187             :       // .) The minus signs are already "baked in" to the residuals and so do not appear below.
     188             : 
     189             :       // Store this stage's non-time residual.  We are calling operator=
     190             :       // here, and that calls close().
     191       46657 :       *_stage_residuals[_stage - 1] = *_Re_non_time;
     192             : 
     193             :       // Build up the residual for this stage.
     194       46657 :       residual.add(1., *_Re_time);
     195      139917 :       for (unsigned int j = 0; j < _stage; ++j)
     196       93260 :         residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
     197       46657 :       residual.close();
     198             :     }
     199             :     else
     200             :     {
     201             :       // The update step is a final solve:
     202             :       //
     203             :       // R := M*(y_{n+1} - y_n)/dt - \sum_{j=1}^s b_j * f(t_n + c_j*dt, Y_j) = 0
     204             :       //
     205             :       // We could potentially fold _b up into an extra row of _a and
     206             :       // just do one more stage, but I think handling it separately like
     207             :       // this is easier to understand and doesn't create too much code
     208             :       // repitition.
     209       39811 :       residual.add(1., *_Re_time);
     210      159244 :       for (unsigned int j = 0; j < 3; ++j)
     211      119433 :         residual.add(_b[j], *_stage_residuals[j]);
     212       39811 :       residual.close();
     213             :     }
     214             :   }
     215       89874 : }