Line data Source code
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 "LStableDirk4.h"
11 : #include "NonlinearSystemBase.h"
12 : #include "FEProblem.h"
13 : #include "PetscSupport.h"
14 :
15 : using namespace libMesh;
16 :
17 : registerMooseObject("MooseApp", LStableDirk4);
18 :
19 : InputParameters
20 14616 : LStableDirk4::validParams()
21 : {
22 14616 : InputParameters params = TimeIntegrator::validParams();
23 14616 : params.addClassDescription(
24 : "Fourth-order diagonally implicit Runge Kutta method (Dirk) with five stages.");
25 14616 : return params;
26 0 : }
27 :
28 : // Initialize static data
29 : const Real LStableDirk4::_c[LStableDirk4::_n_stages] = {.25, 0., .5, 1., 1.};
30 :
31 : const Real LStableDirk4::_a[LStableDirk4::_n_stages][LStableDirk4::_n_stages] = {
32 : {.25, 0, 0, 0, 0},
33 : {-.25, .25, 0, 0, 0},
34 : {.125, .125, .25, 0, 0},
35 : {-1.5, .75, 1.5, .25, 0},
36 : {0, 1. / 6, 2. / 3, -1. / 12, .25}};
37 :
38 204 : LStableDirk4::LStableDirk4(const InputParameters & parameters)
39 204 : : TimeIntegrator(parameters), _stage(1)
40 : {
41 204 : mooseInfo("LStableDirk4 and other multistage TimeIntegrators are known not to work with "
42 : "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
43 :
44 : // Name the stage residuals "residual_stage1", "residual_stage2", etc.
45 1224 : for (unsigned int stage = 0; stage < _n_stages; ++stage)
46 : {
47 1020 : std::ostringstream oss;
48 1020 : oss << "residual_stage" << stage + 1;
49 1020 : _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
50 1020 : }
51 204 : }
52 :
53 : void
54 18558 : LStableDirk4::computeTimeDerivatives()
55 : {
56 : // We are multiplying by the method coefficients in postResidual(), so
57 : // the time derivatives are of the same form at every stage although
58 : // the current solution varies depending on the stage.
59 18558 : if (!_sys.solutionUDot())
60 0 : mooseError("LStableDirk4: Time derivative of solution (`u_dot`) is not stored. Please set "
61 : "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
62 :
63 18558 : NumericVector<Number> & u_dot = *_sys.solutionUDot();
64 18558 : u_dot = *_solution;
65 18558 : computeTimeDerivativeHelper(u_dot, _solution_old);
66 18558 : u_dot.close();
67 18558 : computeDuDotDu();
68 18558 : }
69 :
70 : void
71 0 : LStableDirk4::computeADTimeDerivatives(ADReal & ad_u_dot,
72 : const dof_id_type & dof,
73 : ADReal & /*ad_u_dotdot*/) const
74 : {
75 0 : computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
76 0 : }
77 :
78 : void
79 222 : LStableDirk4::solve()
80 : {
81 : // Time at end of step
82 222 : Real time_old = _fe_problem.timeOld();
83 :
84 : // Reset iteration counts
85 222 : _n_nonlinear_iterations = 0;
86 222 : _n_linear_iterations = 0;
87 :
88 : // A for-loop would increment _stage too far, so we use an extra
89 : // loop counter.
90 1332 : for (unsigned int current_stage = 1; current_stage <= _n_stages; ++current_stage)
91 : {
92 : // Set the current stage value
93 1110 : _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.
98 1110 : _fe_problem.initPetscOutputAndSomeSolverSettings();
99 :
100 1110 : _console << "Stage " << _stage << std::endl;
101 :
102 : // Set the time for this stage
103 1110 : _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
104 :
105 : // If we previously used coloring, destroy the old object so it doesn't leak when we allocate a
106 : // new object in the following lines
107 1110 : _nl->destroyColoring();
108 :
109 : // Potentially setup finite differencing contexts for the solve
110 1110 : _nl->potentiallySetupFiniteDifferencing();
111 :
112 : // Do the solve
113 1110 : _nl->system().solve();
114 :
115 : // Update the iteration counts
116 1110 : _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
117 1110 : _n_linear_iterations += getNumLinearIterationsLastSolve();
118 :
119 : // Abort time step immediately on stage failure - see TimeIntegrator doc page
120 1110 : if (!_fe_problem.converged(_nl->number()))
121 0 : return;
122 : }
123 : }
124 :
125 : void
126 21898 : LStableDirk4::postResidual(NumericVector<Number> & residual)
127 : {
128 : // Error if _stage got messed up somehow.
129 21898 : if (_stage > _n_stages)
130 : // the explicit cast prevents strange compiler weirdness with the static
131 : // const variable and the variadic mooseError function
132 0 : mooseError("LStableDirk4::postResidual(): Member variable _stage can only have values 1-",
133 0 : (unsigned int)_n_stages,
134 : ".");
135 :
136 : // In the standard RK notation, the residual of stage 1 of s is given by:
137 : //
138 : // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
139 : //
140 : // where:
141 : // .) M is the mass matrix
142 : // .) Y_i is the stage solution
143 : // .) dt is the timestep, and is accounted for in the _Re_time residual.
144 : // .) f are the "non-time" residuals evaluated for a given stage solution.
145 : // .) The minus signs are already "baked in" to the residuals and so do not appear below.
146 :
147 : // Store this stage's non-time residual. We are calling operator=
148 : // here, and that calls close().
149 21898 : *_stage_residuals[_stage - 1] = *_Re_non_time;
150 :
151 : // Build up the residual for this stage.
152 21898 : residual.add(1., *_Re_time);
153 86821 : for (unsigned int j = 0; j < _stage; ++j)
154 64923 : residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
155 21898 : residual.close();
156 21898 : }
|
