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 4249 : OpenMCProblemBase::validParams()
39 : {
40 4249 : InputParameters params = CardinalProblem::validParams();
41 8498 : params.addParam<PostprocessorName>(
42 : "power", "Power (Watts) to normalize the OpenMC tallies; only used for k-eigenvalue mode");
43 8498 : params.addParam<PostprocessorName>(
44 : "source_strength",
45 : "Neutrons/second to normalize the OpenMC tallies; only used for fixed source mode");
46 8498 : params.addParam<bool>("verbose", false, "Whether to print diagnostic information");
47 :
48 8498 : params.addParam<MooseEnum>("tally_type", getTallyTypeEnum(), "Type of tally to use in OpenMC");
49 :
50 12747 : params.addRangeCheckedParam<Real>(
51 : "scaling",
52 8498 : 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 8498 : 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 8498 : 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 8498 : params.addParam<PostprocessorName>("particles",
69 : "Number of particles to run in each OpenMC batch; this "
70 : "overrides the setting in the XML files.");
71 8498 : 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 8498 : params.addParam<bool>("reuse_source",
77 8498 : false,
78 : "Whether to take the initial fission source "
79 : "for interation n to be the converged source bank from iteration n-1");
80 8498 : params.addParam<bool>(
81 : "skip_statepoint",
82 8498 : 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 8498 : params.addParam<bool>(
86 : "reset_seed",
87 8498 : false,
88 : "Whether to reset OpenMC's seed to the initial starting seed before each OpenMC solve");
89 :
90 : // Kinetics parameters.
91 8498 : params.addParam<bool>("calc_kinetics_params",
92 8498 : false,
93 : "Whether or not Cardinal should enable the calculation of kinetics "
94 : "parameters (Lambda effective and beta effective).");
95 8498 : 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 8498 : params.addParam<FileName>(
100 : "xml_directory", "./", "The directory in which to look for OpenMC XML files.");
101 4249 : return params;
102 0 : }
103 :
104 2141 : OpenMCProblemBase::OpenMCProblemBase(const InputParameters & params)
105 : : CardinalProblem(params),
106 : PostprocessorInterface(this),
107 2141 : _verbose(getParam<bool>("verbose")),
108 4282 : _reuse_source(getParam<bool>("reuse_source")),
109 2141 : _specified_scaling(params.isParamSetByUser("scaling")),
110 4282 : _scaling(getParam<Real>("scaling")),
111 4282 : _skip_statepoint(getParam<bool>("skip_statepoint")),
112 2141 : _fixed_point_iteration(-1),
113 2141 : _total_n_particles(0),
114 2141 : _has_adaptivity(getMooseApp().actionWarehouse().hasActions("set_adaptivity_options")),
115 2141 : _run_on_adaptivity_cycle(true),
116 4282 : _calc_kinetics_params(getParam<bool>("calc_kinetics_params")),
117 4282 : _reset_seed(getParam<bool>("reset_seed")),
118 2141 : _initial_seed(openmc::openmc_get_seed()),
119 6423 : _xml_directory(getParam<FileName>("xml_directory"))
120 : {
121 4282 : 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 2141 : std::vector<std::string> argv_vec = {"openmc"};
128 8564 : 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 2141 : argv_vec.push_back(_xml_directory);
135 :
136 : std::vector<char *> argv;
137 :
138 6423 : for (const auto & arg : argv_vec)
139 : {
140 4282 : argv.push_back(const_cast<char *>(arg.data()));
141 : }
142 : // Add terminating nullptr
143 2141 : argv.push_back(nullptr);
144 :
145 2141 : openmc_init(argv.size() - 1, argv.data(), &_communicator.get());
146 :
147 : // ensure that any mapped cells have their distribcell indices generated in OpenMC
148 2141 : 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 2141 : _run_mode = openmc::settings::run_mode;
159 2141 : const auto & tally_actions = getMooseApp().actionWarehouse().getActions<AddTallyAction>();
160 2141 : const auto & mgxs_actions = getMooseApp().actionWarehouse().getActions<SetupMGXSAction>();
161 2141 : 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 2019 : if (tally_actions.size() > 0 || mgxs_actions.size() > 0)
167 : {
168 3300 : checkRequiredParam(params, "power", "running in k-eigenvalue mode");
169 1650 : _power = &getPostprocessorValue("power");
170 : }
171 : else
172 738 : checkUnusedParam(params, "power", "no tallies have been added");
173 :
174 4038 : checkUnusedParam(params, "source_strength", "running in k-eigenvalue mode");
175 2019 : 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 2135 : _n_cell_digits = std::to_string(openmc::model::cells.size()).length();
204 :
205 2135 : if (openmc::settings::libmesh_comm)
206 0 : mooseWarning("libMesh communicator already set in OpenMC.");
207 :
208 2135 : openmc::settings::libmesh_comm = &_mesh.comm();
209 :
210 4270 : if (isParamValid("openmc_verbosity"))
211 0 : openmc::settings::verbosity = getParam<unsigned int>("openmc_verbosity");
212 :
213 4270 : if (isParamValid("inactive_batches"))
214 64 : openmc::settings::n_inactive = getParam<unsigned int>("inactive_batches");
215 :
216 4270 : if (isParamValid("particles"))
217 68 : _particles = &getPostprocessorValue("particles");
218 :
219 4270 : 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 2133 : 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 2129 : }
256 :
257 1819 : 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 4970 : OpenMCProblemBase::nParticles() const
282 : {
283 4970 : return openmc::settings::n_particles;
284 : }
285 :
286 : std::string
287 6726 : 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 6726 : if (index == -1)
292 18 : return "VOID";
293 :
294 : const char * name;
295 6708 : int err = openmc_material_get_name(index, &name);
296 6708 : catchOpenMCError(err, "get material name for material with index " + std::to_string(index));
297 :
298 6708 : std::string n = name;
299 :
300 : // if the material does not have a name, just return the ID instead
301 6708 : if (n.empty())
302 7284 : n = std::to_string(materialID(index));
303 :
304 6708 : return n;
305 : }
306 :
307 : int32_t
308 17397064 : OpenMCProblemBase::cellID(const int32_t index) const
309 : {
310 : int32_t id;
311 17397064 : int err = openmc_cell_get_id(index, &id);
312 17397064 : catchOpenMCError(err, "get ID for cell with index " + std::to_string(index));
313 17397064 : return id;
314 : }
315 :
316 : int32_t
317 2581822 : OpenMCProblemBase::materialID(const int32_t index) const
318 : {
319 2581822 : if (index == openmc::MATERIAL_VOID)
320 : return -1;
321 :
322 : int32_t id;
323 2581750 : int err = openmc_material_get_id(index, &id);
324 2581750 : catchOpenMCError(err, "get ID for material with index " + std::to_string(index));
325 2581750 : 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 72 : OpenMCProblemBase::firstSolve() const
348 : {
349 72 : return _fixed_point_iteration < 0;
350 : }
351 :
352 : void
353 2484 : OpenMCProblemBase::externalSolve()
354 : {
355 4968 : 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 2484 : 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 2466 : _console << " Running OpenMC with " << nParticles() << " particles per batch..." << std::endl;
366 :
367 : // apply a new starting fission source
368 2466 : 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 2466 : executeEditors();
377 :
378 2460 : if (_reset_seed)
379 : {
380 56 : openmc_hard_reset();
381 56 : openmc_set_seed(_initial_seed);
382 : }
383 :
384 : int err;
385 2460 : if (_criticality_search)
386 52 : _criticality_search->searchForCriticality();
387 : else
388 : {
389 2408 : err = openmc_run();
390 2408 : if (err)
391 0 : mooseError(openmc_err_msg);
392 : }
393 :
394 2456 : _total_n_particles += nParticles();
395 :
396 2456 : err = openmc_reset_timers();
397 2456 : if (err)
398 0 : mooseError(openmc_err_msg);
399 :
400 2456 : _fixed_point_iteration++;
401 :
402 : // save the latest fission source for re-use in the next iteration
403 2456 : if (_reuse_source)
404 48 : writeSourceBank(sourceBankFileName());
405 2474 : }
406 :
407 : void
408 1920 : OpenMCProblemBase::initialSetup()
409 : {
410 1920 : CardinalProblem::initialSetup();
411 :
412 : // Initialize the IFP parameters tally.
413 1920 : if (_calc_kinetics_params)
414 : {
415 : // For \Lambda_eff, \beta_{eff}, and the denominator of \beta_{eff,i}
416 17 : _ifp_common_tally_index = openmc::model::tallies.size();
417 17 : _ifp_common_tally = openmc::Tally::create();
418 17 : _ifp_common_tally->set_scores({"ifp-time-numerator", "ifp-denominator", "ifp-beta-numerator"});
419 17 : _ifp_common_tally->estimator_ = openmc::TallyEstimator::COLLISION;
420 :
421 : // For \beta_{eff,i}. A separate tally is required when sieving by delayed group to compute
422 : // standard deviations and relative errors correctly for the total \beta_eff (due to covariances
423 : // between delayed groups).
424 17 : _ifp_mg_beta_tally_index = openmc::model::tallies.size();
425 17 : _ifp_mg_beta_tally = openmc::Tally::create();
426 17 : _ifp_mg_beta_tally->set_scores({"ifp-beta-numerator"});
427 17 : _ifp_mg_beta_tally->estimator_ = openmc::TallyEstimator::COLLISION;
428 :
429 : auto dnp_grp_filter =
430 17 : dynamic_cast<openmc::DelayedGroupFilter *>(openmc::Filter::create("delayedgroup"));
431 17 : std::vector<int> grps{1, 2, 3, 4, 5, 6};
432 17 : dnp_grp_filter->set_groups(openmc::span<int>(grps));
433 :
434 17 : std::vector<openmc::Filter *> df{dnp_grp_filter};
435 17 : _ifp_mg_beta_tally->set_filters({df});
436 17 : }
437 :
438 : // Find a criticality search object
439 1920 : TheWarehouse::Query query = theWarehouse().query().condition<AttribSystem>("CriticalitySearch");
440 : std::vector<CriticalitySearchBase *> objs;
441 : query.queryInto(objs);
442 :
443 1920 : if (objs.size() > 1)
444 0 : mooseError("Cannot have more than one CriticalitySearch object");
445 :
446 1920 : if (objs.size())
447 52 : _criticality_search = objs[0];
448 1920 : }
449 :
450 : void
451 4976 : OpenMCProblemBase::syncSolutions(ExternalProblem::Direction direction)
452 : {
453 : // Always run OpenMC on the first timestep in a steady solve with adaptivity. This
454 : // ensures that OpenMC runs at least once during each Picard iteration.
455 4976 : _run_on_adaptivity_cycle |= (timeStep() == 1 && !isTransient());
456 4976 : }
457 :
458 : bool
459 643 : OpenMCProblemBase::adaptMesh()
460 : {
461 643 : _run_on_adaptivity_cycle = CardinalProblem::adaptMesh() || isTransient();
462 643 : return _run_on_adaptivity_cycle;
463 : }
464 :
465 : void
466 24 : OpenMCProblemBase::writeSourceBank(const std::string & filename)
467 : {
468 24 : hid_t file_id = openmc::file_open(filename, 'w', true);
469 : openmc::write_attribute(file_id, "filetype", "source");
470 24 : openmc::write_attribute(file_id, "version", openmc::VERSION_STATEPOINT);
471 24 : openmc::write_source_bank(
472 : file_id, openmc::simulation::source_bank, openmc::simulation::work_index);
473 24 : openmc::file_close(file_id);
474 24 : }
475 :
476 : unsigned int
477 2417 : OpenMCProblemBase::numElemsInSubdomain(const SubdomainID & id) const
478 : {
479 2417 : unsigned int n = 0;
480 7076365 : for (unsigned int e = 0; e < _mesh.nElem(); ++e)
481 : {
482 7073948 : const auto * elem = _mesh.queryElemPtr(e);
483 :
484 7073948 : if (!isLocalElem(elem) || !elem->active())
485 3380064 : continue;
486 :
487 : const auto subdomain_id = elem->subdomain_id();
488 3693884 : if (id == subdomain_id)
489 1533864 : n += 1;
490 : }
491 :
492 2417 : _communicator.sum(n);
493 :
494 2417 : return n;
495 : }
496 :
497 : bool
498 22700252 : OpenMCProblemBase::isLocalElem(const Elem * elem) const
499 : {
500 22700252 : if (!elem)
501 : {
502 : // we should only not be able to find an element if the mesh is distributed
503 : libmesh_assert(!_mesh.getMesh().is_serial());
504 : return false;
505 : }
506 :
507 14484310 : if (elem->processor_id() == _communicator.rank())
508 11716100 : return true;
509 :
510 : return false;
511 : }
512 :
513 : bool
514 4 : OpenMCProblemBase::cellHasZeroInstances(const cellInfo & cell_info) const
515 : {
516 4 : auto n = openmc::model::cells.at(cell_info.first)->n_instances();
517 4 : return !n;
518 : }
519 :
520 : void
521 11393062 : OpenMCProblemBase::setCellTemperature(const int32_t & index,
522 : const int32_t & instance,
523 : const Real & T,
524 : const cellInfo & cell_info) const
525 : {
526 11393062 : int err = openmc_cell_set_temperature(index, T, &instance, false);
527 11393062 : if (err)
528 : {
529 : std::string descriptor =
530 12 : "set cell " + printCell(cell_info) + " to temperature " + Moose::stringify(T) + " (K)";
531 :
532 : // special error message if cell has zero instances
533 4 : if (cellHasZeroInstances(cell_info))
534 0 : mooseError("Failed to set the temperature for cell " + printCell(cell_info) +
535 : " with zero instances.");
536 :
537 12 : mooseError("In attempting to set cell " + printCell(cell_info) + " to temperature " +
538 4 : Moose::stringify(T) + " (K), OpenMC reported:\n\n",
539 4 : std::string(openmc_err_msg) + "\n\n" +
540 : "If you are trying to debug a model setup, you can set 'initial_properties = "
541 : "xml' to use the initial temperature and density in the OpenMC XML files for "
542 : "OpenMC's first run.");
543 : }
544 11393058 : }
545 :
546 : std::vector<int32_t>
547 441711 : OpenMCProblemBase::cellFill(const cellInfo & cell_info, int & fill_type) const
548 : {
549 441711 : int32_t * materials = nullptr;
550 441711 : int n_materials = 0;
551 :
552 441711 : int err = openmc_cell_get_fill(cell_info.first, &fill_type, &materials, &n_materials);
553 883422 : catchOpenMCError(err, "get fill of cell " + printCell(cell_info));
554 :
555 : std::vector<int32_t> material_indices;
556 441711 : material_indices.assign(materials, materials + n_materials);
557 441711 : return material_indices;
558 0 : }
559 :
560 : bool
561 441711 : OpenMCProblemBase::materialFill(const cellInfo & cell_info, int32_t & material_index) const
562 : {
563 : int fill_type;
564 441711 : auto material_indices = cellFill(cell_info, fill_type);
565 :
566 441711 : if (fill_type != static_cast<int>(openmc::Fill::MATERIAL))
567 : return false;
568 :
569 : // The number of materials in a cell is either 1, or equal to the number of instances
570 : // (if distributed materials were used).
571 441709 : if (material_indices.size() == 1)
572 438781 : material_index = material_indices[0];
573 : else
574 2928 : material_index = material_indices[cell_info.second];
575 :
576 : return true;
577 441711 : }
578 :
579 : void
580 1522 : OpenMCProblemBase::setCellDensity(const Real & density, const cellInfo & cell_info) const
581 : {
582 : // OpenMC technically allows a density of >= 0.0, but we can impose a tighter
583 : // check here with a better error message than the Excepts() in material->set_density
584 : // because it could be a very common mistake to forget to set an initial condition
585 : // for density if OpenMC runs first
586 1522 : if (density <= 0.0)
587 4 : mooseError("Densities less than or equal to zero cannot be set in the OpenMC model!\n\n cell " +
588 4 : printCell(cell_info) + " set to density " + Moose::stringify(density) + " (kg/m3)");
589 :
590 : int32_t material_index;
591 1520 : auto is_material_cell = materialFill(cell_info, material_index);
592 :
593 1520 : if (!is_material_cell)
594 0 : mooseError(
595 : "Density transfer does not currently support cells filled with universes or lattices!");
596 :
597 : // throw a special error if the cell is void, because the OpenMC error isn't very
598 : // clear what the mistake is
599 1520 : if (material_index == MATERIAL_VOID)
600 : {
601 12 : mooseWarning("Skipping setting density for cell " + printCell(cell_info) +
602 : " because this cell is void (vacuum)");
603 4 : return;
604 : }
605 :
606 : // Compute the density. We multiply density by 0.001 to convert from kg/m3
607 : // (the units assumed in the 'density' auxvariable as well as the MOOSE fluid
608 : // properties module) to g/cm3
609 1514 : int err = openmc_cell_set_density(
610 1514 : cell_info.first, _density_conversion_factor * density, &cell_info.second, false);
611 :
612 1514 : if (err)
613 : {
614 : // special error message if cell has zero instances
615 0 : if (cellHasZeroInstances(cell_info))
616 0 : mooseError("Failed to set the density for cell " + printCell(cell_info) +
617 : " with zero instances.");
618 :
619 0 : mooseError("In attempting to set cell " + printCell(cell_info) + " to density " +
620 0 : Moose::stringify(density) + " (kg/m3), OpenMC reported:\n\n",
621 0 : std::string(openmc_err_msg) + "\n\n" +
622 : "If you are trying to debug a model setup, you can set 'initial_properties = "
623 : "xml' to use the initial temperature and density in the OpenMC XML files for "
624 : "OpenMC's first run.");
625 : }
626 : }
627 :
628 : std::string
629 4539874 : OpenMCProblemBase::printCell(const cellInfo & cell_info, const bool brief) const
630 : {
631 4539874 : int32_t id = cellID(cell_info.first);
632 :
633 4539874 : std::stringstream msg;
634 4539874 : if (!brief)
635 4526955 : msg << "id ";
636 :
637 9079748 : msg << std::setw(_n_cell_digits) << Moose::stringify(id) << ", instance "
638 9079748 : << std::setw(_n_cell_digits) << Moose::stringify(cell_info.second) << " (of "
639 4539874 : << std::setw(_n_cell_digits)
640 18159496 : << Moose::stringify(openmc::model::cells.at(cell_info.first)->n_instances()) << ")";
641 :
642 4539874 : return msg.str();
643 4539874 : }
644 :
645 : void
646 2 : OpenMCProblemBase::importProperties() const
647 : {
648 2 : _console << "Reading temperature and density from properties.h5" << std::endl;
649 :
650 2 : int err = openmc_properties_import("properties.h5");
651 2 : catchOpenMCError(err, "load temperature and density from a properties.h5 file");
652 0 : }
653 :
654 : OMCTensor
655 4150 : OpenMCProblemBase::relativeError(const OMCTensor & sum,
656 : const OMCTensor & sum_sq,
657 : const int & n_realizations) const
658 : {
659 4150 : auto rel_err = openmc::tensor::zeros<double>({sum.size()});
660 :
661 478000 : for (unsigned int i = 0; i < sum.size(); ++i)
662 : {
663 473850 : auto mean = sum(i) / n_realizations;
664 473850 : auto std_dev = std::sqrt((sum_sq(i) / n_realizations - mean * mean) / (n_realizations - 1));
665 473850 : rel_err[i] = mean != 0.0 ? std_dev / std::abs(mean) : 0.0;
666 : }
667 :
668 4150 : return rel_err;
669 : }
670 :
671 : Real
672 456 : OpenMCProblemBase::relativeError(const Real & sum,
673 : const Real & sum_sq,
674 : const int & n_realizations) const
675 : {
676 456 : auto mean = sum / n_realizations;
677 456 : auto std_dev = std::sqrt((sum_sq / n_realizations - mean * mean) / (n_realizations - 1));
678 456 : return mean != 0.0 ? std_dev / std::abs(mean) : 0.0;
679 : }
680 :
681 : OMCTensor
682 5956 : OpenMCProblemBase::tallySum(const openmc::Tally * tally, const unsigned int & score) const
683 : {
684 5956 : return OMCTensor(tally->results_.slice(
685 11912 : openmc::tensor::all, score, static_cast<int>(openmc::TallyResult::SUM)));
686 : }
687 :
688 : double
689 2186 : OpenMCProblemBase::tallySumAcrossBins(std::vector<const openmc::Tally *> tally,
690 : const unsigned int & score) const
691 : {
692 : double sum = 0.0;
693 :
694 4372 : for (const auto & t : tally)
695 : {
696 2186 : auto mean = tallySum(t, score);
697 2186 : sum += mean.sum();
698 : }
699 :
700 2186 : return sum;
701 : }
702 :
703 : double
704 0 : OpenMCProblemBase::tallyMeanAcrossBins(std::vector<const openmc::Tally *> tally,
705 : const unsigned int & score) const
706 : {
707 : int n = 0;
708 0 : for (const auto & t : tally)
709 0 : n += t->n_realizations_;
710 :
711 0 : return tallySumAcrossBins(tally, score) / n;
712 : }
713 :
714 : std::string
715 2531 : OpenMCProblemBase::enumToTallyScore(const std::string & score) const
716 : {
717 : // the MultiMooseEnum is all caps, but the MooseEnum is already the correct case,
718 : // so we need to treat these as separate
719 2531 : std::string s = score;
720 2531 : if (std::all_of(
721 18731 : s.begin(), s.end(), [](unsigned char c) { return !std::isalpha(c) || std::isupper(c); }))
722 : {
723 20258 : std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
724 :
725 : // we need to revert back to some letters being uppercase for certain scores
726 2029 : if (s == "h3_production")
727 : s = "H3_production";
728 : }
729 :
730 : // MOOSE enums use underscores, OpenMC uses dashes
731 : std::replace(s.begin(), s.end(), '_', '-');
732 2531 : return s;
733 : }
734 :
735 : std::string
736 0 : OpenMCProblemBase::tallyScoreToEnum(const std::string & score) const
737 : {
738 : // MOOSE enums use underscores, OpenMC uses dashes
739 0 : std::string s = score;
740 : std::replace(s.begin(), s.end(), '-', '_');
741 0 : return s;
742 : }
743 :
744 : openmc::TallyEstimator
745 306 : OpenMCProblemBase::tallyEstimator(tally::TallyEstimatorEnum estimator) const
746 : {
747 : switch (estimator)
748 : {
749 : case tally::tracklength:
750 : return openmc::TallyEstimator::TRACKLENGTH;
751 : case tally::collision:
752 : return openmc::TallyEstimator::COLLISION;
753 : case tally::analog:
754 : return openmc::TallyEstimator::ANALOG;
755 0 : default:
756 0 : mooseError("Unhandled TallyEstimatorEnum!");
757 : }
758 : }
759 :
760 : std::string
761 0 : OpenMCProblemBase::estimatorToString(openmc::TallyEstimator estimator) const
762 : {
763 0 : switch (estimator)
764 : {
765 : case openmc::TallyEstimator::TRACKLENGTH:
766 0 : return "tracklength";
767 : case openmc::TallyEstimator::COLLISION:
768 0 : return "collision";
769 : case openmc::TallyEstimator::ANALOG:
770 0 : return "analog";
771 0 : default:
772 0 : mooseError("Unhandled TallyEstimatorEnum!");
773 : }
774 : }
775 :
776 : openmc::TriggerMetric
777 128 : OpenMCProblemBase::triggerMetric(std::string trigger) const
778 : {
779 128 : if (trigger == "variance")
780 : return openmc::TriggerMetric::variance;
781 128 : else if (trigger == "std_dev")
782 : return openmc::TriggerMetric::standard_deviation;
783 128 : else if (trigger == "rel_err")
784 : return openmc::TriggerMetric::relative_error;
785 0 : else if (trigger == "none")
786 : return openmc::TriggerMetric::not_active;
787 : else
788 0 : mooseError("Unhandled TallyTriggerTypeEnum: ", trigger);
789 : }
790 :
791 : openmc::TriggerMetric
792 2109 : OpenMCProblemBase::triggerMetric(trigger::TallyTriggerTypeEnum trigger) const
793 : {
794 : switch (trigger)
795 : {
796 : case trigger::variance:
797 : return openmc::TriggerMetric::variance;
798 : case trigger::std_dev:
799 : return openmc::TriggerMetric::standard_deviation;
800 : case trigger::rel_err:
801 : return openmc::TriggerMetric::relative_error;
802 : case trigger::none:
803 : return openmc::TriggerMetric::not_active;
804 0 : default:
805 0 : mooseError("Unhandled TallyTriggerTypeEnum!");
806 : }
807 : }
808 :
809 : bool
810 0 : OpenMCProblemBase::cellIsVoid(const cellInfo & cell_info) const
811 : {
812 : // material_index will be unchanged if the cell is filled by a universe or lattice.
813 : // Otherwise, this will get set to the material index in the cell.
814 0 : int32_t material_index = 0;
815 0 : materialFill(cell_info, material_index);
816 0 : return material_index == MATERIAL_VOID;
817 : }
818 :
819 : void
820 1083 : OpenMCProblemBase::geometryType(bool & has_csg_universe, bool & has_dag_universe) const
821 : {
822 1083 : has_csg_universe = false;
823 1083 : has_dag_universe = false;
824 :
825 : // Loop over universes and check if type is DAGMC
826 6157 : for (const auto & universe : openmc::model::universes)
827 : {
828 5074 : if (universe->geom_type() == openmc::GeometryType::DAG)
829 48 : has_dag_universe = true;
830 5026 : else if (universe->geom_type() == openmc::GeometryType::CSG)
831 5026 : has_csg_universe = true;
832 : else
833 0 : mooseError("Unhandled GeometryType enum!");
834 : }
835 1083 : }
836 :
837 : long unsigned int
838 2143 : OpenMCProblemBase::numCells() const
839 : {
840 : long unsigned int n_openmc_cells = 0;
841 90371 : for (const auto & c : openmc::model::cells)
842 88228 : n_openmc_cells += c->n_instances();
843 :
844 2143 : return n_openmc_cells;
845 : }
846 :
847 : const openmc::Tally &
848 228 : OpenMCProblemBase::getCommonKineticsTally()
849 : {
850 228 : if (!_ifp_common_tally)
851 0 : mooseError("Internal error: kinetics parameters have not been enabled.");
852 :
853 228 : return *_ifp_common_tally;
854 : }
855 :
856 : const openmc::Tally &
857 198 : OpenMCProblemBase::getMGBetaTally()
858 : {
859 198 : return *_ifp_mg_beta_tally;
860 : }
861 :
862 : bool
863 189054 : OpenMCProblemBase::isReactionRateScore(const std::string & score) const
864 : {
865 : const std::set<std::string> viable_scores = {"H3-production",
866 : "total",
867 : "absorption",
868 : "scatter",
869 : "nu-scatter",
870 : "fission",
871 : "nu-fission",
872 : "prompt-nu-fission",
873 189054 : "delayed-nu-fission"};
874 189054 : return viable_scores.count(score);
875 : }
876 :
877 : bool
878 475637 : OpenMCProblemBase::isHeatingScore(const std::string & score) const
879 : {
880 : const std::set<std::string> viable_scores = {
881 475637 : "heating", "heating-local", "kappa-fission", "fission-q-prompt", "fission-q-recoverable"};
882 475637 : return viable_scores.count(score);
883 : }
884 :
885 : unsigned int
886 9299 : OpenMCProblemBase::addExternalVariable(const std::string & name,
887 : const std::string & system,
888 : const std::vector<SubdomainName> * block)
889 : {
890 9299 : auto var_params = _factory.getValidParams("MooseVariable");
891 18598 : var_params.set<MooseEnum>("family") = "MONOMIAL";
892 18598 : var_params.set<MooseEnum>("order") = "CONSTANT";
893 :
894 9299 : if (block)
895 11950 : var_params.set<std::vector<SubdomainName>>("block") = *block;
896 :
897 9299 : checkDuplicateVariableName(name, system);
898 18590 : addAuxVariable("MooseVariable", name, var_params);
899 18590 : return _aux->getFieldVariable<Real>(0, name).number();
900 9295 : }
901 :
902 : std::string
903 5224 : OpenMCProblemBase::subdomainName(const SubdomainID & id) const
904 : {
905 5224 : std::string name = _mesh.getSubdomainName(id);
906 5224 : if (name.empty())
907 10004 : name = std::to_string(id);
908 5224 : return name;
909 : }
910 :
911 : void
912 1920 : OpenMCProblemBase::getOpenMCUserObjects()
913 : {
914 1920 : _cell_transform_uos.clear();
915 :
916 1920 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("UserObject");
917 : std::vector<UserObject *> userobjs;
918 : uo_query.queryInto(userobjs);
919 :
920 9203 : for (const auto & u : userobjs)
921 : {
922 7283 : OpenMCNuclideDensities * c = dynamic_cast<OpenMCNuclideDensities *>(u);
923 7283 : if (c)
924 44 : _nuclide_densities_uos.push_back(c);
925 :
926 7283 : OpenMCTallyEditor * e = dynamic_cast<OpenMCTallyEditor *>(u);
927 7283 : if (e)
928 60 : _tally_editor_uos.push_back(e);
929 :
930 7283 : OpenMCDomainFilterEditor * f = dynamic_cast<OpenMCDomainFilterEditor *>(u);
931 7283 : if (f)
932 28 : _filter_editor_uos.push_back(f);
933 :
934 7283 : OpenMCCellTransform * t = dynamic_cast<OpenMCCellTransform *>(u);
935 7283 : if (t)
936 54 : _cell_transform_uos.push_back(t);
937 : }
938 :
939 1920 : checkOpenMCUserObjectIDs();
940 1916 : }
941 :
942 : bool
943 1967 : OpenMCProblemBase::hasCellTransform() const
944 : {
945 1967 : return !_cell_transform_uos.empty();
946 : }
947 :
948 : void
949 1920 : OpenMCProblemBase::checkOpenMCUserObjectIDs() const
950 : {
951 : std::set<int32_t> tally_ids;
952 1978 : for (const auto & te : _tally_editor_uos)
953 : {
954 60 : int32_t tally_id = te->tallyId();
955 : if (tally_ids.count(tally_id) != 0)
956 2 : te->duplicateTallyError(tally_id);
957 58 : tally_ids.insert(tally_id);
958 : }
959 :
960 : std::set<int32_t> filter_ids;
961 1944 : for (const auto & fe : _filter_editor_uos)
962 : {
963 28 : int32_t filter_id = fe->filterId();
964 : if (filter_ids.count(filter_id) != 0)
965 2 : fe->duplicateFilterError(filter_id);
966 26 : filter_ids.insert(filter_id);
967 : }
968 1916 : }
969 :
970 : void
971 1750 : OpenMCProblemBase::checkTallyEditorIDs() const
972 : {
973 1750 : std::vector<int32_t> mapped_tally_ids = getMappedTallyIDs();
974 :
975 1804 : for (const auto & te : _tally_editor_uos)
976 : {
977 56 : int32_t tally_id = te->tallyId();
978 :
979 : // ensure that the TallyEditor IDs don't apply to any mapped tally objects
980 56 : if (std::find(mapped_tally_ids.begin(), mapped_tally_ids.end(), tally_id) !=
981 : mapped_tally_ids.end())
982 2 : te->mappedTallyError(tally_id);
983 : }
984 1748 : }
985 :
986 : void
987 2466 : OpenMCProblemBase::executeFilterEditors()
988 : {
989 2466 : executeControls(EXEC_FILTER_EDITORS);
990 :
991 2466 : if (!_filter_editor_uos.size())
992 : return;
993 :
994 24 : _console << "Executing filter editors..." << std::endl;
995 48 : for (const auto & fe : _filter_editor_uos)
996 24 : fe->execute();
997 : }
998 :
999 : void
1000 2466 : OpenMCProblemBase::executeTallyEditors()
1001 : {
1002 2466 : executeControls(EXEC_TALLY_EDITORS);
1003 :
1004 2466 : if (!_tally_editor_uos.size())
1005 : return;
1006 :
1007 54 : _console << "Executing tally editors..." << std::endl;
1008 102 : for (const auto & te : _tally_editor_uos)
1009 54 : te->execute();
1010 : }
1011 :
1012 : void
1013 2466 : OpenMCProblemBase::executeEditors()
1014 : {
1015 2466 : executeFilterEditors();
1016 2466 : executeTallyEditors();
1017 2460 : }
1018 :
1019 : void
1020 2484 : OpenMCProblemBase::sendNuclideDensitiesToOpenMC()
1021 : {
1022 2484 : if (_nuclide_densities_uos.size() == 0)
1023 : return;
1024 :
1025 : // We could probably put this somewhere better, but it's good for now
1026 44 : executeControls(EXEC_SEND_OPENMC_DENSITIES);
1027 :
1028 44 : _console << "Sending nuclide compositions to OpenMC... ";
1029 84 : for (const auto & uo : _nuclide_densities_uos)
1030 44 : uo->setValue();
1031 : }
1032 :
1033 : #endif
|
