https://mooseframework.inl.gov
LStableDirk3.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 "LStableDirk3.h"
11 #include "NonlinearSystem.h"
12 #include "FEProblem.h"
13 #include "PetscSupport.h"
14 using namespace libMesh;
15 
17 
20 {
22  params.addClassDescription(
23  "Third order diagonally implicit Runge Kutta method (Dirk) with three stages.");
24  return params;
25 }
26 
28  : TimeIntegrator(parameters),
29  _stage(1),
30  _gamma(-std::sqrt(2.) * std::cos(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. +
31  std::sqrt(6.) * std::sin(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. + 1.)
32 {
33  mooseInfo("LStableDirk3 and other multistage TimeIntegrators are known not to work with "
34  "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
35 
36  // Name the stage residuals "residual_stage1", "residual_stage2", etc.
37  for (unsigned int stage = 0; stage < 3; ++stage)
38  {
39  std::ostringstream oss;
40  oss << "residual_stage" << stage + 1;
41  _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
42  }
43 
44  // Initialize parameters
45  _c[0] = _gamma;
46  _c[1] = .5 * (1 + _gamma);
47  _c[2] = 1.0;
48 
49  _a[0][0] = _gamma;
50  _a[1][0] = .5 * (1 - _gamma);
51  _a[1][1] = _gamma;
52  _a[2][0] = .25 * (-6 * _gamma * _gamma + 16 * _gamma - 1);
53  _a[2][1] = .25 * (6 * _gamma * _gamma - 20 * _gamma + 5);
54  _a[2][2] = _gamma;
55 }
56 
57 void
59 {
60  // We are multiplying by the method coefficients in postResidual(), so
61  // the time derivatives are of the same form at every stage although
62  // the current solution varies depending on the stage.
63  if (!_sys.solutionUDot())
64  mooseError("LStableDirk3: Time derivative of solution (`u_dot`) is not stored. Please set "
65  "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
66 
68  u_dot = *_solution;
70  u_dot.close();
72 }
73 
74 void
76  const dof_id_type & dof,
77  ADReal & /*ad_u_dotdot*/) const
78 {
80 }
81 
82 void
84 {
85  // Time at end of step
86  Real time_old = _fe_problem.timeOld();
87 
88  // Reset iteration counts
91 
92  // A for-loop would increment _stage too far, so we use an extra
93  // loop counter.
94  for (unsigned int current_stage = 1; current_stage < 4; ++current_stage)
95  {
96  // Set the current stage value
97  _stage = current_stage;
98 
99  // This ensures that all the Output objects in the OutputWarehouse
100  // have had solveSetup() called, and sets the default solver
101  // parameters for PETSc.
103 
104  _console << "Stage " << _stage << std::endl;
105 
106  // Set the time for this stage
107  _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
108 
109  // If we previously used coloring, destroy the old object so it doesn't leak when we allocate a
110  // new object in the following lines
111  _nl->destroyColoring();
112 
113  // Potentially setup finite differencing contexts for the solve
115 
116  // Do the solve
117  _nl->system().solve();
118 
119  // Update the iteration counts
122 
123  // Abort time step immediately on stage failure - see TimeIntegrator doc page
124  if (!_fe_problem.converged(_nl->number()))
125  return;
126  }
127 }
128 
129 void
131 {
132  // Error if _stage got messed up somehow.
133  if (_stage > 3)
134  mooseError(
135  "LStableDirk3::postResidual(): Member variable _stage can only have values 1, 2, or 3.");
136 
137  // In the standard RK notation, the residual of stage 1 of s is given by:
138  //
139  // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
140  //
141  // where:
142  // .) M is the mass matrix
143  // .) Y_i is the stage solution
144  // .) dt is the timestep, and is accounted for in the _Re_time residual.
145  // .) f are the "non-time" residuals evaluated for a given stage solution.
146  // .) The minus signs are already "baked in" to the residuals and so do not appear below.
147 
148  // Store this stage's non-time residual. We are calling operator=
149  // here, and that calls close().
151 
152  // Build up the residual for this stage.
153  residual.add(1., *_Re_time);
154  for (unsigned int j = 0; j < _stage; ++j)
155  residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
156  residual.close();
157 }
virtual void initPetscOutputAndSomeSolverSettings()
Reinitialize PETSc output for proper linear/nonlinear iteration display.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
virtual Real & time() const
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the NonLinearTimeIntegratorInterface called immediately after the residuals are computed ...
Definition: LStableDirk3.C:130
FEProblemBase & _fe_problem
Reference to the problem.
void mooseInfo(Args &&... args) const
NumericVector< Number > * _stage_residuals[3]
Definition: LStableDirk3.h:67
virtual void solve() override
Solves the time step and sets the number of nonlinear and linear iterations.
Definition: LStableDirk3.C:83
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.
Third order diagonally implicit Runge Kutta method (Dirk) with three stages.
Definition: LStableDirk3.h:40
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
static InputParameters validParams()
Definition: LStableDirk3.C:19
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
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 computeTimeDerivativeHelper(T &u_dot, const T2 &u_old) const
Helper function that actually does the math for computing the time derivative.
Definition: LStableDirk3.h:86
Real _a[3][3]
Definition: LStableDirk3.h:81
registerMooseObject("MooseApp", LStableDirk3)
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
GHOSTED
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 _stage
Definition: LStableDirk3.h:64
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
virtual void computeTimeDerivatives() override
Computes the time derivative and the Jacobian of the time derivative.
Definition: LStableDirk3.C:58
void destroyColoring()
Destroy the coloring object if it exists.
virtual void solve()
virtual void close()=0
virtual 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: LStableDirk3.C:75
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.
const Real _gamma
Definition: LStableDirk3.h:70
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
LStableDirk3(const InputParameters &parameters)
Definition: LStableDirk3.C:27
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.
Real _c[3]
Definition: LStableDirk3.h:74
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
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.