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 4761 : OpenMCProblemBase::validParams()
44 : {
45 4761 : InputParameters params = CardinalProblem::validParams();
46 9522 : params.addParam<PostprocessorName>(
47 : "power", "Power (Watts) to normalize the OpenMC tallies; only used for k-eigenvalue mode");
48 9522 : params.addParam<PostprocessorName>(
49 : "source_strength",
50 : "Neutrons/second to normalize the OpenMC tallies; only used for fixed source mode");
51 9522 : params.addParam<bool>("verbose", false, "Whether to print diagnostic information");
52 :
53 9522 : params.addParam<MooseEnum>("tally_type", getTallyTypeEnum(), "Type of tally to use in OpenMC");
54 :
55 14283 : params.addRangeCheckedParam<Real>(
56 : "scaling",
57 9522 : 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 9522 : 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 9522 : 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 9522 : params.addParam<PostprocessorName>("particles",
74 : "Number of particles to run in each OpenMC batch; this "
75 : "overrides the setting in the XML files.");
76 9522 : 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 9522 : params.addParam<bool>("reuse_source",
82 9522 : false,
83 : "Whether to take the initial fission source "
84 : "for interation n to be the converged source bank from iteration n-1");
85 9522 : params.addParam<bool>(
86 : "skip_statepoint",
87 9522 : 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 9522 : params.addParam<bool>(
91 : "reset_seed",
92 9522 : false,
93 : "Whether to reset OpenMC's seed to the initial starting seed before each OpenMC solve");
94 :
95 9522 : params.addParam<FileName>(
96 : "xml_directory", "./", "The directory in which to look for OpenMC XML files.");
97 :
98 : // Kinetics parameters.
99 9522 : params.addParam<bool>("calc_kinetics_params",
100 9522 : false,
101 : "Whether or not Cardinal should enable the calculation of kinetics "
102 : "parameters (Lambda effective and beta effective).");
103 9522 : 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 9522 : 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 9522 : 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 4761 : return params;
120 0 : }
121 :
122 2399 : OpenMCProblemBase::OpenMCProblemBase(const InputParameters & params)
123 : : CardinalProblem(params),
124 : PostprocessorInterface(this),
125 2399 : _verbose(getParam<bool>("verbose")),
126 4798 : _reuse_source(getParam<bool>("reuse_source")),
127 2399 : _specified_scaling(params.isParamSetByUser("scaling")),
128 4798 : _scaling(getParam<Real>("scaling")),
129 4798 : _skip_statepoint(getParam<bool>("skip_statepoint")),
130 2399 : _fixed_point_iteration(-1),
131 2399 : _total_n_particles(0),
132 2399 : _has_adaptivity(getMooseApp().actionWarehouse().hasActions("set_adaptivity_options")),
133 2399 : _run_on_adaptivity_cycle(true),
134 4798 : _calc_kinetics_params(getParam<bool>("calc_kinetics_params")),
135 4798 : _reset_seed(getParam<bool>("reset_seed")),
136 2399 : _initial_seed(openmc::openmc_get_seed()),
137 7197 : _xml_directory(getParam<FileName>("xml_directory"))
138 : {
139 4798 : 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 2399 : 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 2399 : _run_mode = openmc::settings::run_mode;
155 2399 : const auto & tally_actions = getMooseApp().actionWarehouse().getActions<AddTallyAction>();
156 2399 : const auto & mgxs_actions = getMooseApp().actionWarehouse().getActions<SetupMGXSAction>();
157 2399 : 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 2245 : if (tally_actions.size() > 0 || mgxs_actions.size() > 0)
163 : {
164 3472 : checkRequiredParam(params, "power", "running in k-eigenvalue mode");
165 1736 : _power = &getPostprocessorValue("power");
166 : }
167 : else
168 1018 : checkUnusedParam(params, "power", "no tallies have been added");
169 :
170 4490 : checkUnusedParam(params, "source_strength", "running in k-eigenvalue mode");
171 2245 : 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 2393 : _n_cell_digits = std::to_string(openmc::model::cells.size()).length();
202 :
203 2393 : if (openmc::settings::libmesh_comm)
204 0 : mooseWarning("libMesh communicator already set in OpenMC.");
205 :
206 2393 : openmc::settings::libmesh_comm = &_mesh.comm();
207 :
208 4786 : if (isParamValid("openmc_verbosity"))
209 0 : openmc::settings::verbosity = getParam<unsigned int>("openmc_verbosity");
210 :
211 4786 : if (isParamValid("inactive_batches"))
212 250 : openmc::settings::n_inactive = getParam<unsigned int>("inactive_batches");
213 :
214 4786 : if (isParamValid("particles"))
215 161 : _particles = &getPostprocessorValue("particles");
216 :
217 2393 : if (!runRandomRay())
218 : {
219 4524 : checkUnusedParam(params, "inactive_distance", "not running in random ray mode");
220 4524 : checkUnusedParam(params, "active_distance", "not running in random ray mode");
221 : }
222 :
223 4786 : if (isParamValid("inactive_distance"))
224 64 : openmc::RandomRay::distance_inactive_ = getParam<Real>("inactive_distance");
225 :
226 4786 : if (isParamValid("active_distance"))
227 64 : openmc::RandomRay::distance_active_ = getParam<Real>("active_distance");
228 :
229 4786 : 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 2391 : 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 2385 : }
270 :
271 2015 : 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 5496 : OpenMCProblemBase::nParticles() const
296 : {
297 5496 : return openmc::settings::n_particles;
298 : }
299 :
300 : std::string
301 7982 : 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 7982 : if (index == -1)
306 51 : return "VOID";
307 :
308 : const char * name;
309 7931 : int err = openmc_material_get_name(index, &name);
310 7931 : catchOpenMCError(err, "get material name for material with index " + std::to_string(index));
311 :
312 7931 : std::string n = name;
313 :
314 : // if the material does not have a name, just return the ID instead
315 7931 : if (n.empty())
316 7460 : n = std::to_string(materialID(index));
317 :
318 7931 : return n;
319 : }
320 :
321 : int32_t
322 20133831 : OpenMCProblemBase::cellID(const int32_t index) const
323 : {
324 : int32_t id;
325 20133831 : int err = openmc_cell_get_id(index, &id);
326 20133831 : catchOpenMCError(err, "get ID for cell with index " + std::to_string(index));
327 20133831 : return id;
328 : }
329 :
330 : int32_t
331 1303570 : OpenMCProblemBase::materialID(const int32_t index) const
332 : {
333 1303570 : if (index == openmc::MATERIAL_VOID)
334 : return -1;
335 :
336 : int32_t id;
337 1297666 : int err = openmc_material_get_id(index, &id);
338 1297666 : catchOpenMCError(err, "get ID for material with index " + std::to_string(index));
339 1297666 : 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 2795 : OpenMCProblemBase::firstSolve() const
362 : {
363 2795 : return _fixed_point_iteration < 0;
364 : }
365 :
366 : void
367 2747 : OpenMCProblemBase::externalSolve()
368 : {
369 5494 : 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 2747 : 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 2729 : _console << " Running OpenMC with " << nParticles() << " particles per batch..." << std::endl;
380 :
381 : // apply a new starting fission source
382 2729 : 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 2729 : executeEditors();
391 :
392 2723 : if (_reset_seed)
393 : {
394 56 : openmc_hard_reset();
395 56 : openmc_set_seed(_initial_seed);
396 : }
397 :
398 : int err = 0;
399 2723 : if (!firstSolve())
400 : {
401 694 : err = openmc_reset_timers();
402 694 : if (err)
403 0 : mooseError(openmc_err_msg);
404 : }
405 :
406 2723 : if (_criticality_search)
407 960 : _criticality_search->searchForCriticality([&]() { this->critSearchStep(); });
408 : else
409 : {
410 2647 : if (runRandomRay())
411 143 : openmc_run_random_ray();
412 : else
413 2504 : err = openmc_run();
414 :
415 2647 : if (err)
416 0 : mooseError(openmc_err_msg);
417 : }
418 :
419 2719 : _total_n_particles += nParticles();
420 :
421 2719 : _fixed_point_iteration++;
422 :
423 : // save the latest fission source for re-use in the next iteration
424 2719 : if (_reuse_source)
425 48 : writeSourceBank(sourceBankFileName());
426 2737 : }
427 :
428 : void
429 2128 : OpenMCProblemBase::initialSetup()
430 : {
431 2128 : CardinalProblem::initialSetup();
432 :
433 : // Initialize the IFP parameters tally.
434 2128 : 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 2128 : TheWarehouse::Query query = theWarehouse().query().condition<AttribSystem>("CriticalitySearch");
461 : std::vector<CriticalitySearchBase *> objs;
462 : query.queryInto(objs);
463 :
464 2128 : if (objs.size() > 1)
465 0 : mooseError("Cannot have more than one CriticalitySearch object");
466 :
467 2128 : if (objs.size())
468 76 : _criticality_search = objs[0];
469 :
470 : // Find model modifier objects
471 2128 : TheWarehouse::Query mm_query = theWarehouse().query().condition<AttribSystem>("ModelModifiers");
472 : std::vector<ModelModifiersBase *> mm_objs;
473 : mm_query.queryInto(mm_objs);
474 2160 : for (const auto & m : mm_objs)
475 32 : m->modifyOpenMCModel();
476 2128 : }
477 :
478 : void
479 5502 : 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 5502 : _run_on_adaptivity_cycle |= (timeStep() == 1 && !isTransient());
484 5502 : }
485 :
486 : bool
487 704 : OpenMCProblemBase::adaptMesh()
488 : {
489 704 : _run_on_adaptivity_cycle = CardinalProblem::adaptMesh() || isTransient();
490 704 : 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 2751 : OpenMCProblemBase::numElemsInSubdomain(const SubdomainID & id) const
506 : {
507 2751 : unsigned int n = 0;
508 7179057 : for (unsigned int e = 0; e < _mesh.nElem(); ++e)
509 : {
510 7176306 : const auto * elem = _mesh.queryElemPtr(e);
511 :
512 7176306 : if (!isLocalElem(elem) || !elem->active())
513 3390864 : continue;
514 :
515 : const auto subdomain_id = elem->subdomain_id();
516 3785442 : if (id == subdomain_id)
517 1549656 : n += 1;
518 : }
519 :
520 2751 : _communicator.sum(n);
521 :
522 2751 : return n;
523 : }
524 :
525 : bool
526 22962160 : OpenMCProblemBase::isLocalElem(const Elem * elem) const
527 : {
528 22962160 : 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 14714078 : if (elem->processor_id() == _communicator.rank())
536 11868880 : 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 11393840 : OpenMCProblemBase::setCellTemperature(const int32_t & index,
550 : const int32_t & instance,
551 : const Real & T,
552 : const cellInfo & cell_info) const
553 : {
554 11393840 : int err = openmc_cell_set_temperature(index, T, &instance, false);
555 11393840 : 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 11393836 : }
573 :
574 : std::vector<int32_t>
575 1740161 : OpenMCProblemBase::cellFill(const cellInfo & cell_info, int & fill_type) const
576 : {
577 1740161 : int32_t * materials = nullptr;
578 1740161 : int n_materials = 0;
579 :
580 1740161 : int err = openmc_cell_get_fill(cell_info.first, &fill_type, &materials, &n_materials);
581 3480322 : catchOpenMCError(err, "get fill of cell " + printCell(cell_info));
582 :
583 : std::vector<int32_t> material_indices;
584 1740161 : material_indices.assign(materials, materials + n_materials);
585 1740161 : return material_indices;
586 0 : }
587 :
588 : bool
589 1740161 : OpenMCProblemBase::materialFill(const cellInfo & cell_info, int32_t & material_index) const
590 : {
591 : int fill_type;
592 1740161 : auto material_indices = cellFill(cell_info, fill_type);
593 :
594 1740161 : 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 1740161 : if (material_indices.size() == 1)
600 1331009 : material_index = material_indices[0];
601 : else
602 409152 : material_index = material_indices[cell_info.second];
603 :
604 : return true;
605 1740161 : }
606 :
607 : const Real
608 1335860 : 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 1335860 : return openmc::settings::run_CE ? _density_conversion_factor : 1.0;
621 : }
622 :
623 : void
624 2650 : OpenMCProblemBase::setCellDensity(const int32_t & index,
625 : const int32_t & instance,
626 : const Real & density,
627 : const cellInfo & cell_info) const
628 : {
629 : // OpenMC technically allows a density of >= 0.0, but we can impose a tighter
630 : // check here with a better error message than the Excepts() in material->set_density
631 : // because it could be a very common mistake to forget to set an initial condition
632 : // for density if OpenMC runs first
633 2650 : if (density <= 0.0)
634 4 : mooseError("Densities less than or equal to zero cannot be set in the OpenMC model!\n\n cell " +
635 4 : printCell(cell_info) + " set to density " + Moose::stringify(density) + " (kg/m3)");
636 :
637 : // Compute the density. We multiply density by 0.001 to convert from kg/m3
638 : // (the units assumed in the 'density' auxvariable as well as the MOOSE fluid
639 : // properties module) to g/cm3
640 2648 : int err = openmc_cell_set_density(index, densityConversionFactor() * density, &instance, false);
641 :
642 2648 : if (err)
643 : {
644 : // special error message if cell has zero instances
645 0 : if (cellHasZeroInstances(cell_info))
646 0 : mooseError("Failed to set the density for cell " + printCell(cell_info) +
647 : " with zero instances.");
648 :
649 0 : mooseError("In attempting to set cell " + printCell(cell_info) + " to density " +
650 0 : Moose::stringify(density) + " (kg/m3), OpenMC reported:\n\n",
651 0 : std::string(openmc_err_msg) + "\n\n" +
652 : "If you are trying to debug a model setup, you can set 'initial_properties = "
653 : "xml' to use the initial temperature and density in the OpenMC XML files for "
654 : "OpenMC's first run.");
655 : }
656 2648 : }
657 :
658 : std::string
659 7178996 : OpenMCProblemBase::printCell(const cellInfo & cell_info, const bool brief) const
660 : {
661 7178996 : int32_t id = cellID(cell_info.first);
662 :
663 7178996 : std::stringstream msg;
664 7178996 : if (!brief)
665 7164297 : msg << "id ";
666 :
667 14357992 : msg << std::setw(_n_cell_digits) << Moose::stringify(id) << ", instance "
668 14357992 : << std::setw(_n_cell_digits) << Moose::stringify(cell_info.second) << " (of "
669 7178996 : << std::setw(_n_cell_digits)
670 28715984 : << Moose::stringify(openmc::model::cells.at(cell_info.first)->n_instances()) << ")";
671 :
672 7178996 : return msg.str();
673 7178996 : }
674 :
675 : void
676 2 : OpenMCProblemBase::importProperties() const
677 : {
678 2 : _console << "Reading temperature and density from properties.h5" << std::endl;
679 :
680 2 : int err = openmc_properties_import("properties.h5");
681 2 : catchOpenMCError(err, "load temperature and density from a properties.h5 file");
682 0 : }
683 :
684 : OMCTensor
685 4416 : OpenMCProblemBase::relativeError(const OMCTensor & sum,
686 : const OMCTensor & sum_sq,
687 : const int & n_realizations) const
688 : {
689 4416 : auto rel_err = openmc::tensor::zeros<double>({sum.size()});
690 :
691 480270 : for (unsigned int i = 0; i < sum.size(); ++i)
692 : {
693 475854 : auto mean = sum(i) / n_realizations;
694 475854 : auto std_dev = std::sqrt((sum_sq(i) / n_realizations - mean * mean) / (n_realizations - 1));
695 475854 : rel_err[i] = mean != 0.0 ? std_dev / std::abs(mean) : 0.0;
696 : }
697 :
698 4416 : return rel_err;
699 : }
700 :
701 : Real
702 456 : OpenMCProblemBase::relativeError(const Real & sum,
703 : const Real & sum_sq,
704 : const int & n_realizations) const
705 : {
706 456 : auto mean = sum / n_realizations;
707 456 : auto std_dev = std::sqrt((sum_sq / n_realizations - mean * mean) / (n_realizations - 1));
708 456 : return mean != 0.0 ? std_dev / std::abs(mean) : 0.0;
709 : }
710 :
711 : OMCTensor
712 6372 : OpenMCProblemBase::tallySum(const openmc::Tally * tally, const unsigned int & score) const
713 : {
714 6372 : return OMCTensor(tally->results_.slice(
715 12744 : openmc::tensor::all, score, static_cast<int>(openmc::TallyResult::SUM)));
716 : }
717 :
718 : double
719 2338 : OpenMCProblemBase::tallySumAcrossBins(std::vector<const openmc::Tally *> tally,
720 : const unsigned int & score) const
721 : {
722 : double sum = 0.0;
723 :
724 4676 : for (const auto & t : tally)
725 : {
726 2338 : auto mean = tallySum(t, score);
727 2338 : sum += mean.sum();
728 : }
729 :
730 2338 : return sum;
731 : }
732 :
733 : double
734 0 : OpenMCProblemBase::tallyMeanAcrossBins(std::vector<const openmc::Tally *> tally,
735 : const unsigned int & score) const
736 : {
737 : int n = 0;
738 0 : for (const auto & t : tally)
739 0 : n += t->n_realizations_;
740 :
741 0 : return tallySumAcrossBins(tally, score) / n;
742 : }
743 :
744 : bool
745 15594 : OpenMCProblemBase::runRandomRay() const
746 : {
747 15594 : return openmc::settings::solver_type == openmc::SolverType::RANDOM_RAY;
748 : }
749 :
750 : std::string
751 2763 : OpenMCProblemBase::enumToTallyScore(const std::string & score) const
752 : {
753 : // the MultiMooseEnum is all caps, but the MooseEnum is already the correct case,
754 : // so we need to treat these as separate
755 2763 : std::string s = score;
756 2763 : if (std::all_of(
757 20372 : s.begin(), s.end(), [](unsigned char c) { return !std::isalpha(c) || std::isupper(c); }))
758 : {
759 22077 : std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
760 :
761 : // we need to revert back to some letters being uppercase for certain scores
762 2234 : if (s == "h3_production")
763 : s = "H3_production";
764 : }
765 :
766 : // MOOSE enums use underscores, OpenMC uses dashes
767 : std::replace(s.begin(), s.end(), '_', '-');
768 2763 : return s;
769 : }
770 :
771 : std::string
772 0 : OpenMCProblemBase::tallyScoreToEnum(const std::string & score) const
773 : {
774 : // MOOSE enums use underscores, OpenMC uses dashes
775 0 : std::string s = score;
776 : std::replace(s.begin(), s.end(), '-', '_');
777 0 : return s;
778 : }
779 :
780 : openmc::TallyEstimator
781 306 : OpenMCProblemBase::tallyEstimator(tally::TallyEstimatorEnum estimator) const
782 : {
783 : switch (estimator)
784 : {
785 : case tally::tracklength:
786 : return openmc::TallyEstimator::TRACKLENGTH;
787 : case tally::collision:
788 : return openmc::TallyEstimator::COLLISION;
789 : case tally::analog:
790 : return openmc::TallyEstimator::ANALOG;
791 0 : default:
792 0 : mooseError("Unhandled TallyEstimatorEnum!");
793 : }
794 : }
795 :
796 : std::string
797 0 : OpenMCProblemBase::estimatorToString(openmc::TallyEstimator estimator) const
798 : {
799 0 : switch (estimator)
800 : {
801 : case openmc::TallyEstimator::TRACKLENGTH:
802 0 : return "tracklength";
803 : case openmc::TallyEstimator::COLLISION:
804 0 : return "collision";
805 : case openmc::TallyEstimator::ANALOG:
806 0 : return "analog";
807 0 : default:
808 0 : mooseError("Unhandled TallyEstimatorEnum!");
809 : }
810 : }
811 :
812 : openmc::TriggerMetric
813 128 : OpenMCProblemBase::triggerMetric(std::string trigger) const
814 : {
815 128 : if (trigger == "variance")
816 : return openmc::TriggerMetric::variance;
817 128 : else if (trigger == "std_dev")
818 : return openmc::TriggerMetric::standard_deviation;
819 128 : else if (trigger == "rel_err")
820 : return openmc::TriggerMetric::relative_error;
821 0 : else if (trigger == "none")
822 : return openmc::TriggerMetric::not_active;
823 : else
824 0 : mooseError("Unhandled TallyTriggerTypeEnum: ", trigger);
825 : }
826 :
827 : openmc::TriggerMetric
828 3004 : OpenMCProblemBase::triggerMetric(trigger::TallyTriggerTypeEnum trigger) const
829 : {
830 : switch (trigger)
831 : {
832 : case trigger::variance:
833 : return openmc::TriggerMetric::variance;
834 : case trigger::std_dev:
835 : return openmc::TriggerMetric::standard_deviation;
836 : case trigger::rel_err:
837 : return openmc::TriggerMetric::relative_error;
838 : case trigger::none:
839 : return openmc::TriggerMetric::not_active;
840 0 : default:
841 0 : mooseError("Unhandled TallyTriggerTypeEnum!");
842 : }
843 : }
844 :
845 : bool
846 0 : OpenMCProblemBase::cellIsVoid(const cellInfo & cell_info) const
847 : {
848 : // material_index will be unchanged if the cell is filled by a universe or lattice.
849 : // Otherwise, this will get set to the material index in the cell.
850 0 : int32_t material_index = 0;
851 0 : materialFill(cell_info, material_index);
852 0 : return material_index == MATERIAL_VOID;
853 : }
854 :
855 : void
856 1221 : OpenMCProblemBase::geometryType(bool & has_csg_universe, bool & has_dag_universe) const
857 : {
858 1221 : has_csg_universe = false;
859 1221 : has_dag_universe = false;
860 :
861 : // Loop over universes and check if type is DAGMC
862 6504 : for (const auto & universe : openmc::model::universes)
863 : {
864 5283 : if (universe->geom_type() == openmc::GeometryType::DAG)
865 68 : has_dag_universe = true;
866 5215 : else if (universe->geom_type() == openmc::GeometryType::CSG)
867 5215 : has_csg_universe = true;
868 : else
869 0 : mooseError("Unhandled GeometryType enum!");
870 : }
871 1221 : }
872 :
873 : long unsigned int
874 3040 : OpenMCProblemBase::numCells() const
875 : {
876 : long unsigned int n_openmc_cells = 0;
877 96039 : for (const auto & c : openmc::model::cells)
878 92999 : n_openmc_cells += c->n_instances();
879 :
880 3040 : return n_openmc_cells;
881 : }
882 :
883 : const openmc::Tally &
884 228 : OpenMCProblemBase::getCommonKineticsTally()
885 : {
886 228 : if (!_ifp_common_tally)
887 0 : mooseError("Internal error: kinetics parameters have not been enabled.");
888 :
889 228 : return *_ifp_common_tally;
890 : }
891 :
892 : const openmc::Tally &
893 198 : OpenMCProblemBase::getMGBetaTally()
894 : {
895 198 : return *_ifp_mg_beta_tally;
896 : }
897 :
898 : bool
899 192614 : OpenMCProblemBase::isReactionRateScore(const std::string & score) const
900 : {
901 : const std::set<std::string> viable_scores = {"H3-production",
902 : "total",
903 : "absorption",
904 : "scatter",
905 : "nu-scatter",
906 : "fission",
907 : "nu-fission",
908 : "prompt-nu-fission",
909 192614 : "delayed-nu-fission"};
910 192614 : return viable_scores.count(score);
911 : }
912 :
913 : bool
914 479906 : OpenMCProblemBase::isHeatingScore(const std::string & score) const
915 : {
916 : const std::set<std::string> viable_scores = {
917 479906 : "heating", "heating-local", "kappa-fission", "fission-q-prompt", "fission-q-recoverable"};
918 479906 : return viable_scores.count(score);
919 : }
920 :
921 : bool
922 202 : OpenMCProblemBase::validRandomRayScore(const std::string & score) const
923 : {
924 : const std::set<std::string> viable_scores = {
925 202 : "flux", "total", "fission", "nu-fission", "kappa-fission"};
926 202 : return viable_scores.count(score);
927 : }
928 :
929 : unsigned int
930 10290 : OpenMCProblemBase::addExternalVariable(const std::string & name,
931 : const std::string & system,
932 : const std::vector<SubdomainName> * block)
933 : {
934 10290 : auto var_params = _factory.getValidParams("MooseVariable");
935 20580 : var_params.set<MooseEnum>("family") = "MONOMIAL";
936 20580 : var_params.set<MooseEnum>("order") = "CONSTANT";
937 :
938 10290 : if (block)
939 13228 : var_params.set<std::vector<SubdomainName>>("block") = *block;
940 :
941 10290 : checkDuplicateVariableName(name, system);
942 20572 : addAuxVariable("MooseVariable", name, var_params);
943 20572 : return _aux->getFieldVariable<Real>(0, name).number();
944 10286 : }
945 :
946 : std::string
947 6210 : OpenMCProblemBase::subdomainName(const SubdomainID & id) const
948 : {
949 6210 : std::string name = _mesh.getSubdomainName(id);
950 6210 : if (name.empty())
951 11792 : name = std::to_string(id);
952 6210 : return name;
953 : }
954 :
955 : void
956 2128 : OpenMCProblemBase::getOpenMCUserObjects()
957 : {
958 2128 : _cell_transform_uos.clear();
959 :
960 2128 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("UserObject");
961 : std::vector<UserObject *> userobjs;
962 : uo_query.queryInto(userobjs);
963 :
964 11124 : for (const auto & u : userobjs)
965 : {
966 8996 : OpenMCNuclideDensities * c = dynamic_cast<OpenMCNuclideDensities *>(u);
967 8996 : if (c)
968 44 : _nuclide_densities_uos.push_back(c);
969 :
970 8996 : OpenMCTallyEditor * e = dynamic_cast<OpenMCTallyEditor *>(u);
971 8996 : if (e)
972 60 : _tally_editor_uos.push_back(e);
973 :
974 8996 : OpenMCDomainFilterEditor * f = dynamic_cast<OpenMCDomainFilterEditor *>(u);
975 8996 : if (f)
976 28 : _filter_editor_uos.push_back(f);
977 :
978 8996 : OpenMCCellTransform * t = dynamic_cast<OpenMCCellTransform *>(u);
979 8996 : if (t)
980 54 : _cell_transform_uos.push_back(t);
981 : }
982 :
983 2128 : checkOpenMCUserObjectIDs();
984 2124 : }
985 :
986 : bool
987 2175 : OpenMCProblemBase::hasCellTransform() const
988 : {
989 2175 : return !_cell_transform_uos.empty();
990 : }
991 :
992 : void
993 2128 : OpenMCProblemBase::checkOpenMCUserObjectIDs() const
994 : {
995 : std::set<int32_t> tally_ids;
996 2186 : for (const auto & te : _tally_editor_uos)
997 : {
998 60 : int32_t tally_id = te->tallyId();
999 : if (tally_ids.count(tally_id) != 0)
1000 2 : te->duplicateTallyError(tally_id);
1001 58 : tally_ids.insert(tally_id);
1002 : }
1003 :
1004 : std::set<int32_t> filter_ids;
1005 2152 : for (const auto & fe : _filter_editor_uos)
1006 : {
1007 28 : int32_t filter_id = fe->filterId();
1008 : if (filter_ids.count(filter_id) != 0)
1009 2 : fe->duplicateFilterError(filter_id);
1010 26 : filter_ids.insert(filter_id);
1011 : }
1012 2124 : }
1013 :
1014 : void
1015 2189 : OpenMCProblemBase::checkTallyEditorIDs() const
1016 : {
1017 2189 : std::vector<int32_t> mapped_tally_ids = getMappedTallyIDs();
1018 :
1019 2243 : for (const auto & te : _tally_editor_uos)
1020 : {
1021 56 : int32_t tally_id = te->tallyId();
1022 :
1023 : // ensure that the TallyEditor IDs don't apply to any mapped tally objects
1024 56 : if (std::find(mapped_tally_ids.begin(), mapped_tally_ids.end(), tally_id) !=
1025 : mapped_tally_ids.end())
1026 2 : te->mappedTallyError(tally_id);
1027 : }
1028 2187 : }
1029 :
1030 : void
1031 2729 : OpenMCProblemBase::executeFilterEditors()
1032 : {
1033 2729 : executeControls(EXEC_FILTER_EDITORS);
1034 :
1035 2729 : if (!_filter_editor_uos.size())
1036 : return;
1037 :
1038 24 : _console << "Executing filter editors..." << std::endl;
1039 48 : for (const auto & fe : _filter_editor_uos)
1040 24 : fe->execute();
1041 : }
1042 :
1043 : void
1044 2729 : OpenMCProblemBase::executeTallyEditors()
1045 : {
1046 2729 : executeControls(EXEC_TALLY_EDITORS);
1047 :
1048 2729 : if (!_tally_editor_uos.size())
1049 : return;
1050 :
1051 54 : _console << "Executing tally editors..." << std::endl;
1052 102 : for (const auto & te : _tally_editor_uos)
1053 54 : te->execute();
1054 : }
1055 :
1056 : void
1057 2729 : OpenMCProblemBase::executeEditors()
1058 : {
1059 2729 : executeFilterEditors();
1060 2729 : executeTallyEditors();
1061 2723 : }
1062 :
1063 : void
1064 3555 : OpenMCProblemBase::sendNuclideDensitiesToOpenMC()
1065 : {
1066 3555 : if (_nuclide_densities_uos.size() == 0)
1067 : return;
1068 :
1069 : // We could probably put this somewhere better, but it's good for now
1070 44 : executeControls(EXEC_SEND_OPENMC_DENSITIES);
1071 :
1072 44 : _console << "Sending nuclide compositions to OpenMC... ";
1073 84 : for (const auto & uo : _nuclide_densities_uos)
1074 44 : uo->setValue();
1075 : }
1076 :
1077 : Real
1078 8 : OpenMCProblemBase::tallyNormalizationValue() const
1079 : {
1080 8 : return _run_mode == openmc::RunMode::FIXED_SOURCE ? *_source_strength : *_power;
1081 : }
1082 :
1083 : #endif
|
