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