https://mooseframework.inl.gov
LStableDirk2.C
Go to the documentation of this file.
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 "LStableDirk2.h"
11 #include "NonlinearSystem.h"
12 #include "FEProblem.h"
13 #include "PetscSupport.h"
14 
15 using namespace libMesh;
16 
18 
21 {
23  params.addClassDescription(
24  "Second order diagonally implicit Runge Kutta method (Dirk) with two stages.");
25  return params;
26 }
27 
29  : TimeIntegrator(parameters),
30  _stage(1),
31  _residual_stage1(addVector("residual_stage1", false, GHOSTED)),
32  _residual_stage2(addVector("residual_stage2", false, GHOSTED)),
33  _alpha(1. - 0.5 * std::sqrt(2))
34 {
35  mooseInfo("LStableDirk2 and other multistage TimeIntegrators are known not to work with "
36  "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
37 }
38 
39 void
41 {
42  // We are multiplying by the method coefficients in postResidual(), so
43  // the time derivatives are of the same form at every stage although
44  // the current solution varies depending on the stage.
45  if (!_sys.solutionUDot())
46  mooseError("LStableDirk2: Time derivative of solution (`u_dot`) is not stored. Please set "
47  "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
48 
50  u_dot = *_solution;
52  u_dot.close();
54 }
55 
56 void
58  const dof_id_type & dof,
59  ADReal & /*ad_u_dotdot*/) const
60 {
62 }
63 
64 void
66 {
67  // Time at end of step
68  Real time_new = _fe_problem.time();
69 
70  // Time at beginning of step
71  Real time_old = _fe_problem.timeOld();
72 
73  // Time at stage 1
74  Real time_stage1 = time_old + _alpha * _dt;
75 
76  // Reset iteration counts
79 
80  // Compute first stage
82  _console << "1st stage" << std::endl;
83  _stage = 1;
84  _fe_problem.time() = time_stage1;
85  _nl->system().solve();
89 
90  // Abort time step immediately on stage failure - see TimeIntegrator doc page
92  return;
93 
94  // Compute second stage
96  _console << "2nd stage" << std::endl;
97  _stage = 2;
98  _fe_problem.timeOld() = time_stage1;
99  _fe_problem.time() = time_new;
101  _nl->system().solve();
104 
105  // Reset time at beginning of step to its original value
106  _fe_problem.timeOld() = time_old;
107 }
108 
109 void
111 {
112  if (_stage == 1)
113  {
114  // In the standard RK notation, the stage 1 residual is given by:
115  //
116  // R := (Y_1 - y_n)/dt - alpha*f(t_n + alpha*dt, Y_1) = 0
117  //
118  // where:
119  // .) f(t_n + alpha*dt, Y_1) corresponds to the residuals of the
120  // non-time kernels. We'll save this as "_residual_stage1" to use
121  // later.
122  // .) (Y_1 - y_n)/dt corresponds to the residual of the time kernels.
123  // .) The minus sign in front of alpha is already "baked in" to
124  // the non-time residuals, so it does not appear here.
127 
128  residual.add(1., *_Re_time);
129  residual.add(_alpha, *_residual_stage1);
130  residual.close();
131  }
132  else if (_stage == 2)
133  {
134  // In the standard RK notation, the stage 2 residual is given by:
135  //
136  // R := (Y_2 - y_n)/dt - (1-alpha)*f(t_n + alpha*dt, Y_1) - alpha*f(t_n + dt, Y_2) = 0
137  //
138  // where:
139  // .) f(t_n + alpha*dt, Y_1) has already been saved as _residual_stage1.
140  // .) f(t_n + dt, Y_2) will now be saved as "_residual_stage2".
141  // .) (Y_2 - y_n)/dt corresponds to the residual of the time kernels.
142  // .) The minus signs are once again "baked in" to the non-time
143  // residuals, so they do not appear here.
144  //
145  // The solution at the end of stage 2, i.e. Y_2, is also the final solution.
148 
149  residual.add(1., *_Re_time);
150  residual.add(1. - _alpha, *_residual_stage1);
151  residual.add(_alpha, *_residual_stage2);
152  residual.close();
153  }
154  else
155  mooseError("LStableDirk2::postResidual(): Member variable _stage can only have values 1 or 2.");
156 }
virtual void initPetscOutputAndSomeSolverSettings()
Reinitialize PETSc output for proper linear/nonlinear iteration display.
virtual Real & time() const
virtual void computeTimeDerivatives() override
Computes the time derivative and the Jacobian of the time derivative.
Definition: LStableDirk2.C:40
FEProblemBase & _fe_problem
Reference to the problem.
void mooseInfo(Args &&... args) const
LStableDirk2(const InputParameters &parameters)
Definition: LStableDirk2.C:28
NonlinearSystemBase * _nl
Pointer to the nonlinear system, can happen that we dont have any.
const Real _alpha
Definition: LStableDirk2.h:71
static InputParameters validParams()
Definition: LStableDirk2.C:20
SystemBase & _sys
Reference to the system this time integrator operates on.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
virtual void solve() override
Solves the time step and sets the number of nonlinear and linear iterations.
Definition: LStableDirk2.C:65
DualNumber< Real, DNDerivativeType, true > ADReal
Definition: ADRealForward.h:47
Real & _dt
The current time step size.
NumericVector< Number > * _Re_time
residual vector for time contributions
void computeADTimeDerivatives(ADReal &ad_u_dot, const dof_id_type &dof, ADReal &ad_u_dotdot) const override
method for computing local automatic differentiation time derivatives
Definition: LStableDirk2.C:57
virtual NumericVector< Number > * solutionUDot()
Definition: SystemBase.h:252
virtual bool converged(const unsigned int sys_num)
Eventually we want to convert this virtual over to taking a solver system number argument.
Definition: SubProblem.h:113
unsigned int _n_linear_iterations
Total number of linear iterations over all stages of the time step.
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
NumericVector< Number > * _residual_stage1
Buffer to store non-time residual from first stage solve.
Definition: LStableDirk2.h:65
void destroyColoring()
Destroy the coloring object if it exists.
virtual void solve()
virtual void close()=0
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the NonLinearTimeIntegratorInterface called immediately after the residuals are computed ...
Definition: LStableDirk2.C:110
void computeTimeDerivativeHelper(T &u_dot, const T2 &u_old) const
Helper function that actually does the math for computing the time derivative.
Definition: LStableDirk2.h:76
const NumericVector< Number > *const & _solution
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Base class for time integrators.
unsigned int getNumLinearIterationsLastSolve() const
Gets the number of linear iterations in the most recent solve.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
NumericVector< Number > * _residual_stage2
Buffer to store non-time residual from second stage solve.
Definition: LStableDirk2.h:68
unsigned int _stage
Indicates the current stage (1 or 2).
Definition: LStableDirk2.h:62
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
virtual Real & timeOld() const
const NumericVector< Number > & _solution_old
unsigned int _n_nonlinear_iterations
Total number of nonlinear iterations over all stages of the time step.
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
unsigned int getNumNonlinearIterationsLastSolve() const
Gets the number of nonlinear iterations in the most recent solve.
Second order diagonally implicit Runge Kutta method (Dirk) with two stages.
Definition: LStableDirk2.h:38
virtual void add(const numeric_index_type i, const T value)=0
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
registerMooseObject("MooseApp", LStableDirk2)
virtual void potentiallySetupFiniteDifferencing()
Create finite differencing contexts for assembly of the Jacobian and/or approximating the action of t...
static InputParameters validParams()
uint8_t dof_id_type
void computeDuDotDu()
Compute _du_dot_du.
virtual libMesh::System & system() override
Get the reference to the libMesh system.