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 "Assembly.h"
12 : #include "ExplicitMixedOrder.h"
13 : #include "ExplicitTimeIntegrator.h"
14 : #include "Moose.h"
15 : #include "MooseError.h"
16 : #include "MooseTypes.h"
17 : #include "MooseVariableFieldBase.h"
18 : #include "NonlinearSystem.h"
19 : #include "FEProblem.h"
20 : #include "TimeStepper.h"
21 : #include "TransientBase.h"
22 :
23 : // libMesh includes
24 : #include "TransientBase.h"
25 : #include "libmesh/id_types.h"
26 : #include "libmesh/nonlinear_solver.h"
27 : #include "libmesh/sparse_matrix.h"
28 : #include "DirichletBCBase.h"
29 : #include "libmesh/vector_value.h"
30 : #include <algorithm>
31 : #include <iterator>
32 : #include <utility>
33 :
34 : registerMooseObject("SolidMechanicsApp", ExplicitMixedOrder);
35 : registerMooseObjectRenamed("SolidMechanicsApp",
36 : DirectCentralDifference,
37 : "10/14/2025 00:00",
38 : ExplicitMixedOrder);
39 :
40 : InputParameters
41 87 : ExplicitMixedOrder::validParams()
42 : {
43 87 : InputParameters params = ExplicitTimeIntegrator::validParams();
44 :
45 87 : params.addClassDescription(
46 : "Implementation of explicit time integration without invoking any of the nonlinear solver.");
47 :
48 174 : params.addParam<bool>("use_constant_mass",
49 174 : false,
50 : "If set to true, will only compute the mass matrix in the first time step, "
51 : "and keep using it throughout the simulation.");
52 :
53 174 : params.addParam<bool>(
54 : "recompute_mass_matrix_after_mesh_change",
55 174 : false,
56 : "If set to true, the mass matrix will be recomputed when the mesh changes (e.g. through "
57 : "adaptivity). If use_constant_mass is set to true, adadaptivity is used, and this parameter "
58 : "is not set to true, the simulation will error out when the mesh changes.");
59 :
60 174 : params.addParam<TagName>("mass_matrix_tag", "mass", "The tag for the mass matrix");
61 :
62 174 : params.addParam<std::vector<VariableName>>(
63 : "second_order_vars",
64 : {},
65 : "A subset of variables that require second-order integration (velocity and acceleration) to "
66 : "be applied by this time integrator.");
67 :
68 174 : params.addParam<std::vector<VariableName>>(
69 : "first_order_vars",
70 : {},
71 : "A subset of variables that require first-order integration (velocity only) to be applied by "
72 : "this time integrator.");
73 :
74 : // Prevent users from using variables option by accident.
75 87 : params.suppressParameter<std::vector<VariableName>>("variables");
76 :
77 174 : MooseEnum solve_type("consistent lumped lump_preconditioned", "lumped");
78 87 : params.setParameters("solve_type", solve_type);
79 87 : params.ignoreParameter<MooseEnum>("solve_type");
80 87 : return params;
81 87 : }
82 :
83 58 : ExplicitMixedOrder::ExplicitMixedOrder(const InputParameters & parameters)
84 : : ExplicitTimeIntegrator(parameters),
85 58 : _constant_mass(getParam<bool>("use_constant_mass")),
86 58 : _recompute_mass_matrix_on_mesh_change(
87 58 : getParam<bool>("recompute_mass_matrix_after_mesh_change")),
88 58 : _mesh_changed(true),
89 116 : _mass_matrix_name(getParam<TagName>("mass_matrix_tag")),
90 58 : _mass_matrix_lumped(addVector("mass_matrix_lumped", true, GHOSTED)),
91 58 : _solution_older(_sys.solutionState(2)),
92 116 : _vars_first(declareRestartableData<std::unordered_set<unsigned int>>("first_order_vars")),
93 58 : _local_first_order_indices(
94 58 : declareRestartableData<std::vector<dof_id_type>>("first_local_indices")),
95 116 : _vars_second(declareRestartableData<std::unordered_set<unsigned int>>("second_order_vars")),
96 58 : _local_second_order_indices(
97 116 : declareRestartableData<std::vector<dof_id_type>>("second_local_indices"))
98 : {
99 58 : _fe_problem.setUDotRequested(true);
100 58 : _fe_problem.setUDotOldRequested(true);
101 58 : _fe_problem.setUDotDotRequested(true);
102 58 : }
103 :
104 : void
105 2001 : ExplicitMixedOrder::computeTimeDerivatives()
106 : {
107 : /*
108 : Because this is called in NonLinearSystemBase
109 : this should not actually compute the time derivatives.
110 : Calculating time derivatives here will cause issues for the
111 : solution update.
112 : */
113 2001 : return;
114 : }
115 :
116 : void
117 12 : ExplicitMixedOrder::meshChanged()
118 : {
119 12 : if (_constant_mass && !_recompute_mass_matrix_on_mesh_change)
120 0 : paramError("recompute_mass_matrix_after_mesh_change",
121 : "Must be set to true explicitly by the user to support adaptivity with "
122 : "`use_constant_mass`.");
123 :
124 : // after mesh changes we need to recompute the mass matrix, it is not interpolable as it contains
125 : // a volume integrated quantity!
126 12 : _mesh_changed = true;
127 :
128 12 : ExplicitTimeIntegrator::meshChanged();
129 12 : }
130 :
131 : TagID
132 398 : ExplicitMixedOrder::massMatrixTagID() const
133 : {
134 398 : return _sys.subproblem().getMatrixTagID(_mass_matrix_name);
135 : }
136 :
137 : void
138 369 : ExplicitMixedOrder::solve()
139 : {
140 : // Getting the tagID for the mass matrix
141 369 : auto mass_tag = massMatrixTagID();
142 :
143 : // Reset iteration counts
144 369 : _n_nonlinear_iterations = 0;
145 369 : _n_linear_iterations = 0;
146 :
147 369 : _current_time = _fe_problem.time();
148 :
149 369 : auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
150 :
151 369 : if (_mesh_changed)
152 33 : updateDOFIndices();
153 :
154 : // Compute the mass matrix
155 369 : if (_mesh_changed || !_constant_mass)
156 : {
157 : // We only want to compute "inverted" lumped mass matrix once.
158 78 : _fe_problem.computeJacobianTag(
159 78 : *_nonlinear_implicit_system->current_local_solution, mass_matrix, mass_tag);
160 :
161 : // Calculate and record the lumped mass matrix for use in residual calculation
162 78 : mass_matrix.vector_mult(*_mass_matrix_lumped, *_ones);
163 78 : _mass_matrix_lumped->close();
164 :
165 : // "Invert" the diagonal mass matrix
166 78 : *_mass_matrix_diag_inverted = *_mass_matrix_lumped;
167 78 : _mass_matrix_diag_inverted->reciprocal();
168 78 : _mass_matrix_diag_inverted->close();
169 : }
170 :
171 369 : _mesh_changed = false;
172 :
173 : // Set time to the time at which to evaluate the residual
174 369 : _fe_problem.time() = _fe_problem.timeOld();
175 369 : _nonlinear_implicit_system->update();
176 :
177 : // Evaluate residual and move it to the RHS
178 369 : evaluateRHSResidual();
179 :
180 : // Perform the linear solve
181 369 : bool converged = performExplicitSolve(mass_matrix);
182 369 : _nl->overwriteNodeFace(*_nonlinear_implicit_system->solution);
183 :
184 : // Update the solution
185 369 : *_nonlinear_implicit_system->solution = _nl->solutionOld();
186 369 : *_nonlinear_implicit_system->solution += *_solution_update;
187 :
188 369 : _nonlinear_implicit_system->update();
189 :
190 369 : _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
191 369 : _nonlinear_implicit_system->nonlinear_solver->converged = converged;
192 369 : }
193 :
194 : void
195 369 : ExplicitMixedOrder::evaluateRHSResidual()
196 : {
197 : // Compute the residual
198 369 : _fe_problem.computeResidual(
199 369 : *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
200 :
201 : // Move the residual to the RHS
202 369 : *_explicit_residual *= -1.0;
203 369 : }
204 :
205 : void
206 369 : ExplicitMixedOrder::postResidual(NumericVector<Number> & residual)
207 : {
208 369 : residual += *_Re_time;
209 369 : residual += *_Re_non_time;
210 369 : residual.close();
211 :
212 : // Reset time to the time at which to evaluate nodal BCs, which comes next
213 369 : _fe_problem.time() = _current_time;
214 369 : }
215 :
216 : bool
217 369 : ExplicitMixedOrder::performExplicitSolve(SparseMatrix<Number> &)
218 : {
219 : bool converged = false;
220 :
221 : // Grab all the vectors that we will need
222 369 : auto accel = _sys.solutionUDotDot();
223 369 : auto vel = _sys.solutionUDot();
224 :
225 : // Compute Forward Euler
226 : // Split diag mass and residual vectors into correct subvectors
227 : const std::unique_ptr<NumericVector<Number>> mass_inv_first(
228 369 : NumericVector<Number>::build(_communicator, libMesh::default_solver_package(), PARALLEL));
229 : const std::unique_ptr<NumericVector<Real>> exp_res_first(
230 369 : NumericVector<Number>::build(_communicator, libMesh::default_solver_package(), PARALLEL));
231 369 : _mass_matrix_diag_inverted->create_subvector(*mass_inv_first, _local_first_order_indices, false);
232 369 : _explicit_residual->create_subvector(*exp_res_first, _local_first_order_indices, false);
233 :
234 : // Need velocity vector split into subvectors
235 369 : auto vel_first = vel->get_subvector(_local_first_order_indices);
236 :
237 : // Velocity update for foward euler
238 369 : vel_first->pointwise_mult(*mass_inv_first, *exp_res_first);
239 :
240 : // Restore the velocities
241 369 : vel->restore_subvector(std::move(vel_first), _local_first_order_indices);
242 :
243 : // Compute Central Difference
244 : // Split diag mass and residual vectors into correct subvectors
245 : const std::unique_ptr<NumericVector<Real>> mass_inv_second(
246 369 : NumericVector<Number>::build(_communicator, libMesh::default_solver_package(), PARALLEL));
247 : const std::unique_ptr<NumericVector<Real>> exp_res_second(
248 369 : NumericVector<Number>::build(_communicator, libMesh::default_solver_package(), PARALLEL));
249 369 : _mass_matrix_diag_inverted->create_subvector(
250 369 : *mass_inv_second, _local_second_order_indices, false);
251 369 : _explicit_residual->create_subvector(*exp_res_second, _local_second_order_indices, false);
252 :
253 : // Only need acceleration and old velocity vector for central difference
254 369 : auto accel_second = accel->get_subvector(_local_second_order_indices);
255 :
256 369 : auto vel_second = vel->get_subvector(_local_second_order_indices);
257 :
258 : // Compute acceleration for central difference
259 369 : accel_second->pointwise_mult(*mass_inv_second, *exp_res_second);
260 :
261 : // Scaling the acceleration
262 369 : auto accel_scaled = accel_second->clone();
263 369 : accel_scaled->scale((_dt + _dt_old) / 2);
264 :
265 : // Velocity update for central difference
266 369 : *vel_second += *accel_scaled;
267 :
268 : // Restore acceleration
269 369 : accel->restore_subvector(std::move(accel_second), _local_second_order_indices);
270 :
271 369 : vel->restore_subvector(std::move(vel_second), _local_second_order_indices);
272 :
273 : // Same solution update for both methods
274 369 : *_solution_update = *vel;
275 369 : _solution_update->scale(_dt);
276 :
277 : // Check for convergence by seeing if there is a nan or inf
278 369 : auto sum = _solution_update->sum();
279 : converged = std::isfinite(sum);
280 :
281 : // The linear iteration count remains zero
282 369 : _n_linear_iterations = 0;
283 369 : vel->close();
284 369 : accel->close();
285 :
286 369 : return converged;
287 369 : }
288 :
289 : void
290 29 : ExplicitMixedOrder::init()
291 : {
292 29 : ExplicitTimeIntegrator::init();
293 :
294 : // Compute ICs for velocity
295 29 : computeICs();
296 :
297 : // Seperate variables into first and second time integration order and find
298 : // the local indices for each
299 29 : const auto & var_names_first = getParam<std::vector<VariableName>>("first_order_vars");
300 58 : const auto & var_names_second = getParam<std::vector<VariableName>>("second_order_vars");
301 : std::vector<unsigned int> var_num_vec;
302 :
303 29 : auto & lm_sys = _sys.system();
304 29 : lm_sys.get_all_variable_numbers(var_num_vec);
305 29 : std::unordered_set<unsigned int> var_nums(var_num_vec.begin(), var_num_vec.end());
306 :
307 33 : for (const auto & var_name : var_names_first)
308 4 : if (lm_sys.has_variable(var_name))
309 : {
310 4 : const auto var_num = lm_sys.variable_number(var_name);
311 4 : _vars_first.insert(var_num);
312 : var_nums.erase(var_num);
313 : }
314 :
315 90 : for (const auto & var_name : var_names_second)
316 61 : if (lm_sys.has_variable(var_name))
317 : {
318 61 : const auto var_num = lm_sys.variable_number(var_name);
319 61 : _vars_second.insert(var_num);
320 : var_nums.erase(var_num);
321 : }
322 :
323 : // If var_nums is empty then that means the user has specified all the variables in this system
324 29 : if (!var_nums.empty())
325 2 : mooseError("Not all nonlinear variables have their order specified.");
326 27 : }
327 :
328 : void
329 33 : ExplicitMixedOrder::updateDOFIndices()
330 : {
331 33 : auto & lm_sys = _sys.system();
332 :
333 : std::vector<dof_id_type> var_dof_indices, work_vec;
334 36 : for (const auto var_num : _vars_first)
335 : {
336 3 : work_vec = _local_first_order_indices;
337 3 : _local_first_order_indices.clear();
338 3 : lm_sys.get_dof_map().local_variable_indices(var_dof_indices, lm_sys.get_mesh(), var_num);
339 3 : std::merge(work_vec.begin(),
340 : work_vec.end(),
341 : var_dof_indices.begin(),
342 : var_dof_indices.end(),
343 : std::back_inserter(_local_first_order_indices));
344 : }
345 :
346 33 : work_vec.clear();
347 33 : var_dof_indices.clear();
348 :
349 111 : for (const auto var_num : _vars_second)
350 : {
351 78 : work_vec = _local_second_order_indices;
352 78 : _local_second_order_indices.clear();
353 78 : lm_sys.get_dof_map().local_variable_indices(var_dof_indices, lm_sys.get_mesh(), var_num);
354 78 : std::merge(work_vec.begin(),
355 : work_vec.end(),
356 : var_dof_indices.begin(),
357 : var_dof_indices.end(),
358 : std::back_inserter(_local_second_order_indices));
359 : }
360 33 : }
361 :
362 : void
363 29 : ExplicitMixedOrder::computeICs()
364 : {
365 : // Compute the first-order approximation of the velocity at the current time step
366 : // using the Euler scheme, where the velocity is estimated as the difference
367 : // between the current solution and the previous time step, divided by the time
368 29 : auto vel = _sys.solutionUDot();
369 29 : *vel = *_solution;
370 29 : *vel -= _solution_old;
371 29 : *vel /= _dt;
372 29 : vel->close();
373 29 : }
374 :
375 : ExplicitMixedOrder::TimeOrder
376 591 : ExplicitMixedOrder::findVariableTimeOrder(unsigned int var_num) const
377 : {
378 591 : if (_vars_first.empty() && _vars_second.empty())
379 0 : mooseError("Time order sets are both empty.");
380 : if (_vars_first.count(var_num))
381 60 : return FIRST;
382 531 : else if (_vars_second.count(var_num))
383 531 : return SECOND;
384 : else
385 0 : mooseError("Variable " + _sys.system().variable_name(var_num) +
386 : " does not exist in time order sets.");
387 : }
|
