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 : // MOOSE includes
11 : #include "ExplicitSSPRungeKutta.h"
12 : #include "NonlinearSystem.h"
13 : #include "FEProblem.h"
14 :
15 : // libMesh includes
16 : #include "libmesh/nonlinear_solver.h"
17 :
18 : using namespace libMesh;
19 :
20 : registerMooseObject("MooseApp", ExplicitSSPRungeKutta);
21 :
22 : InputParameters
23 14385 : ExplicitSSPRungeKutta::validParams()
24 : {
25 14385 : InputParameters params = ExplicitTimeIntegrator::validParams();
26 :
27 14385 : MooseEnum orders("1=1 2 3");
28 14385 : params.addRequiredParam<MooseEnum>("order", orders, "Order of time integration");
29 14385 : params.addClassDescription("Explicit strong stability preserving Runge-Kutta methods");
30 28770 : return params;
31 14385 : }
32 :
33 80 : ExplicitSSPRungeKutta::ExplicitSSPRungeKutta(const InputParameters & parameters)
34 : : ExplicitTimeIntegrator(parameters),
35 80 : _order(getParam<MooseEnum>("order")),
36 80 : _stage(0),
37 80 : _solution_intermediate_stage(addVector("solution_intermediate_stage", false, GHOSTED)),
38 80 : _tmp_solution(addVector("tmp_solution", false, GHOSTED)),
39 240 : _tmp_mass_solution_product(addVector("tmp_mass_solution_product", false, GHOSTED))
40 : {
41 : // For SSPRK methods up to order 3, the number of stages equals the order
42 80 : _n_stages = _order;
43 :
44 80 : if (_order == 1)
45 : {
46 64 : _a = {{1.0}};
47 0 : _b = {1.0};
48 32 : _c = {0.0};
49 : }
50 48 : else if (_order == 2)
51 : {
52 72 : _a = {{1.0, 0.0}, {0.5, 0.5}};
53 0 : _b = {1.0, 0.5};
54 24 : _c = {0.0, 1.0};
55 : }
56 24 : else if (_order == 3)
57 : {
58 96 : _a = {{1.0, 0.0, 0.0}, {0.75, 0.25, 0.0}, {1.0 / 3.0, 0.0, 2.0 / 3.0}};
59 0 : _b = {1.0, 0.25, 2.0 / 3.0};
60 24 : _c = {0.0, 1.0, 0.5};
61 : }
62 : else
63 0 : mooseError("Invalid time integration order.");
64 :
65 80 : _solution_stage.resize(_n_stages + 1, nullptr);
66 240 : }
67 :
68 : void
69 858 : ExplicitSSPRungeKutta::computeTimeDerivatives()
70 : {
71 : // Only the Jacobian needs to be computed, since the mass matrix needs it
72 858 : computeDuDotDu();
73 858 : }
74 :
75 : void
76 0 : ExplicitSSPRungeKutta::computeADTimeDerivatives(ADReal & ad_u_dot,
77 : const dof_id_type & dof,
78 : ADReal & /*ad_u_dotdot*/) const
79 : {
80 : // Note that if the solution for the current stage is not a nullptr, then neither
81 : // are the previous stages.
82 0 : if (_solution_stage[_stage])
83 : {
84 0 : for (unsigned int k = 0; k <= _stage; k++)
85 0 : ad_u_dot -= _a[_stage][k] * (*(_solution_stage[k]))(dof);
86 0 : ad_u_dot *= 1.0 / (_b[_stage] * _dt);
87 : }
88 : else
89 : {
90 : // We must be outside the solve loop in order to meet this criterion. In that case are we at
91 : // timestep_begin or timestep_end? We don't know, so I don't think it's meaningful to compute
92 : // derivatives here. Let's put in a quiet NaN which will only signal if we try to do something
93 : // meaningful with it (and then we do want to signal because time derivatives may not be
94 : // meaningful right now)
95 0 : ad_u_dot = std::numeric_limits<typename ADReal::value_type>::quiet_NaN();
96 : }
97 0 : }
98 :
99 : void
100 165 : ExplicitSSPRungeKutta::solve()
101 : {
102 : // Reset iteration counts
103 165 : _n_nonlinear_iterations = 0;
104 165 : _n_linear_iterations = 0;
105 :
106 165 : _current_time = _fe_problem.time();
107 165 : const Real time_old = _fe_problem.timeOld();
108 165 : const Real dt = _current_time - time_old;
109 :
110 165 : bool converged = false;
111 :
112 165 : _solution_stage[0] = &_solution_old;
113 495 : for (_stage = 0; _stage < _n_stages; _stage++)
114 : {
115 330 : if (_stage == 0)
116 : {
117 : // Nothing needs to be done
118 : }
119 165 : else if (_stage == _n_stages - 1)
120 : {
121 110 : _solution_stage[_stage] = _solution;
122 : }
123 : else
124 : {
125 : // Else must be the intermediate stage of the 3-stage method
126 55 : *_solution_intermediate_stage = *_solution;
127 55 : _solution_intermediate_stage->close();
128 55 : _solution_stage[_stage] = _solution_intermediate_stage;
129 : }
130 :
131 : // Set stage time for residual evaluation
132 330 : _fe_problem.time() = time_old + _c[_stage] * dt;
133 330 : _nonlinear_implicit_system->update();
134 :
135 330 : converged = solveStage();
136 330 : if (!converged)
137 0 : return;
138 : }
139 :
140 165 : if (_stage == _n_stages)
141 : // We made it to the end of the solve. We may call functions like computeTimeDerivatives later
142 : // for postprocessing purposes in which case we need to ensure we're accessing our data
143 : // correctly (e.g. not out-of-bounds)
144 165 : --_stage;
145 : }
146 :
147 : bool
148 330 : ExplicitSSPRungeKutta::solveStage()
149 : {
150 : // Compute the mass matrix
151 330 : computeTimeDerivatives();
152 330 : auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
153 330 : _fe_problem.computeJacobianTag(
154 330 : *_nonlinear_implicit_system->current_local_solution, mass_matrix, _Ke_time_tag);
155 :
156 : // Compute RHS vector using previous stage solution in steady-state residual
157 330 : _explicit_residual->zero();
158 330 : _fe_problem.computeResidual(
159 330 : *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
160 :
161 : // Move the residual to the RHS
162 330 : *_explicit_residual *= -1.0;
163 :
164 : // Perform the linear solve
165 330 : bool converged = performExplicitSolve(mass_matrix);
166 :
167 : // Update the solution: u^(s) = u^(s-1) + du^(s)
168 330 : (*_nonlinear_implicit_system->solution).zero();
169 330 : *_nonlinear_implicit_system->solution = *(_solution_stage[_stage]);
170 330 : *_nonlinear_implicit_system->solution += *_solution_update;
171 :
172 : // Enforce contraints on the solution
173 330 : DofMap & dof_map = _nonlinear_implicit_system->get_dof_map();
174 330 : dof_map.enforce_constraints_exactly(*_nonlinear_implicit_system,
175 330 : _nonlinear_implicit_system->solution.get());
176 330 : _nonlinear_implicit_system->update();
177 :
178 330 : _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
179 :
180 330 : _nonlinear_implicit_system->nonlinear_solver->converged = converged;
181 :
182 330 : return converged;
183 : }
184 :
185 : void
186 330 : ExplicitSSPRungeKutta::postResidual(NumericVector<Number> & residual)
187 : {
188 : // The time residual is not included in the steady-state residual
189 330 : residual += *_Re_non_time;
190 :
191 : // Compute \sum_{k=0}^{s-1} a_{s,k} u^(k) - u^(s-1)
192 330 : _tmp_solution->zero();
193 880 : for (unsigned int k = 0; k <= _stage; k++)
194 550 : _tmp_solution->add(_a[_stage][k], *(_solution_stage[k]));
195 330 : _tmp_solution->add(-1.0, *(_solution_stage[_stage]));
196 330 : _tmp_solution->close();
197 :
198 : // Perform mass matrix product with the above vector
199 330 : auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
200 330 : mass_matrix.vector_mult(*_tmp_mass_solution_product, *_tmp_solution);
201 :
202 : // Finish computing residual vector (before modification by nodal BCs)
203 330 : residual -= *_tmp_mass_solution_product;
204 :
205 330 : residual.close();
206 :
207 : // Set time at which to evaluate nodal BCs
208 330 : _fe_problem.time() = _current_time;
209 330 : }
210 :
211 : Real
212 858 : ExplicitSSPRungeKutta::duDotDuCoeff() const
213 : {
214 858 : return Real(1) / _b[_stage];
215 : }
|
