www.mooseframework.org
AStableDirk4.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
18 
19 template <>
22 {
24  params.addParam<bool>("safe_start", true, "If true, use LStableDirk4 to bootstrap this method.");
25  return params;
26 }
27 
29  : TimeIntegrator(parameters),
30  _stage(1),
31  _gamma(0.5 + std::sqrt(3) / 3. * std::cos(libMesh::pi / 18.)),
32  _safe_start(getParam<bool>("safe_start"))
33 {
34  mooseInfo("AStableDirk4 and other multistage TimeIntegrators are known not to work with "
35  "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
36 
37  // Name the stage residuals "residual_stage1", "residual_stage2", etc.
38  for (unsigned int stage = 0; stage < 3; ++stage)
39  {
40  std::ostringstream oss;
41  oss << "residual_stage" << stage + 1;
42  _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
43  }
44 
45  // Initialize parameters
46  _c[0] = _gamma;
47  _c[1] = .5;
48  _c[2] = 1.0 - _gamma;
49 
50  _a[0][0] = _gamma;
51  _a[1][0] = .5 - _gamma;
52  _a[1][1] = _gamma;
53  _a[2][0] = 2. * _gamma;
54  _a[2][1] = 1 - 4. * _gamma;
55  _a[2][2] = _gamma;
56 
57  _b[0] = 1. / (24. * (.5 - _gamma) * (.5 - _gamma));
58  _b[1] = 1. - 1. / (12. * (.5 - _gamma) * (.5 - _gamma));
59  _b[2] = _b[0];
60 
61  // If doing a _safe_start, construct the bootstrapping
62  // TimeIntegrator. Note that this method will also add
63  // residual_stage vectors to the system, but since they have the
64  // same name as the vectors we added, they won't be duplicated.
65  if (_safe_start)
66  {
67  Factory & factory = _app.getFactory();
68  InputParameters params = factory.getValidParams("LStableDirk4");
69 
70  // We need to set some parameters that are normally set in
71  // FEProblemBase::addTimeIntegrator() to ensure that the
72  // getCheckedPointerParam() sanity checking is happy. This is why
73  // constructing MOOSE objects "manually" is generally frowned upon.
74  params.set<SystemBase *>("_sys") = &_sys;
75 
76  _bootstrap_method = factory.create<LStableDirk4>("LStableDirk4", name() + "_bootstrap", params);
77  }
78 }
79 
80 void
82 {
83  // We are multiplying by the method coefficients in postResidual(), so
84  // the time derivatives are of the same form at every stage although
85  // the current solution varies depending on the stage.
86  if (!_sys.solutionUDot())
87  mooseError("AStableDirk4: Time derivative of solution (`u_dot`) is not stored. Please set "
88  "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
89 
90  NumericVector<Number> & u_dot = *_sys.solutionUDot();
91  u_dot = *_solution;
93  u_dot.close();
94  _du_dot_du = 1. / _dt;
95 }
96 
97 void
98 AStableDirk4::computeADTimeDerivatives(DualReal & ad_u_dot, const dof_id_type & dof) const
99 {
101 }
102 
103 void
105 {
106  // Reset iteration counts
109 
110  if (_t_step == 1 && _safe_start)
111  {
112  _bootstrap_method->solve();
113  _n_nonlinear_iterations = _bootstrap_method->getNumNonlinearIterations();
114  _n_linear_iterations = _bootstrap_method->getNumLinearIterations();
115  }
116  else
117  {
118  // Time at end of step
119  Real time_old = _fe_problem.timeOld();
120 
121  // A for-loop would increment _stage too far, so we use an extra
122  // loop counter. The method has three stages and an update stage,
123  // which we treat as just one more (special) stage in the implementation.
124  for (unsigned int current_stage = 1; current_stage < 5; ++current_stage)
125  {
126  // Set the current stage value
127  _stage = current_stage;
128 
129  // This ensures that all the Output objects in the OutputWarehouse
130  // have had solveSetup() called, and sets the default solver
131  // parameters for PETSc.
133 
134  if (current_stage < 4)
135  {
136  _console << "Stage " << _stage << "\n";
137  _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
138  }
139  else
140  {
141  _console << "Update Stage.\n";
142  _fe_problem.time() = time_old + _dt;
143  }
144 
145  // Do the solve
147 
148  // Update the iteration counts
151 
152  // Abort time step immediately on stage failure - see TimeIntegrator doc page
153  if (!_fe_problem.converged())
154  return;
155  }
156  }
157 }
158 
159 void
160 AStableDirk4::postResidual(NumericVector<Number> & residual)
161 {
162  if (_t_step == 1 && _safe_start)
163  _bootstrap_method->postResidual(residual);
164 
165  else
166  {
167  // Error if _stage got messed up somehow.
168  if (_stage > 4)
169  mooseError("AStableDirk4::postResidual(): Member variable _stage can only have values 1-4.");
170 
171  if (_stage < 4)
172  {
173  // In the standard RK notation, the residual of stage 1 of s is given by:
174  //
175  // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
176  //
177  // where:
178  // .) M is the mass matrix
179  // .) Y_i is the stage solution
180  // .) dt is the timestep, and is accounted for in the _Re_time residual.
181  // .) f are the "non-time" residuals evaluated for a given stage solution.
182  // .) The minus signs are already "baked in" to the residuals and so do not appear below.
183 
184  // Store this stage's non-time residual. We are calling operator=
185  // here, and that calls close().
187 
188  // Build up the residual for this stage.
189  residual.add(1., _Re_time);
190  for (unsigned int j = 0; j < _stage; ++j)
191  residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
192  residual.close();
193  }
194  else
195  {
196  // The update step is a final solve:
197  //
198  // R := M*(y_{n+1} - y_n)/dt - \sum_{j=1}^s b_j * f(t_n + c_j*dt, Y_j) = 0
199  //
200  // We could potentially fold _b up into an extra row of _a and
201  // just do one more stage, but I think handling it separately like
202  // this is easier to understand and doesn't create too much code
203  // repitition.
204  residual.add(1., _Re_time);
205  for (unsigned int j = 0; j < 3; ++j)
206  residual.add(_b[j], *_stage_residuals[j]);
207  residual.close();
208  }
209  }
210 }
NonlinearSystemBase & _nl
Fourth-order diagonally implicit Runge Kutta method (Dirk) with five stages.
Definition: LStableDirk4.h:56
virtual NumericVector< Number > * solutionUDot()=0
std::shared_ptr< MooseObject > create(const std::string &obj_name, const std::string &name, InputParameters parameters, THREAD_ID tid=0, bool print_deprecated=true)
Build an object (must be registered) - THIS METHOD IS DEPRECATED (Use create<T>()) ...
Definition: Factory.C:87
void initPetscOutput()
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:104
NonlinearSystemBase & getNonlinearSystemBase()
Generic factory class for build all sorts of objects.
Definition: Factory.h:153
virtual NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const ParallelType type)
Adds a solution length vector to the system.
Definition: SystemBase.C:542
InputParameters getValidParams(const std::string &name)
Get valid parameters for the object.
Definition: Factory.C:67
FEProblemBase & _fe_problem
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
NumericVector< Number > & _Re_non_time
residual vector for non-time contributions
SystemBase & _sys
Real _a[3][3]
Definition: AStableDirk4.h:92
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
DualNumber< Real, NumberArray< AD_MAX_DOFS_PER_ELEM, Real > > DualReal
Definition: DualReal.h:29
AStableDirk4(const InputParameters &parameters)
Definition: AStableDirk4.C:28
void computeADTimeDerivatives(DualReal &ad_u_dot, const dof_id_type &dof) const override
method for computing local automatic differentiation time derivatives
Definition: AStableDirk4.C:98
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Base class for a system (of equations)
Definition: SystemBase.h:92
registerMooseObject("MooseApp", AStableDirk4)
Factory & getFactory()
Retrieve the Factory associated with this App.
Definition: MooseApp.h:286
Real _b[3]
Definition: AStableDirk4.h:96
virtual bool converged() override
const Real _gamma
Definition: AStableDirk4.h:81
InputParameters validParams< AStableDirk4 >()
Definition: AStableDirk4.C:21
unsigned int _stage
Definition: AStableDirk4.h:75
NumericVector< Number > * _stage_residuals[3]
Definition: AStableDirk4.h:78
InputParameters validParams< TimeIntegrator >()
unsigned int _n_linear_iterations
Total number of linear iterations over all stages of the time step.
Fourth-order diagonally implicit Runge Kutta method (Dirk) with three stages plus an update...
Definition: AStableDirk4.h:56
Real & _du_dot_du
solution vector for
const NumericVector< Number > *const & _solution
solution vectors
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
Base class for time integrators.
unsigned int getNumLinearIterationsLastSolve() const
Gets the number of linear iterations in the most recent solve.
virtual System & system() override
Get the reference to the libMesh system.
MooseApp & _app
The MooseApp this object is associated with.
Definition: MooseObject.h:177
Real _c[3]
Definition: AStableDirk4.h:85
const std::string & name() const
Get the name of the object.
Definition: MooseObject.h:59
NumericVector< Number > & _Re_time
residual vector for time contributions
virtual Real & timeOld() const
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
const NumericVector< Number > & _solution_old
virtual void computeTimeDerivatives() override
Definition: AStableDirk4.C:81
void mooseInfo(Args &&... args) const
Definition: MooseObject.h:167
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.
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the TimeIntegrator called immediately after the residuals are computed in NonlinearSystem...
Definition: AStableDirk4.C:160