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