www.mooseframework.org
LStableDirk3.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 #include "LStableDirk3.h"
11 #include "NonlinearSystem.h"
12 #include "FEProblem.h"
13 #include "PetscSupport.h"
14 
16 
19 {
21  params.addClassDescription(
22  "Third order diagonally implicit Runge Kutta method (Dirk) with three stages.");
23  return params;
24 }
25 
27  : TimeIntegrator(parameters),
28  _stage(1),
29  _gamma(-std::sqrt(2.) * std::cos(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. +
30  std::sqrt(6.) * std::sin(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. + 1.)
31 {
32  mooseInfo("LStableDirk3 and other multistage TimeIntegrators are known not to work with "
33  "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
34 
35  // Name the stage residuals "residual_stage1", "residual_stage2", etc.
36  for (unsigned int stage = 0; stage < 3; ++stage)
37  {
38  std::ostringstream oss;
39  oss << "residual_stage" << stage + 1;
40  _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
41  }
42 
43  // Initialize parameters
44  _c[0] = _gamma;
45  _c[1] = .5 * (1 + _gamma);
46  _c[2] = 1.0;
47 
48  _a[0][0] = _gamma;
49  _a[1][0] = .5 * (1 - _gamma);
50  _a[1][1] = _gamma;
51  _a[2][0] = .25 * (-6 * _gamma * _gamma + 16 * _gamma - 1);
52  _a[2][1] = .25 * (6 * _gamma * _gamma - 20 * _gamma + 5);
53  _a[2][2] = _gamma;
54 }
55 
56 void
58 {
59  // We are multiplying by the method coefficients in postResidual(), so
60  // the time derivatives are of the same form at every stage although
61  // the current solution varies depending on the stage.
62  if (!_sys.solutionUDot())
63  mooseError("LStableDirk3: Time derivative of solution (`u_dot`) is not stored. Please set "
64  "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
65 
67  u_dot = *_solution;
69  u_dot.close();
70  _du_dot_du = 1. / _dt;
71 }
72 
73 void
75  const dof_id_type & dof,
76  DualReal & /*ad_u_dotdot*/) const
77 {
79 }
80 
81 void
83 {
84  // Time at end of step
85  Real time_old = _fe_problem.timeOld();
86 
87  // Reset iteration counts
90 
91  // A for-loop would increment _stage too far, so we use an extra
92  // loop counter.
93  for (unsigned int current_stage = 1; current_stage < 4; ++current_stage)
94  {
95  // Set the current stage value
96  _stage = current_stage;
97 
98  // This ensures that all the Output objects in the OutputWarehouse
99  // have had solveSetup() called, and sets the default solver
100  // parameters for PETSc.
102 
103  _console << "Stage " << _stage << std::endl;
104 
105  // Set the time for this stage
106  _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
107 
108  // Do the solve
109  _nl.system().solve();
110 
111  // Update the iteration counts
114 
115  // Abort time step immediately on stage failure - see TimeIntegrator doc page
117  return;
118  }
119 }
120 
121 void
123 {
124  // Error if _stage got messed up somehow.
125  if (_stage > 3)
126  mooseError(
127  "LStableDirk3::postResidual(): Member variable _stage can only have values 1, 2, or 3.");
128 
129  // In the standard RK notation, the residual of stage 1 of s is given by:
130  //
131  // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
132  //
133  // where:
134  // .) M is the mass matrix
135  // .) Y_i is the stage solution
136  // .) dt is the timestep, and is accounted for in the _Re_time residual.
137  // .) f are the "non-time" residuals evaluated for a given stage solution.
138  // .) The minus signs are already "baked in" to the residuals and so do not appear below.
139 
140  // Store this stage's non-time residual. We are calling operator=
141  // here, and that calls close().
143 
144  // Build up the residual for this stage.
145  residual.add(1., _Re_time);
146  for (unsigned int j = 0; j < _stage; ++j)
147  residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
148  residual.close();
149 }
NonlinearSystemBase & _nl
virtual void initPetscOutputAndSomeSolverSettings()
Reinitialize PETSc output for proper linear/nonlinear iteration display.
virtual NumericVector< Number > * solutionUDot()=0
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 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
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
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the TimeIntegrator called immediately after the residuals are computed in NonlinearSystem...
Definition: LStableDirk3.C:122
FEProblemBase & _fe_problem
void mooseInfo(Args &&... args) const
NumericVector< Number > * _stage_residuals[3]
Definition: LStableDirk3.h:66
virtual void solve() override
Solves the time step and sets the number of nonlinear and linear iterations.
Definition: LStableDirk3.C:82
DualNumber< Real, DNDerivativeType, true > DualReal
Definition: DualReal.h:49
NumericVector< Number > & _Re_non_time
residual vector for non-time contributions
SystemBase & _sys
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:18
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:85
MetaPhysicL::DualNumber< V, D, asd > sqrt(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: ADReal.h:35
Real _a[3][3]
Definition: LStableDirk3.h:80
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
unsigned int _n_linear_iterations
Total number of linear iterations over all stages of the time step.
unsigned int _stage
Definition: LStableDirk3.h:63
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1125
virtual void computeTimeDerivatives() override
Computes the time derivative and the Jacobian of the time derivative.
Definition: LStableDirk3.C:57
virtual void close()=0
virtual 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: LStableDirk3.C:74
Real & _du_dot_du
Derivative of time derivative with respect to current solution: .
const NumericVector< Number > *const & _solution
solution vectors
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:69
virtual System & system() override
Get the reference to the libMesh system.
LStableDirk3(const InputParameters &parameters)
Definition: LStableDirk3.C:26
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
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:73
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()
uint8_t dof_id_type