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