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 "MooseApp.h"
11 : #include "SystemBase.h"
12 : #include "Factory.h"
13 : #include "SubProblem.h"
14 : #include "MooseVariableFE.h"
15 : #include "MooseVariableFV.h"
16 : #include "MooseVariableScalar.h"
17 : #include "MooseVariableConstMonomial.h"
18 : #include "Conversion.h"
19 : #include "Parser.h"
20 : #include "AllLocalDofIndicesThread.h"
21 : #include "MooseTypes.h"
22 : #include "InitialCondition.h"
23 : #include "ScalarInitialCondition.h"
24 : #include "Assembly.h"
25 : #include "MooseMesh.h"
26 : #include "MooseUtils.h"
27 : #include "FVBoundaryCondition.h"
28 : #include "FEProblemBase.h"
29 : #include "TimeIntegrator.h"
30 :
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/string_to_enum.h"
33 : #include "libmesh/fe_interface.h"
34 : #include "libmesh/static_condensation.h"
35 :
36 : using namespace libMesh;
37 :
38 : /// Free function used for a libMesh callback
39 : void
40 65414 : extraSendList(std::vector<dof_id_type> & send_list, void * context)
41 : {
42 65414 : SystemBase * sys = static_cast<SystemBase *>(context);
43 65414 : sys->augmentSendList(send_list);
44 65414 : }
45 :
46 : /// Free function used for a libMesh callback
47 : void
48 68579 : extraSparsity(SparsityPattern::Graph & sparsity,
49 : std::vector<dof_id_type> & n_nz,
50 : std::vector<dof_id_type> & n_oz,
51 : void * context)
52 : {
53 68579 : SystemBase * sys = static_cast<SystemBase *>(context);
54 68579 : sys->augmentSparsity(sparsity, n_nz, n_oz);
55 68579 : }
56 :
57 131304 : SystemBase::SystemBase(SubProblem & subproblem,
58 : FEProblemBase & fe_problem,
59 : const std::string & name,
60 131304 : Moose::VarKindType var_kind)
61 : : libMesh::ParallelObject(subproblem),
62 : ConsoleStreamInterface(subproblem.getMooseApp()),
63 131304 : _subproblem(subproblem),
64 131304 : _fe_problem(fe_problem),
65 262608 : _app(subproblem.getMooseApp()),
66 131304 : _factory(_app.getFactory()),
67 131304 : _mesh(subproblem.mesh()),
68 131304 : _name(name),
69 262608 : _vars(libMesh::n_threads()),
70 131304 : _var_map(),
71 131304 : _max_var_number(0),
72 131304 : _u_dot(nullptr),
73 131304 : _u_dotdot(nullptr),
74 131304 : _u_dot_old(nullptr),
75 131304 : _u_dotdot_old(nullptr),
76 131304 : _saved_old(nullptr),
77 131304 : _saved_older(nullptr),
78 131304 : _saved_dot_old(nullptr),
79 131304 : _saved_dotdot_old(nullptr),
80 131304 : _var_kind(var_kind),
81 131304 : _max_var_n_dofs_per_elem(0),
82 131304 : _max_var_n_dofs_per_node(0),
83 131304 : _automatic_scaling(false),
84 131304 : _verbose(false),
85 262608 : _solution_states_initialized(false)
86 : {
87 131304 : }
88 :
89 : MooseVariableFieldBase &
90 4885378 : SystemBase::getVariable(THREAD_ID tid, const std::string & var_name) const
91 : {
92 : MooseVariableFieldBase * var =
93 4885378 : dynamic_cast<MooseVariableFieldBase *>(_vars[tid].getVariable(var_name));
94 4885378 : if (!var)
95 4 : mooseError("Variable '", var_name, "' does not exist in this system");
96 4885374 : return *var;
97 : }
98 :
99 : MooseVariableFieldBase &
100 169801494 : SystemBase::getVariable(THREAD_ID tid, unsigned int var_number) const
101 : {
102 169801494 : if (var_number < _numbered_vars[tid].size())
103 169801494 : if (_numbered_vars[tid][var_number])
104 169801494 : return *_numbered_vars[tid][var_number];
105 :
106 0 : mooseError("Variable #", Moose::stringify(var_number), " does not exist in this system");
107 : }
108 :
109 : template <typename T>
110 : MooseVariableFE<T> &
111 60158 : SystemBase::getFieldVariable(THREAD_ID tid, const std::string & var_name)
112 : {
113 60158 : return *_vars[tid].getFieldVariable<T>(var_name);
114 : }
115 :
116 : template <typename T>
117 : MooseVariableField<T> &
118 110896 : SystemBase::getActualFieldVariable(THREAD_ID tid, const std::string & var_name)
119 : {
120 110896 : return *_vars[tid].getActualFieldVariable<T>(var_name);
121 : }
122 :
123 : template <typename T>
124 : MooseVariableFV<T> &
125 690 : SystemBase::getFVVariable(THREAD_ID tid, const std::string & var_name)
126 : {
127 690 : return *_vars[tid].getFVVariable<T>(var_name);
128 : }
129 :
130 : template <typename T>
131 : MooseVariableFE<T> &
132 0 : SystemBase::getFieldVariable(THREAD_ID tid, unsigned int var_number)
133 : {
134 0 : return *_vars[tid].getFieldVariable<T>(var_number);
135 : }
136 :
137 : template <typename T>
138 : MooseVariableField<T> &
139 147053271 : SystemBase::getActualFieldVariable(THREAD_ID tid, unsigned int var_number)
140 : {
141 147053271 : return *_vars[tid].getActualFieldVariable<T>(var_number);
142 : }
143 :
144 : MooseVariableScalar &
145 60040 : SystemBase::getScalarVariable(THREAD_ID tid, const std::string & var_name) const
146 : {
147 60040 : MooseVariableScalar * var = dynamic_cast<MooseVariableScalar *>(_vars[tid].getVariable(var_name));
148 60040 : if (!var)
149 8 : mooseError("Scalar variable '" + var_name + "' does not exist in this system");
150 60032 : return *var;
151 : }
152 :
153 : MooseVariableScalar &
154 207789 : SystemBase::getScalarVariable(THREAD_ID tid, unsigned int var_number) const
155 : {
156 : MooseVariableScalar * var =
157 207789 : dynamic_cast<MooseVariableScalar *>(_vars[tid].getVariable(var_number));
158 207789 : if (!var)
159 0 : mooseError("variable #" + Moose::stringify(var_number) + " does not exist in this system");
160 207789 : return *var;
161 : }
162 :
163 : const std::set<SubdomainID> *
164 344 : SystemBase::getVariableBlocks(unsigned int var_number)
165 : {
166 : mooseAssert(_var_map.find(var_number) != _var_map.end(), "Variable does not exist.");
167 344 : if (_var_map[var_number].empty())
168 344 : return nullptr;
169 : else
170 0 : return &_var_map[var_number];
171 : }
172 :
173 : void
174 509 : SystemBase::addVariableToZeroOnResidual(std::string var_name)
175 : {
176 509 : _vars_to_be_zeroed_on_residual.push_back(var_name);
177 509 : }
178 :
179 : void
180 339 : SystemBase::addVariableToZeroOnJacobian(std::string var_name)
181 : {
182 339 : _vars_to_be_zeroed_on_jacobian.push_back(var_name);
183 339 : }
184 :
185 : void
186 116 : SystemBase::setVariableGlobalDoFs(const std::string & var_name)
187 : {
188 348 : AllLocalDofIndicesThread aldit(_subproblem, {var_name});
189 116 : ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
190 116 : Threads::parallel_reduce(elem_range, aldit);
191 :
192 : // Gather the dof indices across procs to get all the dof indices for var_name
193 116 : aldit.dofIndicesSetUnion();
194 :
195 116 : const auto & all_dof_indices = aldit.getDofIndices();
196 116 : _var_all_dof_indices.assign(all_dof_indices.begin(), all_dof_indices.end());
197 232 : }
198 :
199 : void
200 3826133 : SystemBase::zeroVariables(std::vector<std::string> & vars_to_be_zeroed)
201 : {
202 3826133 : if (vars_to_be_zeroed.size() > 0)
203 : {
204 10142 : NumericVector<Number> & solution = this->solution();
205 :
206 10142 : auto problem = dynamic_cast<FEProblemBase *>(&_subproblem);
207 10142 : if (!problem)
208 0 : mooseError("System needs to be registered in FEProblemBase for using zeroVariables.");
209 :
210 10142 : AllLocalDofIndicesThread aldit(*problem, vars_to_be_zeroed, true);
211 10142 : ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
212 10142 : Threads::parallel_reduce(elem_range, aldit);
213 :
214 10142 : const auto & dof_indices_to_zero = aldit.getDofIndices();
215 :
216 10142 : solution.close();
217 :
218 3082797 : for (const auto & dof : dof_indices_to_zero)
219 3072655 : solution.set(dof, 0);
220 :
221 10142 : solution.close();
222 :
223 : // Call update to update the current_local_solution for this system
224 10142 : system().update();
225 10142 : }
226 3826133 : }
227 :
228 : void
229 3296816 : SystemBase::zeroVariablesForResidual()
230 : {
231 3296816 : zeroVariables(_vars_to_be_zeroed_on_residual);
232 3296816 : }
233 :
234 : void
235 519693 : SystemBase::zeroVariablesForJacobian()
236 : {
237 519693 : zeroVariables(_vars_to_be_zeroed_on_jacobian);
238 519693 : }
239 :
240 : Order
241 62676 : SystemBase::getMinQuadratureOrder()
242 : {
243 62676 : Order order = CONSTANT;
244 62676 : const std::vector<MooseVariableFieldBase *> & vars = _vars[0].fieldVariables();
245 124120 : for (const auto & var : vars)
246 : {
247 61444 : FEType fe_type = var->feType();
248 61444 : if (fe_type.default_quadrature_order() > order)
249 50359 : order = fe_type.default_quadrature_order();
250 : }
251 :
252 62676 : return order;
253 : }
254 :
255 : void
256 879636850 : SystemBase::prepare(THREAD_ID tid)
257 : {
258 879636850 : if (_subproblem.hasActiveElementalMooseVariables(tid))
259 : {
260 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
261 845703393 : _subproblem.getActiveElementalMooseVariables(tid);
262 845703393 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
263 1646426752 : for (const auto & var : vars)
264 800723359 : var->clearDofIndices();
265 :
266 1899911028 : for (const auto & var : active_elemental_moose_variables)
267 1054207635 : if (&(var->sys()) == this)
268 514606187 : var->prepare();
269 : }
270 : else
271 : {
272 33933457 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
273 70436184 : for (const auto & var : vars)
274 36502727 : var->prepare();
275 : }
276 879636850 : }
277 :
278 : void
279 58958 : SystemBase::prepareFace(THREAD_ID tid, bool resize_data)
280 : {
281 : // We only need to do something if the element prepare was restricted
282 58958 : if (_subproblem.hasActiveElementalMooseVariables(tid))
283 : {
284 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
285 36574 : _subproblem.getActiveElementalMooseVariables(tid);
286 :
287 36574 : std::vector<MooseVariableFieldBase *> newly_prepared_vars;
288 :
289 36574 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
290 66881 : for (const auto & var : vars)
291 : {
292 : mooseAssert(&var->sys() == this,
293 : "I will cry if we store variables in our warehouse that don't belong to us");
294 :
295 : // If it wasn't in the active list, we need to prepare it. This has the potential to duplicate
296 : // prepare if we have these conditions:
297 : //
298 : // 1. We have a displaced problem
299 : // 2. We are using AD
300 : // 3. We are not using global AD indexing
301 : //
302 : // But I think I would rather risk duplicate prepare than introduce an additional member set
303 : // variable for tracking prepared variables. Set insertion is slow and some simulations have a
304 : // ton of variables
305 30307 : if (!active_elemental_moose_variables.count(var))
306 : {
307 3646 : var->prepare();
308 3646 : newly_prepared_vars.push_back(var);
309 : }
310 : }
311 :
312 : // Make sure to resize the residual and jacobian datastructures for all the new variables
313 36574 : if (resize_data)
314 19912 : for (const auto var_ptr : newly_prepared_vars)
315 : {
316 1625 : _subproblem.assembly(tid, number()).prepareVariable(var_ptr);
317 1625 : if (_subproblem.checkNonlocalCouplingRequirement())
318 0 : _subproblem.assembly(tid, number()).prepareVariableNonlocal(var_ptr);
319 : }
320 36574 : }
321 58958 : }
322 :
323 : void
324 9248134 : SystemBase::prepareNeighbor(THREAD_ID tid)
325 : {
326 9248134 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
327 23724231 : for (const auto & var : vars)
328 14476097 : var->prepareNeighbor();
329 9248134 : }
330 :
331 : void
332 652032 : SystemBase::prepareLowerD(THREAD_ID tid)
333 : {
334 652032 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
335 1293713 : for (const auto & var : vars)
336 641681 : var->prepareLowerD();
337 652032 : }
338 :
339 : void
340 442178865 : SystemBase::reinitElem(const Elem * const elem, THREAD_ID tid)
341 : {
342 442178865 : if (_subproblem.hasActiveElementalMooseVariables(tid))
343 : {
344 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
345 425442884 : _subproblem.getActiveElementalMooseVariables(tid);
346 957199962 : for (const auto & var : active_elemental_moose_variables)
347 531757078 : if (&(var->sys()) == this)
348 437064780 : var->computeElemValues();
349 : }
350 : else
351 : {
352 16735981 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
353 33085914 : for (const auto & var : vars)
354 16349933 : var->computeElemValues();
355 : }
356 :
357 442178865 : if (system().has_static_condensation())
358 1159995 : for (auto & [tag, matrix] : _active_tagged_matrices)
359 : {
360 255139 : libmesh_ignore(tag);
361 255139 : cast_ptr<StaticCondensation *>(matrix)->set_current_elem(*elem);
362 : }
363 442178865 : }
364 :
365 : void
366 13259811 : SystemBase::reinitElemFace(const Elem * /*elem*/, unsigned int /*side*/, THREAD_ID tid)
367 : {
368 13259811 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
369 35754299 : for (const auto & var : vars)
370 22494488 : var->computeElemValuesFace();
371 13259811 : }
372 :
373 : void
374 9202694 : SystemBase::reinitNeighborFace(const Elem * /*elem*/, unsigned int /*side*/, THREAD_ID tid)
375 : {
376 9202694 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
377 23645511 : for (const auto & var : vars)
378 14442817 : var->computeNeighborValuesFace();
379 9202694 : }
380 :
381 : void
382 45440 : SystemBase::reinitNeighbor(const Elem * /*elem*/, THREAD_ID tid)
383 : {
384 45440 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
385 78720 : for (const auto & var : vars)
386 33280 : var->computeNeighborValues();
387 45440 : }
388 :
389 : void
390 652032 : SystemBase::reinitLowerD(THREAD_ID tid)
391 : {
392 652032 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
393 1293713 : for (const auto & var : vars)
394 641681 : var->computeLowerDValues();
395 652032 : }
396 :
397 : void
398 60356192 : SystemBase::reinitNode(const Node * /*node*/, THREAD_ID tid)
399 : {
400 60356192 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
401 164461072 : for (const auto & var : vars)
402 : {
403 104104880 : var->reinitNode();
404 104104880 : if (var->isNodalDefined())
405 87974034 : var->computeNodalValues();
406 : }
407 60356192 : }
408 :
409 : void
410 169868338 : SystemBase::reinitNodeFace(const Node * /*node*/, BoundaryID /*bnd_id*/, THREAD_ID tid)
411 : {
412 169868338 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
413 413422589 : for (const auto & var : vars)
414 : {
415 243554251 : var->reinitNode();
416 243554251 : if (var->isNodalDefined())
417 216829859 : var->computeNodalValues();
418 : }
419 169868338 : }
420 :
421 : void
422 13626 : SystemBase::reinitNodes(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
423 : {
424 13626 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
425 20604 : for (const auto & var : vars)
426 : {
427 6978 : var->reinitNodes(nodes);
428 6978 : var->computeNodalValues();
429 : }
430 13626 : }
431 :
432 : void
433 1446 : SystemBase::reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
434 : {
435 1446 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
436 2334 : for (const auto & var : vars)
437 : {
438 888 : var->reinitNodesNeighbor(nodes);
439 888 : var->computeNodalNeighborValues();
440 : }
441 1446 : }
442 :
443 : void
444 18527453 : SystemBase::reinitScalars(THREAD_ID tid, bool reinit_for_derivative_reordering /*=false*/)
445 : {
446 18527453 : const std::vector<MooseVariableScalar *> & vars = _vars[tid].scalars();
447 19053260 : for (const auto & var : vars)
448 525807 : var->reinit(reinit_for_derivative_reordering);
449 18527453 : }
450 :
451 : void
452 65414 : SystemBase::augmentSendList(std::vector<dof_id_type> & send_list)
453 : {
454 65414 : std::set<dof_id_type> & ghosted_elems = _subproblem.ghostedElems();
455 :
456 65414 : DofMap & dof_map = dofMap();
457 :
458 65414 : std::vector<dof_id_type> dof_indices;
459 :
460 65414 : System & sys = system();
461 :
462 65414 : unsigned int sys_num = sys.number();
463 :
464 65414 : unsigned int n_vars = sys.n_vars();
465 :
466 156187 : for (const auto & elem_id : ghosted_elems)
467 : {
468 90773 : Elem * elem = _mesh.elemPtr(elem_id);
469 :
470 90773 : if (elem->active())
471 : {
472 90773 : dof_map.dof_indices(elem, dof_indices);
473 :
474 : // Only need to ghost it if it's actually not on this processor
475 1188329 : for (const auto & dof : dof_indices)
476 1097556 : if (dof < dof_map.first_dof() || dof >= dof_map.end_dof())
477 1073406 : send_list.push_back(dof);
478 :
479 : // Now add the DoFs from all of the nodes. This is necessary because of block
480 : // restricted variables. A variable might not live _on_ this element but it
481 : // might live on nodes connected to this element.
482 508467 : for (unsigned int n = 0; n < elem->n_nodes(); n++)
483 : {
484 417694 : Node * node = elem->node_ptr(n);
485 :
486 : // Have to get each variable's dofs
487 1629124 : for (unsigned int v = 0; v < n_vars; v++)
488 : {
489 1211430 : const Variable & var = sys.variable(v);
490 1211430 : unsigned int var_num = var.number();
491 1211430 : unsigned int n_comp = var.n_components();
492 :
493 : // See if this variable has any dofs at this node
494 1211430 : if (node->n_dofs(sys_num, var_num) > 0)
495 : {
496 : // Loop over components of the variable
497 2128568 : for (unsigned int c = 0; c < n_comp; c++)
498 1064284 : send_list.push_back(node->dof_number(sys_num, var_num, c));
499 : }
500 : }
501 : }
502 : }
503 : }
504 65414 : }
505 :
506 : /**
507 : * Save the old and older solutions.
508 : */
509 : void
510 252 : SystemBase::saveOldSolutions()
511 : {
512 : const auto states =
513 252 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Time)].size();
514 252 : if (states > 1)
515 : {
516 230 : _saved_solution_states.resize(states);
517 690 : for (unsigned int i = 1; i <= states - 1; ++i)
518 460 : if (!_saved_solution_states[i])
519 920 : _saved_solution_states[i] =
520 460 : &addVector("save_solution_state_" + std::to_string(i), false, PARALLEL);
521 :
522 690 : for (unsigned int i = 1; i <= states - 1; ++i)
523 460 : *(_saved_solution_states[i]) = solutionState(i);
524 : }
525 :
526 252 : if (!_saved_dot_old && solutionUDotOld())
527 0 : _saved_dot_old = &addVector("save_solution_dot_old", false, PARALLEL);
528 252 : if (!_saved_dotdot_old && solutionUDotDotOld())
529 0 : _saved_dotdot_old = &addVector("save_solution_dotdot_old", false, PARALLEL);
530 :
531 252 : if (solutionUDotOld())
532 0 : *_saved_dot_old = *solutionUDotOld();
533 :
534 252 : if (solutionUDotDotOld())
535 0 : *_saved_dotdot_old = *solutionUDotDotOld();
536 252 : }
537 :
538 : /**
539 : * Restore the old and older solutions when the saved solutions present.
540 : */
541 : void
542 252 : SystemBase::restoreOldSolutions()
543 : {
544 : const auto states =
545 252 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Time)].size();
546 252 : if (states > 1)
547 690 : for (unsigned int i = 1; i <= states - 1; ++i)
548 460 : if (_saved_solution_states[i])
549 : {
550 460 : solutionState(i) = *(_saved_solution_states[i]);
551 460 : removeVector("save_solution_state_" + std::to_string(i));
552 460 : _saved_solution_states[i] = nullptr;
553 : }
554 :
555 252 : if (_saved_dot_old && solutionUDotOld())
556 : {
557 0 : *solutionUDotOld() = *_saved_dot_old;
558 0 : removeVector("save_solution_dot_old");
559 0 : _saved_dot_old = nullptr;
560 : }
561 252 : if (_saved_dotdot_old && solutionUDotDotOld())
562 : {
563 0 : *solutionUDotDotOld() = *_saved_dotdot_old;
564 0 : removeVector("save_solution_dotdot_old");
565 0 : _saved_dotdot_old = nullptr;
566 : }
567 252 : }
568 :
569 : SparseMatrix<Number> &
570 348 : SystemBase::addMatrix(TagID tag)
571 : {
572 348 : if (!_subproblem.matrixTagExists(tag))
573 0 : mooseError("Cannot add tagged matrix with TagID ",
574 : tag,
575 : " in system '",
576 0 : name(),
577 : "' because the tag does not exist in the problem");
578 :
579 348 : if (hasMatrix(tag))
580 0 : return getMatrix(tag);
581 :
582 348 : const auto matrix_name = _subproblem.matrixTagName(tag);
583 348 : SparseMatrix<Number> & mat = system().add_matrix(matrix_name);
584 348 : associateMatrixToTag(mat, tag);
585 :
586 348 : return mat;
587 348 : }
588 :
589 : void
590 0 : SystemBase::removeMatrix(TagID tag_id)
591 : {
592 0 : if (!_subproblem.matrixTagExists(tag_id))
593 0 : mooseError("Cannot remove the matrix with TagID ",
594 : tag_id,
595 : "\nin system '",
596 0 : name(),
597 : "', because that tag does not exist in the problem");
598 :
599 0 : if (hasMatrix(tag_id))
600 : {
601 0 : const auto matrix_name = _subproblem.matrixTagName(tag_id);
602 0 : system().remove_matrix(matrix_name);
603 0 : _tagged_matrices[tag_id] = nullptr;
604 0 : }
605 0 : }
606 :
607 : NumericVector<Number> &
608 69362 : SystemBase::addVector(const std::string & vector_name, const bool project, const ParallelType type)
609 : {
610 69362 : if (hasVector(vector_name))
611 449 : return getVector(vector_name);
612 :
613 68913 : NumericVector<Number> & vec = system().add_vector(vector_name, project, type);
614 68913 : return vec;
615 : }
616 :
617 : NumericVector<Number> &
618 164787 : SystemBase::addVector(TagID tag, const bool project, const ParallelType type)
619 : {
620 164787 : if (!_subproblem.vectorTagExists(tag))
621 0 : mooseError("Cannot add tagged vector with TagID ",
622 : tag,
623 : " in system '",
624 0 : name(),
625 : "' because the tag does not exist in the problem");
626 :
627 164787 : if (hasVector(tag))
628 : {
629 0 : auto & vec = getVector(tag);
630 :
631 0 : if (type != ParallelType::AUTOMATIC && vec.type() != type)
632 0 : mooseError("Cannot add tagged vector '",
633 0 : _subproblem.vectorTagName(tag),
634 : "', in system '",
635 0 : name(),
636 : "' because a vector with the same name was found with a different parallel type");
637 :
638 0 : return vec;
639 : }
640 :
641 164787 : const auto vector_name = _subproblem.vectorTagName(tag);
642 164787 : NumericVector<Number> & vec = system().add_vector(vector_name, project, type);
643 164787 : associateVectorToTag(vec, tag);
644 :
645 164787 : return vec;
646 164787 : }
647 :
648 : void
649 19024195 : SystemBase::closeTaggedVector(const TagID tag)
650 : {
651 19024195 : if (!_subproblem.vectorTagExists(tag))
652 0 : mooseError("Cannot close vector with TagID ",
653 : tag,
654 : " in system '",
655 0 : name(),
656 : "' because that tag does not exist in the problem");
657 19024195 : else if (!hasVector(tag))
658 0 : mooseError("Cannot close vector tag with name '",
659 0 : _subproblem.vectorTagName(tag),
660 : "' in system '",
661 0 : name(),
662 : "' because there is no vector associated with that tag");
663 19024195 : getVector(tag).close();
664 19024195 : }
665 :
666 : void
667 6593749 : SystemBase::closeTaggedVectors(const std::set<TagID> & tags)
668 : {
669 25617944 : for (const auto tag : tags)
670 19024195 : closeTaggedVector(tag);
671 6593749 : }
672 :
673 : void
674 9512434 : SystemBase::zeroTaggedVector(const TagID tag)
675 : {
676 9512434 : if (!_subproblem.vectorTagExists(tag))
677 0 : mooseError("Cannot zero vector with TagID ",
678 : tag,
679 : " in system '",
680 0 : name(),
681 : "' because that tag does not exist in the problem");
682 9512434 : else if (!hasVector(tag))
683 0 : mooseError("Cannot zero vector tag with name '",
684 0 : _subproblem.vectorTagName(tag),
685 : "' in system '",
686 0 : name(),
687 : "' because there is no vector associated with that tag");
688 9512434 : if (!_subproblem.vectorTagNotZeroed(tag))
689 9512412 : getVector(tag).zero();
690 9512434 : }
691 :
692 : void
693 3296947 : SystemBase::zeroTaggedVectors(const std::set<TagID> & tags)
694 : {
695 12809381 : for (const auto tag : tags)
696 9512434 : zeroTaggedVector(tag);
697 3296947 : }
698 :
699 : void
700 0 : SystemBase::removeVector(TagID tag_id)
701 : {
702 0 : if (!_subproblem.vectorTagExists(tag_id))
703 0 : mooseError("Cannot remove the vector with TagID ",
704 : tag_id,
705 : "\nin system '",
706 0 : name(),
707 : "', because that tag does not exist in the problem");
708 :
709 0 : if (hasVector(tag_id))
710 : {
711 0 : auto vector_name = _subproblem.vectorTagName(tag_id);
712 0 : system().remove_vector(vector_name);
713 0 : _tagged_vectors[tag_id] = nullptr;
714 0 : }
715 0 : }
716 :
717 : void
718 174087 : SystemBase::addVariable(const std::string & var_type,
719 : const std::string & name,
720 : InputParameters & parameters)
721 : {
722 174087 : _numbered_vars.resize(libMesh::n_threads());
723 :
724 174087 : const auto components = parameters.get<unsigned int>("components");
725 :
726 : // Convert the std::vector parameter provided by the user into a std::set for use by libMesh's
727 : // System::add_variable method
728 174087 : std::set<SubdomainID> blocks;
729 174087 : const auto & block_param = parameters.get<std::vector<SubdomainName>>("block");
730 185838 : for (const auto & subdomain_name : block_param)
731 : {
732 11751 : SubdomainID blk_id = _mesh.getSubdomainID(subdomain_name);
733 11751 : blocks.insert(blk_id);
734 : }
735 :
736 : const auto fe_type =
737 348174 : FEType(Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")),
738 522261 : Utility::string_to_enum<FEFamily>(parameters.get<MooseEnum>("family")));
739 174087 : const auto fe_field_type = FEInterface::field_type(fe_type);
740 :
741 : unsigned int var_num;
742 :
743 174087 : if (var_type == "ArrayMooseVariable")
744 : {
745 2375 : if (fe_field_type == TYPE_VECTOR)
746 0 : mooseError("Vector family type cannot be used in an array variable");
747 :
748 2375 : std::vector<std::string> array_var_component_names;
749 2375 : const bool has_array_names = parameters.isParamValid("array_var_component_names");
750 2375 : if (has_array_names)
751 : {
752 : array_var_component_names =
753 160 : parameters.get<std::vector<std::string>>("array_var_component_names");
754 160 : if (array_var_component_names.size() != components)
755 8 : parameters.paramError("array_var_component_names",
756 : "Must be the same size as 'components' (size ",
757 : components,
758 : ") for array variable '",
759 : name,
760 : "'");
761 : }
762 :
763 : // Build up the variable names
764 2371 : std::vector<std::string> var_names;
765 7779 : for (unsigned int i = 0; i < components; i++)
766 : {
767 5408 : if (!has_array_names)
768 4849 : array_var_component_names.push_back(std::to_string(i));
769 5408 : var_names.push_back(name + "_" + array_var_component_names[i]);
770 : }
771 :
772 : // makes sure there is always a name, either the provided one or '1 2 3 ...'
773 4742 : parameters.set<std::vector<std::string>>("array_var_component_names") =
774 2371 : array_var_component_names;
775 :
776 : // The number returned by libMesh is the _last_ variable number... we want to hold onto the
777 : // _first_
778 2371 : var_num = system().add_variables(var_names, fe_type, &blocks) - (components - 1);
779 :
780 : // Set as array variable
781 7113 : if (parameters.isParamSetByUser("array") && !parameters.get<bool>("array"))
782 8 : parameters.paramError("array",
783 : "Must be set to true for variable '",
784 : name,
785 : "' because 'components' > 1 (is an array variable)");
786 2367 : parameters.set<bool>("array") = true;
787 2367 : }
788 : else
789 : {
790 343424 : if (parameters.isParamSetByUser("array_var_component_names"))
791 8 : parameters.paramError("array_var_component_names",
792 : "Should not be set because this variable (",
793 : name,
794 : ") is a non-array variable");
795 171708 : var_num = system().add_variable(name, fe_type, &blocks);
796 : }
797 :
798 348150 : parameters.set<unsigned int>("_var_num") = var_num;
799 348150 : parameters.set<SystemBase *>("_system_base") = this;
800 :
801 364229 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
802 : {
803 380316 : parameters.set<THREAD_ID>("tid") = tid;
804 : std::shared_ptr<MooseVariableBase> var =
805 190158 : _factory.create<MooseVariableBase>(var_type, name, parameters, tid);
806 :
807 190154 : _vars[tid].add(name, var);
808 :
809 190154 : if (auto fe_var = dynamic_cast<MooseVariableFieldBase *>(var.get()))
810 : {
811 186614 : auto required_size = var_num + components;
812 186614 : if (required_size > _numbered_vars[tid].size())
813 186614 : _numbered_vars[tid].resize(required_size);
814 376588 : for (MooseIndex(components) component = 0; component < components; ++component)
815 189974 : _numbered_vars[tid][var_num + component] = fe_var;
816 :
817 186614 : if (auto * const functor = dynamic_cast<Moose::FunctorBase<ADReal> *>(fe_var))
818 182326 : _subproblem.addFunctor(name, *functor, tid);
819 4288 : else if (auto * const functor = dynamic_cast<Moose::FunctorBase<ADRealVectorValue> *>(fe_var))
820 1666 : _subproblem.addFunctor(name, *functor, tid);
821 2622 : else if (auto * const functor = dynamic_cast<Moose::FunctorBase<RealEigenVector> *>(fe_var))
822 2622 : _subproblem.addFunctor(name, *functor, tid);
823 : else
824 0 : mooseError("This should be a functor");
825 : }
826 :
827 190154 : if (auto scalar_var = dynamic_cast<MooseVariableScalar *>(var.get()))
828 : {
829 3540 : if (auto * const functor = dynamic_cast<Moose::FunctorBase<ADReal> *>(scalar_var))
830 3540 : _subproblem.addFunctor(name, *functor, tid);
831 : else
832 0 : mooseError("Scalar variables should be functors");
833 : }
834 :
835 190154 : if (var->blockRestricted())
836 22890 : for (const SubdomainID & id : var->blockIDs())
837 25842 : for (MooseIndex(components) component = 0; component < components; ++component)
838 13150 : _var_map[var_num + component].insert(id);
839 : else
840 362814 : for (MooseIndex(components) component = 0; component < components; ++component)
841 182858 : _var_map[var_num + component] = std::set<SubdomainID>();
842 190154 : }
843 :
844 : // getMaxVariableNumber is an API method used in Rattlesnake
845 174071 : if (var_num > _max_var_number)
846 92888 : _max_var_number = var_num;
847 174071 : _du_dot_du.resize(var_num + 1);
848 174071 : }
849 :
850 : bool
851 6450011 : SystemBase::hasVariable(const std::string & var_name) const
852 : {
853 6450011 : auto & names = getVariableNames();
854 6450011 : if (system().has_variable(var_name))
855 5073575 : return system().variable_type(var_name).family != SCALAR;
856 1376436 : if (std::find(names.begin(), names.end(), var_name) != names.end())
857 : // array variable
858 36173 : return true;
859 : else
860 1340263 : return false;
861 : }
862 :
863 : bool
864 0 : SystemBase::isArrayVariable(const std::string & var_name) const
865 : {
866 0 : auto & names = getVariableNames();
867 0 : if (!system().has_variable(var_name) &&
868 0 : std::find(names.begin(), names.end(), var_name) != names.end())
869 : // array variable
870 0 : return true;
871 : else
872 0 : return false;
873 : }
874 :
875 : bool
876 548989 : SystemBase::hasScalarVariable(const std::string & var_name) const
877 : {
878 548989 : if (system().has_variable(var_name))
879 296719 : return system().variable_type(var_name).family == SCALAR;
880 : else
881 252270 : return false;
882 : }
883 :
884 : bool
885 98106 : SystemBase::isScalarVariable(unsigned int var_num) const
886 : {
887 98106 : return (system().variable(var_num).type().family == SCALAR);
888 : }
889 :
890 : unsigned int
891 59564103 : SystemBase::nVariables() const
892 : {
893 59564103 : unsigned int n = nFieldVariables();
894 59564103 : n += _vars[0].scalars().size();
895 :
896 59564103 : return n;
897 : }
898 :
899 : unsigned int
900 59564103 : SystemBase::nFieldVariables() const
901 : {
902 59564103 : unsigned int n = 0;
903 136699864 : for (auto & var : _vars[0].fieldVariables())
904 77135761 : n += var->count();
905 :
906 59564103 : return n;
907 : }
908 :
909 : unsigned int
910 2539914 : SystemBase::nFVVariables() const
911 : {
912 2539914 : unsigned int n = 0;
913 5444926 : for (auto & var : _vars[0].fieldVariables())
914 2905012 : if (var->isFV())
915 1599006 : n += var->count();
916 :
917 2539914 : return n;
918 : }
919 :
920 : /**
921 : * Check if the named vector exists in the system.
922 : */
923 : bool
924 510351 : SystemBase::hasVector(const std::string & name) const
925 : {
926 510351 : return system().have_vector(name);
927 : }
928 :
929 : /**
930 : * Get a raw NumericVector with the given name.
931 : */
932 : NumericVector<Number> &
933 5737 : SystemBase::getVector(const std::string & name)
934 : {
935 5737 : return system().get_vector(name);
936 : }
937 :
938 : const NumericVector<Number> &
939 0 : SystemBase::getVector(const std::string & name) const
940 : {
941 0 : return system().get_vector(name);
942 : }
943 :
944 : NumericVector<Number> &
945 2775230930 : SystemBase::getVector(TagID tag)
946 : {
947 2775230930 : if (!hasVector(tag))
948 : {
949 0 : if (!_subproblem.vectorTagExists(tag))
950 0 : mooseError("Cannot retreive vector with tag ", tag, " because that tag does not exist");
951 : else
952 0 : mooseError("Cannot retreive vector with tag ",
953 : tag,
954 : " in system '",
955 0 : name(),
956 : "'\nbecause a vector has not been associated with that tag.");
957 : }
958 :
959 2775230930 : return *_tagged_vectors[tag];
960 : }
961 :
962 : const NumericVector<Number> &
963 688 : SystemBase::getVector(TagID tag) const
964 : {
965 688 : if (!hasVector(tag))
966 : {
967 0 : if (!_subproblem.vectorTagExists(tag))
968 0 : mooseError("Cannot retreive vector with tag ", tag, " because that tag does not exist");
969 : else
970 0 : mooseError("Cannot retreive vector with tag ",
971 : tag,
972 : " in system '",
973 0 : name(),
974 : "'\nbecause a vector has not been associated with that tag.");
975 : }
976 :
977 688 : return *_tagged_vectors[tag];
978 : }
979 :
980 : void
981 10806209 : SystemBase::associateVectorToTag(NumericVector<Number> & vec, TagID tag)
982 : {
983 10806209 : if (!_subproblem.vectorTagExists(tag))
984 0 : mooseError("Cannot associate vector to tag ", tag, " because that tag does not exist");
985 :
986 10806209 : if (_tagged_vectors.size() < tag + 1)
987 261055 : _tagged_vectors.resize(tag + 1);
988 :
989 10806209 : _tagged_vectors[tag] = &vec;
990 10806209 : }
991 :
992 : void
993 3321733 : SystemBase::disassociateVectorFromTag(NumericVector<Number> & vec, TagID tag)
994 : {
995 3321733 : if (!_subproblem.vectorTagExists(tag))
996 0 : mooseError("Cannot disassociate vector from tag ", tag, " because that tag does not exist");
997 3321733 : if (hasVector(tag) && &getVector(tag) != &vec)
998 0 : mooseError("You can not disassociate a vector from a tag which it was not associated to");
999 :
1000 3321733 : disassociateVectorFromTag(tag);
1001 3321733 : }
1002 :
1003 : void
1004 3449117 : SystemBase::disassociateVectorFromTag(TagID tag)
1005 : {
1006 3449117 : if (!_subproblem.vectorTagExists(tag))
1007 0 : mooseError("Cannot disassociate vector from tag ", tag, " because that tag does not exist");
1008 :
1009 3449117 : if (_tagged_vectors.size() < tag + 1)
1010 0 : _tagged_vectors.resize(tag + 1);
1011 3449117 : _tagged_vectors[tag] = nullptr;
1012 3449117 : }
1013 :
1014 : void
1015 31846 : SystemBase::disassociateDefaultVectorTags()
1016 : {
1017 31846 : const auto tags = defaultVectorTags();
1018 191076 : for (const auto tag : tags)
1019 159230 : if (_subproblem.vectorTagExists(tag))
1020 127384 : disassociateVectorFromTag(tag);
1021 31846 : }
1022 :
1023 : SparseMatrix<Number> &
1024 2019948483 : SystemBase::getMatrix(TagID tag)
1025 : {
1026 2019948483 : if (!hasMatrix(tag))
1027 : {
1028 0 : if (!_subproblem.matrixTagExists(tag))
1029 0 : mooseError("Cannot retreive matrix with tag ", tag, " because that tag does not exist");
1030 : else
1031 0 : mooseError("Cannot retreive matrix with tag ",
1032 : tag,
1033 : " in system '",
1034 0 : name(),
1035 : "'\nbecause a matrix has not been associated with that tag.");
1036 : }
1037 :
1038 2019948483 : return *_tagged_matrices[tag];
1039 : }
1040 :
1041 : const SparseMatrix<Number> &
1042 344 : SystemBase::getMatrix(TagID tag) const
1043 : {
1044 344 : if (!hasMatrix(tag))
1045 : {
1046 0 : if (!_subproblem.matrixTagExists(tag))
1047 0 : mooseError("Cannot retreive matrix with tag ", tag, " because that tag does not exist");
1048 : else
1049 0 : mooseError("Cannot retreive matrix with tag ",
1050 : tag,
1051 : " in system '",
1052 0 : name(),
1053 : "'\nbecause a matrix has not been associated with that tag.");
1054 : }
1055 :
1056 344 : return *_tagged_matrices[tag];
1057 : }
1058 :
1059 : void
1060 959257 : SystemBase::closeTaggedMatrices(const std::set<TagID> & tags)
1061 : {
1062 2857375 : for (auto tag : tags)
1063 1898118 : if (hasMatrix(tag))
1064 962537 : getMatrix(tag).close();
1065 959257 : }
1066 :
1067 : void
1068 0 : SystemBase::flushTaggedMatrices(const std::set<TagID> & tags)
1069 : {
1070 0 : for (auto tag : tags)
1071 0 : if (hasMatrix(tag))
1072 0 : getMatrix(tag).flush();
1073 0 : }
1074 :
1075 : void
1076 542541 : SystemBase::associateMatrixToTag(SparseMatrix<Number> & matrix, TagID tag)
1077 : {
1078 542541 : if (!_subproblem.matrixTagExists(tag))
1079 0 : mooseError("Cannot associate matrix to tag ", tag, " because that tag does not exist");
1080 :
1081 542541 : if (_tagged_matrices.size() < tag + 1)
1082 40026 : _tagged_matrices.resize(tag + 1);
1083 :
1084 542541 : _tagged_matrices[tag] = &matrix;
1085 542541 : }
1086 :
1087 : void
1088 536435 : SystemBase::disassociateMatrixFromTag(SparseMatrix<Number> & matrix, TagID tag)
1089 : {
1090 536435 : if (!_subproblem.matrixTagExists(tag))
1091 0 : mooseError("Cannot disassociate matrix from tag ", tag, " because that tag does not exist");
1092 536435 : if (hasMatrix(tag) && &getMatrix(tag) != &matrix)
1093 0 : mooseError("You can not disassociate a matrix from a tag which it was not associated to");
1094 :
1095 536435 : disassociateMatrixFromTag(tag);
1096 536435 : }
1097 :
1098 : void
1099 548735 : SystemBase::disassociateMatrixFromTag(TagID tag)
1100 : {
1101 548735 : if (!_subproblem.matrixTagExists(tag))
1102 0 : mooseError("Cannot disassociate matrix from tag ", tag, " because that tag does not exist");
1103 :
1104 548735 : if (_tagged_matrices.size() < tag + 1)
1105 1410 : _tagged_matrices.resize(tag + 1);
1106 548735 : _tagged_matrices[tag] = nullptr;
1107 548735 : }
1108 :
1109 : void
1110 4100 : SystemBase::disassociateDefaultMatrixTags()
1111 : {
1112 4100 : const auto tags = defaultMatrixTags();
1113 16400 : for (const auto tag : tags)
1114 12300 : if (_subproblem.matrixTagExists(tag))
1115 12300 : disassociateMatrixFromTag(tag);
1116 4100 : }
1117 :
1118 : void
1119 3607328 : SystemBase::deactivateAllMatrixTags()
1120 : {
1121 3607328 : auto num_matrix_tags = _subproblem.numMatrixTags();
1122 :
1123 3607328 : _matrix_tag_active_flags.resize(num_matrix_tags);
1124 :
1125 10929968 : for (decltype(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
1126 7322640 : _matrix_tag_active_flags[tag] = false;
1127 3607328 : _active_tagged_matrices.clear();
1128 3607328 : }
1129 :
1130 : void
1131 3822280 : SystemBase::activateAllMatrixTags()
1132 : {
1133 3822280 : auto num_matrix_tags = _subproblem.numMatrixTags();
1134 :
1135 3822280 : _matrix_tag_active_flags.resize(num_matrix_tags);
1136 3822280 : _active_tagged_matrices.clear();
1137 :
1138 11578851 : for (const auto tag : make_range(num_matrix_tags))
1139 7756571 : if (hasMatrix(tag))
1140 : {
1141 535698 : _matrix_tag_active_flags[tag] = true;
1142 535698 : _active_tagged_matrices.emplace(tag, &getMatrix(tag));
1143 : }
1144 : else
1145 7220873 : _matrix_tag_active_flags[tag] = false;
1146 3822280 : }
1147 :
1148 : bool
1149 61996 : SystemBase::matrixTagActive(TagID tag) const
1150 : {
1151 : mooseAssert(_subproblem.matrixTagExists(tag), "Matrix tag " << tag << " does not exist");
1152 :
1153 61996 : return tag < _matrix_tag_active_flags.size() && _matrix_tag_active_flags[tag];
1154 : }
1155 :
1156 : unsigned int
1157 1909696021 : SystemBase::number() const
1158 : {
1159 1909696021 : return system().number();
1160 : }
1161 :
1162 : DofMap &
1163 1967340 : SystemBase::dofMap()
1164 : {
1165 1967340 : return system().get_dof_map();
1166 : }
1167 :
1168 : const DofMap &
1169 0 : SystemBase::dofMap() const
1170 : {
1171 0 : return system().get_dof_map();
1172 : }
1173 :
1174 : void
1175 430 : SystemBase::addVariableToCopy(const std::string & dest_name,
1176 : const std::string & source_name,
1177 : const std::string & timestep)
1178 : {
1179 430 : _var_to_copy.push_back(VarCopyInfo(dest_name, source_name, timestep));
1180 430 : }
1181 :
1182 : void
1183 586 : SystemBase::copyVars(ExodusII_IO & io)
1184 : {
1185 586 : int n_steps = io.get_num_time_steps();
1186 :
1187 586 : bool did_copy = false;
1188 995 : for (const auto & vci : _var_to_copy)
1189 : {
1190 413 : int timestep = -1;
1191 :
1192 413 : if (vci._timestep == "LATEST")
1193 : // Use the last time step in the file from which to retrieve the solution
1194 320 : timestep = n_steps;
1195 : else
1196 : {
1197 93 : timestep = MooseUtils::convert<int>(vci._timestep);
1198 93 : if (timestep > n_steps)
1199 4 : mooseError("Invalid value passed as \"initial_from_file_timestep\". Expected \"LATEST\" or "
1200 : "a valid integer between 1 and ",
1201 : n_steps,
1202 : " inclusive, received ",
1203 4 : vci._timestep);
1204 : }
1205 :
1206 409 : did_copy = true;
1207 :
1208 409 : if (hasVariable(vci._dest_name))
1209 : {
1210 361 : const auto & var = getVariable(0, vci._dest_name);
1211 361 : if (var.isArray())
1212 : {
1213 24 : const auto & array_var = getFieldVariable<RealEigenVector>(0, vci._dest_name);
1214 72 : for (MooseIndex(var.count()) i = 0; i < var.count(); ++i)
1215 : {
1216 48 : const auto & exodus_var = var.arrayVariableComponent(i);
1217 48 : const auto & system_var = array_var.componentName(i);
1218 48 : if (var.isNodal())
1219 24 : io.copy_nodal_solution(system(), exodus_var, system_var, timestep);
1220 : else
1221 24 : io.copy_elemental_solution(system(), exodus_var, system_var, timestep);
1222 : }
1223 : }
1224 : else
1225 : {
1226 337 : if (var.isNodal())
1227 289 : io.copy_nodal_solution(system(), vci._dest_name, vci._source_name, timestep);
1228 : else
1229 48 : io.copy_elemental_solution(system(), vci._dest_name, vci._source_name, timestep);
1230 : }
1231 : }
1232 48 : else if (hasScalarVariable(vci._dest_name))
1233 288 : io.copy_scalar_solution(system(), {vci._dest_name}, {vci._source_name}, timestep);
1234 : else
1235 0 : mooseError("Unrecognized variable ", vci._dest_name, " in variables to copy.");
1236 : }
1237 :
1238 582 : if (did_copy)
1239 297 : solution().close();
1240 678 : }
1241 :
1242 : void
1243 4300710 : SystemBase::update()
1244 : {
1245 4300710 : system().update();
1246 4300710 : }
1247 :
1248 : void
1249 0 : SystemBase::solve()
1250 : {
1251 0 : system().solve();
1252 0 : }
1253 :
1254 : /**
1255 : * Copy current solution into old and older
1256 : */
1257 : void
1258 112411 : SystemBase::copySolutionsBackwards()
1259 : {
1260 112411 : system().update();
1261 112411 : copyOldSolutions();
1262 112411 : copyPreviousNonlinearSolutions();
1263 112411 : }
1264 :
1265 : /**
1266 : * Shifts the solutions backwards in nonlinear iteration history
1267 : */
1268 : void
1269 112411 : SystemBase::copyPreviousNonlinearSolutions()
1270 : {
1271 : const auto states =
1272 112411 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Nonlinear)].size();
1273 112411 : if (states > 1)
1274 400 : for (unsigned int i = states - 1; i > 0; --i)
1275 400 : solutionState(i, Moose::SolutionIterationType::Nonlinear) =
1276 200 : solutionState(i - 1, Moose::SolutionIterationType::Nonlinear);
1277 :
1278 112411 : if (solutionPreviousNewton())
1279 200 : *solutionPreviousNewton() = *currentSolution();
1280 112411 : }
1281 :
1282 : /**
1283 : * Shifts the solutions backwards in time
1284 : */
1285 : void
1286 684478 : SystemBase::copyOldSolutions()
1287 : {
1288 : // copy the solutions backward: current->old, old->older
1289 : const auto states =
1290 684478 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Time)].size();
1291 684478 : if (states > 1)
1292 1023938 : for (unsigned int i = states - 1; i > 0; --i)
1293 526256 : solutionState(i) = solutionState(i - 1);
1294 :
1295 684478 : if (solutionUDotOld())
1296 11496 : *solutionUDotOld() = *solutionUDot();
1297 684478 : if (solutionUDotDotOld())
1298 11496 : *solutionUDotDotOld() = *solutionUDotDot();
1299 684478 : }
1300 :
1301 : void
1302 665748 : SystemBase::copyPreviousFixedPointSolutions()
1303 : {
1304 : const auto n_states =
1305 665748 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::FixedPoint)]
1306 665748 : .size();
1307 665748 : if (n_states > 1)
1308 91958 : for (unsigned int i = n_states - 1; i > 0; --i)
1309 91958 : solutionState(i, Moose::SolutionIterationType::FixedPoint) =
1310 45979 : solutionState(i - 1, Moose::SolutionIterationType::FixedPoint);
1311 665748 : }
1312 :
1313 : /**
1314 : * Restore current solutions (call after your solve failed)
1315 : */
1316 : void
1317 10652 : SystemBase::restoreSolutions()
1318 : {
1319 10652 : if (!hasSolutionState(1))
1320 0 : mooseError("Cannot restore solutions without old solution");
1321 :
1322 10652 : *(const_cast<NumericVector<Number> *&>(currentSolution())) = solutionOld();
1323 10652 : solution() = solutionOld();
1324 10652 : if (solutionUDotOld())
1325 22 : *solutionUDot() = *solutionUDotOld();
1326 10652 : if (solutionUDotDotOld())
1327 22 : *solutionUDotDot() = *solutionUDotDotOld();
1328 10652 : if (solutionPreviousNewton())
1329 2 : *solutionPreviousNewton() = solutionOld();
1330 10652 : system().update();
1331 10652 : }
1332 :
1333 : void
1334 460 : SystemBase::removeVector(const std::string & name)
1335 : {
1336 460 : system().remove_vector(name);
1337 460 : }
1338 :
1339 : const std::string &
1340 458133 : SystemBase::name() const
1341 : {
1342 458133 : return system().name();
1343 : }
1344 :
1345 : NumericVector<Number> *
1346 124377 : SystemBase::solutionPreviousNewton()
1347 : {
1348 124377 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
1349 1516 : return &getVector(Moose::PREVIOUS_NL_SOLUTION_TAG);
1350 : else
1351 122861 : return nullptr;
1352 : }
1353 :
1354 : const NumericVector<Number> *
1355 0 : SystemBase::solutionPreviousNewton() const
1356 : {
1357 0 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
1358 0 : return &getVector(Moose::PREVIOUS_NL_SOLUTION_TAG);
1359 : else
1360 0 : return nullptr;
1361 : }
1362 :
1363 : void
1364 126267 : SystemBase::initSolutionState()
1365 : {
1366 : // Default is the current solution
1367 126267 : unsigned int state = 0;
1368 :
1369 : // Add additional states as required by the variable states requested
1370 292288 : for (const auto & var : getVariables(/* tid = */ 0))
1371 166021 : state = std::max(state, var->oldestSolutionStateRequested());
1372 129467 : for (const auto & var : getScalarVariables(/* tid = */ 0))
1373 3200 : state = std::max(state, var->oldestSolutionStateRequested());
1374 :
1375 126267 : needSolutionState(state, Moose::SolutionIterationType::Time);
1376 :
1377 126267 : _solution_states_initialized = true;
1378 126267 : }
1379 :
1380 : TagName
1381 65274 : SystemBase::oldSolutionStateVectorName(const unsigned int state,
1382 : const Moose::SolutionIterationType iteration_type) const
1383 : {
1384 : mooseAssert(state != 0, "Not an old state");
1385 :
1386 65274 : if (iteration_type == Moose::SolutionIterationType::Time)
1387 : {
1388 64222 : if (state == 1)
1389 62736 : return Moose::OLD_SOLUTION_TAG;
1390 1486 : else if (state == 2)
1391 1486 : return Moose::OLDER_SOLUTION_TAG;
1392 : }
1393 1052 : else if (iteration_type == Moose::SolutionIterationType::Nonlinear && state == 1)
1394 214 : return Moose::PREVIOUS_NL_SOLUTION_TAG;
1395 838 : else if (iteration_type == Moose::SolutionIterationType::FixedPoint && state == 1)
1396 838 : return Moose::PREVIOUS_FP_SOLUTION_TAG;
1397 :
1398 0 : return "solution_state_" + std::to_string(state) + "_" + Moose::stringify(iteration_type);
1399 : }
1400 :
1401 : const NumericVector<Number> &
1402 166820198 : SystemBase::solutionState(const unsigned int state,
1403 : const Moose::SolutionIterationType iteration_type) const
1404 : {
1405 166820198 : if (!hasSolutionState(state, iteration_type))
1406 0 : mooseError("For iteration type '",
1407 0 : Moose::stringify(iteration_type),
1408 : "': solution state ",
1409 : state,
1410 : " was requested in ",
1411 0 : name(),
1412 : " but only up to state ",
1413 0 : (_solution_states[static_cast<unsigned short>(iteration_type)].size() == 0)
1414 0 : ? 0
1415 0 : : _solution_states[static_cast<unsigned short>(iteration_type)].size() - 1,
1416 : " is available.");
1417 :
1418 166820198 : const auto & solution_states = _solution_states[static_cast<unsigned short>(iteration_type)];
1419 :
1420 166820198 : if (state == 0)
1421 : mooseAssert(solution_states[0] == &solutionInternal(), "Inconsistent current solution");
1422 : else
1423 : mooseAssert(solution_states[state] ==
1424 : &getVector(oldSolutionStateVectorName(state, iteration_type)),
1425 : "Inconsistent solution state");
1426 :
1427 166820198 : return *solution_states[state];
1428 : }
1429 :
1430 : NumericVector<Number> &
1431 62195316 : SystemBase::solutionState(const unsigned int state,
1432 : const Moose::SolutionIterationType iteration_type)
1433 : {
1434 62195316 : if (!hasSolutionState(state, iteration_type))
1435 71142 : needSolutionState(state, iteration_type);
1436 62195316 : return *_solution_states[static_cast<unsigned short>(iteration_type)][state];
1437 : }
1438 :
1439 : libMesh::ParallelType
1440 1254 : SystemBase::solutionStateParallelType(const unsigned int state,
1441 : const Moose::SolutionIterationType iteration_type) const
1442 : {
1443 1254 : if (!hasSolutionState(state, iteration_type))
1444 0 : mooseError("solutionStateParallelType() may only be called if the solution state exists.");
1445 :
1446 1254 : return _solution_states[static_cast<unsigned short>(iteration_type)][state]->type();
1447 : }
1448 :
1449 : void
1450 215887 : SystemBase::needSolutionState(const unsigned int state,
1451 : const Moose::SolutionIterationType iteration_type,
1452 : const libMesh::ParallelType parallel_type)
1453 : {
1454 : libmesh_parallel_only(this->comm());
1455 : mooseAssert(!Threads::in_threads,
1456 : "This routine is not thread-safe. Request the solution state before using it in "
1457 : "a threaded region.");
1458 :
1459 215887 : if (hasSolutionState(state, iteration_type))
1460 90018 : return;
1461 :
1462 125869 : auto & solution_states = _solution_states[static_cast<unsigned short>(iteration_type)];
1463 125869 : solution_states.resize(state + 1);
1464 :
1465 : // The 0-th (current) solution state is owned by libMesh
1466 125869 : if (!solution_states[0])
1467 124527 : solution_states[0] = &solutionInternal();
1468 : else
1469 : mooseAssert(solution_states[0] == &solutionInternal(), "Inconsistent current solution");
1470 :
1471 : // We will manually add all states past current
1472 192397 : for (unsigned int i = 1; i <= state; ++i)
1473 66528 : if (!solution_states[i])
1474 : {
1475 65274 : auto tag = _subproblem.addVectorTag(oldSolutionStateVectorName(i, iteration_type),
1476 : Moose::VECTOR_TAG_SOLUTION);
1477 65274 : solution_states[i] = &addVector(tag, true, parallel_type);
1478 : }
1479 : else
1480 : {
1481 : // If the existing parallel type is PARALLEL and GHOSTED is now requested,
1482 : // this would require an upgrade, which is risky if anybody has already
1483 : // stored a pointer to the existing vector, since the upgrade would create
1484 : // a new vector and make that pointer null. If the existing parallel type
1485 : // is GHOSTED and PARALLEL is now requested, we don't need to do anything.
1486 1254 : if (parallel_type == GHOSTED && solutionStateParallelType(i, iteration_type) == PARALLEL)
1487 0 : mooseError("The solution state has already been declared as PARALLEL");
1488 :
1489 : mooseAssert(solution_states[i] == &getVector(oldSolutionStateVectorName(i, iteration_type)),
1490 : "Inconsistent solution state");
1491 : }
1492 : }
1493 :
1494 : void
1495 715 : SystemBase::applyScalingFactors(const std::vector<Real> & inverse_scaling_factors)
1496 : {
1497 1500 : for (MooseIndex(_vars) thread = 0; thread < _vars.size(); ++thread)
1498 : {
1499 785 : auto & field_variables = _vars[thread].fieldVariables();
1500 1977 : for (MooseIndex(field_variables) i = 0, p = 0; i < field_variables.size(); ++i)
1501 : {
1502 1192 : auto factors = field_variables[i]->arrayScalingFactor();
1503 2395 : for (unsigned int j = 0; j < field_variables[i]->count(); ++j, ++p)
1504 1203 : factors[j] /= inverse_scaling_factors[p];
1505 :
1506 1192 : field_variables[i]->scalingFactor(factors);
1507 1192 : }
1508 :
1509 785 : auto offset = field_variables.size();
1510 :
1511 785 : auto & scalar_variables = _vars[thread].scalars();
1512 909 : for (MooseIndex(scalar_variables) i = 0; i < scalar_variables.size(); ++i)
1513 372 : scalar_variables[i]->scalingFactor(
1514 124 : {1. / inverse_scaling_factors[offset + i] * scalar_variables[i]->scalingFactor()});
1515 :
1516 785 : if (thread == 0 && _verbose)
1517 : {
1518 286 : _console << "Automatic scaling factors:\n";
1519 286 : auto original_flags = _console.flags();
1520 286 : auto original_precision = _console.precision();
1521 286 : _console.unsetf(std::ios_base::floatfield);
1522 286 : _console.precision(6);
1523 :
1524 738 : for (const auto & field_variable : field_variables)
1525 : {
1526 452 : const auto & factors = field_variable->arrayScalingFactor();
1527 452 : _console << " " << field_variable->name() << ":";
1528 914 : for (const auto i : make_range(field_variable->count()))
1529 462 : _console << " " << factors[i];
1530 452 : _console << "\n";
1531 : }
1532 382 : for (const auto & scalar_variable : scalar_variables)
1533 192 : _console << " " << scalar_variable->name() << ": " << scalar_variable->scalingFactor()
1534 96 : << "\n";
1535 286 : _console << "\n" << std::endl;
1536 :
1537 : // restore state
1538 286 : _console.flags(original_flags);
1539 286 : _console.precision(original_precision);
1540 : }
1541 : }
1542 715 : }
1543 :
1544 : void
1545 363 : SystemBase::addScalingVector()
1546 : {
1547 363 : addVector("scaling_factors", /*project=*/false, libMesh::ParallelType::GHOSTED);
1548 363 : _subproblem.hasScalingVector(number());
1549 363 : }
1550 :
1551 : bool
1552 78662077 : SystemBase::computingScalingJacobian() const
1553 : {
1554 78662077 : return _subproblem.computingScalingJacobian();
1555 : }
1556 :
1557 : void
1558 125633 : SystemBase::initialSetup()
1559 : {
1560 262781 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1561 137148 : _vars[tid].initialSetup();
1562 :
1563 : // If we need raw gradients, we initialize them here.
1564 125633 : bool gradient_storage_initialized = false;
1565 291209 : for (const auto & field_var : _vars[0].fieldVariables())
1566 165576 : if (!gradient_storage_initialized && field_var->needsGradientVectorStorage())
1567 : {
1568 979 : _raw_grad_container.clear();
1569 2542 : for (const auto i : make_range(this->_mesh.dimension()))
1570 : {
1571 1563 : libmesh_ignore(i);
1572 1563 : _raw_grad_container.push_back(currentSolution()->zero_clone());
1573 : }
1574 : }
1575 125633 : }
1576 :
1577 : void
1578 641225 : SystemBase::timestepSetup()
1579 : {
1580 1339155 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1581 697930 : _vars[tid].timestepSetup();
1582 641225 : }
1583 :
1584 : void
1585 4134900 : SystemBase::customSetup(const ExecFlagType & exec_type)
1586 : {
1587 8639522 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1588 4504622 : _vars[tid].customSetup(exec_type);
1589 4134900 : }
1590 :
1591 : void
1592 6188103 : SystemBase::subdomainSetup()
1593 : {
1594 13309515 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1595 7121412 : _vars[tid].subdomainSetup();
1596 6188103 : }
1597 :
1598 : void
1599 6961449 : SystemBase::residualSetup()
1600 : {
1601 14547321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1602 7585872 : _vars[tid].residualSetup();
1603 6961449 : }
1604 :
1605 : void
1606 1107868 : SystemBase::jacobianSetup()
1607 : {
1608 2317352 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1609 1209484 : _vars[tid].jacobianSetup();
1610 1107868 : }
1611 :
1612 : void
1613 772096 : SystemBase::clearAllDofIndices()
1614 : {
1615 1611948 : for (auto & var_warehouse : _vars)
1616 839852 : var_warehouse.clearAllDofIndices();
1617 772096 : }
1618 :
1619 : void
1620 15411218 : SystemBase::setActiveVariableCoupleableVectorTags(const std::set<TagID> & vtags, THREAD_ID tid)
1621 : {
1622 15411218 : _vars[tid].setActiveVariableCoupleableVectorTags(vtags);
1623 15411218 : }
1624 :
1625 : void
1626 112062 : SystemBase::setActiveScalarVariableCoupleableVectorTags(const std::set<TagID> & vtags,
1627 : THREAD_ID tid)
1628 : {
1629 112062 : _vars[tid].setActiveScalarVariableCoupleableVectorTags(vtags);
1630 112062 : }
1631 :
1632 : void
1633 62506 : SystemBase::addDotVectors()
1634 : {
1635 62506 : if (_fe_problem.uDotRequested())
1636 125012 : _u_dot = &addVector("u_dot", true, GHOSTED);
1637 62506 : if (_fe_problem.uDotOldRequested())
1638 656 : _u_dot_old = &addVector("u_dot_old", true, GHOSTED);
1639 62506 : if (_fe_problem.uDotDotRequested())
1640 656 : _u_dotdot = &addVector("u_dotdot", true, GHOSTED);
1641 62506 : if (_fe_problem.uDotDotOldRequested())
1642 656 : _u_dotdot_old = &addVector("u_dotdot_old", true, GHOSTED);
1643 62506 : }
1644 :
1645 : NumericVector<Number> &
1646 1140 : SystemBase::serializedSolution()
1647 : {
1648 1140 : if (!_serialized_solution.get())
1649 : {
1650 113 : _serialized_solution = NumericVector<Number>::build(_communicator);
1651 113 : _serialized_solution->init(system().n_dofs(), false, SERIAL);
1652 : }
1653 :
1654 1140 : return *_serialized_solution;
1655 : }
1656 :
1657 : void
1658 62532 : SystemBase::addTimeIntegrator(const std::string & type,
1659 : const std::string & name,
1660 : InputParameters & parameters)
1661 : {
1662 125064 : parameters.set<SystemBase *>("_sys") = this;
1663 62532 : _time_integrators.push_back(_factory.create<TimeIntegrator>(type, name, parameters));
1664 62532 : }
1665 :
1666 : void
1667 3816 : SystemBase::copyTimeIntegrators(const SystemBase & other_sys)
1668 : {
1669 3816 : _time_integrators = other_sys._time_integrators;
1670 3816 : }
1671 :
1672 : const TimeIntegrator *
1673 739061 : SystemBase::queryTimeIntegrator(const unsigned int var_num) const
1674 : {
1675 740021 : for (auto & ti : _time_integrators)
1676 471242 : if (ti->integratesVar(var_num))
1677 470282 : return ti.get();
1678 :
1679 268779 : return nullptr;
1680 : }
1681 :
1682 : const TimeIntegrator &
1683 92 : SystemBase::getTimeIntegrator(const unsigned int var_num) const
1684 : {
1685 92 : const auto * const ti = queryTimeIntegrator(var_num);
1686 :
1687 92 : if (ti)
1688 92 : return *ti;
1689 : else
1690 0 : mooseError("No time integrator found that integrates variable number ",
1691 0 : std::to_string(var_num));
1692 : }
1693 :
1694 : const std::vector<std::shared_ptr<TimeIntegrator>> &
1695 587636 : SystemBase::getTimeIntegrators()
1696 : {
1697 587636 : return _time_integrators;
1698 : }
1699 :
1700 : const Number &
1701 1275643343 : SystemBase::duDotDu(const unsigned int var_num) const
1702 : {
1703 1275643343 : return _du_dot_du[var_num];
1704 : }
1705 :
1706 : const std::set<SubdomainID> &
1707 363 : SystemBase::getSubdomainsForVar(const std::string & var_name) const
1708 : {
1709 363 : return getSubdomainsForVar(getVariable(0, var_name).number());
1710 : }
1711 :
1712 : std::string
1713 77655 : SystemBase::prefix() const
1714 : {
1715 231181 : return system().prefix_with_name() ? system().prefix() : "";
1716 : }
1717 :
1718 : void
1719 125289 : SystemBase::sizeVariableMatrixData()
1720 : {
1721 262713 : for (const auto & warehouse : _vars)
1722 275818 : for (const auto & [var_num, var_ptr] : warehouse.numberToVariableMap())
1723 138394 : var_ptr->sizeMatrixTagData();
1724 125289 : }
1725 :
1726 : template MooseVariableFE<Real> & SystemBase::getFieldVariable<Real>(THREAD_ID tid,
1727 : const std::string & var_name);
1728 :
1729 : template MooseVariableFE<RealVectorValue> &
1730 : SystemBase::getFieldVariable<RealVectorValue>(THREAD_ID tid, const std::string & var_name);
1731 :
1732 : template MooseVariableFE<RealEigenVector> &
1733 : SystemBase::getFieldVariable<RealEigenVector>(THREAD_ID tid, const std::string & var_name);
1734 :
1735 : template MooseVariableFE<Real> & SystemBase::getFieldVariable<Real>(THREAD_ID tid,
1736 : unsigned int var_number);
1737 :
1738 : template MooseVariableFE<RealVectorValue> &
1739 : SystemBase::getFieldVariable<RealVectorValue>(THREAD_ID tid, unsigned int var_number);
1740 :
1741 : template MooseVariableFE<RealEigenVector> &
1742 : SystemBase::getFieldVariable<RealEigenVector>(THREAD_ID tid, unsigned int var_number);
1743 :
1744 : template MooseVariableField<Real> &
1745 : SystemBase::getActualFieldVariable<Real>(THREAD_ID tid, const std::string & var_name);
1746 :
1747 : template MooseVariableField<RealVectorValue> &
1748 : SystemBase::getActualFieldVariable<RealVectorValue>(THREAD_ID tid, const std::string & var_name);
1749 :
1750 : template MooseVariableField<RealEigenVector> &
1751 : SystemBase::getActualFieldVariable<RealEigenVector>(THREAD_ID tid, const std::string & var_name);
1752 :
1753 : template MooseVariableField<Real> &
1754 : SystemBase::getActualFieldVariable<Real>(THREAD_ID tid, unsigned int var_number);
1755 :
1756 : template MooseVariableField<RealVectorValue> &
1757 : SystemBase::getActualFieldVariable<RealVectorValue>(THREAD_ID tid, unsigned int var_number);
1758 :
1759 : template MooseVariableField<RealEigenVector> &
1760 : SystemBase::getActualFieldVariable<RealEigenVector>(THREAD_ID tid, unsigned int var_number);
1761 :
1762 : template MooseVariableFV<Real> & SystemBase::getFVVariable<Real>(THREAD_ID tid,
1763 : const std::string & var_name);
|
