https://mooseframework.inl.gov
AStableDirk4.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 // 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 
20 
23 {
25  params.addClassDescription("Fourth-order diagonally implicit Runge Kutta method (Dirk) with "
26  "three stages plus an update.");
27  params.addParam<bool>("safe_start", true, "If true, use LStableDirk4 to bootstrap this method.");
28  return params;
29 }
30 
32  : TimeIntegrator(parameters),
33  _stage(1),
34  _gamma(0.5 + std::sqrt(3) / 3. * std::cos(libMesh::pi / 18.)),
35  _safe_start(getParam<bool>("safe_start"))
36 {
37  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  for (unsigned int stage = 0; stage < 3; ++stage)
42  {
43  std::ostringstream oss;
44  oss << "residual_stage" << stage + 1;
45  _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
46  }
47 
48  // Initialize parameters
49  _c[0] = _gamma;
50  _c[1] = .5;
51  _c[2] = 1.0 - _gamma;
52 
53  _a[0][0] = _gamma;
54  _a[1][0] = .5 - _gamma;
55  _a[1][1] = _gamma;
56  _a[2][0] = 2. * _gamma;
57  _a[2][1] = 1 - 4. * _gamma;
58  _a[2][2] = _gamma;
59 
60  _b[0] = 1. / (24. * (.5 - _gamma) * (.5 - _gamma));
61  _b[1] = 1. - 1. / (12. * (.5 - _gamma) * (.5 - _gamma));
62  _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  if (_safe_start)
69  {
70  Factory & factory = _app.getFactory();
71  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  params.set<SystemBase *>("_sys") = &_sys;
78 
79  _bootstrap_method = factory.create<LStableDirk4>("LStableDirk4", name() + "_bootstrap", params);
80  }
81 }
82 
83 void
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  if (!_sys.solutionUDot())
90  mooseError("AStableDirk4: Time derivative of solution (`u_dot`) is not stored. Please set "
91  "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
92 
94  u_dot = *_solution;
96  u_dot.close();
98 }
99 
100 void
102  const dof_id_type & dof,
103  ADReal & /*ad_u_dotdot*/) const
104 {
106 }
107 
108 void
110 {
111  // Reset iteration counts
114 
115  if (_t_step == 1 && _safe_start)
116  {
117  _bootstrap_method->solve();
118  _n_nonlinear_iterations = _bootstrap_method->getNumNonlinearIterations();
119  _n_linear_iterations = _bootstrap_method->getNumLinearIterations();
120  }
121  else
122  {
123  // Time at end of step
124  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  for (unsigned int current_stage = 1; current_stage < 5; ++current_stage)
130  {
131  // Set the current stage value
132  _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.
138 
139  if (current_stage < 4)
140  {
141  _console << "Stage " << _stage << std::endl;
142  _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
143  }
144  else
145  {
146  _console << "Update Stage." << std::endl;
147  _fe_problem.time() = time_old + _dt;
148  }
149 
150  // Do the solve
151  _nl->system().solve();
152 
153  // Update the iteration counts
156 
157  // Abort time step immediately on stage failure - see TimeIntegrator doc page
158  if (!_fe_problem.converged(_nl->number()))
159  return;
160  }
161  }
162 }
163 
164 void
166 {
167  if (_t_step == 1 && _safe_start)
168  _bootstrap_method->postResidual(residual);
169 
170  else
171  {
172  // Error if _stage got messed up somehow.
173  if (_stage > 4)
174  mooseError("AStableDirk4::postResidual(): Member variable _stage can only have values 1-4.");
175 
176  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().
192 
193  // Build up the residual for this stage.
194  residual.add(1., *_Re_time);
195  for (unsigned int j = 0; j < _stage; ++j)
196  residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
197  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  residual.add(1., *_Re_time);
210  for (unsigned int j = 0; j < 3; ++j)
211  residual.add(_b[j], *_stage_residuals[j]);
212  residual.close();
213  }
214  }
215 }
Fourth-order diagonally implicit Runge Kutta method (Dirk) with five stages.
Definition: LStableDirk4.h:51
virtual void initPetscOutputAndSomeSolverSettings()
Reinitialize PETSc output for proper linear/nonlinear iteration display.
std::shared_ptr< LStableDirk4 > _bootstrap_method
Definition: AStableDirk4.h:104
virtual Real & time() const
virtual void solve() override
Solves the time step and sets the number of nonlinear and linear iterations.
Definition: AStableDirk4.C:109
FEProblemBase & _fe_problem
Reference to the problem.
void mooseInfo(Args &&... args) const
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
NonlinearSystemBase * _nl
Pointer to the nonlinear system, can happen that we dont have any.
SystemBase & _sys
Reference to the system this time integrator operates on.
Real _a[3][3]
Definition: AStableDirk4.h:92
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
AStableDirk4(const InputParameters &parameters)
Definition: AStableDirk4.C:31
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Base class for a system (of equations)
Definition: SystemBase.h:84
DualNumber< Real, DNDerivativeType, true > ADReal
Definition: ADRealForward.h:46
Real & _dt
The current time step size.
registerMooseObject("MooseApp", AStableDirk4)
virtual const std::string & name() const
Get the name of the class.
Definition: MooseBase.h:57
NumericVector< Number > * _Re_time
residual vector for time contributions
Factory & getFactory()
Retrieve a writable reference to the Factory associated with this App.
Definition: MooseApp.h:434
Real _b[3]
Definition: AStableDirk4.h:96
static InputParameters validParams()
Definition: AStableDirk4.C:22
const Real _gamma
Definition: AStableDirk4.h:81
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
virtual std::unique_ptr< Base > create()=0
unsigned int _stage
Definition: AStableDirk4.h:75
GHOSTED
virtual NumericVector< Number > * solutionUDot()
Definition: SystemBase.h:252
NumericVector< Number > * _stage_residuals[3]
Definition: AStableDirk4.h:78
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.
MooseApp & _app
The MOOSE application this is associated with.
Definition: MooseBase.h:84
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
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: AStableDirk4.C:101
Fourth-order diagonally implicit Runge Kutta method (Dirk) with three stages plus an update...
Definition: AStableDirk4.h:51
virtual void solve()
virtual void close()=0
const NumericVector< Number > *const & _solution
void computeTimeDerivativeHelper(T &u_dot, const T2 &u_old) const
Helper function that actually does the math for computing the time derivative.
Definition: AStableDirk4.h:109
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
Real _c[3]
Definition: AStableDirk4.h:85
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
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
const NumericVector< Number > & _solution_old
virtual void computeTimeDerivatives() override
Computes the time derivative and the Jacobian of the time derivative.
Definition: AStableDirk4.C:84
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.
NumericVector< Number > * addVector(const std::string &name, const bool project, const libMesh::ParallelType type)
Wrapper around vector addition for nonlinear time integrators.
virtual void add(const numeric_index_type i, const T value)=0
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
int & _t_step
The current time step number.
static InputParameters validParams()
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the NonLinearTimeIntegratorInterface called immediately after the residuals are computed ...
Definition: AStableDirk4.C:165
uint8_t dof_id_type
void computeDuDotDu()
Compute _du_dot_du.
const Real pi
virtual libMesh::System & system() override
Get the reference to the libMesh system.