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