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 102 : ExplicitMixedOrder::validParams()
42 : {
43 102 : InputParameters params = ExplicitTimeIntegrator::validParams();
44 :
45 102 : params.addClassDescription(
46 : "Implementation of explicit time integration without invoking any of the nonlinear solver.");
47 :
48 204 : params.addParam<bool>("use_constant_mass",
49 204 : 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 204 : params.addParam<TagName>("mass_matrix_tag", "mass", "The tag for the mass matrix");
53 :
54 204 : params.addParam<std::vector<VariableName>>(
55 : "second_order_vars",
56 : {},
57 : "A subset of variables that require second-order integration (velocity and acceleration) to "
58 : "be applied by this time integrator.");
59 :
60 204 : params.addParam<std::vector<VariableName>>(
61 : "first_order_vars",
62 : {},
63 : "A subset of variables that require first-order integration (velocity only) to be applied by "
64 : "this time integrator.");
65 :
66 : // Prevent users from using variables option by accident.
67 102 : params.suppressParameter<std::vector<VariableName>>("variables");
68 :
69 204 : MooseEnum solve_type("consistent lumped lump_preconditioned", "lumped");
70 102 : params.setParameters("solve_type", solve_type);
71 102 : params.ignoreParameter<MooseEnum>("solve_type");
72 102 : return params;
73 102 : }
74 :
75 68 : ExplicitMixedOrder::ExplicitMixedOrder(const InputParameters & parameters)
76 : : ExplicitTimeIntegrator(parameters),
77 68 : _constant_mass(getParam<bool>("use_constant_mass")),
78 136 : _mass_matrix(getParam<TagName>("mass_matrix_tag")),
79 68 : _solution_older(_sys.solutionState(2)),
80 136 : _vars_first(declareRestartableData<std::unordered_set<unsigned int>>("first_order_vars")),
81 68 : _local_first_order_indices(
82 68 : declareRestartableData<std::vector<dof_id_type>>("first_local_indices")),
83 136 : _vars_second(declareRestartableData<std::unordered_set<unsigned int>>("second_order_vars")),
84 68 : _local_second_order_indices(
85 136 : declareRestartableData<std::vector<dof_id_type>>("second_local_indices"))
86 : {
87 68 : _fe_problem.setUDotRequested(true);
88 68 : _fe_problem.setUDotOldRequested(true);
89 68 : _fe_problem.setUDotDotRequested(true);
90 68 : }
91 :
92 : void
93 1950 : ExplicitMixedOrder::computeTimeDerivatives()
94 : {
95 : /*
96 : Because this is called in NonLinearSystemBase
97 : this should not actually compute the time derivatives.
98 : Calculating time derivatives here will cause issues for the
99 : solution update.
100 : */
101 1950 : return;
102 : }
103 :
104 : TagID
105 382 : ExplicitMixedOrder::massMatrixTagID() const
106 : {
107 382 : return _sys.subproblem().getMatrixTagID(_mass_matrix);
108 : }
109 :
110 : void
111 348 : ExplicitMixedOrder::solve()
112 : {
113 : // Getting the tagID for the mass matrix
114 348 : auto mass_tag = massMatrixTagID();
115 :
116 : // Reset iteration counts
117 348 : _n_nonlinear_iterations = 0;
118 348 : _n_linear_iterations = 0;
119 :
120 348 : _current_time = _fe_problem.time();
121 :
122 348 : auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
123 :
124 : // Compute the mass matrix
125 348 : if (!_constant_mass || _t_step == 1)
126 : {
127 : // We only want to compute "inverted" lumped mass matrix once.
128 120 : _fe_problem.computeJacobianTag(
129 120 : *_nonlinear_implicit_system->current_local_solution, mass_matrix, mass_tag);
130 :
131 : // Calculating the lumped mass matrix for use in residual calculation
132 120 : mass_matrix.vector_mult(*_mass_matrix_diag_inverted, *_ones);
133 :
134 : // "Invert" the diagonal mass matrix
135 120 : _mass_matrix_diag_inverted->reciprocal();
136 120 : _mass_matrix_diag_inverted->close();
137 : }
138 :
139 : // Set time to the time at which to evaluate the residual
140 348 : _fe_problem.time() = _fe_problem.timeOld();
141 348 : _nonlinear_implicit_system->update();
142 :
143 : // Compute the residual
144 348 : _explicit_residual->zero();
145 348 : _fe_problem.computeResidual(
146 348 : *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
147 :
148 : // Move the residual to the RHS
149 348 : *_explicit_residual *= -1.0;
150 :
151 : // Perform the linear solve
152 348 : bool converged = performExplicitSolve(mass_matrix);
153 348 : _nl->overwriteNodeFace(*_nonlinear_implicit_system->solution);
154 :
155 : // Update the solution
156 348 : *_nonlinear_implicit_system->solution = _nl->solutionOld();
157 348 : *_nonlinear_implicit_system->solution += *_solution_update;
158 :
159 348 : _nonlinear_implicit_system->update();
160 :
161 348 : _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
162 348 : _nonlinear_implicit_system->nonlinear_solver->converged = converged;
163 348 : }
164 :
165 : void
166 348 : ExplicitMixedOrder::postResidual(NumericVector<Number> & residual)
167 : {
168 348 : residual += *_Re_time;
169 348 : residual += *_Re_non_time;
170 348 : residual.close();
171 :
172 : // Reset time to the time at which to evaluate nodal BCs, which comes next
173 348 : _fe_problem.time() = _current_time;
174 348 : }
175 :
176 : bool
177 348 : ExplicitMixedOrder::performExplicitSolve(SparseMatrix<Number> &)
178 : {
179 :
180 : bool converged = false;
181 :
182 : // Grab all the vectors that we will need
183 348 : auto accel = _sys.solutionUDotDot();
184 348 : auto vel = _sys.solutionUDot();
185 :
186 : // Compute Forward Euler
187 : // Split diag mass and residual vectors into correct subvectors
188 : const std::unique_ptr<NumericVector<Number>> mass_first(
189 348 : NumericVector<Number>::build(_communicator));
190 : const std::unique_ptr<NumericVector<Real>> exp_res_first(
191 348 : NumericVector<Number>::build(_communicator));
192 348 : _mass_matrix_diag_inverted->create_subvector(*mass_first, _local_first_order_indices, false);
193 348 : _explicit_residual->create_subvector(*exp_res_first, _local_first_order_indices, false);
194 :
195 : // Need velocity vector split into subvectors
196 348 : auto vel_first = vel->get_subvector(_local_first_order_indices);
197 :
198 : // Velocity update for foward euler
199 348 : vel_first->pointwise_mult(*mass_first, *exp_res_first);
200 :
201 : // Restore the velocities
202 348 : vel->restore_subvector(std::move(vel_first), _local_first_order_indices);
203 :
204 : // Compute Central Difference
205 : // Split diag mass and residual vectors into correct subvectors
206 : const std::unique_ptr<NumericVector<Real>> mass_second(
207 348 : NumericVector<Number>::build(_communicator));
208 : const std::unique_ptr<NumericVector<Real>> exp_res_second(
209 348 : NumericVector<Number>::build(_communicator));
210 348 : _mass_matrix_diag_inverted->create_subvector(*mass_second, _local_second_order_indices, false);
211 348 : _explicit_residual->create_subvector(*exp_res_second, _local_second_order_indices, false);
212 :
213 : // Only need acceleration and old velocity vector for central difference
214 348 : auto accel_second = accel->get_subvector(_local_second_order_indices);
215 :
216 348 : auto vel_second = vel->get_subvector(_local_second_order_indices);
217 :
218 : // Compute acceleration for central difference
219 348 : accel_second->pointwise_mult(*mass_second, *exp_res_second);
220 :
221 : // Scaling the acceleration
222 348 : auto accel_scaled = accel_second->clone();
223 348 : accel_scaled->scale((_dt + _dt_old) / 2);
224 :
225 : // Velocity update for central difference
226 348 : *vel_second += *accel_scaled;
227 :
228 : // Restore acceleration
229 348 : accel->restore_subvector(std::move(accel_second), _local_second_order_indices);
230 :
231 348 : vel->restore_subvector(std::move(vel_second), _local_second_order_indices);
232 :
233 : // Same solution update for both methods
234 348 : *_solution_update = *vel;
235 348 : _solution_update->scale(_dt);
236 :
237 : // Check for convergence by seeing if there is a nan or inf
238 348 : auto sum = _solution_update->sum();
239 : converged = std::isfinite(sum);
240 :
241 : // The linear iteration count remains zero
242 348 : _n_linear_iterations = 0;
243 348 : vel->close();
244 348 : accel->close();
245 :
246 348 : return converged;
247 348 : }
248 :
249 : void
250 34 : ExplicitMixedOrder::init()
251 : {
252 34 : ExplicitTimeIntegrator::init();
253 :
254 : // Compute ICs for velocity
255 34 : computeICs();
256 :
257 : // Seperate variables into first and second time integration order and find
258 : // the local indices for each
259 34 : const auto & var_names_first = getParam<std::vector<VariableName>>("first_order_vars");
260 68 : const auto & var_names_second = getParam<std::vector<VariableName>>("second_order_vars");
261 : std::vector<unsigned int> var_num_vec;
262 :
263 34 : auto & lm_sys = _sys.system();
264 34 : lm_sys.get_all_variable_numbers(var_num_vec);
265 34 : std::unordered_set<unsigned int> var_nums(var_num_vec.begin(), var_num_vec.end());
266 :
267 42 : for (const auto & var_name : var_names_first)
268 8 : if (lm_sys.has_variable(var_name))
269 : {
270 8 : const auto var_num = lm_sys.variable_number(var_name);
271 8 : _vars_first.insert(var_num);
272 : var_nums.erase(var_num);
273 : }
274 :
275 84 : for (const auto & var_name : var_names_second)
276 50 : if (lm_sys.has_variable(var_name))
277 : {
278 50 : const auto var_num = lm_sys.variable_number(var_name);
279 50 : _vars_second.insert(var_num);
280 : var_nums.erase(var_num);
281 : }
282 :
283 : // If var_nums is empty then that means the user has specified all the variables in this system
284 34 : if (!var_nums.empty())
285 4 : mooseError("Not all nonlinear variables have their order specified.");
286 :
287 : std::vector<dof_id_type> var_dof_indices, work_vec;
288 36 : for (const auto var_num : _vars_first)
289 : {
290 6 : work_vec = _local_first_order_indices;
291 6 : _local_first_order_indices.clear();
292 6 : lm_sys.get_dof_map().local_variable_indices(var_dof_indices, lm_sys.get_mesh(), var_num);
293 6 : std::merge(work_vec.begin(),
294 : work_vec.end(),
295 : var_dof_indices.begin(),
296 : var_dof_indices.end(),
297 : std::back_inserter(_local_first_order_indices));
298 : }
299 :
300 : work_vec.clear();
301 : var_dof_indices.clear();
302 :
303 78 : for (const auto var_num : _vars_second)
304 : {
305 48 : work_vec = _local_second_order_indices;
306 48 : _local_second_order_indices.clear();
307 48 : lm_sys.get_dof_map().local_variable_indices(var_dof_indices, lm_sys.get_mesh(), var_num);
308 48 : std::merge(work_vec.begin(),
309 : work_vec.end(),
310 : var_dof_indices.begin(),
311 : var_dof_indices.end(),
312 : std::back_inserter(_local_second_order_indices));
313 : }
314 30 : }
315 :
316 : void
317 34 : ExplicitMixedOrder::computeICs()
318 : {
319 : // Compute the first-order approximation of the velocity at the current time step
320 : // using the Euler scheme, where the velocity is estimated as the difference
321 : // between the current solution and the previous time step, divided by the time
322 34 : auto vel = _sys.solutionUDot();
323 34 : *vel = *_solution;
324 34 : *vel -= _solution_old;
325 34 : *vel /= _dt;
326 34 : vel->close();
327 34 : }
328 :
329 : ExplicitMixedOrder::TimeOrder
330 372 : ExplicitMixedOrder::findVariableTimeOrder(unsigned int var_num) const
331 : {
332 372 : if (_vars_first.empty() && _vars_second.empty())
333 0 : mooseError("Time order sets are both empty.");
334 : if (_vars_first.count(var_num))
335 120 : return FIRST;
336 252 : else if (_vars_second.count(var_num))
337 252 : return SECOND;
338 : else
339 0 : mooseError("Variable " + _sys.system().variable_name(var_num) +
340 : " does not exist in time order sets.");
341 : }
|
