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