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 "Simulation.h"
11 : #include "FEProblemBase.h"
12 : #include "AddVariableAction.h"
13 : #include "MooseObjectAction.h"
14 : #include "Transient.h"
15 : #include "HeatConductionModel.h"
16 : #include "HeatStructureInterface.h"
17 : #include "FlowChannelBase.h"
18 : #include "FlowJunction.h"
19 :
20 : #include "ClosuresBase.h"
21 : #include "FluidProperties.h"
22 : #include "THMControl.h"
23 : #include "TerminateControl.h"
24 : #include "THMMesh.h"
25 : #include "RelationshipManager.h"
26 : #include "NonlinearSystemBase.h"
27 : #include "TimeIntegrator.h"
28 : #include "ExplicitTimeIntegrator.h"
29 : #include "ExplicitEuler.h"
30 : #include "ExplicitRK2.h"
31 : #include "ExplicitTVDRK2.h"
32 :
33 : #include "libmesh/string_to_enum.h"
34 :
35 : using namespace libMesh;
36 :
37 : std::map<VariableName, int> Simulation::_component_variable_order_map;
38 :
39 : void
40 85528 : Simulation::setComponentVariableOrder(const VariableName & var, int index)
41 : {
42 85528 : _component_variable_order_map[var] = index;
43 85528 : }
44 :
45 3904 : Simulation::Simulation(FEProblemBase & fe_problem, const InputParameters & pars)
46 : : ParallelObject(fe_problem.comm()),
47 3904 : LoggingInterface(_log),
48 3904 : _thm_mesh(*static_cast<THMMesh *>(pars.get<MooseMesh *>("mesh"))),
49 3904 : _fe_problem(fe_problem),
50 3904 : _thm_app(static_cast<ThermalHydraulicsApp &>(*pars.get<MooseApp *>("_moose_app"))),
51 3904 : _thm_factory(_thm_app.getFactory()),
52 3904 : _thm_pars(pars),
53 : _flow_fe_type(FEType(CONSTANT, MONOMIAL)),
54 3904 : _implicit_time_integration(true),
55 3904 : _check_jacobian(false),
56 3904 : _output_vector_velocity(true),
57 7808 : _zero(0)
58 : {
59 3904 : bool second_order_mesh = pars.get<bool>("2nd_order_mesh");
60 3904 : HeatConductionModel::_fe_type =
61 3904 : second_order_mesh ? FEType(SECOND, LAGRANGE) : FEType(FIRST, LAGRANGE);
62 3904 : }
63 :
64 3723 : Simulation::~Simulation()
65 : {
66 14167 : for (auto && k : _control_data)
67 10444 : delete k.second;
68 11169 : }
69 :
70 : void
71 18241 : Simulation::augmentSparsity(const dof_id_type & elem_id1, const dof_id_type & elem_id2)
72 : {
73 : auto it = _sparsity_elem_augmentation.find(elem_id1);
74 18241 : if (it == _sparsity_elem_augmentation.end())
75 15501 : it = _sparsity_elem_augmentation.insert(_sparsity_elem_augmentation.begin(),
76 15501 : {elem_id1, std::vector<dof_id_type>()});
77 18241 : it->second.push_back(elem_id2);
78 :
79 : it = _sparsity_elem_augmentation.find(elem_id2);
80 18241 : if (it == _sparsity_elem_augmentation.end())
81 11016 : it = _sparsity_elem_augmentation.insert(_sparsity_elem_augmentation.begin(),
82 11016 : {elem_id2, std::vector<dof_id_type>()});
83 18241 : it->second.push_back(elem_id1);
84 18241 : }
85 :
86 : void
87 3900 : Simulation::buildMesh()
88 : {
89 3900 : if (_components.size() == 0)
90 : return;
91 :
92 : // build mesh
93 19395 : for (auto && comp : _components)
94 15638 : comp->executeSetupMesh();
95 : }
96 :
97 : void
98 3270 : Simulation::setupQuadrature()
99 : {
100 3270 : if (_components.size() == 0)
101 : return;
102 :
103 : Order order = CONSTANT;
104 : unsigned int n_flow_channels = 0;
105 : unsigned int n_heat_structures = 0;
106 :
107 16110 : for (auto && comp : _components)
108 : {
109 12983 : auto flow_channel = dynamic_cast<FlowChannelBase *>(comp.get());
110 12983 : if (flow_channel != nullptr)
111 3323 : n_flow_channels++;
112 :
113 12983 : auto hs_interface = dynamic_cast<HeatStructureInterface *>(comp.get());
114 12983 : if (hs_interface)
115 1762 : n_heat_structures++;
116 : }
117 :
118 3127 : if (n_flow_channels > 0)
119 : {
120 : const FEType & fe_type = getFlowFEType();
121 2045 : if (fe_type.default_quadrature_order() > order)
122 : order = fe_type.default_quadrature_order();
123 : }
124 3127 : if (n_heat_structures > 0)
125 : {
126 : const FEType & fe_type = HeatConductionModel::feType();
127 : if (fe_type.default_quadrature_order() > order)
128 : order = fe_type.default_quadrature_order();
129 : }
130 :
131 3127 : _fe_problem.createQRules(QGAUSS, order, order, order);
132 : }
133 :
134 : void
135 3900 : Simulation::initSimulation()
136 : {
137 : // sort the components using dependency resolver
138 3900 : DependencyResolver<std::shared_ptr<Component>> dependency_resolver;
139 19538 : for (const auto & comp : _components)
140 : {
141 15638 : dependency_resolver.addNode(comp);
142 26424 : for (const auto & dep : comp->getDependencies())
143 10786 : if (hasComponent(dep))
144 10727 : dependency_resolver.addEdge(_comp_by_name[dep], comp);
145 : }
146 :
147 3900 : _components = dependency_resolver.dfs();
148 3900 : }
149 :
150 : void
151 3891 : Simulation::initComponents()
152 : {
153 : // initialize components
154 19431 : for (auto && comp : _components)
155 15542 : comp->executeInit();
156 :
157 : // perform secondary initialization, which relies on init() being called
158 : // already for all components
159 19429 : for (auto && comp : _components)
160 15540 : comp->executeInitSecondary();
161 3889 : }
162 :
163 : void
164 3889 : Simulation::identifyLoops()
165 : {
166 : // loop over junctions and boundaries (non-geometrical components)
167 19429 : for (const auto & component : _components)
168 : {
169 : const auto flow_connection =
170 15540 : MooseSharedNamespace::dynamic_pointer_cast<Component1DConnection>(component);
171 15540 : if (flow_connection)
172 : {
173 : // create vector of names of this component and its connected flow channels, and then sort
174 : // them
175 6451 : std::vector<std::string> names = flow_connection->getConnectedComponentNames();
176 6451 : names.push_back(component->name());
177 6451 : std::sort(names.begin(), names.end());
178 :
179 : // pick first name alphabetically to be the proposed loop name
180 6451 : std::string proposed_loop_name = names[0];
181 :
182 21125 : for (const std::string & name : names)
183 : {
184 : // if the name is not yet in the map
185 14674 : if (_component_name_to_loop_name.find(name) == _component_name_to_loop_name.end())
186 : // just add the new map key; nothing else needs updating
187 10588 : _component_name_to_loop_name[name] = proposed_loop_name;
188 : else
189 : {
190 : // compare to the existing loop name for this component to make sure the
191 : // proposed loop name is first alphabetically
192 4086 : const std::string current_loop_name = _component_name_to_loop_name[name];
193 : // if proposed loop name comes later, change map values of the current
194 : // loop name to be the proposed name, and then update the proposed name
195 : // to be the current name
196 4086 : if (proposed_loop_name > current_loop_name)
197 : {
198 15921 : for (auto && entry : _component_name_to_loop_name)
199 13011 : if (entry.second == proposed_loop_name)
200 : entry.second = current_loop_name;
201 : proposed_loop_name = current_loop_name;
202 : }
203 : // if proposed loop name comes earlier, change map values of the current
204 : // loop name to be the proposed name
205 1176 : else if (proposed_loop_name < current_loop_name)
206 : {
207 2723 : for (auto && entry : _component_name_to_loop_name)
208 2358 : if (entry.second == current_loop_name)
209 : entry.second = proposed_loop_name;
210 : }
211 : }
212 : }
213 6451 : }
214 : }
215 :
216 : // get the list of loops
217 : std::vector<std::string> loops;
218 14477 : for (const auto & entry : _component_name_to_loop_name)
219 10588 : if (std::find(loops.begin(), loops.end(), entry.second) == loops.end())
220 2682 : loops.push_back(entry.second);
221 :
222 : // fill the map of loop name to model ID
223 6571 : for (const auto & loop : loops)
224 : {
225 : // find a flow channel in this loop and get its model ID
226 : THM::FlowModelID model_id;
227 : bool found_model_id = false;
228 3479 : for (const auto & component : _components)
229 : {
230 : const auto flow_chan_base_component =
231 3473 : MooseSharedNamespace::dynamic_pointer_cast<FlowChannelBase>(component);
232 3473 : if (flow_chan_base_component && (_component_name_to_loop_name[component->name()] == loop))
233 : {
234 2676 : model_id = flow_chan_base_component->getFlowModelID();
235 : found_model_id = true;
236 : break;
237 : }
238 : }
239 : // set the value in the map or throw an error
240 : if (found_model_id)
241 2676 : _loop_name_to_model_id[loop] = model_id;
242 : else
243 : logError("No FlowChannelBase-derived components were found in loop '", loop, "'");
244 : }
245 3889 : }
246 :
247 : void
248 16 : Simulation::printComponentLoops() const
249 : {
250 : // get the list of loops
251 : std::vector<std::string> loops;
252 152 : for (auto && entry : _component_name_to_loop_name)
253 136 : if (std::find(loops.begin(), loops.end(), entry.second) == loops.end())
254 32 : loops.push_back(entry.second);
255 :
256 : // for each loop
257 : Moose::out << "\nListing of component loops:" << std::endl;
258 48 : for (unsigned int i = 0; i < loops.size(); i++)
259 : {
260 32 : Moose::out << "\n Loop " << i + 1 << ":" << std::endl;
261 :
262 : // print out each component in the loop
263 304 : for (auto && entry : _component_name_to_loop_name)
264 272 : if (entry.second == loops[i])
265 : Moose::out << " " << entry.first << std::endl;
266 : }
267 : Moose::out << std::endl;
268 16 : }
269 :
270 : void
271 505 : Simulation::addSimVariable(bool nl, const VariableName & name, FEType fe_type, Real scaling_factor)
272 : {
273 505 : checkVariableNameLength(name);
274 :
275 505 : if (fe_type.family != SCALAR)
276 0 : mooseError("This method should only be used for scalar variables.");
277 :
278 505 : if (_vars.find(name) == _vars.end()) // variable is new
279 : {
280 505 : VariableInfo vi;
281 : InputParameters & params = vi._params;
282 :
283 505 : vi._nl = nl;
284 : vi._var_type = "MooseVariableScalar";
285 505 : params = _thm_factory.getValidParams(vi._var_type);
286 :
287 505 : auto family = AddVariableAction::getNonlinearVariableFamilies();
288 505 : family = Utility::enum_to_string(fe_type.family);
289 505 : params.set<MooseEnum>("family") = family;
290 :
291 505 : auto order = AddVariableAction::getNonlinearVariableOrders();
292 505 : order = Utility::enum_to_string<Order>(fe_type.order);
293 505 : params.set<MooseEnum>("order") = order;
294 :
295 505 : if (nl)
296 486 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
297 262 : else if (!MooseUtils::absoluteFuzzyEqual(scaling_factor, 1.0))
298 0 : mooseError("Aux variables cannot be provided a residual scaling factor.");
299 :
300 505 : _vars[name] = vi;
301 505 : }
302 : else
303 : // One of the two cases is true:
304 : // - This variable was previously added as a scalar variable, and scalar
305 : // variables should not be added more than once, since there is no block
306 : // restriction to extend, as there is in the field variable version of this
307 : // method.
308 : // - This variable was previously added as a field variable, and a variable
309 : // may have only one type (this method is used for scalar variables only).
310 0 : mooseError("The variable '", name, "' was already added.");
311 505 : }
312 :
313 : void
314 72044 : Simulation::addSimVariable(bool nl,
315 : const VariableName & name,
316 : FEType fe_type,
317 : const std::vector<SubdomainName> & subdomain_names,
318 : Real scaling_factor)
319 : {
320 72044 : checkVariableNameLength(name);
321 :
322 72044 : if (fe_type.family == SCALAR)
323 : mooseDeprecated(
324 : "The version of Simulation::addSimVariable() with subdomain names can no longer be used "
325 : "with scalar variables since scalar variables cannot be block-restricted. Use the version "
326 : "of Simulation::addSimVariable() without subdomain names instead.");
327 :
328 : #ifdef DEBUG
329 : for (const auto & subdomain_name : subdomain_names)
330 : mooseAssert(subdomain_name != "ANY_BLOCK_ID",
331 : "'ANY_BLOCK_ID' cannot be used for adding field variables in components.");
332 : #endif
333 :
334 72044 : if (_vars.find(name) == _vars.end()) // variable is new
335 : {
336 44140 : VariableInfo vi;
337 : InputParameters & params = vi._params;
338 :
339 44140 : vi._nl = nl;
340 : vi._var_type = "MooseVariable";
341 44140 : params = _thm_factory.getValidParams(vi._var_type);
342 44140 : params.set<std::vector<SubdomainName>>("block") = subdomain_names;
343 :
344 44140 : auto family = AddVariableAction::getNonlinearVariableFamilies();
345 44140 : family = Utility::enum_to_string(fe_type.family);
346 44140 : params.set<MooseEnum>("family") = family;
347 :
348 44140 : auto order = AddVariableAction::getNonlinearVariableOrders();
349 44140 : order = Utility::enum_to_string<Order>(fe_type.order);
350 44140 : params.set<MooseEnum>("order") = order;
351 :
352 44140 : if (nl)
353 25500 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
354 31390 : else if (!MooseUtils::absoluteFuzzyEqual(scaling_factor, 1.0))
355 0 : mooseError("Aux variables cannot be provided a residual scaling factor.");
356 :
357 44140 : _vars[name] = vi;
358 44140 : }
359 : else // variable was previously added
360 : {
361 27904 : VariableInfo & vi = _vars[name];
362 27904 : InputParameters & params = vi._params;
363 :
364 27904 : if (vi._nl != nl)
365 0 : mooseError("The variable '",
366 : name,
367 : "' has already been added in a different system (nonlinear or aux).");
368 :
369 27904 : if (vi._var_type != "MooseVariable")
370 0 : mooseError("The variable '",
371 : name,
372 : "' has already been added with a different type than 'MooseVariable'.");
373 :
374 27904 : auto family = AddVariableAction::getNonlinearVariableFamilies();
375 55808 : family = Utility::enum_to_string(fe_type.family);
376 27904 : if (!params.get<MooseEnum>("family").compareCurrent(family))
377 0 : mooseError("The variable '", name, "' has already been added with a different FE family.");
378 :
379 27904 : auto order = AddVariableAction::getNonlinearVariableOrders();
380 55808 : order = Utility::enum_to_string<Order>(fe_type.order);
381 27904 : if (!params.get<MooseEnum>("order").compareCurrent(order))
382 0 : mooseError("The variable '", name, "' has already been added with a different FE order.");
383 :
384 : // If already block-restricted, extend the block restriction
385 55808 : if (params.isParamValid("block"))
386 : {
387 27904 : auto blocks = params.get<std::vector<SubdomainName>>("block");
388 56160 : for (const auto & subdomain_name : subdomain_names)
389 28256 : if (std::find(blocks.begin(), blocks.end(), subdomain_name) == blocks.end())
390 27346 : blocks.push_back(subdomain_name);
391 27904 : params.set<std::vector<SubdomainName>>("block") = blocks;
392 27904 : }
393 : else
394 0 : params.set<std::vector<SubdomainName>>("block") = subdomain_names;
395 :
396 55808 : if (params.isParamValid("scaling"))
397 6492 : if (!MooseUtils::absoluteFuzzyEqual(params.get<std::vector<Real>>("scaling")[0],
398 : scaling_factor))
399 0 : mooseError(
400 : "The variable '", name, "' has already been added with a different scaling factor.");
401 27904 : }
402 72044 : }
403 :
404 : void
405 0 : Simulation::addSimVariable(bool nl,
406 : const std::string & var_type,
407 : const VariableName & name,
408 : const InputParameters & params)
409 : {
410 0 : checkVariableNameLength(name);
411 :
412 0 : if (_vars.find(name) == _vars.end()) // variable is new
413 : {
414 0 : VariableInfo vi;
415 0 : vi._nl = nl;
416 : vi._var_type = var_type;
417 0 : vi._params = params;
418 :
419 0 : _vars[name] = vi;
420 : }
421 : else // variable was previously added
422 : {
423 0 : VariableInfo & vi = _vars[name];
424 0 : InputParameters & vi_params = vi._params;
425 :
426 0 : if (vi._nl != nl)
427 0 : mooseError("The variable '",
428 : name,
429 : "' has already been added in a different system (nonlinear or aux).");
430 :
431 0 : if (vi._var_type != var_type)
432 0 : mooseError("The variable '",
433 : name,
434 : "' has already been added with a different type than '",
435 : var_type,
436 : "'.");
437 :
438 : // Check that all valid parameters (other than 'block') are consistent
439 0 : for (auto it = params.begin(); it != params.end(); it++)
440 : {
441 0 : const std::string param_name = it->first;
442 0 : if (param_name == "block")
443 : {
444 0 : if (vi_params.isParamValid("block"))
445 : {
446 0 : auto blocks = vi_params.get<std::vector<SubdomainName>>("block");
447 0 : const auto new_blocks = params.get<std::vector<SubdomainName>>("block");
448 0 : for (const auto & subdomain_name : new_blocks)
449 0 : if (std::find(blocks.begin(), blocks.end(), subdomain_name) == blocks.end())
450 0 : blocks.push_back(subdomain_name);
451 0 : vi_params.set<std::vector<SubdomainName>>("block") = blocks;
452 0 : }
453 : else
454 0 : mooseError("The variable '", name, "' was added previously without block restriction.");
455 : }
456 0 : else if (params.isParamValid(param_name))
457 : {
458 0 : if (vi_params.isParamValid(param_name))
459 : {
460 0 : if (params.rawParamVal(param_name) != vi_params.rawParamVal(param_name))
461 0 : mooseError("The variable '",
462 : name,
463 : "' was added previously with a different value for the parameter '",
464 : param_name,
465 : "'.");
466 : }
467 : else
468 0 : mooseError("The variable '",
469 : name,
470 : "' was added previously without the parameter '",
471 : param_name,
472 : "'.");
473 : }
474 : }
475 : }
476 0 : }
477 :
478 : void
479 72549 : Simulation::checkVariableNameLength(const std::string & name) const
480 : {
481 72549 : if (name.size() > THM::MAX_VARIABLE_LENGTH)
482 0 : mooseError(
483 : "Variable name '", name, "' is too long. The limit is ", THM::MAX_VARIABLE_LENGTH, ".");
484 72549 : }
485 :
486 : void
487 0 : Simulation::addControl(const std::string & type, const std::string & name, InputParameters params)
488 : {
489 0 : params.addPrivateParam<FEProblemBase *>("_fe_problem_base", &_fe_problem);
490 0 : std::shared_ptr<Control> control = _thm_factory.create<Control>(type, name, params);
491 0 : _fe_problem.getControlWarehouse().addObject(control);
492 0 : }
493 :
494 : void
495 67289 : Simulation::addSimInitialCondition(const std::string & type,
496 : const std::string & name,
497 : InputParameters params)
498 : {
499 67289 : if (hasInitialConditionsFromFile())
500 : return;
501 :
502 67281 : if (_ics.find(name) == _ics.end())
503 : {
504 67281 : ICInfo ic(type, params);
505 67281 : _ics[name] = ic;
506 : }
507 : else
508 0 : mooseError("Initial condition with name '", name, "' already exists.");
509 : }
510 :
511 : void
512 6458 : Simulation::addConstantIC(const VariableName & var_name,
513 : Real value,
514 : const std::vector<SubdomainName> & block_names)
515 : {
516 6458 : if (hasInitialConditionsFromFile())
517 60 : return;
518 :
519 6398 : std::string blk_str = block_names[0];
520 6416 : for (unsigned int i = 1; i < block_names.size(); i++)
521 36 : blk_str += ":" + block_names[i];
522 :
523 6398 : std::string class_name = "ConstantIC";
524 6398 : InputParameters params = _thm_factory.getValidParams(class_name);
525 6398 : params.set<VariableName>("variable") = var_name;
526 6398 : params.set<Real>("value") = value;
527 6398 : params.set<std::vector<SubdomainName>>("block") = block_names;
528 12796 : addSimInitialCondition(class_name, genName(var_name, blk_str, "ic"), params);
529 6398 : }
530 :
531 : void
532 4734 : Simulation::addFunctionIC(const VariableName & var_name,
533 : const std::string & func_name,
534 : const std::vector<SubdomainName> & block_names)
535 : {
536 4734 : if (hasInitialConditionsFromFile())
537 20 : return;
538 :
539 4714 : std::string blk_str = block_names[0];
540 5453 : for (unsigned int i = 1; i < block_names.size(); i++)
541 1478 : blk_str += ":" + block_names[i];
542 :
543 4714 : std::string class_name = "FunctionIC";
544 4714 : InputParameters params = _thm_factory.getValidParams(class_name);
545 4714 : params.set<VariableName>("variable") = var_name;
546 9428 : params.set<std::vector<SubdomainName>>("block") = block_names;
547 9428 : params.set<FunctionName>("function") = func_name;
548 9428 : addSimInitialCondition(class_name, genName(var_name, blk_str, "ic"), params);
549 4714 : }
550 :
551 : void
552 489 : Simulation::addConstantScalarIC(const VariableName & var_name, Real value)
553 : {
554 489 : if (hasInitialConditionsFromFile())
555 0 : return;
556 :
557 489 : std::string class_name = "ScalarConstantIC";
558 489 : InputParameters params = _thm_factory.getValidParams(class_name);
559 489 : params.set<VariableName>("variable") = var_name;
560 489 : params.set<Real>("value") = value;
561 978 : addSimInitialCondition(class_name, genName(var_name, "ic"), params);
562 489 : }
563 :
564 : void
565 0 : Simulation::addComponentScalarIC(const VariableName & var_name, const std::vector<Real> & value)
566 : {
567 0 : if (hasInitialConditionsFromFile())
568 0 : return;
569 :
570 0 : std::string class_name = "ScalarComponentIC";
571 0 : InputParameters params = _thm_factory.getValidParams(class_name);
572 0 : params.set<VariableName>("variable") = var_name;
573 0 : params.set<std::vector<Real>>("values") = value;
574 0 : addSimInitialCondition(class_name, genName(var_name, "ic"), params);
575 0 : }
576 :
577 : std::vector<VariableName>
578 3592 : Simulation::sortAddedComponentVariables() const
579 : {
580 : // Check that no index in order map is used more than once.
581 : // Also, convert the map to a vector of pairs to be sorted.
582 : std::set<int> indices;
583 : std::vector<std::pair<VariableName, int>> registered_var_index_pairs;
584 32328 : for (const auto & var_and_index : _component_variable_order_map)
585 : {
586 28736 : registered_var_index_pairs.push_back(var_and_index);
587 :
588 28736 : const auto ind = var_and_index.second;
589 28736 : auto insert_return = indices.insert(ind);
590 28736 : if (!insert_return.second)
591 0 : mooseError("The index ", ind, " was used for multiple component variables.");
592 : }
593 :
594 : // Collect all of the added variable names into an unsorted vector.
595 : std::vector<VariableName> vars_unsorted;
596 48237 : for (const auto & var_and_data : _vars)
597 44645 : vars_unsorted.push_back(var_and_data.first);
598 :
599 : // The sorting works as follows. For those variables that are listed in
600 : // _component_variable_order_map, these are ordered before those that are not,
601 : // in the order of their indices in the map. Those not in the map are sorted
602 : // alphabetically.
603 :
604 : // Sort registered_var_index_pairs by value (index)
605 3592 : std::sort(registered_var_index_pairs.begin(),
606 : registered_var_index_pairs.end(),
607 : [](const std::pair<VariableName, int> & a, const std::pair<VariableName, int> & b)
608 71840 : { return a.second < b.second; });
609 :
610 : // Loop over the ordered, registered variable names and add a variable to the
611 : // sorted list if in vars_unsorted. When this happens, delete the element from
612 : // vars_unsorted, leaving only unregistered variable names after the loop.
613 : std::vector<VariableName> vars_sorted;
614 32328 : for (const auto & var_index_pair : registered_var_index_pairs)
615 : {
616 28736 : const auto & var = var_index_pair.first;
617 28736 : if (std::find(vars_unsorted.begin(), vars_unsorted.end(), var) != vars_unsorted.end())
618 : {
619 11388 : vars_sorted.push_back(var);
620 11388 : vars_unsorted.erase(std::remove(vars_unsorted.begin(), vars_unsorted.end(), var),
621 : vars_unsorted.end());
622 : }
623 : }
624 :
625 : // Sort the remaining (unregistered) variables alphabetically and then add
626 : // them to the end of the full list.
627 3592 : std::sort(vars_unsorted.begin(), vars_unsorted.end());
628 3592 : vars_sorted.insert(vars_sorted.end(), vars_unsorted.begin(), vars_unsorted.end());
629 :
630 3592 : return vars_sorted;
631 3592 : }
632 :
633 : void
634 3735 : Simulation::addVariables()
635 : {
636 3735 : TransientBase * trex = dynamic_cast<TransientBase *>(getApp().getExecutioner());
637 3735 : if (trex)
638 : {
639 : Moose::TimeIntegratorType ti_type = trex->getTimeScheme();
640 : // This is only needed for the listed time integrators that are using the original approach to
641 : // explicit integration in MOOSE. Currently, the new time integrators like
642 : // ActuallyExplicitEuler do not need the implicit flag to be set.
643 3663 : if (ti_type == Moose::TI_EXPLICIT_TVD_RK_2 || ti_type == Moose::TI_EXPLICIT_MIDPOINT ||
644 : ti_type == Moose::TI_EXPLICIT_EULER)
645 0 : _implicit_time_integration = false;
646 : }
647 :
648 3735 : if (_components.size() == 0)
649 143 : return;
650 :
651 : // Cache the variables that components request to add
652 18504 : for (auto && comp : _components)
653 14912 : comp->addVariables();
654 :
655 : // Sort the variables for a consistent ordering
656 3592 : const auto var_names = sortAddedComponentVariables();
657 :
658 : // Report the ordering if the executioner is verbose
659 10776 : if (_fe_problem.getParam<MooseEnum>("verbose_setup") != "false")
660 : {
661 0 : std::stringstream ss;
662 0 : ss << "The system ordering of variables added by Components is as follows:\n";
663 0 : for (const auto & var : var_names)
664 0 : ss << " " << var << "\n";
665 0 : mooseInfo(ss.str());
666 0 : }
667 :
668 : // Add the variables to the problem
669 48237 : for (const auto & name : var_names)
670 : {
671 44645 : VariableInfo & vi = _vars[name];
672 :
673 44645 : if (vi._nl)
674 12993 : _fe_problem.addVariable(vi._var_type, name, vi._params);
675 : else
676 31652 : _fe_problem.addAuxVariable(vi._var_type, name, vi._params);
677 : }
678 :
679 3592 : if (hasInitialConditionsFromFile())
680 54 : setupInitialConditionsFromFile();
681 : else
682 3538 : setupInitialConditionObjects();
683 3592 : }
684 :
685 : void
686 54 : Simulation::setupInitialConditionsFromFile()
687 : {
688 108 : const UserObjectName suo_name = genName("thm", "suo");
689 : {
690 54 : const std::string class_name = "SolutionUserObject";
691 54 : InputParameters params = _thm_factory.getValidParams(class_name);
692 108 : params.set<MeshFileName>("mesh") = _thm_pars.get<FileName>("initial_from_file");
693 54 : params.set<std::string>("timestep") = _thm_pars.get<std::string>("initial_from_file_timestep");
694 54 : _fe_problem.addUserObject(class_name, suo_name, params);
695 54 : }
696 :
697 552 : for (auto && v : _vars)
698 : {
699 : const VariableName & var_name = v.first;
700 : const VariableInfo & vi = v.second;
701 :
702 498 : if (vi._var_type == "MooseVariableScalar")
703 : {
704 8 : std::string class_name = "ScalarSolutionIC";
705 8 : InputParameters params = _thm_factory.getValidParams(class_name);
706 8 : params.set<VariableName>("variable") = var_name;
707 8 : params.set<VariableName>("from_variable") = var_name;
708 8 : params.set<UserObjectName>("solution_uo") = suo_name;
709 16 : _fe_problem.addInitialCondition(class_name, genName(var_name, "ic"), params);
710 8 : }
711 : else
712 : {
713 490 : std::string class_name = "SolutionIC";
714 490 : InputParameters params = _thm_factory.getValidParams(class_name);
715 490 : params.set<VariableName>("variable") = var_name;
716 490 : params.set<VariableName>("from_variable") = var_name;
717 490 : params.set<UserObjectName>("solution_uo") = suo_name;
718 980 : if (vi._params.isParamValid("block"))
719 490 : params.set<std::vector<SubdomainName>>("block") =
720 1470 : vi._params.get<std::vector<SubdomainName>>("block");
721 980 : _fe_problem.addInitialCondition(class_name, genName(var_name, "ic"), params);
722 490 : }
723 : }
724 54 : }
725 :
726 : void
727 3538 : Simulation::setupInitialConditionObjects()
728 : {
729 70819 : for (auto && i : _ics)
730 : {
731 67281 : const std::string & name = i.first;
732 : ICInfo & ic = i.second;
733 67281 : _fe_problem.addInitialCondition(ic._type, name, ic._params);
734 : }
735 3538 : }
736 :
737 : void
738 3735 : Simulation::addMooseObjects()
739 : {
740 18647 : for (auto && comp : _components)
741 14912 : comp->addMooseObjects();
742 3735 : }
743 :
744 : void
745 3900 : Simulation::addRelationshipManagers()
746 : {
747 : {
748 3900 : const std::string class_name = "AugmentSparsityBetweenElements";
749 3900 : auto params = _thm_factory.getValidParams(class_name);
750 3900 : params.set<Moose::RelationshipManagerType>("rm_type") =
751 : Moose::RelationshipManagerType::COUPLING | Moose::RelationshipManagerType::ALGEBRAIC |
752 : Moose::RelationshipManagerType::GEOMETRIC;
753 3900 : params.set<std::string>("for_whom") = _fe_problem.name();
754 3900 : params.set<MooseMesh *>("mesh") = &_thm_mesh;
755 3900 : params.set<std::map<dof_id_type, std::vector<dof_id_type>> *>("_elem_map") =
756 3900 : &_sparsity_elem_augmentation;
757 : auto rm =
758 3900 : _thm_factory.create<RelationshipManager>(class_name, "thm:sparsity_btw_elems", params);
759 11700 : if (!_thm_app.addRelationshipManager(rm))
760 0 : _thm_factory.releaseSharedObjects(*rm);
761 3900 : }
762 :
763 19538 : for (auto && comp : _components)
764 15638 : comp->addRelationshipManagers(Moose::RelationshipManagerType::COUPLING |
765 : Moose::RelationshipManagerType::ALGEBRAIC |
766 : Moose::RelationshipManagerType::GEOMETRIC);
767 3900 : }
768 :
769 : void
770 3748 : Simulation::setupCoordinateSystem()
771 : {
772 7496 : MultiMooseEnum coord_types("XYZ RZ RSPHERICAL");
773 : std::vector<SubdomainName> blocks;
774 :
775 19296 : for (auto && comp : _components)
776 : {
777 15548 : if (comp->parent() == nullptr)
778 : {
779 15548 : const auto & subdomains = comp->getSubdomainNames();
780 15548 : const auto & coord_sys = comp->getCoordSysTypes();
781 :
782 23412 : for (unsigned int i = 0; i < subdomains.size(); i++)
783 : {
784 7864 : blocks.push_back(subdomains[i]);
785 : // coord_types.push_back("XYZ");
786 23592 : coord_types.setAdditionalValue(coord_sys[i] == Moose::COORD_RZ ? "RZ" : "XYZ");
787 : }
788 : }
789 : }
790 3748 : _fe_problem.setCoordSystem(blocks, coord_types);
791 :
792 : // RZ geometries are always aligned with x-axis
793 7496 : MooseEnum rz_coord_axis("X=0 Y=1", "X");
794 3748 : _fe_problem.setAxisymmetricCoordAxis(rz_coord_axis);
795 3748 : }
796 :
797 : void
798 3891 : Simulation::setupMesh()
799 : {
800 3891 : if (_components.size() == 0)
801 : return;
802 :
803 3748 : setupCoordinateSystem();
804 : }
805 :
806 : void
807 3729 : Simulation::couplingMatrixIntegrityCheck() const
808 : {
809 3729 : if (!_fe_problem.shouldSolve())
810 : return;
811 :
812 : const TimeIntegrator * ti = nullptr;
813 : const auto & time_integrators =
814 3468 : _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrators();
815 3468 : if (!time_integrators.empty())
816 : ti = time_integrators.front().get();
817 : // Yes, this is horrible. Don't ask why...
818 3412 : if ((dynamic_cast<const ExplicitTimeIntegrator *>(ti) != nullptr) ||
819 3210 : (dynamic_cast<const ExplicitEuler *>(ti) != nullptr) ||
820 6678 : (dynamic_cast<const ExplicitRK2 *>(ti) != nullptr) ||
821 3210 : (dynamic_cast<const ExplicitTVDRK2 *>(ti) != nullptr))
822 : return;
823 :
824 3210 : const CouplingMatrix * cm = _fe_problem.couplingMatrix(/*nl_sys_num=*/0);
825 3210 : if (cm == nullptr)
826 0 : mooseError("Coupling matrix does not exists. Something really bad happened.");
827 :
828 : bool full = true;
829 15519 : for (unsigned int i = 0; i < cm->size(); i++)
830 82230 : for (unsigned int j = 0; j < cm->size(); j++)
831 : full &= (*cm)(i, j);
832 :
833 3210 : if (!full)
834 0 : mooseError(
835 : "Single matrix preconditioning with full coupling is required to run. Please, check that "
836 : "your input file has the following preconditioning block:\n\n"
837 : "[Preconditioning]\n"
838 : " [pc]\n"
839 : " type = SMP\n"
840 : " full = true\n"
841 : " []\n"
842 : "[].\n");
843 : }
844 :
845 : void
846 3889 : Simulation::integrityCheck() const
847 : {
848 3889 : if (_components.size() == 0)
849 223 : return;
850 :
851 3746 : if (_check_jacobian)
852 : return;
853 :
854 : // go over components and put flow channels into one "bucket"
855 : std::vector<Component *> flow_channels;
856 18936 : for (auto && comp : _components)
857 : {
858 15270 : auto flow_channel = dynamic_cast<FlowChannelBase *>(comp.get());
859 15270 : if (flow_channel != nullptr)
860 4050 : flow_channels.push_back(flow_channel);
861 : }
862 :
863 : // initialize number of connected flow channel inlets and outlets to zero
864 : std::map<std::string, unsigned int> flow_channel_inlets;
865 : std::map<std::string, unsigned int> flow_channel_outlets;
866 7716 : for (auto && comp : flow_channels)
867 : {
868 4050 : flow_channel_inlets[comp->name()] = 0;
869 4050 : flow_channel_outlets[comp->name()] = 0;
870 : }
871 :
872 : // mark connections of any Component1DConnection components
873 18936 : for (const auto & comp : _components)
874 : {
875 15270 : auto pc_comp = dynamic_cast<Component1DConnection *>(comp.get());
876 15270 : if (pc_comp != nullptr)
877 : {
878 14454 : for (const auto & connection : pc_comp->getConnections())
879 : {
880 8098 : if (connection._end_type == Component1DConnection::IN)
881 4048 : flow_channel_inlets[connection._component_name]++;
882 4050 : else if (connection._end_type == Component1DConnection::OUT)
883 4050 : flow_channel_outlets[connection._component_name]++;
884 : }
885 : }
886 : }
887 :
888 : // finally, check that each flow channel has exactly one input and one output
889 7716 : for (auto && comp : flow_channels)
890 : {
891 4050 : if (flow_channel_inlets[comp->name()] == 0)
892 6 : logError("Component '", comp->name(), "' does not have connected inlet.");
893 4044 : else if (flow_channel_inlets[comp->name()] > 1)
894 2 : logError("Multiple inlets specified for component '", comp->name(), "'.");
895 :
896 4050 : if (flow_channel_outlets[comp->name()] == 0)
897 0 : logError("Component '", comp->name(), "' does not have connected outlet.");
898 4050 : else if (flow_channel_outlets[comp->name()] > 1)
899 0 : logError("Multiple outlets specified for component '", comp->name(), "'.");
900 : }
901 :
902 : // let components check themselves
903 18936 : for (auto && comp : _components)
904 15270 : comp->executeCheck();
905 :
906 3666 : _log.emitLoggedWarnings();
907 3662 : _log.emitLoggedErrors();
908 : }
909 :
910 : void
911 3729 : Simulation::controlDataIntegrityCheck()
912 : {
913 3729 : if (_check_jacobian)
914 : return;
915 :
916 : // check that control data are consistent
917 14095 : for (auto && i : _control_data)
918 : {
919 10446 : if (!i.second->getDeclared())
920 : logError("Control data '",
921 2 : i.first,
922 : "' was requested, but was not declared by any active control object.");
923 : }
924 :
925 3649 : _log.emitLoggedErrors();
926 :
927 7294 : auto & ctrl_wh = _fe_problem.getControlWarehouse()[EXEC_TIMESTEP_BEGIN];
928 :
929 : // initialize THM control objects
930 14655 : for (auto && i : ctrl_wh.getObjects())
931 : {
932 11008 : THMControl * ctrl = dynamic_cast<THMControl *>(i.get());
933 11008 : if (ctrl != nullptr)
934 11008 : ctrl->init();
935 : }
936 :
937 14655 : for (auto && i : ctrl_wh.getObjects())
938 : {
939 11008 : THMControl * ctrl = dynamic_cast<THMControl *>(i.get());
940 : // if it is a THM control
941 11008 : if (ctrl != nullptr)
942 : {
943 : // get its dependencies on control data
944 : auto & cd_deps = ctrl->getControlDataDependencies();
945 11417 : for (auto && cd_name : cd_deps)
946 : {
947 409 : ControlDataValue * cdv = _control_data[cd_name];
948 : // find out which control object built the control data
949 409 : std::string dep_name = cdv->getControl()->name();
950 : auto & deps = ctrl->getDependencies();
951 : // and if it is not in its dependency list, add it
952 409 : auto it = std::find(deps.begin(), deps.end(), dep_name);
953 409 : if (it == deps.end())
954 409 : deps.push_back(dep_name);
955 : }
956 : }
957 : }
958 :
959 : // Find all `TerminateControl`s and all their dependencies. Then add those
960 : // objects into TIMESTEP_END control warehouse
961 : MooseObjectWarehouse<Control> & ctrl_wh_tse =
962 7294 : _fe_problem.getControlWarehouse()[EXEC_TIMESTEP_END];
963 14655 : for (auto && i : ctrl_wh.getObjects())
964 : {
965 11008 : if (TerminateControl * ctrl = dynamic_cast<TerminateControl *>(i.get()))
966 : {
967 : std::list<const THMControl *> l;
968 10 : l.push_back(ctrl);
969 40 : while (l.size() > 0)
970 : {
971 30 : const THMControl * ctrl = l.front();
972 : auto & cd_deps = ctrl->getControlDataDependencies();
973 50 : for (auto && cd_name : cd_deps)
974 : {
975 20 : ControlDataValue * cdv = _control_data[cd_name];
976 20 : l.push_back(cdv->getControl());
977 : }
978 60 : ctrl_wh_tse.addObject(ctrl_wh.getObject(ctrl->name()));
979 : l.pop_front();
980 : }
981 : }
982 : }
983 : }
984 :
985 : void
986 0 : Simulation::run()
987 : {
988 0 : }
989 :
990 : void
991 15640 : Simulation::addComponent(const std::string & type, const std::string & name, InputParameters params)
992 : {
993 15640 : std::shared_ptr<Component> comp = _thm_factory.create<Component>(type, name, params);
994 15638 : if (_comp_by_name.find(name) == _comp_by_name.end())
995 15638 : _comp_by_name[name] = comp;
996 : else
997 : logError("Component with name '", name, "' already exists");
998 15638 : _components.push_back(comp);
999 15638 : }
1000 :
1001 : bool
1002 10822 : Simulation::hasComponent(const std::string & name) const
1003 : {
1004 : auto it = _comp_by_name.find(name);
1005 10822 : return (it != _comp_by_name.end());
1006 : }
1007 :
1008 : void
1009 2628 : Simulation::addClosures(const std::string & type, const std::string & name, InputParameters params)
1010 : {
1011 2628 : std::shared_ptr<ClosuresBase> obj_ptr = _thm_factory.create<ClosuresBase>(type, name, params);
1012 2628 : if (_closures_by_name.find(name) == _closures_by_name.end())
1013 2628 : _closures_by_name[name] = obj_ptr;
1014 : else
1015 : logError("A closures object with the name '", name, "' already exists.");
1016 2628 : }
1017 :
1018 : bool
1019 0 : Simulation::hasClosures(const std::string & name) const
1020 : {
1021 0 : return _closures_by_name.find(name) != _closures_by_name.end();
1022 : }
1023 :
1024 : std::shared_ptr<ClosuresBase>
1025 4179 : Simulation::getClosures(const std::string & name) const
1026 : {
1027 : auto it = _closures_by_name.find(name);
1028 4179 : if (it != _closures_by_name.end())
1029 4179 : return it->second;
1030 : else
1031 0 : mooseError("The requested closures object '", name, "' does not exist.");
1032 : }
1033 :
1034 : void
1035 3640 : Simulation::addFileOutputter(const std::string & name)
1036 : {
1037 3640 : _outputters_all.push_back(name);
1038 3640 : _outputters_file.push_back(name);
1039 3640 : }
1040 :
1041 : void
1042 3729 : Simulation::addScreenOutputter(const std::string & name)
1043 : {
1044 3729 : _outputters_all.push_back(name);
1045 3729 : _outputters_screen.push_back(name);
1046 3729 : }
1047 :
1048 : std::vector<OutputName>
1049 0 : Simulation::getOutputsVector(const std::string & key) const
1050 : {
1051 0 : std::string key_lowercase = key;
1052 : std::transform(key_lowercase.begin(), key_lowercase.end(), key_lowercase.begin(), ::tolower);
1053 :
1054 : std::vector<OutputName> outputs;
1055 0 : if (key_lowercase == "none")
1056 0 : outputs.push_back("none"); // provide non-existent name, so it does not get printed out
1057 0 : else if (key_lowercase == "screen")
1058 0 : outputs = _outputters_screen;
1059 0 : else if (key_lowercase == "file")
1060 0 : outputs = _outputters_file;
1061 0 : else if (key_lowercase == "both")
1062 0 : outputs = _outputters_all;
1063 : else
1064 0 : mooseError("The outputs vector key '" + key_lowercase + "' is invalid");
1065 :
1066 0 : return outputs;
1067 0 : }
1068 :
1069 : bool
1070 87145 : Simulation::hasInitialConditionsFromFile() const
1071 : {
1072 174290 : return _thm_pars.isParamValid("initial_from_file");
1073 : }
1074 :
1075 : void
1076 31899 : Simulation::advanceState()
1077 : {
1078 108118 : for (auto && i : _control_data)
1079 76219 : i.second->copyValuesBack();
1080 31899 : }
|
