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 62554 : extraSendList(std::vector<dof_id_type> & send_list, void * context)
40 : {
41 62554 : SystemBase * sys = static_cast<SystemBase *>(context);
42 62554 : sys->augmentSendList(send_list);
43 62554 : }
44 :
45 : /// Free function used for a libMesh callback
46 : void
47 61920 : 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 61920 : SystemBase * sys = static_cast<SystemBase *>(context);
53 61920 : sys->augmentSparsity(sparsity, n_nz, n_oz);
54 61920 : }
55 :
56 120366 : SystemBase::SystemBase(SubProblem & subproblem,
57 : FEProblemBase & fe_problem,
58 : const std::string & name,
59 120366 : Moose::VarKindType var_kind)
60 : : libMesh::ParallelObject(subproblem),
61 : ConsoleStreamInterface(subproblem.getMooseApp()),
62 120366 : _subproblem(subproblem),
63 120366 : _fe_problem(fe_problem),
64 240732 : _app(subproblem.getMooseApp()),
65 120366 : _factory(_app.getFactory()),
66 120366 : _mesh(subproblem.mesh()),
67 120366 : _name(name),
68 120366 : _vars(libMesh::n_threads()),
69 120366 : _var_map(),
70 120366 : _max_var_number(0),
71 120366 : _u_dot(nullptr),
72 120366 : _u_dotdot(nullptr),
73 120366 : _u_dot_old(nullptr),
74 120366 : _u_dotdot_old(nullptr),
75 120366 : _saved_old(nullptr),
76 120366 : _saved_older(nullptr),
77 120366 : _saved_dot_old(nullptr),
78 120366 : _saved_dotdot_old(nullptr),
79 120366 : _var_kind(var_kind),
80 120366 : _max_var_n_dofs_per_elem(0),
81 120366 : _max_var_n_dofs_per_node(0),
82 120366 : _automatic_scaling(false),
83 120366 : _verbose(false),
84 240732 : _solution_states_initialized(false)
85 : {
86 120366 : }
87 :
88 : MooseVariableFieldBase &
89 4441737 : SystemBase::getVariable(THREAD_ID tid, const std::string & var_name) const
90 : {
91 : MooseVariableFieldBase * var =
92 4441737 : dynamic_cast<MooseVariableFieldBase *>(_vars[tid].getVariable(var_name));
93 4441737 : if (!var)
94 4 : mooseError("Variable '", var_name, "' does not exist in this system");
95 4441733 : return *var;
96 : }
97 :
98 : MooseVariableFieldBase &
99 152247346 : SystemBase::getVariable(THREAD_ID tid, unsigned int var_number) const
100 : {
101 152247346 : if (var_number < _numbered_vars[tid].size())
102 152247346 : if (_numbered_vars[tid][var_number])
103 152247346 : 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 53913 : SystemBase::getFieldVariable(THREAD_ID tid, const std::string & var_name)
111 : {
112 53913 : return *_vars[tid].getFieldVariable<T>(var_name);
113 : }
114 :
115 : template <typename T>
116 : MooseVariableField<T> &
117 102839 : SystemBase::getActualFieldVariable(THREAD_ID tid, const std::string & var_name)
118 : {
119 102839 : return *_vars[tid].getActualFieldVariable<T>(var_name);
120 : }
121 :
122 : template <typename T>
123 : MooseVariableFV<T> &
124 662 : SystemBase::getFVVariable(THREAD_ID tid, const std::string & var_name)
125 : {
126 662 : 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 131740896 : SystemBase::getActualFieldVariable(THREAD_ID tid, unsigned int var_number)
139 : {
140 131740896 : return *_vars[tid].getActualFieldVariable<T>(var_number);
141 : }
142 :
143 : MooseVariableScalar &
144 57895 : SystemBase::getScalarVariable(THREAD_ID tid, const std::string & var_name) const
145 : {
146 57895 : MooseVariableScalar * var = dynamic_cast<MooseVariableScalar *>(_vars[tid].getVariable(var_name));
147 57895 : if (!var)
148 8 : mooseError("Scalar variable '" + var_name + "' does not exist in this system");
149 57887 : return *var;
150 : }
151 :
152 : MooseVariableScalar &
153 185323 : SystemBase::getScalarVariable(THREAD_ID tid, unsigned int var_number) const
154 : {
155 : MooseVariableScalar * var =
156 185323 : dynamic_cast<MooseVariableScalar *>(_vars[tid].getVariable(var_number));
157 185323 : if (!var)
158 0 : mooseError("variable #" + Moose::stringify(var_number) + " does not exist in this system");
159 185323 : return *var;
160 : }
161 :
162 : const std::set<SubdomainID> *
163 328 : SystemBase::getVariableBlocks(unsigned int var_number)
164 : {
165 : mooseAssert(_var_map.find(var_number) != _var_map.end(), "Variable does not exist.");
166 328 : if (_var_map[var_number].empty())
167 328 : return nullptr;
168 : else
169 0 : return &_var_map[var_number];
170 : }
171 :
172 : void
173 475 : SystemBase::addVariableToZeroOnResidual(std::string var_name)
174 : {
175 475 : _vars_to_be_zeroed_on_residual.push_back(var_name);
176 475 : }
177 :
178 : void
179 316 : SystemBase::addVariableToZeroOnJacobian(std::string var_name)
180 : {
181 316 : _vars_to_be_zeroed_on_jacobian.push_back(var_name);
182 316 : }
183 :
184 : void
185 112 : SystemBase::setVariableGlobalDoFs(const std::string & var_name)
186 : {
187 224 : AllLocalDofIndicesThread aldit(_subproblem, {var_name});
188 112 : ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
189 112 : Threads::parallel_reduce(elem_range, aldit);
190 :
191 : // Gather the dof indices across procs to get all the dof indices for var_name
192 112 : aldit.dofIndicesSetUnion();
193 :
194 112 : const auto & all_dof_indices = aldit.getDofIndices();
195 112 : _var_all_dof_indices.assign(all_dof_indices.begin(), all_dof_indices.end());
196 224 : }
197 :
198 : void
199 3529535 : SystemBase::zeroVariables(std::vector<std::string> & vars_to_be_zeroed)
200 : {
201 3529535 : if (vars_to_be_zeroed.size() > 0)
202 : {
203 9265 : NumericVector<Number> & solution = this->solution();
204 :
205 9265 : auto problem = dynamic_cast<FEProblemBase *>(&_subproblem);
206 9265 : if (!problem)
207 0 : mooseError("System needs to be registered in FEProblemBase for using zeroVariables.");
208 :
209 9265 : AllLocalDofIndicesThread aldit(*problem, vars_to_be_zeroed, true);
210 9265 : ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
211 9265 : Threads::parallel_reduce(elem_range, aldit);
212 :
213 9265 : const auto & dof_indices_to_zero = aldit.getDofIndices();
214 :
215 9265 : solution.close();
216 :
217 2794663 : for (const auto & dof : dof_indices_to_zero)
218 2785398 : solution.set(dof, 0);
219 :
220 9265 : solution.close();
221 :
222 : // Call update to update the current_local_solution for this system
223 9265 : system().update();
224 9265 : }
225 3529535 : }
226 :
227 : void
228 3043777 : SystemBase::zeroVariablesForResidual()
229 : {
230 3043777 : zeroVariables(_vars_to_be_zeroed_on_residual);
231 3043777 : }
232 :
233 : void
234 476969 : SystemBase::zeroVariablesForJacobian()
235 : {
236 476969 : zeroVariables(_vars_to_be_zeroed_on_jacobian);
237 476969 : }
238 :
239 : Order
240 57397 : SystemBase::getMinQuadratureOrder()
241 : {
242 57397 : Order order = CONSTANT;
243 57397 : const std::vector<MooseVariableFieldBase *> & vars = _vars[0].fieldVariables();
244 113899 : for (const auto & var : vars)
245 : {
246 56502 : FEType fe_type = var->feType();
247 56502 : if (fe_type.default_quadrature_order() > order)
248 46362 : order = fe_type.default_quadrature_order();
249 : }
250 :
251 57397 : return order;
252 : }
253 :
254 : void
255 786744495 : SystemBase::prepare(THREAD_ID tid)
256 : {
257 786744495 : if (_subproblem.hasActiveElementalMooseVariables(tid))
258 : {
259 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
260 756185527 : _subproblem.getActiveElementalMooseVariables(tid);
261 756185527 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
262 1470565411 : for (const auto & var : vars)
263 714379884 : var->clearDofIndices();
264 :
265 1702680110 : for (const auto & var : active_elemental_moose_variables)
266 946494583 : if (&(var->sys()) == this)
267 462056741 : var->prepare();
268 : }
269 : else
270 : {
271 30558968 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
272 63331450 : for (const auto & var : vars)
273 32772482 : var->prepare();
274 : }
275 786744495 : }
276 :
277 : void
278 53182 : SystemBase::prepareFace(THREAD_ID tid, bool resize_data)
279 : {
280 : // We only need to do something if the element prepare was restricted
281 53182 : if (_subproblem.hasActiveElementalMooseVariables(tid))
282 : {
283 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
284 33042 : _subproblem.getActiveElementalMooseVariables(tid);
285 :
286 33042 : std::vector<MooseVariableFieldBase *> newly_prepared_vars;
287 :
288 33042 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
289 61097 : 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 28055 : if (!active_elemental_moose_variables.count(var))
305 : {
306 3313 : var->prepare();
307 3313 : 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 33042 : if (resize_data)
313 18009 : for (const auto var_ptr : newly_prepared_vars)
314 : {
315 1488 : _subproblem.assembly(tid, number()).prepareVariable(var_ptr);
316 1488 : if (_subproblem.checkNonlocalCouplingRequirement())
317 0 : _subproblem.assembly(tid, number()).prepareVariableNonlocal(var_ptr);
318 : }
319 33042 : }
320 53182 : }
321 :
322 : void
323 8655063 : SystemBase::prepareNeighbor(THREAD_ID tid)
324 : {
325 8655063 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
326 22118861 : for (const auto & var : vars)
327 13463798 : var->prepareNeighbor();
328 8655063 : }
329 :
330 : void
331 599590 : SystemBase::prepareLowerD(THREAD_ID tid)
332 : {
333 599590 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
334 1188469 : for (const auto & var : vars)
335 588879 : var->prepareLowerD();
336 599590 : }
337 :
338 : void
339 395446762 : SystemBase::reinitElem(const Elem * /*elem*/, THREAD_ID tid)
340 : {
341 395446762 : if (_subproblem.hasActiveElementalMooseVariables(tid))
342 : {
343 : const std::set<MooseVariableFieldBase *> & active_elemental_moose_variables =
344 380384185 : _subproblem.getActiveElementalMooseVariables(tid);
345 857808427 : for (const auto & var : active_elemental_moose_variables)
346 477424242 : if (&(var->sys()) == this)
347 391659611 : var->computeElemValues();
348 : }
349 : else
350 : {
351 15062577 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
352 29869107 : for (const auto & var : vars)
353 14806530 : var->computeElemValues();
354 : }
355 395446762 : }
356 :
357 : void
358 12425634 : SystemBase::reinitElemFace(const Elem * /*elem*/, unsigned int /*side*/, THREAD_ID tid)
359 : {
360 12425634 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
361 33918349 : for (const auto & var : vars)
362 21492715 : var->computeElemValuesFace();
363 12425634 : }
364 :
365 : void
366 8615303 : SystemBase::reinitNeighborFace(const Elem * /*elem*/, unsigned int /*side*/, THREAD_ID tid)
367 : {
368 8615303 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
369 22049981 : for (const auto & var : vars)
370 13434678 : var->computeNeighborValuesFace();
371 8615303 : }
372 :
373 : void
374 39760 : SystemBase::reinitNeighbor(const Elem * /*elem*/, THREAD_ID tid)
375 : {
376 39760 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
377 68880 : for (const auto & var : vars)
378 29120 : var->computeNeighborValues();
379 39760 : }
380 :
381 : void
382 599590 : SystemBase::reinitLowerD(THREAD_ID tid)
383 : {
384 599590 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
385 1188469 : for (const auto & var : vars)
386 588879 : var->computeLowerDValues();
387 599590 : }
388 :
389 : void
390 53498192 : SystemBase::reinitNode(const Node * /*node*/, THREAD_ID tid)
391 : {
392 53498192 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
393 145953116 : for (const auto & var : vars)
394 : {
395 92454924 : var->reinitNode();
396 92454924 : if (var->isNodalDefined())
397 78143188 : var->computeNodalValues();
398 : }
399 53498192 : }
400 :
401 : void
402 151195034 : SystemBase::reinitNodeFace(const Node * /*node*/, BoundaryID /*bnd_id*/, THREAD_ID tid)
403 : {
404 151195034 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
405 367116012 : for (const auto & var : vars)
406 : {
407 215920978 : var->reinitNode();
408 215920978 : if (var->isNodalDefined())
409 192645179 : var->computeNodalValues();
410 : }
411 151195034 : }
412 :
413 : void
414 12090 : SystemBase::reinitNodes(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
415 : {
416 12090 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
417 18280 : for (const auto & var : vars)
418 : {
419 6190 : var->reinitNodes(nodes);
420 6190 : var->computeNodalValues();
421 : }
422 12090 : }
423 :
424 : void
425 1310 : SystemBase::reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
426 : {
427 1310 : const std::vector<MooseVariableFieldBase *> & vars = _vars[tid].fieldVariables();
428 2110 : for (const auto & var : vars)
429 : {
430 800 : var->reinitNodesNeighbor(nodes);
431 800 : var->computeNodalNeighborValues();
432 : }
433 1310 : }
434 :
435 : void
436 17201914 : SystemBase::reinitScalars(THREAD_ID tid, bool reinit_for_derivative_reordering /*=false*/)
437 : {
438 17201914 : const std::vector<MooseVariableScalar *> & vars = _vars[tid].scalars();
439 17679321 : for (const auto & var : vars)
440 477407 : var->reinit(reinit_for_derivative_reordering);
441 17201914 : }
442 :
443 : void
444 62554 : SystemBase::augmentSendList(std::vector<dof_id_type> & send_list)
445 : {
446 62554 : std::set<dof_id_type> & ghosted_elems = _subproblem.ghostedElems();
447 :
448 62554 : DofMap & dof_map = dofMap();
449 :
450 62554 : std::vector<dof_id_type> dof_indices;
451 :
452 62554 : System & sys = system();
453 :
454 62554 : unsigned int sys_num = sys.number();
455 :
456 62554 : unsigned int n_vars = sys.n_vars();
457 :
458 153239 : for (const auto & elem_id : ghosted_elems)
459 : {
460 90685 : Elem * elem = _mesh.elemPtr(elem_id);
461 :
462 90685 : if (elem->active())
463 : {
464 90685 : dof_map.dof_indices(elem, dof_indices);
465 :
466 : // Only need to ghost it if it's actually not on this processor
467 1187889 : for (const auto & dof : dof_indices)
468 1097204 : if (dof < dof_map.first_dof() || dof >= dof_map.end_dof())
469 1073062 : 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 508027 : for (unsigned int n = 0; n < elem->n_nodes(); n++)
475 : {
476 417342 : Node * node = elem->node_ptr(n);
477 :
478 : // Have to get each variable's dofs
479 1628420 : for (unsigned int v = 0; v < n_vars; v++)
480 : {
481 1211078 : const Variable & var = sys.variable(v);
482 1211078 : unsigned int var_num = var.number();
483 1211078 : unsigned int n_comp = var.n_components();
484 :
485 : // See if this variable has any dofs at this node
486 1211078 : if (node->n_dofs(sys_num, var_num) > 0)
487 : {
488 : // Loop over components of the variable
489 2127864 : for (unsigned int c = 0; c < n_comp; c++)
490 1063932 : send_list.push_back(node->dof_number(sys_num, var_num, c));
491 : }
492 : }
493 : }
494 : }
495 : }
496 62554 : }
497 :
498 : /**
499 : * Save the old and older solutions.
500 : */
501 : void
502 230 : SystemBase::saveOldSolutions()
503 : {
504 : const auto states =
505 230 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Time)].size();
506 230 : if (states > 1)
507 : {
508 210 : _saved_solution_states.resize(states);
509 630 : for (unsigned int i = 1; i <= states - 1; ++i)
510 420 : if (!_saved_solution_states[i])
511 840 : _saved_solution_states[i] =
512 420 : &addVector("save_solution_state_" + std::to_string(i), false, PARALLEL);
513 :
514 630 : for (unsigned int i = 1; i <= states - 1; ++i)
515 420 : *(_saved_solution_states[i]) = solutionState(i);
516 : }
517 :
518 230 : if (!_saved_dot_old && solutionUDotOld())
519 0 : _saved_dot_old = &addVector("save_solution_dot_old", false, PARALLEL);
520 230 : if (!_saved_dotdot_old && solutionUDotDotOld())
521 0 : _saved_dotdot_old = &addVector("save_solution_dotdot_old", false, PARALLEL);
522 :
523 230 : if (solutionUDotOld())
524 0 : *_saved_dot_old = *solutionUDotOld();
525 :
526 230 : if (solutionUDotDotOld())
527 0 : *_saved_dotdot_old = *solutionUDotDotOld();
528 230 : }
529 :
530 : /**
531 : * Restore the old and older solutions when the saved solutions present.
532 : */
533 : void
534 230 : SystemBase::restoreOldSolutions()
535 : {
536 : const auto states =
537 230 : _solution_states[static_cast<unsigned short>(Moose::SolutionIterationType::Time)].size();
538 230 : if (states > 1)
539 630 : for (unsigned int i = 1; i <= states - 1; ++i)
540 420 : if (_saved_solution_states[i])
541 : {
542 420 : solutionState(i) = *(_saved_solution_states[i]);
543 420 : removeVector("save_solution_state_" + std::to_string(i));
544 420 : _saved_solution_states[i] = nullptr;
545 : }
546 :
547 230 : 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 230 : 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 230 : }
560 :
561 : SparseMatrix<Number> &
562 320 : SystemBase::addMatrix(TagID tag)
563 : {
564 320 : 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 320 : if (hasMatrix(tag))
572 0 : return getMatrix(tag);
573 :
574 320 : const auto matrix_name = _subproblem.matrixTagName(tag);
575 320 : SparseMatrix<Number> & mat = system().add_matrix(matrix_name);
576 320 : associateMatrixToTag(mat, tag);
577 :
578 320 : return mat;
579 320 : }
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 63393 : SystemBase::addVector(const std::string & vector_name, const bool project, const ParallelType type)
601 : {
602 63393 : if (hasVector(vector_name))
603 434 : return getVector(vector_name);
604 :
605 62959 : NumericVector<Number> & vec = system().add_vector(vector_name, project, type);
606 62959 : return vec;
607 : }
608 :
609 : NumericVector<Number> &
610 150642 : SystemBase::addVector(TagID tag, const bool project, const ParallelType type)
611 : {
612 150642 : 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 150642 : 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 150642 : const auto vector_name = _subproblem.vectorTagName(tag);
634 150642 : NumericVector<Number> & vec = system().add_vector(vector_name, project, type);
635 150642 : associateVectorToTag(vec, tag);
636 :
637 150642 : return vec;
638 150642 : }
639 :
640 : void
641 17548323 : SystemBase::closeTaggedVector(const TagID tag)
642 : {
643 17548323 : 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 17548323 : 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 17548323 : getVector(tag).close();
656 17548323 : }
657 :
658 : void
659 6087665 : SystemBase::closeTaggedVectors(const std::set<TagID> & tags)
660 : {
661 23635988 : for (const auto tag : tags)
662 17548323 : closeTaggedVector(tag);
663 6087665 : }
664 :
665 : void
666 8774488 : SystemBase::zeroTaggedVector(const TagID tag)
667 : {
668 8774488 : 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 8774488 : 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 8774488 : if (!_subproblem.vectorTagNotZeroed(tag))
681 8774466 : getVector(tag).zero();
682 8774488 : }
683 :
684 : void
685 3043901 : SystemBase::zeroTaggedVectors(const std::set<TagID> & tags)
686 : {
687 11818389 : for (const auto tag : tags)
688 8774488 : zeroTaggedVector(tag);
689 3043901 : }
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 160283 : SystemBase::addVariable(const std::string & var_type,
711 : const std::string & name,
712 : InputParameters & parameters)
713 : {
714 160283 : _numbered_vars.resize(libMesh::n_threads());
715 :
716 160283 : 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 160283 : std::set<SubdomainID> blocks;
721 160283 : const auto & block_param = parameters.get<std::vector<SubdomainName>>("block");
722 169973 : for (const auto & subdomain_name : block_param)
723 : {
724 9690 : SubdomainID blk_id = _mesh.getSubdomainID(subdomain_name);
725 9690 : blocks.insert(blk_id);
726 : }
727 :
728 : const auto fe_type =
729 320566 : FEType(Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")),
730 480849 : Utility::string_to_enum<FEFamily>(parameters.get<MooseEnum>("family")));
731 160283 : const auto fe_field_type = FEInterface::field_type(fe_type);
732 :
733 : unsigned int var_num;
734 :
735 160283 : if (var_type == "ArrayMooseVariable")
736 : {
737 2225 : if (fe_field_type == TYPE_VECTOR)
738 0 : mooseError("Vector family type cannot be used in an array variable");
739 :
740 2225 : std::vector<std::string> array_var_component_names;
741 2225 : const bool has_array_names = parameters.isParamValid("array_var_component_names");
742 2225 : if (has_array_names)
743 : {
744 : array_var_component_names =
745 148 : parameters.get<std::vector<std::string>>("array_var_component_names");
746 148 : 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 2221 : std::vector<std::string> var_names;
758 7287 : for (unsigned int i = 0; i < components; i++)
759 : {
760 5066 : if (!has_array_names)
761 4550 : array_var_component_names.push_back(std::to_string(i));
762 5066 : 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 4442 : parameters.set<std::vector<std::string>>("array_var_component_names") =
767 2221 : 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 2221 : var_num = system().add_variables(var_names, fe_type, &blocks) - (components - 1);
772 :
773 : // Set as array variable
774 2221 : 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 2217 : parameters.set<bool>("array") = true;
779 2217 : }
780 : else
781 : {
782 158058 : 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 158054 : var_num = system().add_variable(name, fe_type, &blocks);
787 : }
788 :
789 160271 : parameters.set<unsigned int>("_var_num") = var_num;
790 160271 : parameters.set<SystemBase *>("_system_base") = this;
791 :
792 336421 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
793 : {
794 176154 : parameters.set<THREAD_ID>("tid") = tid;
795 : std::shared_ptr<MooseVariableBase> var =
796 176154 : _factory.create<MooseVariableBase>(var_type, name, parameters, tid);
797 :
798 176150 : _vars[tid].add(name, var);
799 :
800 176150 : if (auto fe_var = dynamic_cast<MooseVariableFieldBase *>(var.get()))
801 : {
802 172860 : auto required_size = var_num + components;
803 172860 : if (required_size > _numbered_vars[tid].size())
804 172860 : _numbered_vars[tid].resize(required_size);
805 348888 : for (MooseIndex(components) component = 0; component < components; ++component)
806 176028 : _numbered_vars[tid][var_num + component] = fe_var;
807 :
808 172860 : if (auto * const functor = dynamic_cast<Moose::FunctorBase<ADReal> *>(fe_var))
809 168817 : _subproblem.addFunctor(name, *functor, tid);
810 4043 : else if (auto * const functor = dynamic_cast<Moose::FunctorBase<ADRealVectorValue> *>(fe_var))
811 1571 : _subproblem.addFunctor(name, *functor, tid);
812 2472 : else if (auto * const functor = dynamic_cast<Moose::FunctorBase<RealEigenVector> *>(fe_var))
813 2472 : _subproblem.addFunctor(name, *functor, tid);
814 : else
815 0 : mooseError("This should be a functor");
816 : }
817 :
818 176150 : if (var->blockRestricted())
819 19093 : for (const SubdomainID & id : var->blockIDs())
820 21412 : for (MooseIndex(components) component = 0; component < components; ++component)
821 10919 : _var_map[var_num + component].insert(id);
822 : else
823 337842 : for (MooseIndex(components) component = 0; component < components; ++component)
824 170292 : _var_map[var_num + component] = std::set<SubdomainID>();
825 176150 : }
826 :
827 : // getMaxVariableNumber is an API method used in Rattlesnake
828 160267 : if (var_num > _max_var_number)
829 85850 : _max_var_number = var_num;
830 160267 : _du_dot_du.resize(var_num + 1);
831 160267 : }
832 :
833 : bool
834 5862993 : SystemBase::hasVariable(const std::string & var_name) const
835 : {
836 5862993 : auto & names = getVariableNames();
837 5862993 : if (system().has_variable(var_name))
838 4601552 : return system().variable_type(var_name).family != SCALAR;
839 1261441 : if (std::find(names.begin(), names.end(), var_name) != names.end())
840 : // array variable
841 33548 : return true;
842 : else
843 1227893 : 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 506543 : SystemBase::hasScalarVariable(const std::string & var_name) const
860 : {
861 506543 : if (system().has_variable(var_name))
862 272646 : return system().variable_type(var_name).family == SCALAR;
863 : else
864 233897 : return false;
865 : }
866 :
867 : bool
868 91025 : SystemBase::isScalarVariable(unsigned int var_num) const
869 : {
870 91025 : return (system().variable(var_num).type().family == SCALAR);
871 : }
872 :
873 : unsigned int
874 54797502 : SystemBase::nVariables() const
875 : {
876 54797502 : unsigned int n = nFieldVariables();
877 54797502 : n += _vars[0].scalars().size();
878 :
879 54797502 : return n;
880 : }
881 :
882 : unsigned int
883 54797502 : SystemBase::nFieldVariables() const
884 : {
885 54797502 : unsigned int n = 0;
886 126188836 : for (auto & var : _vars[0].fieldVariables())
887 71391334 : n += var->count();
888 :
889 54797502 : return n;
890 : }
891 :
892 : unsigned int
893 2452323 : SystemBase::nFVVariables() const
894 : {
895 2452323 : unsigned int n = 0;
896 5191115 : for (auto & var : _vars[0].fieldVariables())
897 2738792 : if (var->isFV())
898 1528425 : n += var->count();
899 :
900 2452323 : return n;
901 : }
902 :
903 : /**
904 : * Check if the named vector exists in the system.
905 : */
906 : bool
907 464143 : SystemBase::hasVector(const std::string & name) const
908 : {
909 464143 : return system().have_vector(name);
910 : }
911 :
912 : /**
913 : * Get a raw NumericVector with the given name.
914 : */
915 : NumericVector<Number> &
916 5259 : SystemBase::getVector(const std::string & name)
917 : {
918 5259 : 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 2486813279 : SystemBase::getVector(TagID tag)
929 : {
930 2486813279 : 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 2486813279 : return *_tagged_vectors[tag];
943 : }
944 :
945 : const NumericVector<Number> &
946 600 : SystemBase::getVector(TagID tag) const
947 : {
948 600 : 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 600 : return *_tagged_vectors[tag];
961 : }
962 :
963 : void
964 9969592 : SystemBase::associateVectorToTag(NumericVector<Number> & vec, TagID tag)
965 : {
966 9969592 : if (!_subproblem.vectorTagExists(tag))
967 0 : mooseError("Cannot associate vector to tag ", tag, " because that tag does not exist");
968 :
969 9969592 : if (_tagged_vectors.size() < tag + 1)
970 238928 : _tagged_vectors.resize(tag + 1);
971 :
972 9969592 : _tagged_vectors[tag] = &vec;
973 9969592 : }
974 :
975 : void
976 3067180 : SystemBase::disassociateVectorFromTag(NumericVector<Number> & vec, TagID tag)
977 : {
978 3067180 : if (!_subproblem.vectorTagExists(tag))
979 0 : mooseError("Cannot disassociate vector from tag ", tag, " because that tag does not exist");
980 3067180 : 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 3067180 : disassociateVectorFromTag(tag);
984 3067180 : }
985 :
986 : void
987 3188684 : SystemBase::disassociateVectorFromTag(TagID tag)
988 : {
989 3188684 : if (!_subproblem.vectorTagExists(tag))
990 0 : mooseError("Cannot disassociate vector from tag ", tag, " because that tag does not exist");
991 :
992 3188684 : if (_tagged_vectors.size() < tag + 1)
993 0 : _tagged_vectors.resize(tag + 1);
994 3188684 : _tagged_vectors[tag] = nullptr;
995 3188684 : }
996 :
997 : void
998 30376 : SystemBase::disassociateDefaultVectorTags()
999 : {
1000 30376 : const auto tags = defaultVectorTags();
1001 182256 : for (const auto tag : tags)
1002 151880 : if (_subproblem.vectorTagExists(tag))
1003 121504 : disassociateVectorFromTag(tag);
1004 30376 : }
1005 :
1006 : SparseMatrix<Number> &
1007 1831297793 : SystemBase::getMatrix(TagID tag)
1008 : {
1009 1831297793 : 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 1831297793 : return *_tagged_matrices[tag];
1022 : }
1023 :
1024 : const SparseMatrix<Number> &
1025 300 : SystemBase::getMatrix(TagID tag) const
1026 : {
1027 300 : 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 300 : return *_tagged_matrices[tag];
1040 : }
1041 :
1042 : void
1043 879422 : SystemBase::closeTaggedMatrices(const std::set<TagID> & tags)
1044 : {
1045 2619696 : for (auto tag : tags)
1046 1740274 : if (hasMatrix(tag))
1047 882423 : getMatrix(tag).close();
1048 879422 : }
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 497920 : SystemBase::associateMatrixToTag(SparseMatrix<Number> & matrix, TagID tag)
1060 : {
1061 497920 : if (!_subproblem.matrixTagExists(tag))
1062 0 : mooseError("Cannot associate matrix to tag ", tag, " because that tag does not exist");
1063 :
1064 497920 : if (_tagged_matrices.size() < tag + 1)
1065 36640 : _tagged_matrices.resize(tag + 1);
1066 :
1067 497920 : _tagged_matrices[tag] = &matrix;
1068 497920 : }
1069 :
1070 : void
1071 492623 : SystemBase::disassociateMatrixFromTag(SparseMatrix<Number> & matrix, TagID tag)
1072 : {
1073 492623 : if (!_subproblem.matrixTagExists(tag))
1074 0 : mooseError("Cannot disassociate matrix from tag ", tag, " because that tag does not exist");
1075 492623 : 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 492623 : disassociateMatrixFromTag(tag);
1079 492623 : }
1080 :
1081 : void
1082 504182 : SystemBase::disassociateMatrixFromTag(TagID tag)
1083 : {
1084 504182 : if (!_subproblem.matrixTagExists(tag))
1085 0 : mooseError("Cannot disassociate matrix from tag ", tag, " because that tag does not exist");
1086 :
1087 504182 : if (_tagged_matrices.size() < tag + 1)
1088 1335 : _tagged_matrices.resize(tag + 1);
1089 504182 : _tagged_matrices[tag] = nullptr;
1090 504182 : }
1091 :
1092 : void
1093 3853 : SystemBase::disassociateDefaultMatrixTags()
1094 : {
1095 3853 : const auto tags = defaultMatrixTags();
1096 15412 : for (const auto tag : tags)
1097 11559 : if (_subproblem.matrixTagExists(tag))
1098 11559 : disassociateMatrixFromTag(tag);
1099 3853 : }
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 3327025 : SystemBase::deactiveAllMatrixTags()
1127 : {
1128 3327025 : auto num_matrix_tags = _subproblem.numMatrixTags();
1129 :
1130 3327025 : _matrix_tag_active_flags.resize(num_matrix_tags);
1131 :
1132 10083712 : for (decltype(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
1133 6756687 : _matrix_tag_active_flags[tag] = false;
1134 3327025 : }
1135 :
1136 : void
1137 3525731 : SystemBase::activeAllMatrixTags()
1138 : {
1139 3525731 : auto num_matrix_tags = _subproblem.numMatrixTags();
1140 :
1141 3525731 : _matrix_tag_active_flags.resize(num_matrix_tags);
1142 :
1143 10683645 : for (decltype(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
1144 7157914 : if (hasMatrix(tag))
1145 491514 : _matrix_tag_active_flags[tag] = true;
1146 : else
1147 6666400 : _matrix_tag_active_flags[tag] = false;
1148 3525731 : }
1149 :
1150 : bool
1151 54032 : SystemBase::matrixTagActive(TagID tag) const
1152 : {
1153 : mooseAssert(_subproblem.matrixTagExists(tag), "Matrix tag " << tag << " does not exist");
1154 :
1155 54032 : return tag < _matrix_tag_active_flags.size() && _matrix_tag_active_flags[tag];
1156 : }
1157 :
1158 : unsigned int
1159 1762527614 : SystemBase::number() const
1160 : {
1161 1762527614 : return system().number();
1162 : }
1163 :
1164 : DofMap &
1165 1799556 : SystemBase::dofMap()
1166 : {
1167 1799556 : 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 395 : SystemBase::addVariableToCopy(const std::string & dest_name,
1178 : const std::string & source_name,
1179 : const std::string & timestep)
1180 : {
1181 395 : _var_to_copy.push_back(VarCopyInfo(dest_name, source_name, timestep));
1182 395 : }
1183 :
1184 : void
1185 537 : SystemBase::copyVars(ExodusII_IO & io)
1186 : {
1187 537 : int n_steps = io.get_num_time_steps();
1188 :
1189 537 : bool did_copy = false;
1190 911 : for (const auto & vci : _var_to_copy)
1191 : {
1192 378 : int timestep = -1;
1193 :
1194 378 : if (vci._timestep == "LATEST")
1195 : // Use the last time step in the file from which to retrieve the solution
1196 292 : timestep = n_steps;
1197 : else
1198 : {
1199 86 : timestep = MooseUtils::convert<int>(vci._timestep);
1200 86 : 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 374 : did_copy = true;
1209 :
1210 374 : if (hasVariable(vci._dest_name))
1211 : {
1212 330 : const auto & var = getVariable(0, vci._dest_name);
1213 330 : if (var.isArray())
1214 : {
1215 22 : const auto & array_var = getFieldVariable<RealEigenVector>(0, vci._dest_name);
1216 66 : for (MooseIndex(var.count()) i = 0; i < var.count(); ++i)
1217 : {
1218 44 : const auto & exodus_var = var.arrayVariableComponent(i);
1219 44 : const auto & system_var = array_var.componentName(i);
1220 44 : if (var.isNodal())
1221 22 : io.copy_nodal_solution(system(), exodus_var, system_var, timestep);
1222 : else
1223 22 : io.copy_elemental_solution(system(), exodus_var, system_var, timestep);
1224 : }
1225 : }
1226 : else
1227 : {
1228 308 : if (var.isNodal())
1229 264 : io.copy_nodal_solution(system(), vci._dest_name, vci._source_name, timestep);
1230 : else
1231 44 : io.copy_elemental_solution(system(), vci._dest_name, vci._source_name, timestep);
1232 : }
1233 : }
1234 44 : else if (hasScalarVariable(vci._dest_name))
1235 176 : 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 533 : if (did_copy)
1241 272 : solution().close();
1242 621 : }
1243 :
1244 : void
1245 3804672 : SystemBase::update()
1246 : {
1247 3804672 : system().update();
1248 3804672 : }
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 102117 : SystemBase::copySolutionsBackwards()
1261 : {
1262 102117 : system().update();
1263 102117 : copyOldSolutions();
1264 102117 : copyPreviousNonlinearSolutions();
1265 102117 : }
1266 :
1267 : /**
1268 : * Shifts the solutions backwards in nonlinear iteration history
1269 : */
1270 : void
1271 102117 : SystemBase::copyPreviousNonlinearSolutions()
1272 : {
1273 : // 1 is for nonlinear, 0 is for time, we do this for nonlinear only here
1274 102117 : const auto states = _solution_states[1].size();
1275 102117 : if (states > 1)
1276 332 : for (unsigned int i = states - 1; i > 0; --i)
1277 332 : solutionState(i, Moose::SolutionIterationType::Nonlinear) =
1278 166 : solutionState(i - 1, Moose::SolutionIterationType::Nonlinear);
1279 :
1280 102117 : if (solutionPreviousNewton())
1281 166 : *solutionPreviousNewton() = *currentSolution();
1282 102117 : }
1283 :
1284 : /**
1285 : * Shifts the solutions backwards in time
1286 : */
1287 : void
1288 623655 : 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 623655 : const auto states = _solution_states[0].size();
1293 623655 : if (states > 1)
1294 933880 : for (unsigned int i = states - 1; i > 0; --i)
1295 480792 : solutionState(i) = solutionState(i - 1);
1296 :
1297 623655 : if (solutionUDotOld())
1298 10538 : *solutionUDotOld() = *solutionUDot();
1299 623655 : if (solutionUDotDotOld())
1300 10538 : *solutionUDotDotOld() = *solutionUDotDot();
1301 623655 : }
1302 :
1303 : /**
1304 : * Restore current solutions (call after your solve failed)
1305 : */
1306 : void
1307 8428 : SystemBase::restoreSolutions()
1308 : {
1309 8428 : if (!hasSolutionState(1))
1310 0 : mooseError("Cannot restore solutions without old solution");
1311 :
1312 8428 : *(const_cast<NumericVector<Number> *&>(currentSolution())) = solutionOld();
1313 8428 : solution() = solutionOld();
1314 8428 : if (solutionUDotOld())
1315 0 : *solutionUDot() = *solutionUDotOld();
1316 8428 : if (solutionUDotDotOld())
1317 0 : *solutionUDotDot() = *solutionUDotDotOld();
1318 8428 : if (solutionPreviousNewton())
1319 0 : *solutionPreviousNewton() = solutionOld();
1320 8428 : system().update();
1321 8428 : }
1322 :
1323 : void
1324 420 : SystemBase::removeVector(const std::string & name)
1325 : {
1326 420 : system().remove_vector(name);
1327 420 : }
1328 :
1329 : const std::string &
1330 415733 : SystemBase::name() const
1331 : {
1332 415733 : return system().name();
1333 : }
1334 :
1335 : NumericVector<Number> *
1336 111723 : SystemBase::solutionPreviousNewton()
1337 : {
1338 111723 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
1339 1344 : return &getVector(Moose::PREVIOUS_NL_SOLUTION_TAG);
1340 : else
1341 110379 : 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 115356 : SystemBase::initSolutionState()
1355 : {
1356 : // Default is the current solution
1357 115356 : unsigned int state = 0;
1358 :
1359 : // Add additional states as required by the variable states requested
1360 267819 : for (const auto & var : getVariables(/* tid = */ 0))
1361 152463 : state = std::max(state, var->oldestSolutionStateRequested());
1362 118310 : for (const auto & var : getScalarVariables(/* tid = */ 0))
1363 2954 : state = std::max(state, var->oldestSolutionStateRequested());
1364 :
1365 115356 : needSolutionState(state, Moose::SolutionIterationType::Time);
1366 :
1367 115356 : _solution_states_initialized = true;
1368 115356 : }
1369 :
1370 : TagName
1371 58725 : SystemBase::oldSolutionStateVectorName(const unsigned int state,
1372 : const Moose::SolutionIterationType iteration_type) const
1373 : {
1374 : mooseAssert(state != 0, "Not an old state");
1375 :
1376 58725 : if (iteration_type == Moose::SolutionIterationType::Time)
1377 : {
1378 58547 : if (state == 1)
1379 57131 : return Moose::OLD_SOLUTION_TAG;
1380 1416 : else if (state == 2)
1381 1416 : return Moose::OLDER_SOLUTION_TAG;
1382 : }
1383 178 : else if (iteration_type == Moose::SolutionIterationType::Nonlinear && state == 1)
1384 178 : 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 145782782 : SystemBase::solutionState(const unsigned int state,
1391 : const Moose::SolutionIterationType iteration_type) const
1392 : {
1393 145782782 : 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 145782782 : const auto & solution_states = _solution_states[static_cast<unsigned short>(iteration_type)];
1407 :
1408 145782782 : 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 145782782 : return *solution_states[state];
1416 : }
1417 :
1418 : NumericVector<Number> &
1419 55980496 : SystemBase::solutionState(const unsigned int state,
1420 : const Moose::SolutionIterationType iteration_type)
1421 : {
1422 55980496 : if (!hasSolutionState(state, iteration_type))
1423 64894 : needSolutionState(state, iteration_type);
1424 55980496 : return *_solution_states[static_cast<unsigned short>(iteration_type)][state];
1425 : }
1426 :
1427 : void
1428 196717 : 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 196717 : if (hasSolutionState(state, iteration_type))
1437 82146 : return;
1438 :
1439 114571 : auto & solution_states = _solution_states[static_cast<unsigned short>(iteration_type)];
1440 114571 : solution_states.resize(state + 1);
1441 :
1442 : // The 0-th (current) solution state is owned by libMesh
1443 114571 : if (!solution_states[0])
1444 113108 : 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 174500 : for (unsigned int i = 1; i <= state; ++i)
1450 59929 : if (!solution_states[i])
1451 : {
1452 58725 : auto tag = _subproblem.addVectorTag(oldSolutionStateVectorName(i, iteration_type),
1453 : Moose::VECTOR_TAG_SOLUTION);
1454 58725 : 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 706 : SystemBase::applyScalingFactors(const std::vector<Real> & inverse_scaling_factors)
1463 : {
1464 1483 : for (MooseIndex(_vars) thread = 0; thread < _vars.size(); ++thread)
1465 : {
1466 777 : auto & field_variables = _vars[thread].fieldVariables();
1467 1958 : for (MooseIndex(field_variables) i = 0, p = 0; i < field_variables.size(); ++i)
1468 : {
1469 1181 : auto factors = field_variables[i]->arrayScalingFactor();
1470 2373 : for (unsigned int j = 0; j < field_variables[i]->count(); ++j, ++p)
1471 1192 : factors[j] /= inverse_scaling_factors[p];
1472 :
1473 1181 : field_variables[i]->scalingFactor(factors);
1474 1181 : }
1475 :
1476 777 : auto offset = field_variables.size();
1477 :
1478 777 : auto & scalar_variables = _vars[thread].scalars();
1479 898 : for (MooseIndex(scalar_variables) i = 0; i < scalar_variables.size(); ++i)
1480 242 : scalar_variables[i]->scalingFactor(
1481 121 : {1. / inverse_scaling_factors[offset + i] * scalar_variables[i]->scalingFactor()});
1482 :
1483 777 : if (thread == 0 && _verbose)
1484 : {
1485 279 : _console << "Automatic scaling factors:\n";
1486 279 : auto original_flags = _console.flags();
1487 279 : auto original_precision = _console.precision();
1488 279 : _console.unsetf(std::ios_base::floatfield);
1489 279 : _console.precision(6);
1490 :
1491 721 : for (const auto & field_variable : field_variables)
1492 : {
1493 442 : const auto & factors = field_variable->arrayScalingFactor();
1494 442 : _console << " " << field_variable->name() << ":";
1495 894 : for (const auto i : make_range(field_variable->count()))
1496 452 : _console << " " << factors[i];
1497 442 : _console << "\n";
1498 : }
1499 372 : for (const auto & scalar_variable : scalar_variables)
1500 186 : _console << " " << scalar_variable->name() << ": " << scalar_variable->scalingFactor()
1501 93 : << "\n";
1502 279 : _console << "\n" << std::endl;
1503 :
1504 : // restore state
1505 279 : _console.flags(original_flags);
1506 279 : _console.precision(original_precision);
1507 : }
1508 : }
1509 706 : }
1510 :
1511 : void
1512 352 : SystemBase::addScalingVector()
1513 : {
1514 352 : addVector("scaling_factors", /*project=*/false, libMesh::ParallelType::GHOSTED);
1515 352 : _subproblem.hasScalingVector(number());
1516 352 : }
1517 :
1518 : bool
1519 71028219 : SystemBase::computingScalingJacobian() const
1520 : {
1521 71028219 : return _subproblem.computingScalingJacobian();
1522 : }
1523 :
1524 : void
1525 114764 : SystemBase::initialSetup()
1526 : {
1527 240813 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1528 126049 : _vars[tid].initialSetup();
1529 :
1530 : // If we need raw gradients, we initialize them here.
1531 114764 : bool gradient_storage_initialized = false;
1532 266843 : for (const auto & field_var : _vars[0].fieldVariables())
1533 152079 : if (!gradient_storage_initialized && field_var->needsGradientVectorStorage())
1534 : {
1535 642 : _raw_grad_container.clear();
1536 1663 : for (const auto i : make_range(this->_mesh.dimension()))
1537 : {
1538 1021 : libmesh_ignore(i);
1539 1021 : _raw_grad_container.push_back(currentSolution()->zero_clone());
1540 : }
1541 : }
1542 114764 : }
1543 :
1544 : void
1545 582995 : SystemBase::timestepSetup()
1546 : {
1547 1222161 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1548 639166 : _vars[tid].timestepSetup();
1549 582995 : }
1550 :
1551 : void
1552 3768531 : SystemBase::customSetup(const ExecFlagType & exec_type)
1553 : {
1554 7904408 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1555 4135877 : _vars[tid].customSetup(exec_type);
1556 3768531 : }
1557 :
1558 : void
1559 5540300 : SystemBase::subdomainSetup()
1560 : {
1561 11988097 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1562 6447797 : _vars[tid].subdomainSetup();
1563 5540300 : }
1564 :
1565 : void
1566 6427558 : SystemBase::residualSetup()
1567 : {
1568 13477985 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1569 7050427 : _vars[tid].residualSetup();
1570 6427558 : }
1571 :
1572 : void
1573 1016050 : SystemBase::jacobianSetup()
1574 : {
1575 2133578 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1576 1117528 : _vars[tid].jacobianSetup();
1577 1016050 : }
1578 :
1579 : void
1580 704104 : SystemBase::clearAllDofIndices()
1581 : {
1582 1475650 : for (auto & var_warehouse : _vars)
1583 771546 : var_warehouse.clearAllDofIndices();
1584 704104 : }
1585 :
1586 : void
1587 13746547 : SystemBase::setActiveVariableCoupleableVectorTags(const std::set<TagID> & vtags, THREAD_ID tid)
1588 : {
1589 13746547 : _vars[tid].setActiveVariableCoupleableVectorTags(vtags);
1590 13746547 : }
1591 :
1592 : void
1593 99768 : SystemBase::setActiveScalarVariableCoupleableVectorTags(const std::set<TagID> & vtags,
1594 : THREAD_ID tid)
1595 : {
1596 99768 : _vars[tid].setActiveScalarVariableCoupleableVectorTags(vtags);
1597 99768 : }
1598 :
1599 : void
1600 56928 : SystemBase::addDotVectors()
1601 : {
1602 56928 : if (_fe_problem.uDotRequested())
1603 56928 : _u_dot = &addVector("u_dot", true, GHOSTED);
1604 56928 : if (_fe_problem.uDotOldRequested())
1605 304 : _u_dot_old = &addVector("u_dot_old", true, GHOSTED);
1606 56928 : if (_fe_problem.uDotDotRequested())
1607 304 : _u_dotdot = &addVector("u_dotdot", true, GHOSTED);
1608 56928 : if (_fe_problem.uDotDotOldRequested())
1609 304 : _u_dotdot_old = &addVector("u_dotdot_old", true, GHOSTED);
1610 56928 : }
1611 :
1612 : NumericVector<Number> &
1613 1044 : SystemBase::serializedSolution()
1614 : {
1615 1044 : if (!_serialized_solution.get())
1616 : {
1617 104 : _serialized_solution = NumericVector<Number>::build(_communicator);
1618 104 : _serialized_solution->init(system().n_dofs(), false, SERIAL);
1619 : }
1620 :
1621 1044 : return *_serialized_solution;
1622 : }
1623 :
1624 : void
1625 56952 : SystemBase::addTimeIntegrator(const std::string & type,
1626 : const std::string & name,
1627 : InputParameters & parameters)
1628 : {
1629 56952 : parameters.set<SystemBase *>("_sys") = this;
1630 56952 : _time_integrators.push_back(_factory.create<TimeIntegrator>(type, name, parameters));
1631 56952 : }
1632 :
1633 : void
1634 3492 : SystemBase::copyTimeIntegrators(const SystemBase & other_sys)
1635 : {
1636 3492 : _time_integrators = other_sys._time_integrators;
1637 3492 : }
1638 :
1639 : const TimeIntegrator *
1640 684624 : SystemBase::queryTimeIntegrator(const unsigned int var_num) const
1641 : {
1642 685568 : for (auto & ti : _time_integrators)
1643 435892 : if (ti->integratesVar(var_num))
1644 434948 : return ti.get();
1645 :
1646 249676 : return nullptr;
1647 : }
1648 :
1649 : const TimeIntegrator &
1650 90 : SystemBase::getTimeIntegrator(const unsigned int var_num) const
1651 : {
1652 90 : const auto * const ti = queryTimeIntegrator(var_num);
1653 :
1654 90 : if (ti)
1655 90 : 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 535678 : SystemBase::getTimeIntegrators()
1663 : {
1664 535678 : return _time_integrators;
1665 : }
1666 :
1667 : const Number &
1668 1130712840 : SystemBase::duDotDu(const unsigned int var_num) const
1669 : {
1670 1130712840 : return _du_dot_du[var_num];
1671 : }
1672 :
1673 : const std::set<SubdomainID> &
1674 336 : SystemBase::getSubdomainsForVar(const std::string & var_name) const
1675 : {
1676 336 : return getSubdomainsForVar(getVariable(0, var_name).number());
1677 : }
1678 :
1679 : std::string
1680 71052 : SystemBase::prefix() const
1681 : {
1682 142104 : 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);
|