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 "ExplicitTVDRK2.h"
      11             : #include "NonlinearSystemBase.h"
      12             : #include "FEProblem.h"
      13             : #include "PetscSupport.h"
      14             : 
      15             : registerMooseObject("MooseApp", ExplicitTVDRK2);
      16             : 
      17             : InputParameters
      18       14424 : ExplicitTVDRK2::validParams()
      19             : {
      20       14424 :   InputParameters params = TimeIntegrator::validParams();
      21       14424 :   params.addClassDescription("Explicit TVD (total-variation-diminishing) second-order Runge-Kutta "
      22             :                              "time integration method.");
      23       14424 :   return params;
      24           0 : }
      25             : 
      26         106 : ExplicitTVDRK2::ExplicitTVDRK2(const InputParameters & parameters)
      27             :   : TimeIntegrator(parameters),
      28         106 :     _stage(1),
      29         106 :     _residual_old(addVector("residual_old", false, libMesh::GHOSTED)),
      30         212 :     _solution_older(_sys.solutionState(2))
      31             : {
      32         106 :   mooseInfo("ExplicitTVDRK2 and other multistage TimeIntegrators are known not to work with "
      33             :             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
      34         106 : }
      35             : 
      36             : void
      37         436 : ExplicitTVDRK2::preSolve()
      38             : {
      39         436 :   if (_dt == _dt_old)
      40         436 :     _fe_problem.setConstJacobian(true);
      41             :   else
      42           0 :     _fe_problem.setConstJacobian(false);
      43         436 : }
      44             : 
      45             : void
      46        1462 : ExplicitTVDRK2::computeTimeDerivatives()
      47             : {
      48             :   // Since advanceState() is called in between stages 2 and 3, this
      49             :   // changes the meaning of "_solution_old".  In the second stage,
      50             :   // "_solution_older" is actually the original _solution_old.
      51        1462 :   if (!_sys.solutionUDot())
      52           0 :     mooseError("ExplicitTVDRK2: Time derivative of solution (`u_dot`) is not stored. Please set "
      53             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      54             : 
      55        1462 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      56        1462 :   u_dot = *_solution;
      57        1462 :   computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older);
      58             : 
      59        1462 :   computeDuDotDu();
      60        1462 :   u_dot.close();
      61        1462 : }
      62             : 
      63             : void
      64           0 : ExplicitTVDRK2::computeADTimeDerivatives(ADReal & ad_u_dot,
      65             :                                          const dof_id_type & dof,
      66             :                                          ADReal & /*ad_u_dotdot*/) const
      67             : {
      68           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof), _solution_older(dof));
      69           0 : }
      70             : 
      71             : void
      72         436 : ExplicitTVDRK2::solve()
      73             : {
      74         436 :   Real time_new = _fe_problem.time();
      75         436 :   Real time_old = _fe_problem.timeOld();
      76         436 :   Real time_stage2 = time_old + _dt;
      77             : 
      78             :   // Reset numbers of iterations
      79         436 :   _n_nonlinear_iterations = 0;
      80         436 :   _n_linear_iterations = 0;
      81             : 
      82             :   // There is no work to do for the first stage (Y_1 = y_n).  The
      83             :   // first solve therefore happens in the second stage.  Note that the
      84             :   // non-time Kernels (which should be marked implicit=false) are
      85             :   // evaluated at the old solution during this stage.
      86         436 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
      87         436 :   _console << "1st solve" << std::endl;
      88         436 :   _stage = 2;
      89         436 :   _fe_problem.timeOld() = time_old;
      90         436 :   _fe_problem.time() = time_stage2;
      91         436 :   _nl->system().solve();
      92         436 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
      93         436 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
      94             : 
      95             :   // Abort time step immediately on stage failure - see TimeIntegrator doc page
      96         436 :   if (!_fe_problem.converged(_nl->number()))
      97           0 :     return;
      98             : 
      99             :   // Advance solutions old->older, current->old.  Also moves Material
     100             :   // properties and other associated state forward in time.
     101         436 :   _fe_problem.advanceState();
     102             : 
     103             :   // The "update" stage (which we call stage 3) requires an additional
     104             :   // solve with the mass matrix.
     105         436 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
     106         436 :   _console << "2nd solve" << std::endl;
     107         436 :   _stage = 3;
     108         436 :   _fe_problem.timeOld() = time_stage2;
     109         436 :   _fe_problem.time() = time_new;
     110         436 :   _nl->system().solve();
     111         436 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     112         436 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
     113             : 
     114             :   // Reset time at beginning of step to its original value
     115         436 :   _fe_problem.timeOld() = time_old;
     116             : }
     117             : 
     118             : void
     119        1000 : ExplicitTVDRK2::postResidual(NumericVector<Number> & residual)
     120             : {
     121        1000 :   if (_stage == 1)
     122             :   {
     123             :     // If postResidual() is called before solve(), _stage==1 and we don't
     124             :     // need to do anything.
     125             :   }
     126        1000 :   else if (_stage == 2)
     127             :   {
     128             :     // In the standard RK notation, the stage 2 residual is given by:
     129             :     //
     130             :     // R := M*(Y_2 - y_n)/dt - f(t_n, Y_1) = 0
     131             :     //
     132             :     // where:
     133             :     // .) M is the mass matrix.
     134             :     // .) f(t_n, Y_1) is the residual we are currently computing,
     135             :     //    since this method is intended to be used with "implicit=false"
     136             :     //    kernels.
     137             :     // .) M*(Y_2 - y_n)/dt corresponds to the residual of the time kernels.
     138             :     // .) The minus signs are "baked in" to the non-time residuals, so
     139             :     //    they do not appear here.
     140             :     // .) The current non-time residual is saved for the next stage.
     141         564 :     *_residual_old = *_Re_non_time;
     142         564 :     _residual_old->close();
     143             : 
     144         564 :     residual.add(1.0, *_Re_time);
     145         564 :     residual.add(1.0, *_residual_old);
     146         564 :     residual.close();
     147             :   }
     148         436 :   else if (_stage == 3)
     149             :   {
     150             :     // In the standard RK notation, the update step residual is given by:
     151             :     //
     152             :     // R := M*(2*y_{n+1} - Y_2 - y_n)/(2*dt) - (1/2)*f(t_n+dt/2, Y_2) = 0
     153             :     //
     154             :     // where:
     155             :     // .) M is the mass matrix.
     156             :     // .) f(t_n+dt/2, Y_2) is the residual from stage 2.
     157             :     // .) The minus sign is already baked in to the non-time
     158             :     //    residuals, so it does not appear here.
     159             :     // .) Although this is an update step, we have to do a "solve"
     160             :     //    using the mass matrix.
     161         436 :     residual.add(1.0, *_Re_time);
     162         436 :     residual.add(0.5, *_Re_non_time);
     163         436 :     residual.close();
     164             :   }
     165             :   else
     166           0 :     mooseError(
     167           0 :         "ExplicitTVDRK2::postResidual(): _stage = ", _stage, ", only _stage = 1-3 is allowed.");
     168        1000 : }