Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #ifdef ENABLE_OPENMC_COUPLING
20 :
21 : #include "OpenMCProblemBase.h"
22 :
23 : #include "CardinalAppTypes.h"
24 : #include "AddTallyAction.h"
25 : #include "SetupMGXSAction.h"
26 :
27 : #include "OpenMCNuclideDensities.h"
28 : #include "OpenMCDomainFilterEditor.h"
29 : #include "OpenMCTallyEditor.h"
30 :
31 : #include "openmc/random_lcg.h"
32 :
33 : InputParameters
34 3918 : OpenMCProblemBase::validParams()
35 : {
36 3918 : InputParameters params = CardinalProblem::validParams();
37 7836 : params.addParam<PostprocessorName>(
38 : "power", "Power (Watts) to normalize the OpenMC tallies; only used for k-eigenvalue mode");
39 7836 : params.addParam<PostprocessorName>(
40 : "source_strength",
41 : "Neutrons/second to normalize the OpenMC tallies; only used for fixed source mode");
42 7836 : params.addParam<bool>("verbose", false, "Whether to print diagnostic information");
43 :
44 7836 : params.addParam<MooseEnum>("tally_type", getTallyTypeEnum(), "Type of tally to use in OpenMC");
45 :
46 11754 : params.addRangeCheckedParam<Real>(
47 : "scaling",
48 7836 : 1.0,
49 : "scaling > 0.0",
50 : "Scaling factor to apply to [Mesh] to get to units of centimeters that OpenMC expects; "
51 : "setting 'scaling = 100.0', for instance, indicates that the [Mesh] is in units of meters");
52 :
53 : // interfaces to directly set some OpenMC parameters
54 7836 : params.addRangeCheckedParam<unsigned int>(
55 : "openmc_verbosity",
56 : "openmc_verbosity >= 1 & openmc_verbosity <= 10",
57 : "OpenMC verbosity level; this overrides the setting in the XML files. Note that we cannot "
58 : "influence the verbosity of OpenMC's initialization routines, since these are run before "
59 : "Cardinal is initialized.");
60 7836 : params.addRangeCheckedParam<unsigned int>(
61 : "inactive_batches",
62 : "inactive_batches >= 0",
63 : "Number of inactive batches to run in OpenMC; this overrides the setting in the XML files.");
64 7836 : params.addRangeCheckedParam<unsigned int>("particles",
65 : "particles > 0 ",
66 : "Number of particles to run in each OpenMC batch; this "
67 : "overrides the setting in the XML files.");
68 7836 : params.addRangeCheckedParam<unsigned int>(
69 : "batches",
70 : "batches > 0",
71 : "Number of batches to run in OpenMC; this overrides the setting in the XML files.");
72 :
73 7836 : params.addParam<bool>("reuse_source",
74 7836 : false,
75 : "Whether to take the initial fission source "
76 : "for interation n to be the converged source bank from iteration n-1");
77 7836 : params.addParam<bool>(
78 : "skip_statepoint",
79 7836 : false,
80 : "Whether to skip writing any statepoint files from OpenMC; this is a performance "
81 : "optimization for scenarios where you may not want the statepoint files anyways");
82 7836 : params.addParam<bool>(
83 : "reset_seed",
84 7836 : false,
85 : "Whether to reset OpenMC's seed to the initial starting seed before each OpenMC solve");
86 :
87 : // Kinetics parameters.
88 7836 : params.addParam<bool>("calc_kinetics_params",
89 7836 : false,
90 : "Whether or not Cardinal should enable the calculation of kinetics "
91 : "parameters (Lambda effective and beta effective).");
92 7836 : params.addParam<unsigned int>(
93 : "ifp_generations",
94 : openmc::DEFAULT_IFP_N_GENERATION,
95 : "The number of generations to use with the method of iterated fission probabilities.");
96 7836 : params.addParam<FileName>(
97 : "xml_directory", "./", "The directory in which to look for OpenMC XML files.");
98 3918 : return params;
99 0 : }
100 :
101 1976 : OpenMCProblemBase::OpenMCProblemBase(const InputParameters & params)
102 : : CardinalProblem(params),
103 : PostprocessorInterface(this),
104 1976 : _verbose(getParam<bool>("verbose")),
105 3952 : _reuse_source(getParam<bool>("reuse_source")),
106 1976 : _specified_scaling(params.isParamSetByUser("scaling")),
107 3952 : _scaling(getParam<Real>("scaling")),
108 3952 : _skip_statepoint(getParam<bool>("skip_statepoint")),
109 1976 : _fixed_point_iteration(-1),
110 1976 : _total_n_particles(0),
111 1976 : _has_adaptivity(getMooseApp().actionWarehouse().hasActions("set_adaptivity_options")),
112 1976 : _run_on_adaptivity_cycle(true),
113 3952 : _calc_kinetics_params(getParam<bool>("calc_kinetics_params")),
114 3952 : _reset_seed(getParam<bool>("reset_seed")),
115 1976 : _initial_seed(openmc::openmc_get_seed()),
116 5928 : _xml_directory(getParam<FileName>("xml_directory"))
117 : {
118 3952 : if (isParamValid("tally_type"))
119 0 : mooseError("The tally system used by OpenMCProblemBase derived classes has been deprecated. "
120 : "Please add tallies with the [Tallies] block instead.");
121 :
122 5928 : std::vector<std::string> argv_vec = {"openmc", _xml_directory};
123 :
124 : std::vector<char *> argv;
125 :
126 5928 : for (const auto & arg : argv_vec)
127 : {
128 3952 : argv.push_back(const_cast<char *>(arg.data()));
129 : }
130 : // Add terminating nullptr
131 1976 : argv.push_back(nullptr);
132 :
133 1976 : openmc_init(argv.size() - 1, argv.data(), &_communicator.get());
134 :
135 : // ensure that any mapped cells have their distribcell indices generated in OpenMC
136 1976 : if (!openmc::settings::material_cell_offsets)
137 : {
138 0 : mooseWarning("Distributed properties for material cells are disabled "
139 : "in the OpenMC settings. Enabling...");
140 0 : openmc::settings::material_cell_offsets = true;
141 0 : openmc::prepare_distribcell();
142 : }
143 :
144 : // ensure that unsupported run modes are not used, while also checking for
145 : // necessary/unused input parameters for the valid run modes
146 1976 : _run_mode = openmc::settings::run_mode;
147 1976 : const auto & tally_actions = getMooseApp().actionWarehouse().getActions<AddTallyAction>();
148 1976 : const auto & mgxs_actions = getMooseApp().actionWarehouse().getActions<SetupMGXSAction>();
149 1976 : switch (_run_mode)
150 : {
151 : case openmc::RunMode::EIGENVALUE:
152 : {
153 : // Jumping through hoops to see if we're going to add tallies down the line.
154 1870 : if (tally_actions.size() > 0 || mgxs_actions.size() > 0)
155 : {
156 3206 : checkRequiredParam(params, "power", "running in k-eigenvalue mode");
157 1603 : _power = &getPostprocessorValue("power");
158 : }
159 : else
160 534 : checkUnusedParam(params, "power", "no tallies have been added");
161 :
162 3740 : checkUnusedParam(params, "source_strength", "running in k-eigenvalue mode");
163 1870 : break;
164 : }
165 : case openmc::RunMode::FIXED_SOURCE:
166 : {
167 100 : if (tally_actions.size() > 0 || mgxs_actions.size() > 0)
168 : {
169 184 : checkRequiredParam(params, "source_strength", "running in fixed source mode");
170 92 : _source_strength = &getPostprocessorValue("source_strength");
171 : }
172 : else
173 16 : checkUnusedParam(params, "source_strength", "no tallies have been added");
174 :
175 200 : checkUnusedParam(params, "inactive_batches", "running in fixed source mode");
176 200 : checkUnusedParam(params, "reuse_source", "running in fixed source mode");
177 200 : checkUnusedParam(params, "power", "running in fixed source mode");
178 100 : _reuse_source = false;
179 100 : break;
180 : }
181 6 : case openmc::RunMode::PLOTTING:
182 : case openmc::RunMode::PARTICLE:
183 : case openmc::RunMode::VOLUME:
184 6 : mooseError("Running OpenMC in plotting, particle, and volume modes is not supported through "
185 : "Cardinal! Please just run using the OpenMC executable (e.g., openmc --plot for "
186 : "plot mode).");
187 0 : default:
188 0 : mooseError("Unhandled openmc::RunMode enum in OpenMCInitAction!");
189 : }
190 :
191 1970 : _n_cell_digits = std::to_string(openmc::model::cells.size()).length();
192 :
193 1970 : if (openmc::settings::libmesh_comm)
194 0 : mooseWarning("libMesh communicator already set in OpenMC.");
195 :
196 1970 : openmc::settings::libmesh_comm = &_mesh.comm();
197 :
198 3940 : if (isParamValid("openmc_verbosity"))
199 0 : openmc::settings::verbosity = getParam<unsigned int>("openmc_verbosity");
200 :
201 3940 : if (isParamValid("inactive_batches"))
202 200 : openmc::settings::n_inactive = getParam<unsigned int>("inactive_batches");
203 :
204 3940 : if (isParamValid("particles"))
205 248 : openmc::settings::n_particles = getParam<unsigned int>("particles");
206 :
207 3940 : if (isParamValid("batches"))
208 : {
209 130 : auto xml_n_batches = openmc::settings::n_batches;
210 :
211 260 : int err = openmc_set_n_batches(getParam<unsigned int>("batches"),
212 : true /* set the max batches */,
213 130 : true /* add the last batch for statepoint writing */);
214 258 : catchOpenMCError(err, "set the number of batches");
215 :
216 : // if we set the batches from Cardinal, remove whatever statepoint file was
217 : // created for the #batches set in the XML files; this is just to reduce the
218 : // number of statepoint files by removing an unnecessary point
219 : openmc::settings::statepoint_batch.erase(xml_n_batches);
220 : }
221 :
222 : // The OpenMC wrapping doesn't require material properties itself, but we might
223 : // define them on some blocks of the domain for other auxiliary kernel purposes
224 : setMaterialCoverageCheck(false);
225 :
226 : // If the user requests kinetics parameters, make sure it's enabled in OpenMC.
227 1968 : if (_calc_kinetics_params)
228 : {
229 12 : if (_run_mode != openmc::RunMode::EIGENVALUE)
230 2 : paramError("calc_kinetics_params",
231 : "Kinetic parameters can only be calculated in k-eigenvalue mode!");
232 :
233 10 : openmc::settings::ifp_on = true;
234 10 : openmc::settings::ifp_parameter = openmc::IFPParameter::Both;
235 :
236 20 : openmc::settings::ifp_n_generation = getParam<unsigned int>("ifp_generations");
237 10 : if (openmc::settings::ifp_n_generation > openmc::settings::n_inactive)
238 2 : paramError("ifp_generations",
239 : "'ifp_generations' must be less than or equal to the number of inactive batches!");
240 : }
241 5916 : }
242 :
243 1709 : OpenMCProblemBase::~OpenMCProblemBase() { openmc_finalize(); }
244 :
245 : void
246 227834083 : OpenMCProblemBase::catchOpenMCError(const int & err, const std::string descriptor) const
247 : {
248 227834083 : if (err)
249 8 : mooseError(
250 8 : "In attempting to ", descriptor, ", OpenMC reported:\n\n", std::string(openmc_err_msg));
251 227834075 : }
252 :
253 : void
254 0 : OpenMCProblemBase::fillElementalAuxVariable(const unsigned int & var_num,
255 : const std::vector<unsigned int> & elem_ids,
256 : const Real & value)
257 : {
258 0 : auto & solution = _aux->solution();
259 0 : auto sys_number = _aux->number();
260 :
261 : // loop over all the elements and set the specified variable to the specified value
262 0 : for (const auto & e : elem_ids)
263 : {
264 0 : auto elem_ptr = _mesh.queryElemPtr(e);
265 :
266 0 : if (!isLocalElem(elem_ptr))
267 0 : continue;
268 :
269 0 : auto dof_idx = elem_ptr->dof_number(sys_number, var_num, 0);
270 0 : solution.set(dof_idx, value);
271 : }
272 0 : }
273 :
274 : int
275 6591 : OpenMCProblemBase::nParticles() const
276 : {
277 6591 : return openmc::settings::n_particles;
278 : }
279 :
280 : std::string
281 6855 : OpenMCProblemBase::materialName(const int32_t index) const
282 : {
283 : // OpenMC uses -1 to indicate void materials, which don't have a name. So we return
284 : // one ourselves, or else openmc_material_get_name will throw an error.
285 6855 : if (index == -1)
286 16 : return "VOID";
287 :
288 : const char * name;
289 6839 : int err = openmc_material_get_name(index, &name);
290 6839 : catchOpenMCError(err, "get material name for material with index " + std::to_string(index));
291 :
292 6839 : std::string n = name;
293 :
294 : // if the material does not have a name, just return the ID instead
295 6839 : if (n.empty())
296 6974 : n = std::to_string(materialID(index));
297 :
298 6839 : return n;
299 : }
300 :
301 : int32_t
302 130215152 : OpenMCProblemBase::cellID(const int32_t index) const
303 : {
304 : int32_t id;
305 130215152 : int err = openmc_cell_get_id(index, &id);
306 130215152 : catchOpenMCError(err, "get ID for cell with index " + std::to_string(index));
307 130215152 : return id;
308 : }
309 :
310 : int32_t
311 2850909 : OpenMCProblemBase::materialID(const int32_t index) const
312 : {
313 2850909 : if (index == openmc::MATERIAL_VOID)
314 : return -1;
315 :
316 : int32_t id;
317 2850837 : int err = openmc_material_get_id(index, &id);
318 2850837 : catchOpenMCError(err, "get ID for material with index " + std::to_string(index));
319 2850837 : return id;
320 : }
321 :
322 : std::string
323 0 : OpenMCProblemBase::printMaterial(const int32_t & index) const
324 : {
325 0 : int32_t id = materialID(index);
326 0 : std::stringstream msg;
327 0 : msg << "material " << id;
328 0 : return msg.str();
329 0 : }
330 :
331 : std::string
332 10 : OpenMCProblemBase::printPoint(const Point & p) const
333 : {
334 10 : std::stringstream msg;
335 10 : msg << "(" << std::setprecision(6) << std::setw(7) << p(0) << ", " << std::setprecision(6)
336 10 : << std::setw(7) << p(1) << ", " << std::setprecision(6) << std::setw(7) << p(2) << ")";
337 10 : return msg.str();
338 10 : }
339 :
340 : bool
341 72 : OpenMCProblemBase::firstSolve() const
342 : {
343 72 : return _fixed_point_iteration < 0;
344 : }
345 :
346 : void
347 2332 : OpenMCProblemBase::externalSolve()
348 : {
349 4664 : TIME_SECTION("solveOpenMC", 1, "Solving OpenMC", false);
350 :
351 : // Check to see if this is a steady solve. If so, we can skip extra OpenMC runs
352 : // once the mesh stops getting adapted.
353 2332 : if (_has_adaptivity && !_run_on_adaptivity_cycle)
354 : {
355 30 : _console << " Skipping running OpenMC as the mesh has not changed!" << std::endl;
356 : return;
357 : }
358 :
359 2302 : _console << " Running OpenMC with " << nParticles() << " particles per batch..." << std::endl;
360 :
361 : // apply a new starting fission source
362 2302 : if (_reuse_source && !firstSolve())
363 : {
364 16 : openmc::free_memory_source();
365 16 : openmc::model::external_sources.push_back(
366 48 : std::make_unique<openmc::FileSource>(sourceBankFileName()));
367 : }
368 :
369 : // update tallies as needed before starting the OpenMC run
370 2302 : executeEditors();
371 :
372 2296 : if (_reset_seed)
373 : {
374 24 : openmc_hard_reset();
375 24 : openmc_set_seed(_initial_seed);
376 : }
377 :
378 2296 : int err = openmc_run();
379 2296 : if (err)
380 0 : mooseError(openmc_err_msg);
381 :
382 2296 : _total_n_particles += nParticles();
383 :
384 2296 : err = openmc_reset_timers();
385 2296 : if (err)
386 0 : mooseError(openmc_err_msg);
387 :
388 2296 : _fixed_point_iteration++;
389 :
390 : // save the latest fission source for re-use in the next iteration
391 2296 : if (_reuse_source)
392 48 : writeSourceBank(sourceBankFileName());
393 : }
394 :
395 : void
396 1797 : OpenMCProblemBase::initialSetup()
397 : {
398 1797 : ExternalProblem::initialSetup();
399 :
400 : // Initialize the IFP parameters tally.
401 1797 : if (_calc_kinetics_params)
402 : {
403 8 : _ifp_tally_index = openmc::model::tallies.size();
404 8 : _ifp_tally = openmc::Tally::create();
405 8 : _ifp_tally->set_scores({"ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"});
406 8 : _ifp_tally->estimator_ = openmc::TallyEstimator::COLLISION;
407 : }
408 1797 : }
409 :
410 : void
411 4674 : OpenMCProblemBase::syncSolutions(ExternalProblem::Direction direction)
412 : {
413 : // Always run OpenMC on the first timestep in a steady solve with adaptivity. This
414 : // ensures that OpenMC runs at least once during each Picard iteration.
415 4674 : _run_on_adaptivity_cycle |= (timeStep() == 1 && !isTransient());
416 4674 : }
417 :
418 : bool
419 605 : OpenMCProblemBase::adaptMesh()
420 : {
421 605 : _run_on_adaptivity_cycle = CardinalProblem::adaptMesh() || isTransient();
422 605 : return _run_on_adaptivity_cycle;
423 : }
424 :
425 : void
426 24 : OpenMCProblemBase::writeSourceBank(const std::string & filename)
427 : {
428 24 : hid_t file_id = openmc::file_open(filename, 'w', true);
429 : openmc::write_attribute(file_id, "filetype", "source");
430 24 : openmc::write_source_bank(
431 : file_id, openmc::simulation::source_bank, openmc::simulation::work_index);
432 24 : openmc::file_close(file_id);
433 24 : }
434 :
435 : unsigned int
436 2507 : OpenMCProblemBase::numElemsInSubdomain(const SubdomainID & id) const
437 : {
438 2507 : unsigned int n = 0;
439 9217023 : for (unsigned int e = 0; e < _mesh.nElem(); ++e)
440 : {
441 9214516 : const auto * elem = _mesh.queryElemPtr(e);
442 :
443 9214516 : if (!isLocalElem(elem) || !elem->active())
444 4442804 : continue;
445 :
446 : const auto subdomain_id = elem->subdomain_id();
447 4771712 : if (id == subdomain_id)
448 1863372 : n += 1;
449 : }
450 :
451 2507 : _communicator.sum(n);
452 :
453 2507 : return n;
454 : }
455 :
456 : bool
457 136448189 : OpenMCProblemBase::isLocalElem(const Elem * elem) const
458 : {
459 136448189 : if (!elem)
460 : {
461 : // we should only not be able to find an element if the mesh is distributed
462 : libmesh_assert(!_mesh.getMesh().is_serial());
463 : return false;
464 : }
465 :
466 96952347 : if (elem->processor_id() == _communicator.rank())
467 57794289 : return true;
468 :
469 : return false;
470 : }
471 :
472 : bool
473 4 : OpenMCProblemBase::cellHasZeroInstances(const cellInfo & cell_info) const
474 : {
475 4 : auto n = openmc::model::cells.at(cell_info.first)->n_instances();
476 4 : return !n;
477 : }
478 :
479 : void
480 396634508 : OpenMCProblemBase::setCellTemperature(const int32_t & index,
481 : const int32_t & instance,
482 : const Real & T,
483 : const cellInfo & cell_info) const
484 : {
485 396634508 : int err = openmc_cell_set_temperature(index, T, &instance, false);
486 396634508 : if (err)
487 : {
488 : std::string descriptor =
489 12 : "set cell " + printCell(cell_info) + " to temperature " + Moose::stringify(T) + " (K)";
490 :
491 : // special error message if cell has zero instances
492 4 : if (cellHasZeroInstances(cell_info))
493 0 : mooseError("Failed to set the temperature for cell " + printCell(cell_info) +
494 : " with zero instances.");
495 :
496 12 : mooseError("In attempting to set cell " + printCell(cell_info) + " to temperature " +
497 4 : Moose::stringify(T) + " (K), OpenMC reported:\n\n",
498 4 : std::string(openmc_err_msg) + "\n\n" +
499 : "If you are trying to debug a model setup, you can set 'initial_properties = "
500 : "xml' to use the initial temperature and density in the OpenMC XML files for "
501 : "OpenMC's first run.");
502 : }
503 396634504 : }
504 :
505 : std::vector<int32_t>
506 86681791 : OpenMCProblemBase::cellFill(const cellInfo & cell_info, int & fill_type) const
507 : {
508 86681791 : int32_t * materials = nullptr;
509 86681791 : int n_materials = 0;
510 :
511 86681791 : int err = openmc_cell_get_fill(cell_info.first, &fill_type, &materials, &n_materials);
512 173363582 : catchOpenMCError(err, "get fill of cell " + printCell(cell_info));
513 :
514 : std::vector<int32_t> material_indices;
515 86681791 : material_indices.assign(materials, materials + n_materials);
516 86681791 : return material_indices;
517 0 : }
518 :
519 : bool
520 86681791 : OpenMCProblemBase::materialFill(const cellInfo & cell_info, int32_t & material_index) const
521 : {
522 : int fill_type;
523 86681791 : auto material_indices = cellFill(cell_info, fill_type);
524 :
525 86681791 : if (fill_type != static_cast<int>(openmc::Fill::MATERIAL))
526 : return false;
527 :
528 : // The number of materials in a cell is either 1, or equal to the number of instances
529 : // (if distributed materials were used).
530 86681789 : if (material_indices.size() == 1)
531 86678861 : material_index = material_indices[0];
532 : else
533 2928 : material_index = material_indices[cell_info.second];
534 :
535 : return true;
536 86681791 : }
537 :
538 : void
539 1964 : OpenMCProblemBase::setCellDensity(const Real & density, const cellInfo & cell_info) const
540 : {
541 : // OpenMC technically allows a density of >= 0.0, but we can impose a tighter
542 : // check here with a better error message than the Excepts() in material->set_density
543 : // because it could be a very common mistake to forget to set an initial condition
544 : // for density if OpenMC runs first
545 1964 : if (density <= 0.0)
546 4 : mooseError("Densities less than or equal to zero cannot be set in the OpenMC model!\n\n cell " +
547 4 : printCell(cell_info) + " set to density " + Moose::stringify(density) + " (kg/m3)");
548 :
549 : int32_t material_index;
550 1962 : auto is_material_cell = materialFill(cell_info, material_index);
551 :
552 1962 : if (!is_material_cell)
553 0 : mooseError(
554 : "Density transfer does not currently support cells filled with universes or lattices!");
555 :
556 : // throw a special error if the cell is void, because the OpenMC error isn't very
557 : // clear what the mistake is
558 1962 : if (material_index == MATERIAL_VOID)
559 : {
560 12 : mooseWarning("Skipping setting density for cell " + printCell(cell_info) +
561 : " because this cell is void (vacuum)");
562 4 : return;
563 : }
564 :
565 : // Compute the density. We multiply density by 0.001 to convert from kg/m3
566 : // (the units assumed in the 'density' auxvariable as well as the MOOSE fluid
567 : // properties module) to g/cm3
568 1956 : int err = openmc_cell_set_density(
569 1956 : cell_info.first, _density_conversion_factor * density, &cell_info.second, false);
570 :
571 1956 : if (err)
572 : {
573 : // special error message if cell has zero instances
574 0 : if (cellHasZeroInstances(cell_info))
575 0 : mooseError("Failed to set the density for cell " + printCell(cell_info) +
576 : " with zero instances.");
577 :
578 0 : mooseError("In attempting to set cell " + printCell(cell_info) + " to density " +
579 0 : Moose::stringify(density) + " (kg/m3), OpenMC reported:\n\n",
580 0 : std::string(openmc_err_msg) + "\n\n" +
581 : "If you are trying to debug a model setup, you can set 'initial_properties = "
582 : "xml' to use the initial temperature and density in the OpenMC XML files for "
583 : "OpenMC's first run.");
584 : }
585 : }
586 :
587 : std::string
588 94811146 : OpenMCProblemBase::printCell(const cellInfo & cell_info, const bool brief) const
589 : {
590 94811146 : int32_t id = cellID(cell_info.first);
591 :
592 94811146 : std::stringstream msg;
593 94811146 : if (!brief)
594 94772141 : msg << "id ";
595 :
596 189622292 : msg << std::setw(_n_cell_digits) << Moose::stringify(id) << ", instance "
597 189622292 : << std::setw(_n_cell_digits) << Moose::stringify(cell_info.second) << " (of "
598 94811146 : << std::setw(_n_cell_digits)
599 379244584 : << Moose::stringify(openmc::model::cells.at(cell_info.first)->n_instances()) << ")";
600 :
601 94811146 : return msg.str();
602 94811146 : }
603 :
604 : void
605 2 : OpenMCProblemBase::importProperties() const
606 : {
607 2 : _console << "Reading temperature and density from properties.h5" << std::endl;
608 :
609 2 : int err = openmc_properties_import("properties.h5");
610 2 : catchOpenMCError(err, "load temperature and density from a properties.h5 file");
611 0 : }
612 :
613 : xt::xtensor<double, 1>
614 4228 : OpenMCProblemBase::relativeError(const xt::xtensor<double, 1> & sum,
615 : const xt::xtensor<double, 1> & sum_sq,
616 : const int & n_realizations) const
617 : {
618 4228 : xt::xtensor<double, 1> rel_err = xt::zeros<double>({sum.size()});
619 :
620 1705846 : for (unsigned int i = 0; i < sum.size(); ++i)
621 : {
622 1701618 : auto mean = sum(i) / n_realizations;
623 1701618 : auto std_dev = std::sqrt((sum_sq(i) / n_realizations - mean * mean) / (n_realizations - 1));
624 1701618 : rel_err[i] = mean != 0.0 ? std_dev / std::abs(mean) : 0.0;
625 : }
626 :
627 4228 : return rel_err;
628 : }
629 :
630 : xt::xtensor<double, 1>
631 5482 : OpenMCProblemBase::tallySum(openmc::Tally * tally, const unsigned int & score) const
632 : {
633 10964 : return xt::view(tally->results_, xt::all(), score, static_cast<int>(openmc::TallyResult::SUM));
634 : }
635 :
636 : double
637 1876 : OpenMCProblemBase::tallySumAcrossBins(std::vector<openmc::Tally *> tally,
638 : const unsigned int & score) const
639 : {
640 : double sum = 0.0;
641 :
642 3752 : for (const auto & t : tally)
643 : {
644 1876 : auto mean = tallySum(t, score);
645 3752 : sum += xt::sum(mean)();
646 1876 : }
647 :
648 1876 : return sum;
649 : }
650 :
651 : double
652 0 : OpenMCProblemBase::tallyMeanAcrossBins(std::vector<openmc::Tally *> tally,
653 : const unsigned int & score) const
654 : {
655 : int n = 0;
656 0 : for (const auto & t : tally)
657 0 : n += t->n_realizations_;
658 :
659 0 : return tallySumAcrossBins(tally, score) / n;
660 : }
661 :
662 : std::string
663 2624 : OpenMCProblemBase::enumToTallyScore(const std::string & score) const
664 : {
665 : // the MultiMooseEnum is all caps, but the MooseEnum is already the correct case,
666 : // so we need to treat these as separate
667 2624 : std::string s = score;
668 2624 : if (std::all_of(
669 17468 : s.begin(), s.end(), [](unsigned char c) { return !std::isalpha(c) || std::isupper(c); }))
670 : {
671 18544 : std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
672 :
673 : // we need to revert back to some letters being uppercase for certain scores
674 1874 : if (s == "h3_production")
675 : s = "H3_production";
676 : }
677 :
678 : // MOOSE enums use underscores, OpenMC uses dashes
679 : std::replace(s.begin(), s.end(), '_', '-');
680 2624 : return s;
681 : }
682 :
683 : std::string
684 0 : OpenMCProblemBase::tallyScoreToEnum(const std::string & score) const
685 : {
686 : // MOOSE enums use underscores, OpenMC uses dashes
687 0 : std::string s = score;
688 : std::replace(s.begin(), s.end(), '-', '_');
689 0 : return s;
690 : }
691 :
692 : openmc::TallyEstimator
693 346 : OpenMCProblemBase::tallyEstimator(tally::TallyEstimatorEnum estimator) const
694 : {
695 : switch (estimator)
696 : {
697 : case tally::tracklength:
698 : return openmc::TallyEstimator::TRACKLENGTH;
699 : case tally::collision:
700 : return openmc::TallyEstimator::COLLISION;
701 : case tally::analog:
702 : return openmc::TallyEstimator::ANALOG;
703 0 : default:
704 0 : mooseError("Unhandled TallyEstimatorEnum!");
705 : }
706 : }
707 :
708 : std::string
709 0 : OpenMCProblemBase::estimatorToString(openmc::TallyEstimator estimator) const
710 : {
711 0 : switch (estimator)
712 : {
713 : case openmc::TallyEstimator::TRACKLENGTH:
714 0 : return "tracklength";
715 : case openmc::TallyEstimator::COLLISION:
716 0 : return "collision";
717 : case openmc::TallyEstimator::ANALOG:
718 0 : return "analog";
719 0 : default:
720 0 : mooseError("Unhandled TallyEstimatorEnum!");
721 : }
722 : }
723 :
724 : openmc::TriggerMetric
725 138 : OpenMCProblemBase::triggerMetric(std::string trigger) const
726 : {
727 138 : if (trigger == "variance")
728 : return openmc::TriggerMetric::variance;
729 138 : else if (trigger == "std_dev")
730 : return openmc::TriggerMetric::standard_deviation;
731 138 : else if (trigger == "rel_err")
732 : return openmc::TriggerMetric::relative_error;
733 8 : else if (trigger == "none")
734 : return openmc::TriggerMetric::not_active;
735 : else
736 0 : mooseError("Unhandled TallyTriggerTypeEnum: ", trigger);
737 : }
738 :
739 : openmc::TriggerMetric
740 1858 : OpenMCProblemBase::triggerMetric(trigger::TallyTriggerTypeEnum trigger) const
741 : {
742 : switch (trigger)
743 : {
744 : case trigger::variance:
745 : return openmc::TriggerMetric::variance;
746 : case trigger::std_dev:
747 : return openmc::TriggerMetric::standard_deviation;
748 : case trigger::rel_err:
749 : return openmc::TriggerMetric::relative_error;
750 : case trigger::none:
751 : return openmc::TriggerMetric::not_active;
752 0 : default:
753 0 : mooseError("Unhandled TallyTriggerTypeEnum!");
754 : }
755 : }
756 :
757 : bool
758 0 : OpenMCProblemBase::cellIsVoid(const cellInfo & cell_info) const
759 : {
760 : // material_index will be unchanged if the cell is filled by a universe or lattice.
761 : // Otherwise, this will get set to the material index in the cell.
762 0 : int32_t material_index = 0;
763 0 : materialFill(cell_info, material_index);
764 0 : return material_index == MATERIAL_VOID;
765 : }
766 :
767 : void
768 1001 : OpenMCProblemBase::geometryType(bool & has_csg_universe, bool & has_dag_universe) const
769 : {
770 1001 : has_csg_universe = false;
771 1001 : has_dag_universe = false;
772 :
773 : // Loop over universes and check if type is DAGMC
774 13900 : for (const auto & universe : openmc::model::universes)
775 : {
776 12899 : if (universe->geom_type() == openmc::GeometryType::DAG)
777 50 : has_dag_universe = true;
778 12849 : else if (universe->geom_type() == openmc::GeometryType::CSG)
779 12849 : has_csg_universe = true;
780 : else
781 0 : mooseError("Unhandled GeometryType enum!");
782 : }
783 1001 : }
784 :
785 : long unsigned int
786 1893 : OpenMCProblemBase::numCells() const
787 : {
788 : long unsigned int n_openmc_cells = 0;
789 289568 : for (const auto & c : openmc::model::cells)
790 287675 : n_openmc_cells += c->n_instances();
791 :
792 1893 : return n_openmc_cells;
793 : }
794 :
795 : const openmc::Tally &
796 48 : OpenMCProblemBase::getKineticsParamTally()
797 : {
798 48 : if (!_ifp_tally)
799 0 : mooseError("Internal error: kinetics parameters have not been enabled.");
800 :
801 48 : return *_ifp_tally;
802 : }
803 :
804 : bool
805 628150 : OpenMCProblemBase::isReactionRateScore(const std::string & score) const
806 : {
807 : const std::set<std::string> viable_scores = {"H3-production",
808 : "total",
809 : "absorption",
810 : "scatter",
811 : "nu-scatter",
812 : "fission",
813 : "nu-fission",
814 : "prompt-nu-fission",
815 628150 : "delayed-nu-fission"};
816 628150 : return viable_scores.count(score);
817 : }
818 :
819 : bool
820 1159348 : OpenMCProblemBase::isHeatingScore(const std::string & score) const
821 : {
822 : const std::set<std::string> viable_scores = {
823 1159348 : "heating", "heating-local", "kappa-fission", "fission-q-prompt", "fission-q-recoverable"};
824 1159348 : return viable_scores.count(score);
825 : }
826 :
827 : unsigned int
828 9091 : OpenMCProblemBase::addExternalVariable(const std::string & name,
829 : const std::vector<SubdomainName> * block)
830 : {
831 9091 : auto var_params = _factory.getValidParams("MooseVariable");
832 18182 : var_params.set<MooseEnum>("family") = "MONOMIAL";
833 18182 : var_params.set<MooseEnum>("order") = "CONSTANT";
834 :
835 9091 : if (block)
836 11890 : var_params.set<std::vector<SubdomainName>>("block") = *block;
837 :
838 9091 : checkDuplicateVariableName(name);
839 18178 : addAuxVariable("MooseVariable", name, var_params);
840 18178 : return _aux->getFieldVariable<Real>(0, name).number();
841 9089 : }
842 :
843 : std::string
844 5328 : OpenMCProblemBase::subdomainName(const SubdomainID & id) const
845 : {
846 5328 : std::string name = _mesh.getSubdomainName(id);
847 5328 : if (name.empty())
848 9620 : name = std::to_string(id);
849 5328 : return name;
850 : }
851 :
852 : void
853 1797 : OpenMCProblemBase::getOpenMCUserObjects()
854 : {
855 1797 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("UserObject");
856 : std::vector<UserObject *> userobjs;
857 : uo_query.queryInto(userobjs);
858 :
859 8695 : for (const auto & u : userobjs)
860 : {
861 6898 : OpenMCNuclideDensities * c = dynamic_cast<OpenMCNuclideDensities *>(u);
862 6898 : if (c)
863 44 : _nuclide_densities_uos.push_back(c);
864 :
865 6898 : OpenMCTallyEditor * e = dynamic_cast<OpenMCTallyEditor *>(u);
866 6898 : if (e)
867 60 : _tally_editor_uos.push_back(e);
868 :
869 6898 : OpenMCDomainFilterEditor * f = dynamic_cast<OpenMCDomainFilterEditor *>(u);
870 6898 : if (f)
871 28 : _filter_editor_uos.push_back(f);
872 : }
873 :
874 1797 : checkOpenMCUserObjectIDs();
875 1793 : }
876 :
877 : void
878 1797 : OpenMCProblemBase::checkOpenMCUserObjectIDs() const
879 : {
880 : std::set<int32_t> tally_ids;
881 1855 : for (const auto & te : _tally_editor_uos)
882 : {
883 60 : int32_t tally_id = te->tallyId();
884 : if (tally_ids.count(tally_id) != 0)
885 2 : te->duplicateTallyError(tally_id);
886 58 : tally_ids.insert(tally_id);
887 : }
888 :
889 : std::set<int32_t> filter_ids;
890 1821 : for (const auto & fe : _filter_editor_uos)
891 : {
892 28 : int32_t filter_id = fe->filterId();
893 : if (filter_ids.count(filter_id) != 0)
894 2 : fe->duplicateFilterError(filter_id);
895 26 : filter_ids.insert(filter_id);
896 : }
897 1793 : }
898 :
899 : void
900 1848 : OpenMCProblemBase::checkTallyEditorIDs() const
901 : {
902 1848 : std::vector<int32_t> mapped_tally_ids = getMappedTallyIDs();
903 :
904 1902 : for (const auto & te : _tally_editor_uos)
905 : {
906 56 : int32_t tally_id = te->tallyId();
907 :
908 : // ensure that the TallyEditor IDs don't apply to any mapped tally objects
909 56 : if (std::find(mapped_tally_ids.begin(), mapped_tally_ids.end(), tally_id) !=
910 : mapped_tally_ids.end())
911 2 : te->mappedTallyError(tally_id);
912 : }
913 1846 : }
914 :
915 : void
916 2302 : OpenMCProblemBase::executeFilterEditors()
917 : {
918 2302 : executeControls(EXEC_FILTER_EDITORS);
919 2302 : _console << "Executing filter editors..." << std::endl;
920 2326 : for (const auto & fe : _filter_editor_uos)
921 24 : fe->execute();
922 2302 : }
923 :
924 : void
925 2302 : OpenMCProblemBase::executeTallyEditors()
926 : {
927 2302 : executeControls(EXEC_TALLY_EDITORS);
928 2302 : _console << "Executing tally editors..." << std::endl;
929 2350 : for (const auto & te : _tally_editor_uos)
930 54 : te->execute();
931 2296 : }
932 :
933 : void
934 2302 : OpenMCProblemBase::executeEditors()
935 : {
936 2302 : executeFilterEditors();
937 2302 : executeTallyEditors();
938 2296 : }
939 :
940 : void
941 2318 : OpenMCProblemBase::sendNuclideDensitiesToOpenMC()
942 : {
943 2318 : if (_nuclide_densities_uos.size() == 0)
944 : return;
945 :
946 : // We could probably put this somewhere better, but it's good for now
947 44 : executeControls(EXEC_SEND_OPENMC_DENSITIES);
948 :
949 44 : _console << "Sending nuclide compositions to OpenMC... ";
950 84 : for (const auto & uo : _nuclide_densities_uos)
951 44 : uo->setValue();
952 : }
953 :
954 : #endif
|
