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