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 "ImplicitMidpoint.h"
      11             : #include "NonlinearSystem.h"
      12             : #include "FEProblem.h"
      13             : #include "PetscSupport.h"
      14             : 
      15             : registerMooseObject("MooseApp", ImplicitMidpoint);
      16             : 
      17             : InputParameters
      18       14400 : ImplicitMidpoint::validParams()
      19             : {
      20       14400 :   InputParameters params = TimeIntegrator::validParams();
      21       14400 :   params.addClassDescription("Second-order Runge-Kutta (implicit midpoint) time integration.");
      22       14400 :   return params;
      23           0 : }
      24             : 
      25          90 : ImplicitMidpoint::ImplicitMidpoint(const InputParameters & parameters)
      26             :   : TimeIntegrator(parameters),
      27          90 :     _stage(1),
      28          90 :     _residual_stage1(addVector("residual_stage1", false, libMesh::GHOSTED))
      29             : {
      30          90 :   mooseInfo("ImplicitMidpoint and other multistage TimeIntegrators are known not to work with "
      31             :             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
      32          90 : }
      33             : 
      34             : void
      35       56339 : ImplicitMidpoint::computeTimeDerivatives()
      36             : {
      37             :   // We are multiplying by the method coefficients in postResidual(), so
      38             :   // the time derivatives are of the same form at every stage although
      39             :   // the current solution varies depending on the stage.
      40       56339 :   if (!_sys.solutionUDot())
      41           0 :     mooseError("ImplicitMidpoint: Time derivative of solution (`u_dot`) is not stored. Please set "
      42             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      43             : 
      44       56339 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      45       56339 :   u_dot = *_solution;
      46       56339 :   computeTimeDerivativeHelper(u_dot, _solution_old);
      47       56339 :   u_dot.close();
      48       56339 :   computeDuDotDu();
      49       56339 : }
      50             : 
      51             : void
      52           0 : ImplicitMidpoint::computeADTimeDerivatives(ADReal & ad_u_dot,
      53             :                                            const dof_id_type & dof,
      54             :                                            ADReal & /*ad_u_dotdot*/) const
      55             : {
      56           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
      57           0 : }
      58             : 
      59             : void
      60        1248 : ImplicitMidpoint::solve()
      61             : {
      62        1248 :   Real time_new = _fe_problem.time();
      63        1248 :   Real time_old = _fe_problem.timeOld();
      64        1248 :   Real time_half = (time_new + time_old) / 2.;
      65             : 
      66             :   // Reset iteration counts
      67        1248 :   _n_nonlinear_iterations = 0;
      68        1248 :   _n_linear_iterations = 0;
      69             : 
      70             :   // Compute first stage
      71        1248 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
      72        1248 :   _console << "1st stage" << std::endl;
      73        1248 :   _stage = 1;
      74        1248 :   _fe_problem.time() = time_half;
      75        1248 :   _nl->system().solve();
      76        1248 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
      77        1248 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
      78             : 
      79             :   // Abort time step immediately on stage failure - see TimeIntegrator doc page
      80        1248 :   if (!_fe_problem.converged(_nl->number()))
      81           0 :     return;
      82             : 
      83             :   // Compute second stage
      84        1248 :   _fe_problem.initPetscOutputAndSomeSolverSettings();
      85        1248 :   _console << "2nd stage" << std::endl;
      86        1248 :   _stage = 2;
      87        1248 :   _fe_problem.time() = time_new;
      88        1248 :   _nl->system().solve();
      89        1248 :   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
      90        1248 :   _n_linear_iterations += getNumLinearIterationsLastSolve();
      91             : }
      92             : 
      93             : void
      94       56339 : ImplicitMidpoint::postResidual(NumericVector<Number> & residual)
      95             : {
      96       56339 :   if (_stage == 1)
      97             :   {
      98             :     // In the standard RK notation, the stage 1 residual is given by:
      99             :     //
     100             :     // R := M*(Y_1 - y_n)/dt - (1/2)*f(t_n + dt/2, Y_1) = 0
     101             :     //
     102             :     // where:
     103             :     // .) M is the mass matrix
     104             :     // .) f(t_n + dt/2, Y_1) is saved in _residual_stage1
     105             :     // .) The minus sign is baked in to the non-time residuals, so it does not appear here.
     106       20396 :     *_residual_stage1 = *_Re_non_time;
     107       20396 :     _residual_stage1->close();
     108             : 
     109       20396 :     residual.add(1., *_Re_time);
     110       20396 :     residual.add(0.5, *_Re_non_time);
     111       20396 :     residual.close();
     112             :   }
     113       35943 :   else if (_stage == 2)
     114             :   {
     115             :     // The update step.  In the standard RK notation, the update
     116             :     // residual is given by:
     117             :     //
     118             :     // R := M*(y_{n+1} - y_n)/dt - f(t_n + dt/2, Y_1) = 0
     119             :     //
     120             :     // where:
     121             :     // .) M is the mass matrix.
     122             :     // .) f(t_n + dt/2, Y_1) is the residual from stage 1, it has already
     123             :     //    been saved as _residual_stage1.
     124             :     // .) The minus signs are "baked in" to the non-time residuals, so
     125             :     //    they do not appear here.
     126       35943 :     residual.add(1., *_Re_time);
     127       35943 :     residual.add(1., *_residual_stage1);
     128       35943 :     residual.close();
     129             :   }
     130             :   else
     131           0 :     mooseError(
     132           0 :         "ImplicitMidpoint::postResidual(): _stage = ", _stage, ", only _stage = 1, 2 is allowed.");
     133       56339 : }