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