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 46288 : Simulation::setComponentVariableOrder(const VariableName & var, int index)
41 : {
42 46288 : _component_variable_order_map[var] = index;
43 46288 : }
44 :
45 2117 : Simulation::Simulation(FEProblemBase & fe_problem, const InputParameters & pars)
46 : : ParallelObject(fe_problem.comm()),
47 2117 : LoggingInterface(_log),
48 2117 : _thm_mesh(*static_cast<THMMesh *>(pars.get<MooseMesh *>("mesh"))),
49 2117 : _fe_problem(fe_problem),
50 2117 : _thm_app(static_cast<ThermalHydraulicsApp &>(*pars.get<MooseApp *>(MooseBase::app_param))),
51 2117 : _thm_factory(_thm_app.getFactory()),
52 2117 : _thm_pars(pars),
53 : _flow_fe_type(FEType(CONSTANT, MONOMIAL)),
54 2117 : _implicit_time_integration(true),
55 2117 : _check_jacobian(false),
56 2117 : _output_vector_velocity(true),
57 4234 : _zero(0)
58 : {
59 2117 : bool second_order_mesh = pars.get<bool>("2nd_order_mesh");
60 2117 : HeatConductionModel::_fe_type =
61 2117 : second_order_mesh ? FEType(SECOND, LAGRANGE) : FEType(FIRST, LAGRANGE);
62 2117 : }
63 :
64 1951 : Simulation::~Simulation()
65 : {
66 7115 : for (auto && k : _control_data)
67 5164 : delete k.second;
68 5853 : }
69 :
70 : void
71 9932 : 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 9932 : if (it == _sparsity_elem_augmentation.end())
75 8158 : it = _sparsity_elem_augmentation.insert(_sparsity_elem_augmentation.begin(),
76 8158 : {elem_id1, std::vector<dof_id_type>()});
77 9932 : it->second.push_back(elem_id2);
78 :
79 : it = _sparsity_elem_augmentation.find(elem_id2);
80 9932 : if (it == _sparsity_elem_augmentation.end())
81 5937 : it = _sparsity_elem_augmentation.insert(_sparsity_elem_augmentation.begin(),
82 5937 : {elem_id2, std::vector<dof_id_type>()});
83 9932 : it->second.push_back(elem_id1);
84 9932 : }
85 :
86 : void
87 2113 : Simulation::buildMesh()
88 : {
89 2113 : if (_components.size() == 0)
90 : return;
91 :
92 : // build mesh
93 10336 : for (auto && comp : _components)
94 8295 : comp->executeSetupMesh();
95 : }
96 :
97 : void
98 1777 : Simulation::setupQuadrature()
99 : {
100 1777 : 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 8600 : for (auto && comp : _components)
108 : {
109 6895 : auto flow_channel = dynamic_cast<FlowChannelBase *>(comp.get());
110 6895 : if (flow_channel != nullptr)
111 1732 : n_flow_channels++;
112 :
113 6895 : auto hs_interface = dynamic_cast<HeatStructureInterface *>(comp.get());
114 6895 : if (hs_interface)
115 986 : n_heat_structures++;
116 : }
117 :
118 1705 : if (n_flow_channels > 0)
119 : {
120 : const FEType & fe_type = getFlowFEType();
121 1091 : if (fe_type.default_quadrature_order() > order)
122 : order = fe_type.default_quadrature_order();
123 : }
124 1705 : if (n_heat_structures > 0)
125 : {
126 : const FEType & fe_type = HeatConductionModel::feType();
127 767 : if (fe_type.default_quadrature_order() > order)
128 : order = fe_type.default_quadrature_order();
129 : }
130 :
131 1705 : _fe_problem.createQRules(QGAUSS, order, order, order);
132 : }
133 :
134 : void
135 2113 : Simulation::initSimulation()
136 : {
137 : // sort the components using dependency resolver
138 2113 : DependencyResolver<std::shared_ptr<Component>> dependency_resolver;
139 10408 : for (const auto & comp : _components)
140 : {
141 8295 : dependency_resolver.addNode(comp);
142 13922 : for (const auto & dep : comp->getDependencies())
143 5627 : if (hasComponent(dep))
144 5576 : dependency_resolver.addEdge(_comp_by_name[dep], comp);
145 : }
146 :
147 2113 : _components = dependency_resolver.dfs();
148 2113 : }
149 :
150 : void
151 2106 : Simulation::initComponents()
152 : {
153 : // initialize components
154 10323 : for (auto && comp : _components)
155 8219 : comp->executeInit();
156 :
157 : // perform secondary initialization, which relies on init() being called
158 : // already for all components
159 10321 : for (auto && comp : _components)
160 8217 : comp->executeInitSecondary();
161 2104 : }
162 :
163 : void
164 2104 : Simulation::identifyLoops()
165 : {
166 : // loop over junctions and boundaries (non-geometrical components)
167 10321 : for (const auto & component : _components)
168 : {
169 : const auto flow_connection =
170 8217 : MooseSharedNamespace::dynamic_pointer_cast<Component1DConnection>(component);
171 8217 : if (flow_connection)
172 : {
173 : // create vector of names of this component and its connected flow channels, and then sort
174 : // them
175 3390 : std::vector<std::string> names = flow_connection->getConnectedComponentNames();
176 3390 : names.push_back(component->name());
177 3390 : std::sort(names.begin(), names.end());
178 :
179 : // pick first name alphabetically to be the proposed loop name
180 3390 : std::string proposed_loop_name = names[0];
181 :
182 11012 : for (const std::string & name : names)
183 : {
184 : // if the name is not yet in the map
185 7622 : 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 5520 : _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 2102 : 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 2102 : if (proposed_loop_name > current_loop_name)
197 : {
198 7823 : for (auto && entry : _component_name_to_loop_name)
199 6318 : 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 597 : else if (proposed_loop_name < current_loop_name)
206 : {
207 1253 : for (auto && entry : _component_name_to_loop_name)
208 1057 : if (entry.second == current_loop_name)
209 : entry.second = proposed_loop_name;
210 : }
211 : }
212 : }
213 3390 : }
214 : }
215 :
216 : // get the list of loops
217 : std::vector<std::string> loops;
218 7624 : for (const auto & entry : _component_name_to_loop_name)
219 5520 : if (std::find(loops.begin(), loops.end(), entry.second) == loops.end())
220 1437 : loops.push_back(entry.second);
221 :
222 : // fill the map of loop name to model ID
223 3541 : 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 1890 : for (const auto & component : _components)
229 : {
230 : const auto flow_chan_base_component =
231 1884 : MooseSharedNamespace::dynamic_pointer_cast<FlowChannelBase>(component);
232 1884 : if (flow_chan_base_component && (_component_name_to_loop_name[component->name()] == loop))
233 : {
234 1431 : 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 1431 : _loop_name_to_model_id[loop] = model_id;
242 : else
243 : logError("No FlowChannelBase-derived components were found in loop '", loop, "'");
244 : }
245 2104 : }
246 :
247 : void
248 12 : Simulation::printComponentLoops() const
249 : {
250 : // get the list of loops
251 : std::vector<std::string> loops;
252 114 : for (auto && entry : _component_name_to_loop_name)
253 102 : if (std::find(loops.begin(), loops.end(), entry.second) == loops.end())
254 24 : loops.push_back(entry.second);
255 :
256 : // for each loop
257 : Moose::out << "\nListing of component loops:" << std::endl;
258 36 : for (unsigned int i = 0; i < loops.size(); i++)
259 : {
260 24 : Moose::out << "\n Loop " << i + 1 << ":" << std::endl;
261 :
262 : // print out each component in the loop
263 228 : for (auto && entry : _component_name_to_loop_name)
264 204 : if (entry.second == loops[i])
265 : Moose::out << " " << entry.first << std::endl;
266 : }
267 : Moose::out << std::endl;
268 12 : }
269 :
270 : void
271 230 : Simulation::addSimVariable(bool nl, const VariableName & name, FEType fe_type, Real scaling_factor)
272 : {
273 230 : checkVariableNameLength(name);
274 :
275 230 : if (fe_type.family != SCALAR)
276 0 : mooseError("This method should only be used for scalar variables.");
277 :
278 230 : if (_vars.find(name) == _vars.end()) // variable is new
279 : {
280 230 : VariableInfo vi;
281 : InputParameters & params = vi._params;
282 :
283 230 : vi._nl = nl;
284 : vi._var_type = "MooseVariableScalar";
285 230 : params = _thm_factory.getValidParams(vi._var_type);
286 :
287 230 : auto family = AddVariableAction::getNonlinearVariableFamilies();
288 230 : family = Utility::enum_to_string(fe_type.family);
289 230 : params.set<MooseEnum>("family") = family;
290 :
291 230 : auto order = AddVariableAction::getNonlinearVariableOrders();
292 230 : order = Utility::enum_to_string<Order>(fe_type.order);
293 230 : params.set<MooseEnum>("order") = order;
294 :
295 230 : if (nl)
296 240 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
297 110 : else if (!MooseUtils::absoluteFuzzyEqual(scaling_factor, 1.0))
298 0 : mooseError("Aux variables cannot be provided a residual scaling factor.");
299 :
300 230 : _vars[name] = vi;
301 230 : }
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 230 : }
312 :
313 : void
314 36580 : 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 36580 : checkVariableNameLength(name);
321 :
322 36580 : if (fe_type.family == SCALAR)
323 0 : 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 36580 : if (_vars.find(name) == _vars.end()) // variable is new
335 : {
336 23048 : VariableInfo vi;
337 : InputParameters & params = vi._params;
338 :
339 23048 : vi._nl = nl;
340 : vi._var_type = "MooseVariable";
341 23048 : params = _thm_factory.getValidParams(vi._var_type);
342 23048 : params.set<std::vector<SubdomainName>>("block") = subdomain_names;
343 :
344 23048 : auto family = AddVariableAction::getNonlinearVariableFamilies();
345 23048 : family = Utility::enum_to_string(fe_type.family);
346 23048 : params.set<MooseEnum>("family") = family;
347 :
348 23048 : auto order = AddVariableAction::getNonlinearVariableOrders();
349 23048 : order = Utility::enum_to_string<Order>(fe_type.order);
350 23048 : params.set<MooseEnum>("order") = order;
351 :
352 23048 : if (nl)
353 13514 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
354 16291 : else if (!MooseUtils::absoluteFuzzyEqual(scaling_factor, 1.0))
355 0 : mooseError("Aux variables cannot be provided a residual scaling factor.");
356 :
357 23048 : _vars[name] = vi;
358 23048 : }
359 : else // variable was previously added
360 : {
361 13532 : VariableInfo & vi = _vars[name];
362 13532 : InputParameters & params = vi._params;
363 :
364 13532 : 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 13532 : 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 13532 : auto family = AddVariableAction::getNonlinearVariableFamilies();
375 27064 : family = Utility::enum_to_string(fe_type.family);
376 13532 : if (!params.get<MooseEnum>("family").compareCurrent(family))
377 0 : mooseError("The variable '", name, "' has already been added with a different FE family.");
378 :
379 13532 : auto order = AddVariableAction::getNonlinearVariableOrders();
380 27064 : order = Utility::enum_to_string<Order>(fe_type.order);
381 13532 : 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 27064 : if (params.isParamValid("block"))
386 : {
387 13532 : auto blocks = params.get<std::vector<SubdomainName>>("block");
388 27244 : for (const auto & subdomain_name : subdomain_names)
389 13712 : if (std::find(blocks.begin(), blocks.end(), subdomain_name) == blocks.end())
390 13260 : blocks.push_back(subdomain_name);
391 13532 : params.set<std::vector<SubdomainName>>("block") = blocks;
392 13532 : }
393 : else
394 0 : params.set<std::vector<SubdomainName>>("block") = subdomain_names;
395 :
396 27064 : if (params.isParamValid("scaling"))
397 3136 : 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 13532 : }
402 36580 : }
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 36810 : Simulation::checkVariableNameLength(const std::string & name) const
480 : {
481 36810 : 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 36810 : }
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 33888 : Simulation::addSimInitialCondition(const std::string & type,
496 : const std::string & name,
497 : InputParameters params)
498 : {
499 33888 : if (hasInitialConditionsFromFile())
500 : return;
501 :
502 33883 : if (_ics.find(name) == _ics.end())
503 : {
504 33883 : ICInfo ic(type, params);
505 33883 : _ics[name] = ic;
506 : }
507 : else
508 0 : mooseError("Initial condition with name '", name, "' already exists.");
509 : }
510 :
511 : void
512 3279 : Simulation::addConstantIC(const VariableName & var_name,
513 : Real value,
514 : const std::vector<SubdomainName> & block_names)
515 : {
516 3279 : if (hasInitialConditionsFromFile())
517 46 : return;
518 :
519 3233 : std::string blk_str = block_names[0];
520 3247 : for (unsigned int i = 1; i < block_names.size(); i++)
521 28 : blk_str += ":" + block_names[i];
522 :
523 3233 : std::string class_name = "ConstantIC";
524 3233 : InputParameters params = _thm_factory.getValidParams(class_name);
525 3233 : params.set<VariableName>("variable") = var_name;
526 3233 : params.set<Real>("value") = value;
527 3233 : params.set<std::vector<SubdomainName>>("block") = block_names;
528 6466 : addSimInitialCondition(class_name, genName(var_name, blk_str, "ic"), params);
529 3233 : }
530 :
531 : void
532 2344 : Simulation::addFunctionIC(const VariableName & var_name,
533 : const std::string & func_name,
534 : const std::vector<SubdomainName> & block_names)
535 : {
536 2344 : if (hasInitialConditionsFromFile())
537 15 : return;
538 :
539 2329 : std::string blk_str = block_names[0];
540 2756 : for (unsigned int i = 1; i < block_names.size(); i++)
541 854 : blk_str += ":" + block_names[i];
542 :
543 2329 : std::string class_name = "FunctionIC";
544 2329 : InputParameters params = _thm_factory.getValidParams(class_name);
545 2329 : params.set<VariableName>("variable") = var_name;
546 4658 : params.set<std::vector<SubdomainName>>("block") = block_names;
547 4658 : params.set<FunctionName>("function") = func_name;
548 4658 : addSimInitialCondition(class_name, genName(var_name, blk_str, "ic"), params);
549 2329 : }
550 :
551 : void
552 218 : Simulation::addConstantScalarIC(const VariableName & var_name, Real value)
553 : {
554 218 : if (hasInitialConditionsFromFile())
555 0 : return;
556 :
557 218 : std::string class_name = "ScalarConstantIC";
558 218 : InputParameters params = _thm_factory.getValidParams(class_name);
559 218 : params.set<VariableName>("variable") = var_name;
560 218 : params.set<Real>("value") = value;
561 436 : addSimInitialCondition(class_name, genName(var_name, "ic"), params);
562 218 : }
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 1894 : 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 17046 : for (const auto & var_and_index : _component_variable_order_map)
585 : {
586 15152 : registered_var_index_pairs.push_back(var_and_index);
587 :
588 15152 : const auto ind = var_and_index.second;
589 15152 : auto insert_return = indices.insert(ind);
590 15152 : 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 25172 : for (const auto & var_and_data : _vars)
597 23278 : 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 1894 : 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 37880 : { 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 17046 : for (const auto & var_index_pair : registered_var_index_pairs)
615 : {
616 15152 : const auto & var = var_index_pair.first;
617 15152 : if (std::find(vars_unsorted.begin(), vars_unsorted.end(), var) != vars_unsorted.end())
618 : {
619 5997 : vars_sorted.push_back(var);
620 5997 : 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 1894 : std::sort(vars_unsorted.begin(), vars_unsorted.end());
628 1894 : vars_sorted.insert(vars_sorted.end(), vars_unsorted.begin(), vars_unsorted.end());
629 :
630 1894 : return vars_sorted;
631 1894 : }
632 :
633 : void
634 1966 : Simulation::addVariables()
635 : {
636 1966 : TransientBase * trex = dynamic_cast<TransientBase *>(getApp().getExecutioner());
637 1966 : 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 1922 : 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 1966 : if (_components.size() == 0)
649 72 : return;
650 :
651 : // Cache the variables that components request to add
652 9535 : for (auto && comp : _components)
653 7641 : comp->addVariables();
654 :
655 : // Sort the variables for a consistent ordering
656 1894 : const auto var_names = sortAddedComponentVariables();
657 :
658 : // Report the ordering if the executioner is verbose
659 5682 : 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 25172 : for (const auto & name : var_names)
670 : {
671 23278 : VariableInfo & vi = _vars[name];
672 :
673 23278 : if (vi._nl)
674 6877 : _fe_problem.addVariable(vi._var_type, name, vi._params);
675 : else
676 16401 : _fe_problem.addAuxVariable(vi._var_type, name, vi._params);
677 : }
678 :
679 1894 : if (hasInitialConditionsFromFile())
680 41 : setupInitialConditionsFromFile();
681 : else
682 1853 : setupInitialConditionObjects();
683 1894 : }
684 :
685 : void
686 41 : Simulation::setupInitialConditionsFromFile()
687 : {
688 82 : const UserObjectName suo_name = genName("thm", "suo");
689 : {
690 41 : const std::string class_name = "SolutionUserObject";
691 41 : InputParameters params = _thm_factory.getValidParams(class_name);
692 82 : params.set<MeshFileName>("mesh") = _thm_pars.get<FileName>("initial_from_file");
693 41 : params.set<std::string>("timestep") = _thm_pars.get<std::string>("initial_from_file_timestep");
694 41 : _fe_problem.addUserObject(class_name, suo_name, params);
695 41 : }
696 :
697 421 : for (auto && v : _vars)
698 : {
699 : const VariableName & var_name = v.first;
700 : const VariableInfo & vi = v.second;
701 :
702 380 : if (vi._var_type == "MooseVariableScalar")
703 : {
704 6 : std::string class_name = "ScalarSolutionIC";
705 6 : InputParameters params = _thm_factory.getValidParams(class_name);
706 6 : params.set<VariableName>("variable") = var_name;
707 6 : params.set<VariableName>("from_variable") = var_name;
708 6 : params.set<UserObjectName>("solution_uo") = suo_name;
709 12 : _fe_problem.addInitialCondition(class_name, genName(var_name, "ic"), params);
710 6 : }
711 : else
712 : {
713 374 : std::string class_name = "SolutionIC";
714 374 : InputParameters params = _thm_factory.getValidParams(class_name);
715 374 : params.set<VariableName>("variable") = var_name;
716 374 : params.set<VariableName>("from_variable") = var_name;
717 374 : params.set<UserObjectName>("solution_uo") = suo_name;
718 748 : if (vi._params.isParamValid("block"))
719 374 : params.set<std::vector<SubdomainName>>("block") =
720 1122 : vi._params.get<std::vector<SubdomainName>>("block");
721 748 : _fe_problem.addInitialCondition(class_name, genName(var_name, "ic"), params);
722 374 : }
723 : }
724 41 : }
725 :
726 : void
727 1853 : Simulation::setupInitialConditionObjects()
728 : {
729 35736 : for (auto && i : _ics)
730 : {
731 33883 : const std::string & name = i.first;
732 : ICInfo & ic = i.second;
733 33883 : _fe_problem.addInitialCondition(ic._type, name, ic._params);
734 : }
735 1853 : }
736 :
737 : void
738 1966 : Simulation::addMooseObjects()
739 : {
740 9607 : for (auto && comp : _components)
741 7641 : comp->addMooseObjects();
742 1966 : }
743 :
744 : void
745 2113 : Simulation::addRelationshipManagers()
746 : {
747 : {
748 2113 : const std::string class_name = "AugmentSparsityBetweenElements";
749 2113 : auto params = _thm_factory.getValidParams(class_name);
750 2113 : params.set<Moose::RelationshipManagerType>("rm_type") =
751 : Moose::RelationshipManagerType::COUPLING | Moose::RelationshipManagerType::ALGEBRAIC |
752 : Moose::RelationshipManagerType::GEOMETRIC;
753 2113 : params.set<std::string>("for_whom") = _fe_problem.name();
754 2113 : params.set<MooseMesh *>("mesh") = &_thm_mesh;
755 2113 : params.set<std::map<dof_id_type, std::vector<dof_id_type>> *>("_elem_map") =
756 2113 : &_sparsity_elem_augmentation;
757 : auto rm =
758 2113 : _thm_factory.create<RelationshipManager>(class_name, "thm:sparsity_btw_elems", params);
759 6339 : if (!_thm_app.addRelationshipManager(rm))
760 0 : _thm_factory.releaseSharedObjects(*rm);
761 2113 : }
762 :
763 10408 : for (auto && comp : _components)
764 8295 : comp->addRelationshipManagers(Moose::RelationshipManagerType::COUPLING |
765 : Moose::RelationshipManagerType::ALGEBRAIC |
766 : Moose::RelationshipManagerType::GEOMETRIC);
767 2113 : }
768 :
769 : void
770 2034 : Simulation::setupCoordinateSystem()
771 : {
772 4068 : MultiMooseEnum coord_types("XYZ RZ RSPHERICAL");
773 : std::vector<SubdomainName> blocks;
774 :
775 10259 : for (auto && comp : _components)
776 : {
777 8225 : if (comp->parent() == nullptr)
778 : {
779 8225 : const auto & subdomains = comp->getSubdomainNames();
780 8225 : const auto & coord_sys = comp->getCoordSysTypes();
781 :
782 12448 : for (unsigned int i = 0; i < subdomains.size(); i++)
783 : {
784 4223 : blocks.push_back(subdomains[i]);
785 : // coord_types.push_back("XYZ");
786 12669 : coord_types.setAdditionalValue(coord_sys[i] == Moose::COORD_RZ ? "RZ" : "XYZ");
787 : }
788 : }
789 : }
790 2034 : _fe_problem.setCoordSystem(blocks, coord_types);
791 :
792 : // RZ geometries are always aligned with x-axis
793 4068 : MooseEnum rz_coord_axis("X=0 Y=1", "X");
794 2034 : _fe_problem.setAxisymmetricCoordAxis(rz_coord_axis);
795 2034 : }
796 :
797 : void
798 2106 : Simulation::setupMesh()
799 : {
800 2106 : if (_components.size() == 0)
801 : return;
802 :
803 2034 : setupCoordinateSystem();
804 : }
805 :
806 : void
807 1960 : Simulation::couplingMatrixIntegrityCheck() const
808 : {
809 1960 : if (!_fe_problem.shouldSolve())
810 : return;
811 :
812 : const TimeIntegrator * ti = nullptr;
813 : const auto & time_integrators =
814 1787 : _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrators();
815 1787 : if (!time_integrators.empty())
816 : ti = time_integrators.front().get();
817 : // Yes, this is horrible. Don't ask why...
818 1755 : if ((dynamic_cast<const ExplicitTimeIntegrator *>(ti) != nullptr) ||
819 1645 : (dynamic_cast<const ExplicitEuler *>(ti) != nullptr) ||
820 3432 : (dynamic_cast<const ExplicitRK2 *>(ti) != nullptr) ||
821 1645 : (dynamic_cast<const ExplicitTVDRK2 *>(ti) != nullptr))
822 : return;
823 :
824 1645 : const CouplingMatrix * cm = _fe_problem.couplingMatrix(/*nl_sys_num=*/0);
825 1645 : if (cm == nullptr)
826 0 : mooseError("Coupling matrix does not exists. Something really bad happened.");
827 :
828 : bool full = true;
829 8003 : for (unsigned int i = 0; i < cm->size(); i++)
830 42498 : for (unsigned int j = 0; j < cm->size(); j++)
831 : full &= (*cm)(i, j);
832 :
833 1645 : 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 2104 : Simulation::integrityCheck() const
847 : {
848 2104 : if (_components.size() == 0)
849 113 : return;
850 :
851 2032 : if (_check_jacobian)
852 : return;
853 :
854 : // go over components and put flow channels into one "bucket"
855 : std::vector<Component *> flow_channels;
856 10073 : for (auto && comp : _components)
857 : {
858 8082 : auto flow_channel = dynamic_cast<FlowChannelBase *>(comp.get());
859 8082 : if (flow_channel != nullptr)
860 2088 : 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 4079 : for (auto && comp : flow_channels)
867 : {
868 2088 : flow_channel_inlets[comp->name()] = 0;
869 2088 : flow_channel_outlets[comp->name()] = 0;
870 : }
871 :
872 : // mark connections of any Component1DConnection components
873 10073 : for (const auto & comp : _components)
874 : {
875 8082 : auto pc_comp = dynamic_cast<Component1DConnection *>(comp.get());
876 8082 : if (pc_comp != nullptr)
877 : {
878 7519 : for (const auto & connection : pc_comp->getConnections())
879 : {
880 4174 : if (connection._end_type == Component1DConnection::IN)
881 2086 : flow_channel_inlets[connection._component_name]++;
882 2088 : else if (connection._end_type == Component1DConnection::OUT)
883 2088 : 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 4079 : for (auto && comp : flow_channels)
890 : {
891 2088 : if (flow_channel_inlets[comp->name()] == 0)
892 6 : logError("Component '", comp->name(), "' does not have connected inlet.");
893 2082 : else if (flow_channel_inlets[comp->name()] > 1)
894 2 : logError("Multiple inlets specified for component '", comp->name(), "'.");
895 :
896 2088 : if (flow_channel_outlets[comp->name()] == 0)
897 0 : logError("Component '", comp->name(), "' does not have connected outlet.");
898 2088 : 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 10073 : for (auto && comp : _components)
904 8082 : comp->executeCheck();
905 :
906 1991 : _log.emitLoggedWarnings();
907 1987 : _log.emitLoggedErrors();
908 1853 : }
909 :
910 : void
911 1960 : Simulation::controlDataIntegrityCheck()
912 : {
913 1960 : if (_check_jacobian)
914 : return;
915 :
916 : // check that control data are consistent
917 7097 : for (auto && i : _control_data)
918 : {
919 5178 : 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 1919 : _log.emitLoggedErrors();
926 :
927 3834 : auto & ctrl_wh = _fe_problem.getControlWarehouse()[EXEC_TIMESTEP_BEGIN];
928 :
929 : // initialize THM control objects
930 7360 : for (auto && i : ctrl_wh.getObjects())
931 : {
932 5443 : THMControl * ctrl = dynamic_cast<THMControl *>(i.get());
933 5443 : if (ctrl != nullptr)
934 5443 : ctrl->init();
935 : }
936 :
937 7360 : for (auto && i : ctrl_wh.getObjects())
938 : {
939 5443 : THMControl * ctrl = dynamic_cast<THMControl *>(i.get());
940 : // if it is a THM control
941 5443 : if (ctrl != nullptr)
942 : {
943 : // get its dependencies on control data
944 : auto & cd_deps = ctrl->getControlDataDependencies();
945 5665 : for (auto && cd_name : cd_deps)
946 : {
947 222 : ControlDataValue * cdv = _control_data[cd_name];
948 : // find out which control object built the control data
949 222 : 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 222 : auto it = std::find(deps.begin(), deps.end(), dep_name);
953 222 : if (it == deps.end())
954 222 : 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 3834 : _fe_problem.getControlWarehouse()[EXEC_TIMESTEP_END];
963 7360 : for (auto && i : ctrl_wh.getObjects())
964 : {
965 5443 : if (TerminateControl * ctrl = dynamic_cast<TerminateControl *>(i.get()))
966 : {
967 : std::list<const THMControl *> l;
968 8 : l.push_back(ctrl);
969 32 : while (l.size() > 0)
970 : {
971 24 : const THMControl * ctrl = l.front();
972 : auto & cd_deps = ctrl->getControlDataDependencies();
973 40 : for (auto && cd_name : cd_deps)
974 : {
975 16 : ControlDataValue * cdv = _control_data[cd_name];
976 16 : l.push_back(cdv->getControl());
977 : }
978 48 : 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 8297 : Simulation::addComponent(const std::string & type, const std::string & name, InputParameters params)
992 : {
993 8297 : std::shared_ptr<Component> comp = _thm_factory.create<Component>(type, name, params);
994 8295 : if (_comp_by_name.find(name) == _comp_by_name.end())
995 8295 : _comp_by_name[name] = comp;
996 : else
997 : logError("Component with name '", name, "' already exists");
998 8295 : _components.push_back(comp);
999 8295 : }
1000 :
1001 : bool
1002 5663 : Simulation::hasComponent(const std::string & name) const
1003 : {
1004 : auto it = _comp_by_name.find(name);
1005 5663 : return (it != _comp_by_name.end());
1006 : }
1007 :
1008 : void
1009 1398 : Simulation::addClosures(const std::string & type, const std::string & name, InputParameters params)
1010 : {
1011 1398 : std::shared_ptr<ClosuresBase> obj_ptr = _thm_factory.create<ClosuresBase>(type, name, params);
1012 1398 : if (_closures_by_name.find(name) == _closures_by_name.end())
1013 1398 : _closures_by_name[name] = obj_ptr;
1014 : else
1015 : logError("A closures object with the name '", name, "' already exists.");
1016 1398 : }
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 2151 : Simulation::getClosures(const std::string & name) const
1026 : {
1027 : auto it = _closures_by_name.find(name);
1028 2151 : if (it != _closures_by_name.end())
1029 2151 : return it->second;
1030 : else
1031 0 : mooseError("The requested closures object '", name, "' does not exist.");
1032 : }
1033 :
1034 : void
1035 1849 : Simulation::addFileOutputter(const std::string & name)
1036 : {
1037 1849 : _outputters_all.push_back(name);
1038 1849 : _outputters_file.push_back(name);
1039 1849 : }
1040 :
1041 : void
1042 1960 : Simulation::addScreenOutputter(const std::string & name)
1043 : {
1044 1960 : _outputters_all.push_back(name);
1045 1960 : _outputters_screen.push_back(name);
1046 1960 : }
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 44148 : Simulation::hasInitialConditionsFromFile() const
1071 : {
1072 88296 : return _thm_pars.isParamValid("initial_from_file");
1073 : }
1074 :
1075 : void
1076 19049 : Simulation::advanceState()
1077 : {
1078 59133 : for (auto && i : _control_data)
1079 40084 : i.second->copyValuesBack();
1080 19049 : }
|
