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