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 "OpenMCCellAverageProblem.h"
22 :
23 : #include "DelimitedFileReader.h"
24 : #include "DisplacedProblem.h"
25 : #include "TallyBase.h"
26 : #include "CellTally.h"
27 : #include "AddTallyAction.h"
28 : #include "SetupMGXSAction.h"
29 : #include "OpenMCVolumeCalculation.h"
30 : #include "CreateDisplacedProblemAction.h"
31 : #include "CriticalitySearchBase.h"
32 : #include "OpenMCCellMaterialFill.h"
33 :
34 : #include "openmc/constants.h"
35 : #include "openmc/cross_sections.h"
36 : #include "openmc/dagmc.h"
37 : #include "openmc/error.h"
38 : #include "openmc/lattice.h"
39 : #include "openmc/particle.h"
40 : #include "openmc/photon.h"
41 : #include "openmc/message_passing.h"
42 : #include "openmc/mgxs_interface.h"
43 : #include "openmc/nuclide.h"
44 : #include "openmc/random_lcg.h"
45 : #include "openmc/settings.h"
46 : #include "openmc/summary.h"
47 : #include "openmc/tallies/trigger.h"
48 : #include "openmc/volume_calc.h"
49 : #include "openmc/universe.h"
50 :
51 : registerMooseObject("CardinalApp", OpenMCCellAverageProblem);
52 :
53 : bool OpenMCCellAverageProblem::_first_transfer = true;
54 : bool OpenMCCellAverageProblem::_printed_initial = false;
55 : bool OpenMCCellAverageProblem::_printed_triso_warning = false;
56 :
57 : InputParameters
58 4803 : OpenMCCellAverageProblem::validParams()
59 : {
60 4803 : InputParameters params = OpenMCProblemBase::validParams();
61 9606 : params.addParam<bool>("output_cell_mapping",
62 9606 : true,
63 : "Whether to automatically output the mapping from OpenMC cells to the "
64 : "[Mesh], usually for diagnostic purposes");
65 :
66 9606 : params.addParam<MooseEnum>(
67 : "initial_properties",
68 9606 : getInitialPropertiesEnum(),
69 : "Where to read the temperature and density initial conditions for OpenMC");
70 :
71 9606 : params.addParam<bool>("export_properties",
72 9606 : false,
73 : "Whether to export OpenMC's temperature and density properties to an HDF5 "
74 : "file after updating them from MOOSE.");
75 9606 : params.addParam<bool>(
76 : "normalize_by_global_tally",
77 9606 : true,
78 : "Whether to normalize local tallies by a global tally (true) or else by the sum "
79 : "of the local tally (false)");
80 9606 : params.addParam<bool>("assume_separate_tallies",
81 9606 : false,
82 : "Whether to assume that all tallies added in the XML files or by Cardinal "
83 : "are spatially separate. This is a performance optimization");
84 :
85 : MooseEnum scores_heat(
86 9606 : "heating heating_local kappa_fission fission_q_prompt fission_q_recoverable");
87 9606 : params.addParam<MooseEnum>(
88 : "source_rate_normalization",
89 : scores_heat,
90 : "Score to use for computing the "
91 : "particle source rate (source/sec) for a certain tallies in "
92 : "eigenvalue mode. In other words, the "
93 : "source/sec is computed as (power divided by the global value of this tally)");
94 9606 : params.addParam<std::string>(
95 : "normalization_tally",
96 : "The name of a tally added in [Talliies] to be used when normalizing results in "
97 : "eigenvalue calculations. This tally object must contain the score specified in "
98 : "'source_rate_normalization'.");
99 :
100 9606 : params.addParam<MooseEnum>(
101 : "k_trigger",
102 9606 : getTallyTriggerEnum(),
103 : "Trigger criterion to determine when OpenMC simulation is complete based on k");
104 9606 : params.addRangeCheckedParam<Real>(
105 : "k_trigger_threshold", "k_trigger_threshold > 0", "Threshold for the k trigger");
106 9606 : params.addRangeCheckedParam<unsigned int>(
107 : "max_batches", "max_batches > 0", "Maximum number of batches, when using triggers");
108 14409 : params.addRangeCheckedParam<unsigned int>(
109 9606 : "batch_interval", 1, "batch_interval > 0", "Trigger batch interval");
110 :
111 9606 : params.addParam<std::vector<std::vector<std::string>>>(
112 : "temperature_variables",
113 : "Vector of variable names corresponding to the temperatures sent into OpenMC. Each entry "
114 : "maps to "
115 : "the corresponding entry in 'temperature_blocks.' If not specified, each entry defaults to "
116 : "'temp'");
117 9606 : params.addParam<std::vector<std::vector<SubdomainName>>>(
118 : "temperature_blocks",
119 : "Blocks corresponding to each of the 'temperature_variables'. If not specified, "
120 : "there will be no temperature feedback to OpenMC.");
121 :
122 9606 : params.addParam<std::vector<std::vector<std::string>>>(
123 : "density_variables",
124 : "Vector of variable names corresponding to the densities sent into OpenMC. Each entry maps "
125 : "to the corresponding entry in 'density_blocks.' If not specified, each entry defaults to "
126 : "'density'");
127 9606 : params.addParam<std::vector<std::vector<SubdomainName>>>(
128 : "density_blocks",
129 : "Blocks corresponding to each of the 'density_variables'. If not specified, "
130 : "there will be no density feedback to OpenMC.");
131 9606 : params.addRangeCheckedParam<std::vector<Real>>(
132 : "mgxs_reference_densities_by_block",
133 : "mgxs_reference_densities_by_block > 0.0",
134 : "Reference density values to use when applying density feedback (only used in multi-group "
135 : "mode). These densities represent the initial densities used when generated the multigroup "
136 : "library. Each entry maps to the corresponding row in 'density_blocks.' Units are "
137 : "expected to be kg/m3.");
138 :
139 9606 : params.addParam<unsigned int>("cell_level",
140 : "Coordinate level in OpenMC (across the entire geometry) to use "
141 : "for identifying cells");
142 9606 : params.addParam<unsigned int>(
143 : "lowest_cell_level",
144 : "Lowest coordinate level in OpenMC to use for identifying cells. The cell level for coupling "
145 : "will use the value set with this parameter unless the geometry does not have that many "
146 : "layers of geometry nesting, in which case the locally lowest depth is used");
147 :
148 9606 : params.addParam<std::vector<SubdomainName>>(
149 : "identical_cell_fills",
150 : "Blocks on which the OpenMC cells have identical fill universes; this is an optimization to "
151 : "speed up initialization for TRISO problems while also reducing memory usage. It is assumed "
152 : "that any cell which maps to one of these subdomains has exactly the same universe filling "
153 : "it as all other cells which map to these subdomains. We HIGHLY recommend that the first "
154 : "time you try using this, that you also set 'check_identical_cell_fills = true' to catch "
155 : "any possible user errors which would exclude you from using this option safely.");
156 9606 : params.addParam<bool>(
157 : "check_identical_cell_fills",
158 9606 : false,
159 : "Whether to check that your model does indeed have identical cell fills, allowing "
160 : "you to set 'identical_cell_fills' to speed up initialization");
161 :
162 9606 : params.addParam<MooseEnum>(
163 9606 : "relaxation", getRelaxationEnum(), "Type of relaxation to apply to the OpenMC solution");
164 14409 : params.addRangeCheckedParam<Real>("relaxation_factor",
165 9606 : 0.5,
166 : "relaxation_factor > 0.0 & relaxation_factor < 2.0",
167 : "Relaxation factor for use with constant relaxation");
168 9606 : params.addParam<int>("first_iteration_particles",
169 : "Number of particles to use for first iteration "
170 : "when using Dufek-Gudowski relaxation");
171 :
172 9606 : params.addParam<UserObjectName>(
173 : "symmetry_mapper",
174 : "User object (of type SymmetryPointGenerator) "
175 : "to map from a symmetric OpenMC model to a full-domain [Mesh]. For example, you can use this "
176 : "to map from a quarter-symmetric OpenMC model to a whole-domain [Mesh].");
177 :
178 9606 : params.addParam<UserObjectName>(
179 : "volume_calculation",
180 : "User object that will perform a stochastic volume calculation to get the OpenMC "
181 : "cell volumes. This can be used to check that the MOOSE regions to which the cells map are "
182 : "of approximately the same volume as the true cells.");
183 9606 : params.addParam<UserObjectName>("skinner",
184 : "When using DAGMC geometries, an optional skinner that will "
185 : "regenerate the OpenMC geometry on-the-fly according to "
186 : "iso-contours of temperature and density");
187 4803 : params.addClassDescription(
188 : "Couple OpenMC to MOOSE through cell-averaged temperature, density, and tallies.");
189 :
190 4803 : return params;
191 4803 : }
192 :
193 2410 : OpenMCCellAverageProblem::OpenMCCellAverageProblem(const InputParameters & params)
194 : : OpenMCProblemBase(params),
195 2396 : _serialized_solution(_aux->serializedSolution()),
196 4792 : _output_cell_mapping(getParam<bool>("output_cell_mapping")),
197 2396 : _initial_condition(
198 2396 : getParam<MooseEnum>("initial_properties").getEnum<coupling::OpenMCInitialCondition>()),
199 4792 : _relaxation(getParam<MooseEnum>("relaxation").getEnum<relaxation::RelaxationEnum>()),
200 4792 : _k_trigger(getParam<MooseEnum>("k_trigger").getEnum<trigger::TallyTriggerTypeEnum>()),
201 4792 : _export_properties(getParam<bool>("export_properties")),
202 4792 : _using_skinner(isParamValid("skinner")),
203 : // 'used_displaced' is added to '_need_to_reinit_coupling' later in the ctor.
204 2396 : _need_to_reinit_coupling(_has_adaptivity || _using_skinner),
205 2396 : _has_identical_cell_fills(params.isParamSetByUser("identical_cell_fills")),
206 4792 : _check_identical_cell_fills(getParam<bool>("check_identical_cell_fills")),
207 4792 : _assume_separate_tallies(getParam<bool>("assume_separate_tallies")),
208 2396 : _specified_density_feedback(params.isParamSetByUser("density_blocks")),
209 2396 : _specified_temperature_feedback(params.isParamSetByUser("temperature_blocks")),
210 2396 : _needs_to_map_cells(_specified_density_feedback || _specified_temperature_feedback),
211 : _volume_calc(nullptr),
212 2396 : _symmetry(nullptr),
213 6623 : _initial_num_openmc_surfaces(openmc::model::surfaces.size())
214 : {
215 2396 : const auto & subdomains = mesh().meshSubdomains();
216 6183 : for (const auto & s : subdomains)
217 3789 : if (mesh().getCoordSystem(s) == Moose::COORD_RZ)
218 2 : mooseError(
219 : "OpenMC coupling to axisymmetric meshes is not yet supported! Please convert your mesh "
220 : "block to a 3-D mesh (you may still use axisymmetric meshes for your other physics "
221 : "coupled to OpenMC and transfer data between those apps and a 3-D OpenMC model. You just "
222 : "cannot use an axisymmetric mesh from which OpenMC reads/writes data).");
223 :
224 2394 : if (_specified_temperature_feedback && openmc::settings::temperature_range[1] == 0.0)
225 10 : mooseWarning("For multiphysics simulations, we recommend setting the 'temperature_range' in "
226 : "OpenMC's settings.xml file. This will pre-load nuclear data over a range of "
227 : "temperatures, instead of only the temperatures defined in the XML file.\n\nFor "
228 : "efficiency purposes, OpenMC only checks that cell temperatures are within the "
229 : "global min/max of loaded data, which can be different from data loaded for each "
230 : "nuclide. Run may abort suddenly if requested nuclear data is not available.");
231 :
232 : // Check to see if a displaced problem is being initialized
233 : const auto & dis_actions =
234 2392 : getMooseApp().actionWarehouse().getActions<CreateDisplacedProblemAction>();
235 4784 : for (const auto & act : dis_actions)
236 : {
237 2392 : auto displacements = act->isParamValid("displacements");
238 4784 : auto use = act->getParam<bool>("use_displaced_mesh");
239 2392 : _use_displaced = displacements && use;
240 :
241 : // print a warning if the user added displacements, but are not using them
242 2392 : if (!use && displacements)
243 0 : mooseWarning("When 'use_displaced_mesh' is false, the 'displacements' are unused!");
244 :
245 7176 : if (act->isParamSetByUser("use_displaced_mesh") && use && !displacements)
246 0 : mooseWarning("When 'use_displaced_mesh' is true, but no 'displacements' are provided, then "
247 : "the displaced mesh will not be used.");
248 :
249 2392 : _need_to_reinit_coupling |= _use_displaced;
250 : }
251 :
252 : // Look through the list of AddTallyActions to see if we have a CellTally. If so, we need to map
253 : // cells.
254 2392 : const auto & tally_actions = getMooseApp().actionWarehouse().getActions<AddTallyAction>();
255 4768 : for (const auto & act : tally_actions)
256 2376 : _has_cell_tallies |= act->getMooseObjectType() == "CellTally";
257 :
258 : // Repeat the same check for SetUpMGXSActions.
259 2392 : const auto & mgxs_actions = getMooseApp().actionWarehouse().getActions<SetupMGXSAction>();
260 2450 : for (const auto & act : mgxs_actions)
261 58 : _has_cell_tallies |= act->addingCellTallies();
262 2392 : _needs_to_map_cells |= _has_cell_tallies;
263 :
264 2392 : if (!_needs_to_map_cells)
265 552 : checkUnusedParam(params,
266 : "output_cell_mapping",
267 : "'temperature_blocks', 'density_blocks', and 'tally_blocks' are empty");
268 :
269 2392 : if (!_specified_temperature_feedback && !_specified_density_feedback)
270 1470 : checkUnusedParam(
271 : params, "initial_properties", "'temperature_blocks' and 'density_blocks' are unused");
272 :
273 : // We need to clear and re-initialize OpenMC problem in the cases of:
274 : // - the [Mesh] is being adaptively refined
275 : // - the [Mesh] is deforming in space
276 : //
277 : // If the [Mesh] is changing, then we certainly know that the mesh tallies
278 : // need to be re-initialized because (a) for file-based mesh tallies, we need
279 : // to enforce that the mesh is identical to the [Mesh] and (b) for directly
280 : // tallying on the [Mesh], we need to pass that mesh info into OpenMC. For good
281 : // measure, we also need to re-initialize cell tallies because it's possible
282 : // that as the [Mesh] changes, the mapping from OpenMC cells to the [Mesh]
283 : // also changes, which could open the door to new cell IDs/instances being added
284 : // to the cell instance filter. If we need to re-init tallies, then we can't
285 : // guarantee that the tallies from iteration to iteration correspond to exactly
286 : // the same number of bins or to exactly the same regions of space, so we must
287 : // disable relaxation.
288 2392 : if ((_use_displaced || _has_adaptivity) && _relaxation != relaxation::none)
289 4 : paramError(
290 : "relaxation",
291 : "When adaptivity is requested or a displaced problem is used, the mapping from the "
292 : "OpenMC model to the [Mesh] may vary in time. This means that we have no guarantee that "
293 : "the "
294 : "number of tally bins (or even the regions of space corresponding to each bin) are fixed. "
295 : "Therefore, it is not possible to apply relaxation to the OpenMC tallies because you might "
296 : "end up trying to add vectors of different length (and possibly spatial mapping).");
297 :
298 2388 : if (_run_mode == openmc::RunMode::FIXED_SOURCE)
299 290 : checkUnusedParam(params, "normalize_by_global_tally", "running OpenMC in fixed source mode");
300 :
301 2386 : if (_run_mode != openmc::RunMode::EIGENVALUE && _k_trigger != trigger::none)
302 2 : paramError("k_trigger",
303 : "Cannot specify a 'k_trigger' for OpenMC runs that are not eigenvalue mode!");
304 :
305 : // determine the number of particles set either through XML or the wrapping
306 2384 : if (_relaxation == relaxation::dufek_gudowski)
307 : {
308 32 : checkUnusedParam(params, "particles", "using Dufek-Gudowski relaxation");
309 32 : checkRequiredParam(params, "first_iteration_particles", "using Dufek-Gudowski relaxation");
310 32 : openmc::settings::n_particles = getParam<int>("first_iteration_particles");
311 32 : _n_particles_1 = getParam<int>("first_iteration_particles");
312 : }
313 : else
314 4736 : checkUnusedParam(params, "first_iteration_particles", "not using Dufek-Gudowski relaxation");
315 :
316 : // OpenMC will throw an error if the geometry contains DAG universes but OpenMC wasn't compiled
317 : // with DAGMC. So we can assume that if we have a DAGMC geometry, that we will also by this
318 : // point have DAGMC enabled.
319 : #ifdef ENABLE_DAGMC
320 : bool has_csg;
321 : bool has_dag;
322 1227 : geometryType(has_csg, has_dag);
323 :
324 1227 : if (!has_dag)
325 2319 : checkUnusedParam(
326 : params, "skinner", "the OpenMC model does not contain any DagMC universes", true);
327 67 : else if (_using_skinner)
328 : {
329 : // Loop over all universes to find the DAGMC universe and to check and make sure we only have
330 : // the one.
331 : unsigned int num_dag_universes = 0;
332 116 : for (const auto & universe : openmc::model::universes)
333 : {
334 64 : if (universe->geom_type() == openmc::GeometryType::DAG)
335 : {
336 53 : _dagmc_universe_id = universe->id_;
337 53 : num_dag_universes++;
338 : }
339 : }
340 :
341 52 : if (num_dag_universes != 1)
342 2 : mooseError("The 'skinner' can only be used when the OpenMC geometry contains a single DAGMC "
343 : "universe.\n"
344 1 : "Your geometry contains " +
345 0 : Moose::stringify(num_dag_universes) + " DAGMC universes.");
346 :
347 : // Loop over each element of each lattice to make sure that it doesn't contain the DAGMC
348 : // universe.
349 51 : for (const auto & lattice : openmc::model::lattices)
350 : {
351 3 : for (openmc::LatticeIter it = lattice->begin(); it != lattice->end(); ++it)
352 2 : if (openmc::model::universes[*it]->id_ == _dagmc_universe_id)
353 1 : mooseError("The 'skinner' cannot be used when the DAGMC universe is contained in lattice "
354 : "geometry.");
355 :
356 1 : if (lattice->outer_ != openmc::NO_OUTER_UNIVERSE &&
357 1 : openmc::model::universes[lattice->outer_]->id_ == _dagmc_universe_id)
358 1 : mooseError("The 'skinner' cannot be used when the DAGMC universe is used as the outer "
359 : "universe of a lattice.");
360 : }
361 :
362 : // Need to make sure that there is only a single cell which uses the DAGMC universe as it's
363 : // fill. The root universe must contain that cell, otherwise the DAGMC universe may be
364 : // replicated across the problem.
365 : unsigned int num_dag_instances = 0;
366 223 : for (const auto & cell : openmc::model::cells)
367 : {
368 174 : if (cell->type_ == openmc::Fill::UNIVERSE &&
369 8 : cell->fill_ == openmc::model::universe_map.at(_dagmc_universe_id))
370 : {
371 7 : _dagmc_root_universe = false;
372 7 : num_dag_instances++;
373 7 : _cell_using_dagmc_universe_id = cell->id_;
374 : }
375 : }
376 :
377 49 : if (num_dag_instances > 1)
378 2 : mooseError("The 'skinner' can only be used when the DAGMC universe in the OpenMC geometry is "
379 : "used as a cell "
380 1 : "fill at most once.\n Your geometry contains " +
381 0 : Moose::stringify(num_dag_instances) +
382 : " cells which "
383 : "use the DAGMC universe as their fill.");
384 :
385 48 : if (!_dagmc_root_universe &&
386 5 : openmc::model::cells[openmc::model::cell_map.at(_cell_using_dagmc_universe_id)]
387 5 : ->universe_ != openmc::model::root_universe)
388 1 : mooseError("The 'skinner' can only be used when the cell using the DAGMC universe as a fill "
389 : "is contained in the "
390 : "root universe.");
391 : }
392 : #else
393 2314 : checkUnusedParam(
394 : params, "skinner", "DAGMC geometries in OpenMC are not enabled in this build of Cardinal");
395 : #endif
396 :
397 2378 : if (_relaxation != relaxation::constant)
398 4408 : checkUnusedParam(params, "relaxation_factor", "not using constant relaxation");
399 :
400 2378 : readBlockParameters("identical_cell_fills", _identical_cell_fill_blocks);
401 :
402 2378 : if (!_has_identical_cell_fills)
403 4704 : checkUnusedParam(
404 : params, "check_identical_cell_fills", "'identical_cell_fills' is not specified");
405 :
406 4746 : readBlockVariables("temperature", "temp", _temp_vars_to_blocks, _temp_blocks);
407 4728 : readBlockVariables("density", "density", _density_vars_to_blocks, _density_blocks);
408 :
409 : // When running in multi-group mode, the user needs to provide a reference density if density
410 : // feedback is specified (to convert to the dimensionless MGXS density). In the future, it would
411 : // be nice if OpenMC materials could store their own reference densities (in multi-group mode)
412 : // as this is rather error prone.
413 2360 : if (!openmc::settings::run_CE && _specified_density_feedback)
414 : {
415 56 : checkRequiredParam(params,
416 : "mgxs_reference_densities_by_block",
417 : "running in multi-group mode and using density feedback");
418 28 : const auto & density_scales = getParam<std::vector<Real>>("mgxs_reference_densities_by_block");
419 :
420 : const auto & density_blocks =
421 56 : getParam<std::vector<std::vector<SubdomainName>>>("density_blocks");
422 :
423 28 : if (density_scales.size() != density_blocks.size())
424 2 : paramError(
425 : "mgxs_reference_densities_by_block",
426 : "'mgxs_reference_densities_by_block' must have the same number of entries as rows in "
427 : "'density_blocks'!");
428 :
429 60 : for (unsigned int i = 0; i < density_blocks.size(); ++i)
430 76 : for (const auto & subdomain_name : density_blocks[i])
431 42 : _subdomain_to_ref_density[mesh().getSubdomainID(subdomain_name)] = density_scales[i];
432 : }
433 : else
434 4664 : checkUnusedParam(params,
435 : "mgxs_reference_densities_by_block",
436 : "not running in multi-group mode and using density feedback");
437 :
438 2392 : for (const auto & i : _identical_cell_fill_blocks)
439 36 : if (std::find(_density_blocks.begin(), _density_blocks.end(), i) != _density_blocks.end())
440 2 : paramError(
441 : "identical_cell_fills",
442 : "Entries in 'identical_cell_fills' cannot be contained in 'density_blocks'; the\n"
443 : "identical fill universe optimization is not yet implemented for density feedback.");
444 :
445 2356 : if (_needs_to_map_cells)
446 : {
447 6246 : if (isParamValid("cell_level") == isParamValid("lowest_cell_level"))
448 1 : mooseError("Either 'cell_level' or 'lowest_cell_level' must be specified. You have given "
449 : "either both or none.");
450 :
451 : std::string selected_param;
452 4162 : if (isParamValid("cell_level"))
453 : {
454 4074 : _cell_level = getParam<unsigned int>("cell_level");
455 : selected_param = "cell_level";
456 :
457 2037 : if (_cell_level >= openmc::model::n_coord_levels)
458 4 : paramError(selected_param,
459 : "Coordinate level for finding cells cannot be greater than total number "
460 2 : "of coordinate levels: " +
461 0 : Moose::stringify(openmc::model::n_coord_levels) + "!");
462 : }
463 : else
464 : {
465 88 : _cell_level = getParam<unsigned int>("lowest_cell_level");
466 : selected_param = "lowest_cell_level";
467 : }
468 : }
469 : else
470 : {
471 548 : checkUnusedParam(params,
472 : "cell_level",
473 : "'temperature_blocks', 'density_blocks', and 'tally_blocks' are empty");
474 548 : checkUnusedParam(params,
475 : "lowest_cell_level",
476 : "'temperature_blocks', 'density_blocks', and 'tally_blocks' are empty");
477 : }
478 2353 : }
479 :
480 : const MooseMesh &
481 957745 : OpenMCCellAverageProblem::getMooseMesh() const
482 : {
483 957745 : return mesh(_use_displaced);
484 : }
485 :
486 : MooseMesh &
487 30294033 : OpenMCCellAverageProblem::getMooseMesh()
488 : {
489 : // TODO: this could go into MOOSE framework directly
490 30294033 : if (_use_displaced && !_displaced_problem)
491 0 : mooseWarning("Displaced mesh was requested but the displaced problem does not exist. "
492 : "Regular mesh will be returned");
493 :
494 30294033 : MooseMesh & m = ((_use_displaced && _displaced_problem) ? _displaced_problem->mesh() : mesh());
495 30294033 : return m;
496 : }
497 :
498 : void
499 4746 : OpenMCCellAverageProblem::readBlockVariables(
500 : const std::string & param,
501 : const std::string & default_name,
502 : std::map<std::string, std::vector<SubdomainName>> & vars_to_specified_blocks,
503 : std::vector<SubdomainID> & specified_blocks)
504 : {
505 4746 : std::string b = param + "_blocks";
506 4746 : std::string v = param + "_variables";
507 :
508 4746 : if (!isParamValid(b))
509 : {
510 7818 : checkUnusedParam(parameters(), v, "not setting '" + b + "'");
511 : return;
512 : }
513 :
514 : std::vector<std::vector<SubdomainName>> blocks;
515 4270 : read2DBlockParameters(b, blocks, specified_blocks);
516 :
517 : // now, get the names of those variables
518 : std::vector<std::vector<std::string>> vars;
519 2130 : if (isParamValid(v))
520 : {
521 74 : vars = getParam<std::vector<std::vector<std::string>>>(v);
522 :
523 222 : checkEmptyVector(vars, "'" + v + "");
524 252 : for (const auto & t : vars)
525 534 : checkEmptyVector(t, "Entries in '" + v + "'");
526 :
527 74 : if (vars.size() != blocks.size())
528 24 : mooseError("'" + v + "' and '" + b + "' must be the same length!\n'" + v + "' is of length " +
529 12 : std::to_string(vars.size()) + " and '" + b + "' is of length " +
530 4 : std::to_string(blocks.size()));
531 :
532 : // TODO: for now, we restrict each set of blocks to map to a single variable
533 232 : for (std::size_t i = 0; i < vars.size(); ++i)
534 166 : if (vars[i].size() > 1)
535 12 : mooseError("Each entry in '" + v + "' must be of length 1. Entry " + std::to_string(i) +
536 8 : " is of length " + std::to_string(vars[i].size()));
537 : }
538 : else
539 : {
540 : // set a reasonable default, if not specified
541 2056 : vars.resize(blocks.size(), std::vector<std::string>(1));
542 4112 : for (std::size_t i = 0; i < blocks.size(); ++i)
543 : vars[i][0] = default_name;
544 : }
545 :
546 4340 : for (std::size_t i = 0; i < vars.size(); ++i)
547 5392 : for (std::size_t j = 0; j < blocks[i].size(); ++j)
548 3174 : vars_to_specified_blocks[vars[i][0]].push_back(blocks[i][j]);
549 2122 : }
550 :
551 : void
552 2139 : OpenMCCellAverageProblem::initialSetup()
553 : {
554 2139 : OpenMCProblemBase::initialSetup();
555 :
556 : // Find ModelModifier objects and store them in _cell_material_modifiers map
557 2139 : TheWarehouse::Query mm_query = theWarehouse().query().condition<AttribSystem>("ModelModifiers");
558 : std::vector<ModelModifiersBase *> mm_objs;
559 : mm_query.queryInto(mm_objs);
560 :
561 : // loop through all queried ModelModifiers and add any OpenMCCellMaterialFill
562 : // to the _cell_material_modifiers map
563 2171 : for (const auto & m : mm_objs)
564 : {
565 32 : auto * modifier = dynamic_cast<OpenMCCellMaterialFill *>(m);
566 32 : if (modifier)
567 32 : _cell_material_modifiers[modifier->getCellIndex()] = modifier;
568 : }
569 :
570 2139 : getOpenMCUserObjects();
571 :
572 2135 : if (_use_displaced && !_using_skinner && !hasCellTransform())
573 14 : mooseWarning("Your problem has a moving mesh, but you have not provided a 'skinner' or an "
574 : "OpenMCCellTransform user object (both of which move the OpenMC geometry). The "
575 : "[Mesh] will move, but the underlying OpenMC geometry will remain unchanged. "
576 : "Unexpected behavior may occur.");
577 :
578 : // Coupling re-initialization should be triggered if we have cell transforms which can happen
579 : // even if the mesh isn't moving or adaptive.
580 2134 : _need_to_reinit_coupling |= hasCellTransform();
581 : // The criticality search may modify the geometry.
582 2134 : if (_criticality_search)
583 76 : _need_to_reinit_coupling |= _criticality_search->changingGeometry();
584 :
585 2134 : if (!_needs_to_map_cells)
586 458 : checkUnusedParam(parameters(),
587 : "volume_calculation",
588 : "'temperature_blocks', 'density_blocks', and 'tally_blocks' are empty");
589 3810 : else if (isParamValid("volume_calculation"))
590 : {
591 104 : const auto & name = getParam<UserObjectName>("volume_calculation");
592 104 : auto * base = &getUserObject<UserObject>(name);
593 :
594 104 : _volume_calc = dynamic_cast<OpenMCVolumeCalculation *>(base);
595 :
596 104 : if (!_volume_calc)
597 0 : paramError("volume_calculation",
598 : "The 'volume_calculation' user object must be of type "
599 : "OpenMCVolumeCalculation!");
600 : }
601 :
602 4268 : if (isParamValid("symmetry_mapper"))
603 : {
604 43 : const auto & name = getParam<UserObjectName>("symmetry_mapper");
605 43 : auto base = &getUserObject<UserObject>(name);
606 :
607 43 : _symmetry = dynamic_cast<SymmetryPointGenerator *>(base);
608 :
609 43 : if (!_symmetry)
610 2 : paramError("symmetry_mapper",
611 : "The 'symmetry_mapper' user object has to be of type SymmetryPointGenerator!");
612 : }
613 :
614 : // Get triggers.
615 2132 : getTallyTriggerParameters(_pars);
616 :
617 2130 : setupProblem();
618 :
619 : #ifdef ENABLE_DAGMC
620 1071 : if (_using_skinner)
621 : {
622 46 : std::set<SubdomainID> t(_temp_blocks.begin(), _temp_blocks.end());
623 46 : std::set<SubdomainID> d(_density_blocks.begin(), _density_blocks.end());
624 :
625 46 : if (t != getMooseMesh().meshSubdomains())
626 0 : paramError("temperature_blocks",
627 : "The 'skinner' requires temperature feedback to be applied over the entire mesh. "
628 : "Please update `temperature_blocks` to include all blocks.");
629 :
630 46 : if (d != getMooseMesh().meshSubdomains() && _specified_density_feedback)
631 0 : paramError("density_blocks",
632 : "The 'skinner' requires density feedback to be applied over the entire mesh. "
633 : "Please update `density_blocks` to include all blocks.");
634 :
635 46 : if (t != d && _specified_density_feedback)
636 0 : mooseError("The 'skinner' will apply skinning over the entire domain, and requires that the "
637 : "entire problem uses identical settings for feedback. Please update "
638 : "'temperature_blocks' and 'density_blocks' to include all blocks.");
639 :
640 46 : if (_symmetry)
641 1 : mooseError("Cannot combine the 'skinner' with 'symmetry_mapper'!\n\nWhen using a skinner, "
642 : "the [Mesh] must exactly match the underlying OpenMC model, so there is\n"
643 : "no need to transform spatial coordinates to map between OpenMC and the [Mesh].");
644 :
645 : // Rudimentary error checking to make sure all non-void DAGMC cells are mapped. This helps catch
646 : // errors where the skinned MOOSE mesh deletes DAGMC geometry. Also error if the user is
647 : // attempting to use a skinner when mapping both CSG cells and DAGMC geometry to the MOOSE mesh.
648 : // The skinner is currently not set up to ignore elements that map to cells and will generate
649 : // DAGMC geometry that overlaps with pre-existing CSG cells.
650 : // TODO: This would be nice to fix, but would require a rework of the skinner.
651 : std::set<int32_t> mapped_dag_cells;
652 199 : for (const auto & c : openmc::model::cells)
653 : {
654 832 : for (const auto & [c_info, elem] : _cell_to_elem)
655 : {
656 678 : if (c->geom_type() == openmc::GeometryType::DAG &&
657 671 : c_info.first == openmc::model::cell_map.at(c->id_))
658 102 : mapped_dag_cells.insert(c->id_);
659 576 : else if (c->geom_type() == openmc::GeometryType::CSG &&
660 7 : c_info.first == openmc::model::cell_map.at(c->id_))
661 1 : mooseError("At present, the 'skinner' can only be used when the only OpenMC geometry "
662 : "which maps to the MOOSE mesh is DAGMC geometry. Your geometry contains CSG "
663 : "cells which map to the MOOSE mesh.");
664 : }
665 : }
666 :
667 : unsigned int num_unmapped = 0;
668 : unsigned int num_dag_cells = 0;
669 198 : for (const auto & c : openmc::model::cells)
670 : {
671 : auto no_void =
672 154 : std::find(c->material_.begin(), c->material_.end(), MATERIAL_VOID) == c->material_.end();
673 154 : if (mapped_dag_cells.count(c->id_) == 0 && c->geom_type() == openmc::GeometryType::DAG &&
674 : no_void)
675 1 : num_unmapped++;
676 154 : if (c->geom_type() == openmc::GeometryType::DAG)
677 151 : num_dag_cells++;
678 : }
679 :
680 44 : if (num_unmapped > 0)
681 2 : mooseWarning("Your DAGMC geometry contains unmapped cells! The skinner assumes that "
682 : "the DAG geometry used in the OpenMC model maps one to one to the mesh "
683 : "mirror; if that is not the case the skinner may delete some parts of "
684 1 : "your OpenMC model when the underlying geometry is regenerated. You have " +
685 1 : Moose::stringify(num_unmapped) + " unmapped DAGMC cells out of " +
686 0 : Moose::stringify(num_dag_cells) + " DAGMC cells.");
687 :
688 43 : const auto & name = getParam<UserObjectName>("skinner");
689 43 : auto base = &getUserObject<UserObject>(name);
690 :
691 43 : _skinner = dynamic_cast<MoabSkinner *>(base);
692 :
693 43 : if (!_skinner)
694 1 : paramError("skinner", "The 'skinner' user object must be of type MoabSkinner!");
695 :
696 42 : if (_skinner->hasDensitySkinning() != _specified_density_feedback)
697 1 : mooseError(
698 : "Detected inconsistent settings for density skinning and 'density_blocks'. If applying "
699 : "density feedback with 'density_blocks', then you must apply density skinning in the '",
700 : name,
701 : "' user object (and vice versa)");
702 :
703 41 : if (_initial_condition == coupling::hdf5)
704 1 : paramError("initial_properties",
705 : "Cannot load initial temperature and density properties from "
706 : "HDF5 files because there is no guarantee that the geometry (which is adaptively "
707 : "changing) matches "
708 : "that used to write the HDF5 file.");
709 :
710 : // If the DAGMC universe is the root universe the geometry contains no CSG cells. We need
711 : // to force the skinner to add a graveyard as the problem will contain no boundary contitions
712 : // after skinning is performed. If there are CSG cells in the geometry, this is not the case
713 : // as the DAGMC universe is embedded in a cell (which applies boundary conditions).
714 40 : if (_dagmc_root_universe)
715 37 : _skinner->setGraveyard(true);
716 :
717 39 : _skinner->setScaling(_scaling);
718 39 : _skinner->setVerbosity(_verbose);
719 39 : _skinner->makeDependentOnExternalAction();
720 39 : _skinner->setUseDisplacedMesh(_use_displaced);
721 :
722 : // the skinner expects that there is one OpenMC material per subdomain (otherwise this
723 : // indicates that our [Mesh] doesn't match the .h5m model, because DAGMC itself imposes
724 : // the one-material-per-cell case. In the future, if we generate DAGMC models directly
725 : // from the [Mesh] (bypassing the .h5m), we would not need this error check.
726 39 : _skinner->setMaterialNames(getMaterialInEachSubdomain());
727 38 : _skinner->initialize();
728 : }
729 : #endif
730 2070 : }
731 :
732 : std::vector<std::string>
733 39 : OpenMCCellAverageProblem::getMaterialInEachSubdomain() const
734 : {
735 : std::vector<std::string> mats;
736 133 : for (const auto & s : _subdomain_to_material)
737 : {
738 95 : if (s.second.size() > 1)
739 : {
740 1 : std::stringstream msg;
741 : msg << "The 'skinner' expects to find one OpenMC material mapped to each [Mesh] subdomain, "
742 : "but "
743 3 : << Moose::stringify(s.second.size()) << " materials\nmapped to subdomain " << s.first
744 : << ". This indicates your [Mesh] is not "
745 1 : << "consistent with the .h5m model.\n\nThe materials which mapped to subdomain "
746 1 : << s.first << " are:\n";
747 :
748 3 : for (const auto & m : s.second)
749 4 : msg << "\n" << materialName(m);
750 :
751 1 : mooseError(msg.str());
752 0 : }
753 :
754 188 : mats.push_back(materialName(*(s.second.begin())));
755 : }
756 :
757 38 : return mats;
758 0 : }
759 :
760 : void
761 3051 : OpenMCCellAverageProblem::setupProblem()
762 : {
763 : // establish the local -> global element mapping for convenience
764 3051 : _local_to_global_elem.clear();
765 4051261 : for (unsigned int e = 0; e < getMooseMesh().nElem(); ++e)
766 : {
767 4048210 : const auto * elem = getMooseMesh().queryElemPtr(e);
768 4048210 : if (!isLocalElem(elem) || !elem->active())
769 1939676 : continue;
770 :
771 2108534 : _local_to_global_elem.push_back(e);
772 : }
773 :
774 3051 : _n_openmc_cells = numCells();
775 :
776 3051 : initializeElementToCellMapping();
777 :
778 : // we do this last so that we can at least hit any other errors first before
779 : // spending time on the costly filled cell caching
780 3025 : cacheContainedCells();
781 :
782 : // save the number of contained cells for printing in every transfer if verbose
783 : _cell_to_n_contained.clear();
784 17524 : for (const auto & c : _cell_to_elem)
785 14507 : _cell_to_n_contained[c.first] = numContainedMaterialCells(c.first);
786 :
787 : // the _subdomain_to_material member variable is only used if printing out verbose
788 : // information or if the skinner is used (TODO: though, that statement about the skinner
789 : // will probably go away once we update the skinner to more modern approach of density
790 : // multipliers rather than unique materials per bin). For large problems, determining
791 : // how the subdomains map to OpenMC materials can be costly, so we should make this optional
792 3017 : if (_verbose || _using_skinner)
793 1809 : subdomainsToMaterials();
794 :
795 3015 : initializeTallies();
796 3003 : }
797 :
798 : void
799 2132 : OpenMCCellAverageProblem::getTallyTriggerParameters(const InputParameters & parameters)
800 : {
801 : // parameters needed for k triggers
802 : bool has_tally_trigger = false;
803 2132 : if (_k_trigger != trigger::none)
804 : {
805 68 : checkRequiredParam(parameters, "k_trigger_threshold", "using a k trigger");
806 68 : openmc::settings::keff_trigger.threshold = getParam<Real>("k_trigger_threshold");
807 : has_tally_trigger = true;
808 : }
809 : else
810 4196 : checkUnusedParam(parameters, "k_trigger_threshold", "not using a k trigger");
811 :
812 : // Check to see if any of the local tallies have triggers.
813 4808 : for (const auto & local_tally : _local_tallies)
814 2676 : has_tally_trigger = has_tally_trigger || local_tally->hasTrigger();
815 :
816 2132 : if (has_tally_trigger) // at least one trigger
817 : {
818 106 : openmc::settings::trigger_on = true;
819 212 : checkRequiredParam(parameters, "max_batches", "using triggers");
820 :
821 106 : if (_skip_statepoint)
822 0 : checkUnusedParam(parameters, "skip_statepoint", "using a trigger");
823 :
824 212 : int err = openmc_set_n_batches(getParam<unsigned int>("max_batches"),
825 : true /* set the max batches */,
826 106 : true /* add the last batch for statepoint writing */);
827 106 : catchOpenMCError(err, "set the maximum number of batches");
828 :
829 208 : openmc::settings::trigger_batch_interval = getParam<unsigned int>("batch_interval");
830 : }
831 : else
832 : {
833 4052 : checkUnusedParam(parameters, "max_batches", "not using triggers");
834 4052 : checkUnusedParam(parameters, "batch_interval", "not using triggers");
835 :
836 2026 : if (_skip_statepoint)
837 : openmc::settings::statepoint_batch.clear();
838 : }
839 2130 : }
840 :
841 : const TallyBase *
842 84 : OpenMCCellAverageProblem::getTally(const std::string & name)
843 : {
844 192 : for (const auto & t : _local_tallies)
845 190 : if (t->name() == name)
846 : return t.get();
847 : return nullptr;
848 : }
849 :
850 : std::vector<const MooseVariableFE<Real> *>
851 70 : OpenMCCellAverageProblem::getTallyScoreVariables(const std::string & score,
852 : const std::string & tally_name,
853 : THREAD_ID tid,
854 : const std::string & output,
855 : bool skip_func_exp)
856 : {
857 : std::vector<const MooseVariableFE<Real> *> score_vars;
858 210 : for (const auto & t : _local_tallies)
859 : {
860 140 : if (t->hasScore(score) && t->name() == tally_name)
861 : {
862 70 : auto vars = t->getScoreVars(score);
863 140 : for (unsigned int ext_bin = 0; ext_bin < vars.size(); ++ext_bin)
864 : {
865 70 : if (skip_func_exp && t->extBinSkipped(ext_bin))
866 0 : continue;
867 70 : score_vars.emplace_back(
868 140 : dynamic_cast<const MooseVariableFE<Real> *>(&getVariable(tid, vars[ext_bin] + output)));
869 : }
870 70 : }
871 : }
872 :
873 70 : if (score_vars.size() == 0)
874 0 : mooseError("No tallies contain the requested score " + score + "!");
875 :
876 70 : return score_vars;
877 0 : }
878 :
879 : std::vector<const VariableValue *>
880 70 : OpenMCCellAverageProblem::getTallyScoreVariableValues(const std::string & score,
881 : const std::string & tally_name,
882 : THREAD_ID tid,
883 : const std::string & output,
884 : bool skip_func_exp)
885 : {
886 : std::vector<const VariableValue *> score_vars;
887 210 : for (const auto & t : _local_tallies)
888 : {
889 140 : if (t->hasScore(score) && t->name() == tally_name)
890 : {
891 70 : auto vars = t->getScoreVars(score);
892 140 : for (unsigned int ext_bin = 0; ext_bin < vars.size(); ++ext_bin)
893 : {
894 70 : if (skip_func_exp && t->extBinSkipped(ext_bin))
895 0 : continue;
896 70 : score_vars.emplace_back(
897 210 : &(dynamic_cast<MooseVariableFE<Real> *>(&getVariable(tid, vars[ext_bin] + output))
898 70 : ->sln()));
899 : }
900 70 : }
901 : }
902 :
903 70 : if (score_vars.size() == 0)
904 0 : mooseError("No tallies contain the requested score " + score + "!");
905 :
906 70 : return score_vars;
907 0 : }
908 :
909 : std::vector<const VariableValue *>
910 12 : OpenMCCellAverageProblem::getTallyScoreNeighborVariableValues(const std::string & score,
911 : const std::string & tally_name,
912 : THREAD_ID tid,
913 : const std::string & output,
914 : bool skip_func_exp)
915 : {
916 : std::vector<const VariableValue *> score_vars;
917 40 : for (const auto & t : _local_tallies)
918 : {
919 28 : if (t->hasScore(score) && t->name() == tally_name)
920 : {
921 12 : auto vars = t->getScoreVars(score);
922 24 : for (unsigned int ext_bin = 0; ext_bin < vars.size(); ++ext_bin)
923 : {
924 12 : if (skip_func_exp && t->extBinSkipped(ext_bin))
925 0 : continue;
926 12 : score_vars.emplace_back(
927 36 : &(dynamic_cast<MooseVariableFE<Real> *>(&getVariable(tid, vars[ext_bin] + output))
928 12 : ->slnNeighbor()));
929 : }
930 12 : }
931 : }
932 :
933 12 : if (score_vars.size() == 0)
934 0 : mooseError("No tallies contain the requested score " + score + "!");
935 :
936 12 : return score_vars;
937 0 : }
938 :
939 : bool
940 16 : OpenMCCellAverageProblem::hasOutput(const std::string & score, const std::string & output) const
941 : {
942 18 : for (const auto & t : _local_tallies)
943 30 : if (std::find(t->getOutputs().begin(), t->getOutputs().end(), output) !=
944 16 : t->getOutputs().end() &&
945 14 : t->hasScore(score))
946 : return true;
947 : return false;
948 : }
949 :
950 : void
951 2378 : OpenMCCellAverageProblem::readBlockParameters(const std::string name,
952 : std::unordered_set<SubdomainID> & blocks)
953 : {
954 2378 : if (isParamValid(name))
955 : {
956 26 : auto names = getParam<std::vector<SubdomainName>>(name);
957 52 : checkEmptyVector(names, "'" + name + "'");
958 :
959 : // here, we do not use the displaced mesh because we need to call this during initial
960 : // setup when the displaced problem does not yet exist. However, displacing the mesh
961 : // should not influence the subdomain IDs anyways
962 26 : auto b_ids = mesh().getSubdomainIDs(names);
963 26 : std::copy(b_ids.begin(), b_ids.end(), std::inserter(blocks, blocks.end()));
964 26 : checkBlocksInMesh(name, b_ids, names);
965 26 : }
966 2378 : }
967 :
968 : void
969 2160 : OpenMCCellAverageProblem::checkBlocksInMesh(const std::string name,
970 : const std::vector<SubdomainID> & ids,
971 : const std::vector<SubdomainName> & names) const
972 : {
973 : // here, we do not use the displaced mesh because we need to call this during initial
974 : // setup when the displaced problem does not yet exist. However, displacing the mesh
975 : // should not influence the subdomain IDs anyways
976 2160 : const auto & subdomains = mesh().meshSubdomains();
977 5402 : for (std::size_t b = 0; b < names.size(); ++b)
978 3242 : if (subdomains.find(ids[b]) == subdomains.end())
979 0 : mooseError("Block '" + names[b] + "' specified in '" + name + "' " + "not found in mesh!");
980 2160 : }
981 :
982 : void
983 2140 : OpenMCCellAverageProblem::read2DBlockParameters(const std::string name,
984 : std::vector<std::vector<SubdomainName>> & names,
985 : std::vector<SubdomainID> & flattened_ids)
986 : {
987 2140 : if (isParamValid(name))
988 : {
989 2140 : names = getParam<std::vector<std::vector<SubdomainName>>>(name);
990 :
991 : // check that entire vector is not empty
992 6418 : checkEmptyVector(names, "'" + name + "'");
993 :
994 : // check that each entry in vector is not empty
995 4384 : for (const auto & n : names)
996 6746 : checkEmptyVector(n, "Entries in '" + name + "'");
997 :
998 : // flatten the 2-d set of names into a 1-d vector
999 : std::vector<SubdomainName> flattened_names;
1000 4380 : for (const auto & slice : names)
1001 5452 : for (const auto & i : slice)
1002 3206 : flattened_names.push_back(i);
1003 :
1004 : // here, we do not use the displaced mesh because we need to call this during initial
1005 : // setup when the displaced problem does not yet exist. However, displacing the mesh
1006 : // should not influence the subdomain IDs anyways
1007 4268 : flattened_ids = mesh().getSubdomainIDs(flattened_names);
1008 4268 : checkBlocksInMesh(name, flattened_ids, flattened_names);
1009 :
1010 : // should not be any duplicate blocks
1011 : std::set<SubdomainName> n;
1012 5332 : for (const auto & b : flattened_names)
1013 : {
1014 : if (n.count(b))
1015 4 : mooseError(
1016 4 : "Subdomains cannot be repeated in '" + name + "'! Subdomain '", b, "' is duplicated.");
1017 3198 : n.insert(b);
1018 : }
1019 2130 : }
1020 2130 : }
1021 :
1022 : coupling::CouplingFields
1023 4151710 : OpenMCCellAverageProblem::elemFeedback(const Elem * elem) const
1024 : {
1025 4151710 : const auto & id = elem->subdomain_id();
1026 : bool has_density =
1027 4151710 : std::find(_density_blocks.begin(), _density_blocks.end(), id) != _density_blocks.end();
1028 4151710 : bool has_temp = std::find(_temp_blocks.begin(), _temp_blocks.end(), id) != _temp_blocks.end();
1029 :
1030 4151710 : if (has_density && has_temp)
1031 : return coupling::density_and_temperature;
1032 2911544 : else if (!has_density && has_temp)
1033 : return coupling::temperature;
1034 1066998 : else if (has_density && !has_temp)
1035 : return coupling::density;
1036 : else
1037 1048266 : return coupling::none;
1038 : }
1039 :
1040 : void
1041 3051 : OpenMCCellAverageProblem::storeElementPhase()
1042 : {
1043 : std::set<SubdomainID> excl_temp_blocks;
1044 : std::set<SubdomainID> excl_density_blocks;
1045 : std::set<SubdomainID> intersect;
1046 :
1047 3051 : std::set<SubdomainID> t(_temp_blocks.begin(), _temp_blocks.end());
1048 3051 : std::set<SubdomainID> d(_density_blocks.begin(), _density_blocks.end());
1049 :
1050 3051 : std::set_difference(t.begin(),
1051 : t.end(),
1052 : d.begin(),
1053 : d.end(),
1054 : std::inserter(excl_temp_blocks, excl_temp_blocks.end()));
1055 :
1056 3051 : std::set_difference(d.begin(),
1057 : d.end(),
1058 : t.begin(),
1059 : t.end(),
1060 : std::inserter(excl_density_blocks, excl_density_blocks.end()));
1061 :
1062 3051 : std::set_intersection(
1063 : t.begin(), t.end(), d.begin(), d.end(), std::inserter(intersect, intersect.begin()));
1064 :
1065 3051 : _n_moose_temp_density_elems = 0;
1066 3619 : for (const auto & s : intersect)
1067 568 : _n_moose_temp_density_elems += numElemsInSubdomain(s);
1068 :
1069 3051 : _n_moose_temp_elems = 0;
1070 5080 : for (const auto & s : excl_temp_blocks)
1071 2029 : _n_moose_temp_elems += numElemsInSubdomain(s);
1072 :
1073 3051 : _n_moose_density_elems = 0;
1074 3215 : for (const auto & s : excl_density_blocks)
1075 164 : _n_moose_density_elems += numElemsInSubdomain(s);
1076 :
1077 3051 : _n_moose_none_elems = getMooseMesh().getMesh().n_active_elem() - _n_moose_temp_density_elems -
1078 3051 : _n_moose_temp_elems - _n_moose_density_elems;
1079 3051 : }
1080 :
1081 : void
1082 3047 : OpenMCCellAverageProblem::computeCellMappedVolumes()
1083 : {
1084 : std::vector<Real> volumes;
1085 :
1086 14672 : for (const auto & c : _local_cell_to_elem)
1087 : {
1088 11625 : Real vol = 0.0;
1089 1991021 : for (const auto & e : c.second)
1090 : {
1091 : // we are looping over local elements, so no need to check for nullptr
1092 1979396 : const auto * elem = getMooseMesh().queryElemPtr(globalElemID(e));
1093 1979396 : vol += elem->volume();
1094 : }
1095 :
1096 11625 : volumes.push_back(vol);
1097 : }
1098 :
1099 3047 : gatherCellSum(volumes, _cell_to_elem_volume);
1100 3047 : }
1101 :
1102 : template <typename T>
1103 : void
1104 16740 : OpenMCCellAverageProblem::gatherCellSum(std::vector<T> & local,
1105 : std::map<cellInfo, T> & global) const
1106 : {
1107 : global.clear();
1108 16740 : _communicator.allgather(local);
1109 :
1110 179283 : for (unsigned int i = 0; i < _flattened_ids.size(); ++i)
1111 : {
1112 : cellInfo cell_info = {_flattened_ids[i], _flattened_instances[i]};
1113 :
1114 : if (global.count(cell_info))
1115 74946 : global[cell_info] += local[i];
1116 : else
1117 87597 : global[cell_info] = local[i];
1118 : }
1119 16740 : }
1120 :
1121 : template <typename T>
1122 : void
1123 6561 : OpenMCCellAverageProblem::gatherCellVector(std::vector<T> & local,
1124 : std::vector<unsigned int> & n_local,
1125 : std::map<cellInfo, std::vector<T>> & global)
1126 : {
1127 : global.clear();
1128 6561 : _communicator.allgather(n_local);
1129 6561 : _communicator.allgather(local);
1130 :
1131 : int e = 0;
1132 63217 : for (unsigned int i = 0; i < _flattened_ids.size(); ++i)
1133 : {
1134 : cellInfo cell_info = {_flattened_ids[i], _flattened_instances[i]};
1135 :
1136 7874760 : for (unsigned int j = e; j < e + n_local[i]; ++j)
1137 7818104 : global[cell_info].push_back(local[j]);
1138 :
1139 56656 : e += n_local[i];
1140 : }
1141 6561 : }
1142 :
1143 : coupling::CouplingFields
1144 9410202 : OpenMCCellAverageProblem::cellFeedback(const cellInfo & cell_info) const
1145 : {
1146 : // _cell_to_elem only holds cells that are coupled by feedback to the [Mesh] (for sake of
1147 : // efficiency in cell-based loops for updating temperatures, densities and
1148 : // extracting the tally). But in some auxiliary kernels, we figure out
1149 : // an element's phase in terms of the cell that it maps to. For these cells that
1150 : // do *map* spatially, but just don't participate in coupling, _cell_to_elem doesn't
1151 : // have any notion of those elements
1152 : if (!_cell_phase.count(cell_info))
1153 0 : return coupling::none;
1154 : else
1155 9410202 : return _cell_phase.at(cell_info);
1156 : }
1157 :
1158 : void
1159 3047 : OpenMCCellAverageProblem::getCellMappedPhase()
1160 : {
1161 : std::vector<int> cells_n_temp;
1162 : std::vector<int> cells_n_temp_rho;
1163 : std::vector<int> cells_n_rho;
1164 : std::vector<int> cells_n_none;
1165 :
1166 : // whether each cell maps to a single phase
1167 14672 : for (const auto & c : _local_cell_to_elem)
1168 : {
1169 11625 : std::vector<int> f(4 /* number of coupling options */, 0);
1170 :
1171 : // we are looping over local elements, so no need to check for nullptr
1172 1991021 : for (const auto & e : c.second)
1173 1979396 : f[elemFeedback(getMooseMesh().queryElemPtr(globalElemID(e)))]++;
1174 :
1175 11625 : cells_n_temp.push_back(f[coupling::temperature]);
1176 11625 : cells_n_temp_rho.push_back(f[coupling::density_and_temperature]);
1177 11625 : cells_n_rho.push_back(f[coupling::density]);
1178 11625 : cells_n_none.push_back(f[coupling::none]);
1179 11625 : }
1180 :
1181 3047 : gatherCellSum(cells_n_temp, _n_temp);
1182 3047 : gatherCellSum(cells_n_temp_rho, _n_temp_rho);
1183 3047 : gatherCellSum(cells_n_rho, _n_rho);
1184 3047 : gatherCellSum(cells_n_none, _n_none);
1185 3047 : }
1186 :
1187 : Real
1188 609312 : OpenMCCellAverageProblem::cellVolume(const cellInfo & cell_info) const
1189 : {
1190 : if (_cell_volume.count(cell_info))
1191 609312 : return _cell_volume.at(cell_info);
1192 : else
1193 0 : return 0.0;
1194 : }
1195 :
1196 : void
1197 3033 : OpenMCCellAverageProblem::checkCellMappedPhase()
1198 : {
1199 3033 : if (_volume_calc)
1200 : {
1201 111 : _volume_calc->initializeVolumeCalculation();
1202 109 : _volume_calc->computeVolumes();
1203 : }
1204 :
1205 : VariadicTable<std::string, int, int, int, int, std::string, std::string> vt(
1206 3031 : {"Cell", " T ", " rho ", "T+rho", "Other", "Mapped Vol", "Actual Vol"});
1207 :
1208 : bool has_mapping = false;
1209 :
1210 : std::vector<Real> cv;
1211 : _cell_phase.clear();
1212 17752 : for (const auto & c : _cell_to_elem)
1213 : {
1214 14725 : auto cell_info = c.first;
1215 14725 : int n_temp = _n_temp[cell_info];
1216 14725 : int n_rho = _n_rho[cell_info];
1217 14725 : int n_temp_rho = _n_temp_rho[cell_info];
1218 14725 : int n_none = _n_none[cell_info];
1219 :
1220 14725 : std::ostringstream vol;
1221 14725 : vol << std::setprecision(3) << std::scientific << "";
1222 14725 : if (_volume_calc)
1223 : {
1224 : Real v, std_dev;
1225 878 : _volume_calc->cellVolume(c.first.first, v, std_dev);
1226 876 : cv.push_back(v);
1227 876 : vol << v << " +/- " << std_dev;
1228 : }
1229 :
1230 14723 : std::ostringstream map;
1231 14723 : map << std::setprecision(3) << std::scientific << _cell_to_elem_volume[cell_info];
1232 :
1233 : // okay to print vol.str() here because only rank 0 is printing (which is the only one
1234 : // with meaningful volume data from OpenMC)
1235 29446 : vt.addRow(printCell(cell_info, true), n_temp, n_rho, n_temp_rho, n_none, map.str(), vol.str());
1236 :
1237 : // cells can only map to a single type of feedback
1238 14723 : std::vector<bool> conditions = {n_temp_rho > 0, n_temp > 0, n_rho > 0, n_none > 0};
1239 14723 : if (std::count(conditions.begin(), conditions.end(), true) > 1)
1240 : {
1241 2 : std::stringstream msg;
1242 2 : std::vector<int> conds = {n_temp, n_rho, n_temp_rho, n_none};
1243 2 : int size = std::to_string(*std::max_element(conds.begin(), conds.end())).length();
1244 4 : msg << "Cell " << printCell(cell_info) << " mapped to:\n\n " << std::setw(size) << n_temp
1245 2 : << " elements with temperature feedback\n " << std::setw(size) << n_rho
1246 2 : << " elements with density feedback\n " << std::setw(size) << n_temp_rho
1247 2 : << " elements with both temperature and density feedback\n " << std::setw(size)
1248 : << n_none
1249 : << " uncoupled elements\n\n"
1250 : "Each OpenMC cell (ID, instance) pair must map to elements of the same coupling "
1251 2 : "settings.";
1252 2 : mooseError(msg.str());
1253 0 : }
1254 :
1255 14721 : if (n_temp)
1256 : {
1257 : has_mapping = true;
1258 10205 : _cell_phase[cell_info] = coupling::temperature;
1259 : }
1260 4516 : else if (n_rho)
1261 : {
1262 : has_mapping = true;
1263 510 : _cell_phase[cell_info] = coupling::density;
1264 : }
1265 4006 : else if (n_temp_rho)
1266 : {
1267 : has_mapping = true;
1268 2000 : _cell_phase[cell_info] = coupling::density_and_temperature;
1269 : }
1270 : else
1271 2006 : _cell_phase[cell_info] = coupling::none;
1272 14721 : }
1273 :
1274 : // collect values from rank 0 onto all other ranks, then populate cell_volume
1275 : // (this is necessary because in OpenMC, the stochastic volume calculation only
1276 : // gets meaningful results on rank 0
1277 3027 : if (_volume_calc)
1278 : {
1279 : _cell_volume.clear();
1280 107 : MPI_Bcast(cv.data(), cv.size(), MPI_DOUBLE, 0, _communicator.get());
1281 : int i = 0;
1282 983 : for (const auto & c : _cell_to_elem)
1283 876 : _cell_volume[c.first] = cv[i++];
1284 : }
1285 :
1286 3027 : if (_specified_density_feedback || _specified_temperature_feedback)
1287 1622 : if (!has_mapping)
1288 2 : mooseError("Feedback was specified using 'temperature_blocks' and/or 'density_blocks', but "
1289 : "no MOOSE elements mapped to OpenMC cells!");
1290 :
1291 3025 : if (_verbose && _cell_to_elem.size())
1292 : {
1293 : _console
1294 1745 : << "\n ===================> MAPPING FROM OPENMC TO MOOSE <===================\n"
1295 1745 : << std::endl;
1296 1745 : _console << " T: # elems providing temperature-only feedback" << std::endl;
1297 1745 : _console << " rho: # elems providing density-only feedback" << std::endl;
1298 1745 : _console << " T+rho: # elems providing temperature and density feedback" << std::endl;
1299 1745 : _console << " Other: # elems which do not provide feedback to OpenMC" << std::endl;
1300 1745 : _console << " (but receives a cell tally from OpenMC)" << std::endl;
1301 1745 : _console << " Mapped Vol: volume of MOOSE elems each cell maps to" << std::endl;
1302 1745 : _console << " Actual Vol: OpenMC cell volume (computed with 'volume_calculation')\n"
1303 1745 : << std::endl;
1304 1745 : vt.print(_console);
1305 : }
1306 :
1307 3025 : printAuxVariableIO();
1308 3025 : _printed_initial = true;
1309 3025 : }
1310 :
1311 : void
1312 3025 : OpenMCCellAverageProblem::printAuxVariableIO()
1313 : {
1314 3025 : if (_printed_initial)
1315 : return;
1316 :
1317 2104 : if (!(_specified_density_feedback || _specified_temperature_feedback ||
1318 : _local_tallies.size() > 0))
1319 : return;
1320 :
1321 1896 : _console << "\n ===================> AUXVARIABLES FOR OPENMC I/O <===================\n"
1322 1896 : << std::endl;
1323 :
1324 1896 : if (_specified_density_feedback || _specified_temperature_feedback)
1325 : {
1326 1497 : _console << " Subdomain: subdomain name/ID" << std::endl;
1327 1497 : _console << " Temperature: variable OpenMC reads temperature from (empty if no feedback)"
1328 1497 : << std::endl;
1329 1497 : _console << " Density: variable OpenMC reads density from (empty if no feedback)\n"
1330 1497 : << std::endl;
1331 :
1332 : VariadicTable<std::string, std::string, std::string> aux(
1333 1497 : {"Subdomain", "Temperature", "Density"});
1334 :
1335 4164 : for (const auto & s : getMooseMesh().meshSubdomains())
1336 : {
1337 2667 : std::string temp = _subdomain_to_temp_vars.count(s) ? _subdomain_to_temp_vars[s].second : "";
1338 : std::string rho =
1339 2667 : _subdomain_to_density_vars.count(s) ? _subdomain_to_density_vars[s].second : "";
1340 :
1341 2667 : if (temp == "" && rho == "")
1342 : continue;
1343 :
1344 4862 : aux.addRow(subdomainName(s), temp, rho);
1345 : }
1346 :
1347 1497 : aux.print(_console);
1348 1497 : _console << std::endl;
1349 1497 : }
1350 :
1351 1896 : if (_local_tallies.size() > 0)
1352 : {
1353 1657 : _console << " Tally Name: Cardinal tally object name" << std::endl;
1354 1657 : _console << " Tally Score: OpenMC tally score" << std::endl;
1355 1657 : _console << " AuxVariable: variable where this score is written\n" << std::endl;
1356 :
1357 : VariadicTable<std::string, std::string, std::string> tallies(
1358 1657 : {"Tally Name", "Tally Score", "AuxVariable(s)"});
1359 4317 : for (unsigned int i = 0; i < _local_tallies.size(); ++i)
1360 : {
1361 : const auto & scores = _local_tallies[i]->getScores();
1362 : const auto & names = _local_tallies[i]->getAuxVarNames();
1363 : const auto bins = _local_tallies[i]->numExtFilterBins();
1364 5794 : for (unsigned int j = 0; j < scores.size(); ++j)
1365 : {
1366 3134 : if (names.size() == 0)
1367 388 : continue;
1368 :
1369 6658 : for (unsigned int k = bins * j; k < (j + 1) * bins; ++k)
1370 : {
1371 3912 : const auto l = j == 0 && k == bins * j ? _local_tallies[i]->name() : "";
1372 3912 : const auto c = k == bins * j ? scores[j] : "";
1373 3912 : const auto r = names[k];
1374 7824 : tallies.addRow(l, c, r);
1375 : }
1376 : }
1377 : }
1378 :
1379 1657 : tallies.print(_console);
1380 1657 : }
1381 : }
1382 :
1383 : void
1384 3047 : OpenMCCellAverageProblem::getCellMappedSubdomains()
1385 : {
1386 : std::vector<unsigned int> n_elems;
1387 : std::vector<unsigned int> elem_ids;
1388 :
1389 14672 : for (const auto & c : _local_cell_to_elem)
1390 : {
1391 11625 : n_elems.push_back(c.second.size());
1392 1991021 : for (const auto & e : c.second)
1393 : {
1394 : // we are looping over local elements, so no need to check for nullptr
1395 1979396 : const auto * elem = getMooseMesh().queryElemPtr(globalElemID(e));
1396 1979396 : elem_ids.push_back(elem->subdomain_id());
1397 : }
1398 : }
1399 :
1400 : std::map<cellInfo, std::vector<unsigned int>> cell_to_subdomain_vec;
1401 3047 : gatherCellVector(elem_ids, n_elems, cell_to_subdomain_vec);
1402 :
1403 : // convert to a set
1404 : _cell_to_elem_subdomain.clear();
1405 17866 : for (const auto & c : cell_to_subdomain_vec)
1406 3833059 : for (const auto & s : c.second)
1407 3818240 : _cell_to_elem_subdomain[c.first].insert(s);
1408 :
1409 : // each cell must map to a consistent setting for identical_cell_fills
1410 : // (all of the blocks it maps to must either _all_ be in the identical blocks,
1411 : // or all excluded)
1412 3047 : if (_has_identical_cell_fills)
1413 : {
1414 1140 : for (const auto & c : _cell_to_elem)
1415 : {
1416 1118 : auto cell_info = c.first;
1417 : bool at_least_one_in = false;
1418 : bool at_least_one_out = false;
1419 : SubdomainID in;
1420 : SubdomainID out;
1421 1118 : auto subdomains = _cell_to_elem_subdomain[cell_info];
1422 2238 : for (const auto & s : subdomains)
1423 : {
1424 1120 : if (_identical_cell_fill_blocks.find(s) == _identical_cell_fill_blocks.end())
1425 : {
1426 : at_least_one_out = true;
1427 606 : out = s;
1428 : }
1429 : else
1430 : {
1431 : at_least_one_in = true;
1432 514 : in = s;
1433 : }
1434 : }
1435 :
1436 1118 : if (at_least_one_in && at_least_one_out)
1437 : {
1438 2 : std::stringstream msg;
1439 2 : msg << "Cell " << printCell(cell_info)
1440 : << " mapped to inconsistent 'identical_cell_fills' settings.\n"
1441 6 : << "Subdomain " << in << " is in 'identical_cell_fills', but " << out << " is not.\n\n"
1442 : << "All subdomains to which this cell maps must either ALL be in "
1443 2 : "'identical_cell_fills' or ALL excluded.";
1444 2 : mooseError(msg.str());
1445 0 : }
1446 : }
1447 : }
1448 3045 : }
1449 :
1450 : std::set<SubdomainID>
1451 1807 : OpenMCCellAverageProblem::coupledSubdomains() const
1452 : {
1453 : std::set<SubdomainID> subdomains;
1454 9380 : for (const auto & c : _cell_to_elem)
1455 : {
1456 7573 : const auto & subdomains_spanning_cell = _cell_to_elem_subdomain.at(c.first);
1457 15247 : for (const auto & s : subdomains_spanning_cell)
1458 7674 : subdomains.insert(s);
1459 : }
1460 :
1461 1807 : return subdomains;
1462 : }
1463 :
1464 : void
1465 1809 : OpenMCCellAverageProblem::subdomainsToMaterials()
1466 : {
1467 1809 : const auto time_start = std::chrono::high_resolution_clock::now();
1468 :
1469 3618 : TIME_SECTION("subdomainsToMaterials", 3, "Mapping OpenMC Materials to Mesh", true);
1470 :
1471 : _subdomain_to_material.clear();
1472 :
1473 9386 : for (const auto & c : _cell_to_elem)
1474 : {
1475 7577 : printTrisoHelp(time_start);
1476 :
1477 7577 : const auto & mats = cellHasIdenticalFill(c.first)
1478 7577 : ? _first_identical_cell_materials
1479 7577 : : materialsInCells(_cell_to_contained_material_cells.at(c.first));
1480 :
1481 15255 : for (const auto & s : _cell_to_elem_subdomain.at(c.first))
1482 17216 : for (const auto & m : mats)
1483 9538 : _subdomain_to_material[s].insert(m);
1484 7577 : }
1485 :
1486 : // Warn the user if a reference density is applied to multiple materials.
1487 4732 : for (const auto & [sub, sub_materials] : _subdomain_to_material)
1488 : {
1489 34 : if (_subdomain_to_ref_density.count(sub) && sub_materials.size() > 1)
1490 : {
1491 : std::string materials;
1492 6 : for (auto mat : sub_materials)
1493 8 : materials += materialName(mat) + ", ";
1494 :
1495 6 : mooseWarning("Reference density " + Moose::stringify(_subdomain_to_ref_density.at(sub)) +
1496 4 : " is being applied to a subdomain (" + subdomainName(sub) +
1497 2 : ") which maps to multiple OpenMC materials: " +
1498 2 : materials.substr(0, materials.size() - 2) +
1499 0 : ". If these multiple materials had different densities during " +
1500 0 : "the MGXS generation stage, your model is not consistently " +
1501 0 : "applying density feedback. The solution is to create a " +
1502 : "separate mesh subdomain for each OpenMC material.");
1503 : }
1504 : }
1505 :
1506 1807 : VariadicTable<std::string, std::string> vt_ce({"Subdomain", "Material"});
1507 : VariadicTable<std::string, std::string, std::string> vt_mg(
1508 1807 : {"Subdomain", "Reference Density", "Material"});
1509 1807 : auto subdomains = coupledSubdomains();
1510 4730 : for (const auto & i : subdomains)
1511 : {
1512 : std::map<std::string, int> mat_to_num;
1513 :
1514 8537 : for (const auto & m : _subdomain_to_material[i])
1515 : {
1516 5614 : auto name = materialName(m);
1517 : if (mat_to_num.count(name))
1518 0 : mat_to_num[name] += 1;
1519 : else
1520 5614 : mat_to_num[name] = 1;
1521 : }
1522 :
1523 2923 : std::string mats = "";
1524 8537 : for (const auto & m : mat_to_num)
1525 : {
1526 5614 : std::string extra = m.second > 1 ? " (" + std::to_string(m.second) + ")" : "";
1527 11228 : mats += " " + m.first + extra + ",";
1528 : }
1529 :
1530 : auto ref_density_str = _subdomain_to_ref_density.count(i)
1531 32 : ? Moose::stringify(_subdomain_to_ref_density.at(i))
1532 2923 : : std::string("");
1533 2923 : mats.pop_back();
1534 2923 : if (openmc::settings::run_CE)
1535 5282 : vt_ce.addRow(subdomainName(i), mats);
1536 : else
1537 564 : vt_mg.addRow(subdomainName(i), ref_density_str, mats);
1538 : }
1539 :
1540 1807 : if (_cell_to_elem.size())
1541 : {
1542 : _console
1543 1750 : << "\n ===================> OPENMC SUBDOMAIN MATERIAL MAPPING <====================\n"
1544 1750 : << std::endl;
1545 1750 : _console << " Subdomain: Subdomain name; if unnamed, we show the ID" << std::endl;
1546 1750 : if (!openmc::settings::run_CE)
1547 96 : _console << " Reference Density: Reference density (kg/m3) applied to the subdomain"
1548 96 : << std::endl;
1549 : _console
1550 1750 : << " Material: OpenMC material name(s) in this subdomain; if unnamed, we\n"
1551 1750 : << " show the ID. If N duplicate material names, we show the\n"
1552 1750 : << " number in ( ).\n"
1553 1750 : << std::endl;
1554 1750 : if (openmc::settings::run_CE)
1555 1654 : vt_ce.print(_console);
1556 : else
1557 96 : vt_mg.print(_console);
1558 1750 : _console << std::endl;
1559 : }
1560 1807 : }
1561 :
1562 : void
1563 3051 : OpenMCCellAverageProblem::initializeElementToCellMapping()
1564 : {
1565 : /* We consider five different cases here based on how the MOOSE and OpenMC
1566 : * domains might overlap in space:
1567 : *
1568 : * 1: Perfect overlap, every MOOSE element maps to an OpenMC cell and every
1569 : * OpenMC cell maps to MOOSE element(s)
1570 : *
1571 : * 2: MOOSE domain fully encloses the OpenMC domain, so that not every MOOSE
1572 : * element maps to an OpenMC cell, but every OpenMC cell maps to a MOOSE element
1573 : *
1574 : * 3: OpenMC domain fully encloses the MOOSE domain, so that not every OpenMC
1575 : * cell maps to MOOSE element(s), but every MOOSE element maps to an OpenMC cell
1576 : *
1577 : * 4: MOOSE and OpenMC domains only partially overlap, so that not every MOOSE
1578 : * element maps to an OpenMC and not every OpenMC cell maps to MOOSE element(s)
1579 : *
1580 : * 5: The MOOSE and OpenMC domains do not overlap at all, so no MOOSE elements
1581 : * map to OpenMC cells and no OpenMC cells map to MOOSE elements.
1582 : *
1583 : * We consider situation #5 to be an error, while the others are technically allowed.
1584 : * We need to error here before getting to OpenMC where we don't map to any cells but
1585 : * would still try to set a cell filter based on no cells.
1586 : */
1587 :
1588 : // First, figure out the phase of each element according to the blocks defined by the user
1589 3051 : storeElementPhase();
1590 :
1591 : // perform element to cell mapping
1592 3051 : mapElemsToCells();
1593 :
1594 3047 : if (!_material_cells_only)
1595 : {
1596 : // gather all cell indices from the initial mapping
1597 : std::vector<int32_t> mapped_cells;
1598 212911 : for (const auto & item : _elem_to_cell)
1599 212444 : mapped_cells.push_back(item.first);
1600 :
1601 467 : std::sort(mapped_cells.begin(), mapped_cells.end());
1602 467 : auto new_end = std::unique(mapped_cells.begin(), mapped_cells.end());
1603 : mapped_cells.erase(new_end, mapped_cells.end());
1604 467 : openmc::prepare_distribcell(&mapped_cells);
1605 :
1606 : // perform element to cell mapping again to get correct instances
1607 467 : mapElemsToCells();
1608 467 : }
1609 :
1610 : // For each cell, get one point inside it to speed up the particle search
1611 3047 : getPointInCell();
1612 :
1613 : // Compute the volume that each OpenMC cell maps to in the MOOSE mesh
1614 3047 : computeCellMappedVolumes();
1615 :
1616 : // Get the number of elements of each phase within the cells
1617 3047 : getCellMappedPhase();
1618 :
1619 : // Get the element subdomains within each cell
1620 3047 : getCellMappedSubdomains();
1621 :
1622 3045 : if (_cell_to_elem.size() == 0 && _has_cell_tallies)
1623 2 : mooseError("Did not find any overlap between MOOSE elements and OpenMC cells for "
1624 : "the specified blocks!");
1625 :
1626 9129 : _console << "\nMapping between " + Moose::stringify(getMooseMesh().getMesh().n_active_elem()) +
1627 9129 : " MOOSE elements and " + Moose::stringify(_n_openmc_cells) +
1628 3043 : " OpenMC cells (on " + Moose::stringify(openmc::model::n_coord_levels) +
1629 6086 : " coordinate levels):"
1630 3043 : << std::endl;
1631 :
1632 : VariadicTable<std::string, int, int, int, int> vt(
1633 3043 : {"", "# T Elems", "# rho Elems", "# T+rho Elems", "# Uncoupled Elems"});
1634 3043 : vt.addRow("MOOSE mesh",
1635 : _n_moose_temp_elems,
1636 : _n_moose_density_elems,
1637 : _n_moose_temp_density_elems,
1638 : _n_moose_none_elems);
1639 3043 : vt.addRow("OpenMC cells",
1640 : _n_mapped_temp_elems,
1641 : _n_mapped_density_elems,
1642 : _n_mapped_temp_density_elems,
1643 : _n_mapped_none_elems);
1644 3043 : vt.print(_console);
1645 3043 : _console << std::endl;
1646 :
1647 3043 : if (_needs_to_map_cells)
1648 : {
1649 2474 : if (_n_moose_temp_elems && (_n_mapped_temp_elems != _n_moose_temp_elems))
1650 90 : mooseWarning("The [Mesh] has " + Moose::stringify(_n_moose_temp_elems) +
1651 : " elements providing temperature feedback (the elements in "
1652 32 : "'temperature_blocks'), but only " +
1653 26 : Moose::stringify(_n_mapped_temp_elems) + " got mapped to OpenMC cells.");
1654 :
1655 2468 : if (_n_moose_temp_elems && (_n_mapped_density_elems != _n_moose_density_elems))
1656 4 : mooseWarning("The [Mesh] has " + Moose::stringify(_n_moose_density_elems) +
1657 : " elements providing density feedback (the elements in "
1658 2 : "'density_blocks'), but only " +
1659 0 : Moose::stringify(_n_mapped_density_elems) + " got mapped to OpenMC cells.");
1660 :
1661 2466 : if (_n_moose_temp_density_elems &&
1662 489 : (_n_mapped_temp_density_elems != _n_moose_temp_density_elems))
1663 28 : mooseWarning("The [Mesh] has " + Moose::stringify(_n_moose_temp_density_elems) +
1664 : " elements providing temperature and density feedback (the elements in the "
1665 10 : "intersection of 'temperature_blocks' and 'density_blocks'), but only " +
1666 8 : Moose::stringify(_n_mapped_temp_density_elems) + " got mapped to OpenMC cells.");
1667 :
1668 2464 : if (_n_mapped_none_elems && (_specified_temperature_feedback || _specified_density_feedback))
1669 524 : mooseWarning("Skipping OpenMC multiphysics feedback from " +
1670 524 : Moose::stringify(_n_mapped_none_elems) +
1671 262 : " [Mesh] elements, which occupy a volume of: " +
1672 524 : Moose::stringify(_uncoupled_volume * _scaling * _scaling * _scaling) + " cm3");
1673 :
1674 2464 : if (_n_openmc_cells < _cell_to_elem.size())
1675 0 : mooseError("Internal error: _cell_to_elem has length ",
1676 0 : _cell_to_elem.size(),
1677 : " which should\n"
1678 : "not exceed the number of OpenMC cells, ",
1679 0 : _n_openmc_cells);
1680 : }
1681 :
1682 : // Check that each cell maps to a single phase
1683 3033 : checkCellMappedPhase();
1684 3025 : }
1685 :
1686 : void
1687 15503 : OpenMCCellAverageProblem::setContainedCells(const cellInfo & cell_info,
1688 : const Point & hint,
1689 : std::map<cellInfo, containedCells> & map)
1690 : {
1691 : containedCells contained_cells;
1692 :
1693 15503 : openmc::Position p{hint(0), hint(1), hint(2)};
1694 :
1695 : // we include all material-fill and void cells within the requested cell, because
1696 : // we may want to add cell tallies in void regions, even if no feedback is applied
1697 : // to voids
1698 15503 : const auto & cell = openmc::model::cells[cell_info.first];
1699 15503 : if (cell->type_ == openmc::Fill::MATERIAL)
1700 : {
1701 14984 : std::vector<int32_t> instances = {cell_info.second};
1702 14984 : contained_cells[cell_info.first] = instances;
1703 14984 : }
1704 : else
1705 1038 : contained_cells = cell->get_contained_cells(cell_info.second, &p);
1706 :
1707 15503 : map[cell_info] = contained_cells;
1708 15503 : }
1709 :
1710 : void
1711 22292 : OpenMCCellAverageProblem::printTrisoHelp(
1712 : const std::chrono::time_point<std::chrono::high_resolution_clock> & start) const
1713 : {
1714 22292 : if (!_printed_triso_warning)
1715 : {
1716 22292 : auto stop = std::chrono::high_resolution_clock::now();
1717 22292 : auto elapsed = std::chrono::duration<double, std::milli>(stop - start).count() / 1e3;
1718 22292 : if (elapsed > 120.0)
1719 : {
1720 0 : _printed_triso_warning = true;
1721 0 : _console << "\nThis is taking a long time. Does your problem have TRISOs/other "
1722 : << "highly heterogeneous geometry?\nIf you are repeating the same TRISO/etc. "
1723 0 : "universe many times "
1724 : << "through your OpenMC model, setting\n'identical_cell_fills' will give you a big "
1725 0 : "speedup.\n\n"
1726 : << "For more information, consult the Cardinal documentation: "
1727 0 : "https://tinyurl.com/54kz9aw8"
1728 0 : << std::endl;
1729 : }
1730 : }
1731 22292 : }
1732 :
1733 : void
1734 3025 : OpenMCCellAverageProblem::cacheContainedCells()
1735 : {
1736 6050 : TIME_SECTION("cacheContainedCells", 3, "Caching Contained Cells", true);
1737 :
1738 : bool first_cell = true;
1739 : bool second_cell = false;
1740 : containedCells first_cell_cc;
1741 : containedCells second_cell_cc;
1742 : bool used_cache_shortcut = false;
1743 :
1744 : _cell_to_contained_material_cells.clear();
1745 3025 : _first_identical_cell_materials.clear();
1746 : _instance_offsets.clear();
1747 : _n_offset.clear();
1748 :
1749 : int n = -1;
1750 3025 : const auto time_start = std::chrono::high_resolution_clock::now();
1751 17740 : for (const auto & c : _cell_to_elem)
1752 : {
1753 14715 : auto cell_info = c.first;
1754 14715 : Point hint = transformPointToOpenMC(_cell_to_point[cell_info]);
1755 :
1756 14715 : printTrisoHelp(time_start);
1757 :
1758 : // default to the normal behavior
1759 14715 : if (!cellHasIdenticalFill(cell_info))
1760 14403 : setContainedCells(cell_info, hint, _cell_to_contained_material_cells);
1761 : else
1762 : {
1763 : used_cache_shortcut = true;
1764 312 : _n_offset[cell_info] = ++n;
1765 :
1766 312 : if (first_cell)
1767 : {
1768 12 : setContainedCells(cell_info, hint, _cell_to_contained_material_cells);
1769 12 : first_cell_cc = _cell_to_contained_material_cells[cell_info];
1770 : _first_identical_cell = cell_info;
1771 24 : _first_identical_cell_materials = materialsInCells(first_cell_cc);
1772 : first_cell = false;
1773 : second_cell = true;
1774 : }
1775 300 : else if (second_cell)
1776 : {
1777 12 : setContainedCells(cell_info, hint, _cell_to_contained_material_cells);
1778 12 : second_cell_cc = _cell_to_contained_material_cells[cell_info];
1779 : second_cell = false;
1780 :
1781 : // we will check for equivalence in the end mapping later; but here we still need
1782 : // some checks to make sure the structure is compatible
1783 12 : checkContainedCellsStructure(cell_info, first_cell_cc, second_cell_cc);
1784 :
1785 : // get the offset for each instance for each contained cell
1786 6882 : for (const auto & f : first_cell_cc)
1787 : {
1788 6870 : const auto id = f.first;
1789 : const auto & instances = f.second;
1790 : const auto & new_instances = second_cell_cc[id];
1791 :
1792 : std::vector<int32_t> offsets;
1793 432036 : for (unsigned int i = 0; i < instances.size(); ++i)
1794 425166 : offsets.push_back(new_instances[i] - instances[i]);
1795 :
1796 6870 : _instance_offsets[id] = offsets;
1797 6870 : }
1798 : }
1799 : }
1800 : }
1801 :
1802 : // only need to check if we were attempting the shortcut
1803 3025 : if (_check_identical_cell_fills)
1804 : {
1805 40 : TIME_SECTION("verifyCacheContainedCells", 4, "Verifying Cached Contained Cells", true);
1806 :
1807 : std::map<cellInfo, containedCells> checking_cell_fills;
1808 1096 : for (const auto & c : _cell_to_elem)
1809 1076 : setContainedCells(
1810 2152 : c.first, transformPointToOpenMC(_cell_to_point[c.first]), checking_cell_fills);
1811 :
1812 : std::map<cellInfo, containedCells> current_cell_fills;
1813 1096 : for (const auto & c : _cell_to_elem)
1814 : {
1815 : // Shift the cell instances in-place.
1816 1076 : if (cellHasIdenticalFill(c.first))
1817 : {
1818 312 : current_cell_fills[c.first] = _cell_to_contained_material_cells.at(_first_identical_cell);
1819 187044 : for (auto & [cc_idx, cc_instances] : current_cell_fills[c.first])
1820 11676744 : for (unsigned int instance_idx = 0; instance_idx < cc_instances.size(); instance_idx++)
1821 11490012 : cc_instances[instance_idx] = containedCellInstanceShift(c.first, cc_idx, instance_idx);
1822 : }
1823 : else
1824 764 : current_cell_fills[c.first] = _cell_to_contained_material_cells.at(c.first);
1825 : }
1826 :
1827 : std::map<cellInfo, containedCells> ordered_reference(checking_cell_fills.begin(),
1828 20 : checking_cell_fills.end());
1829 : std::map<cellInfo, containedCells> ordered(current_cell_fills.begin(),
1830 20 : current_cell_fills.end());
1831 20 : compareContainedCells(ordered_reference, ordered);
1832 16 : }
1833 :
1834 3021 : if (_has_identical_cell_fills && !used_cache_shortcut)
1835 10 : mooseWarning("You specified 'identical_cell_fills', but all cells which mapped to these "
1836 : "subdomains were filled \n"
1837 : "by a material (as opposed to a universe/lattice), so the 'identical_cell_fills' "
1838 : "parameter is unused.");
1839 :
1840 : // Check for duplicate contained cells to ensure we don't set the cell temperature or
1841 : // density multiple times erroneously. This occurs if Cardinal maps to multiple cells
1842 : // that are each filled with the same lattice, as OpenMC doesn't add cell instances in
1843 : // that case (lattices must be placed in universes, then those universes may be placed
1844 : // in cells). This check is memory-intensive as it builds a map of all contained
1845 : // cell index + instance pairs, so we ensure its only performed on a single MPI rank
1846 : // and outside of any point where the auxsystem is serialized to reduce peak
1847 : // memory consumption.
1848 3019 : if (_communicator.rank() == 0)
1849 : {
1850 : std::unordered_set<cellInfo> cells_already_set;
1851 8668 : for (const auto & [cell_info, elements] : _cell_to_elem)
1852 : {
1853 : // Skip checking the identical cell fills outside of _first_identical_cell.
1854 : // These mapping errors are caught when verifying contained cells above (if
1855 : // 'check_identical_cell_fills' is true)
1856 7311 : const bool identical_fill = cellHasIdenticalFill(cell_info);
1857 7311 : if (identical_fill && cell_info != _first_identical_cell)
1858 104 : continue;
1859 :
1860 7207 : const auto & contained_cells = _cell_to_contained_material_cells.at(cell_info);
1861 17702 : for (auto & [cc_idx, cc_instances] : contained_cells)
1862 : {
1863 228994 : for (auto cc_instance : cc_instances)
1864 : {
1865 218497 : if (cells_already_set.count({cc_idx, cc_instance}))
1866 6 : mooseError("Cell " + std::to_string(cellID(cc_idx)) + ", instance " +
1867 2 : std::to_string(cc_instance) +
1868 : " has already had its properties (temperature and/or density) set by "
1869 : "Cardinal! This indicates a problem with how you have built your "
1870 : "geometry, because this cell is trying to receive a distribution of "
1871 : "temperatures or densities in space, but each successive set-property "
1872 : "operation is only overwriting the previous value.\n\nThis "
1873 : "error most often appears when you are filling a LATTICE into multiple "
1874 : "cells. One fix is to first place that lattice into a universe, and then "
1875 : "fill that UNIVERSE into multiple cells.\n\nFor more information, "
1876 : "please consult https://github.com/neams-th-coe/cardinal/pull/918.");
1877 :
1878 218497 : cells_already_set.insert({cc_idx, cc_instance});
1879 : }
1880 : }
1881 : }
1882 : }
1883 : // All MPI ranks need to wait until rank zero has finished performing the mapping check.
1884 3017 : _communicator.barrier();
1885 3017 : }
1886 :
1887 : void
1888 1002 : OpenMCCellAverageProblem::checkContainedCellsStructure(const cellInfo & cell_info,
1889 : containedCells & reference,
1890 : containedCells & compare) const
1891 : {
1892 : // make sure the number of keys is the same
1893 1002 : if (reference.size() != compare.size())
1894 0 : mooseError("The cell caching failed to identify identical number of cell IDs filling cell " +
1895 0 : printCell(cell_info) + "\nYou must unset 'identical_cell_fills'");
1896 :
1897 193190 : for (const auto & entry : reference)
1898 : {
1899 192190 : const auto & key = entry.first;
1900 :
1901 : // check that each key exists
1902 : if (!compare.count(key))
1903 6 : mooseError("Not all cells contain cell ID " + Moose::stringify(cellID(key)) +
1904 4 : ". The offender is: cell " + printCell(cell_info) +
1905 : ".\nYou must unset 'identical_cell_fills'!");
1906 :
1907 : // for each int32_t key, compare the std::vector<int32_t> map
1908 : const auto & reference_instances = entry.second;
1909 : const auto & compare_instances = compare[key];
1910 :
1911 : // they should have the same number of instances
1912 192188 : if (reference_instances.size() != compare_instances.size())
1913 0 : mooseError("The cell caching should have identified " +
1914 0 : Moose::stringify(reference_instances.size()) + "cell instances in cell ID " +
1915 0 : Moose::stringify(cellID(key)) + ", but instead found " +
1916 : Moose::stringify(compare_instances.size()) +
1917 : "\nYou must unset 'identical_cell_fills'");
1918 : }
1919 1000 : }
1920 :
1921 : void
1922 20 : OpenMCCellAverageProblem::compareContainedCells(std::map<cellInfo, containedCells> & reference,
1923 : std::map<cellInfo, containedCells> & compare) const
1924 : {
1925 : // check that the number of keys matches
1926 20 : if (reference.size() != compare.size())
1927 0 : mooseError("The cell caching should have identified " + Moose::stringify(reference.size()) +
1928 : " cells, but instead "
1929 0 : "found " +
1930 : Moose::stringify(compare.size()));
1931 :
1932 : // loop over each cellInfo
1933 1006 : for (const auto & entry : reference)
1934 : {
1935 990 : auto cell_info = entry.first;
1936 :
1937 : // make sure the key exists
1938 : if (!compare.count(cell_info))
1939 0 : mooseError("The cell caching failed to map cell " + printCell(cell_info));
1940 :
1941 : // for each cellInfo key, compare the contained cells map
1942 990 : auto reference_map = reference[cell_info];
1943 990 : auto compare_map = compare[cell_info];
1944 :
1945 990 : checkContainedCellsStructure(cell_info, reference_map, compare_map);
1946 :
1947 : // loop over each contained cell
1948 185886 : for (const auto & nested_entry : reference_map)
1949 : {
1950 : // for each int32_t key, compare the std::vector<int32_t> map
1951 184900 : auto reference_instances = nested_entry.second;
1952 184900 : auto compare_instances = compare_map[nested_entry.first];
1953 :
1954 184900 : std::sort(reference_instances.begin(), reference_instances.end());
1955 184900 : std::sort(compare_instances.begin(), compare_instances.end());
1956 :
1957 : // and the instances should exactly match
1958 184900 : if (reference_instances != compare_instances)
1959 2 : mooseError(
1960 2 : "The cell caching failed to get correct instances for material cell ID " +
1961 4 : Moose::stringify(cellID(nested_entry.first)) + " within cell " + printCell(cell_info) +
1962 0 : ". You must unset 'identical_cell_fills'!" + "\n\nThis error might appear if:\n" +
1963 : " - There is a mismatch between your OpenMC model and the [Mesh]\n"
1964 : " - There are additional OpenMC cells filled with this repeatable universe/lattice, "
1965 : "but which are not mapping to the blocks in 'identical_cell_fills'");
1966 184898 : }
1967 : }
1968 16 : }
1969 :
1970 : std::vector<int32_t>
1971 2198 : OpenMCCellAverageProblem::getMappedTallyIDs() const
1972 : {
1973 : std::vector<int32_t> tally_ids;
1974 :
1975 : // local mapped tallies
1976 5466 : for (const auto & t : _local_tallies)
1977 3268 : tally_ids.push_back(t->getTallyID());
1978 : // global normalization tallies
1979 5466 : for (const auto & t : _local_tallies)
1980 3268 : if (t->addingGlobalTally())
1981 1532 : tally_ids.push_back(t->getGlobalTallyID());
1982 :
1983 2198 : return tally_ids;
1984 0 : }
1985 :
1986 : unsigned int
1987 2082896 : OpenMCCellAverageProblem::getCellLevel(const Point & c) const
1988 : {
1989 2082896 : unsigned int level = _cell_level;
1990 2082896 : if (_cell_level > _particle.n_coord() - 1)
1991 : {
1992 8388 : if (isParamValid("lowest_cell_level"))
1993 4192 : level = _particle.n_coord() - 1;
1994 : else
1995 : {
1996 2 : std::string l = Moose::stringify(_cell_level);
1997 4 : mooseError("Requested coordinate level of " + l +
1998 6 : " exceeds number of nested coordinate levels at " + printPoint(c) + ": " +
1999 2 : Moose::stringify(_particle.n_coord()) +
2000 : ".\n\nYou can either change how the OpenMC model is built by nesting universes "
2001 2 : "into deeper levels, or you can try setting 'lowest_cell_level = " +
2002 0 : l +
2003 : "', which will couple on the lowest level found in the geometry at any given x, "
2004 2 : "y, z point, up to and including level " +
2005 0 : l + ".");
2006 : }
2007 : }
2008 :
2009 2082894 : return level;
2010 : }
2011 :
2012 : void
2013 3518 : OpenMCCellAverageProblem::mapElemsToCells()
2014 : {
2015 : // reset counters, flags
2016 3518 : _n_mapped_temp_elems = 0;
2017 3518 : _n_mapped_density_elems = 0;
2018 3518 : _n_mapped_temp_density_elems = 0;
2019 3518 : _n_mapped_none_elems = 0;
2020 3518 : _uncoupled_volume = 0.0;
2021 3518 : _material_cells_only = true;
2022 :
2023 : // reset data structures
2024 3518 : _elem_to_cell.clear();
2025 : _cell_to_elem.clear();
2026 3518 : _flattened_ids.clear();
2027 3518 : _flattened_instances.clear();
2028 :
2029 : int local_elem = -1;
2030 4241200 : for (unsigned int e = 0; e < getMooseMesh().nElem(); ++e)
2031 : {
2032 4237686 : const auto * elem = getMooseMesh().queryElemPtr(e);
2033 :
2034 4237686 : if (!isLocalElem(elem) || !elem->active())
2035 2065372 : continue;
2036 :
2037 2198826 : local_elem++;
2038 :
2039 2198826 : auto id = elem->subdomain_id();
2040 2198826 : const Point & c = elem->vertex_average();
2041 2198826 : Real element_volume = elem->volume();
2042 :
2043 : // find the OpenMC cell at the location 'c' (if any)
2044 2198826 : bool error = findCell(c);
2045 :
2046 : // if we didn't find an OpenMC cell here, then we certainly have an uncoupled region
2047 2198826 : if (error)
2048 : {
2049 26512 : _uncoupled_volume += element_volume;
2050 26512 : _n_mapped_none_elems++;
2051 26512 : continue;
2052 : }
2053 :
2054 : // next, see what type of data is to be sent into OpenMC (to further classify
2055 : // the type of couling)
2056 2172314 : auto phase = elemFeedback(elem);
2057 :
2058 : // Loop over the tallies to check if any CellTally objects map to this element.
2059 : bool elem_mapped_to_cell_tally = false;
2060 4797944 : for (const auto & tally : _local_tallies)
2061 : {
2062 2625630 : auto cell_tally = dynamic_cast<const CellTally *>(tally.get());
2063 2625630 : if (cell_tally)
2064 2346396 : elem_mapped_to_cell_tally |=
2065 : cell_tally->getBlocks().find(id) != cell_tally->getBlocks().end();
2066 : }
2067 :
2068 2172314 : bool requires_mapping = phase != coupling::none || elem_mapped_to_cell_tally;
2069 :
2070 : // get the level in the OpenMC model to fetch mapped cell information. For
2071 : // uncoupled regions, we know we will be successful in finding a cell (because
2072 : // we already screened out uncoupled cells), and the id and instance are unused
2073 : // (so we can just set zero).
2074 2172314 : auto level = requires_mapping ? getCellLevel(c) : 0;
2075 :
2076 : // ensure the mapped cell isn't in a unvierse being used as the "outer"
2077 : // universe of a lattice in the OpenMC model
2078 : if (requires_mapping)
2079 2082894 : latticeOuterCheck(c, level);
2080 :
2081 2082892 : switch (phase)
2082 : {
2083 627054 : case coupling::density_and_temperature:
2084 : {
2085 627054 : _n_mapped_temp_density_elems++;
2086 627054 : break;
2087 : }
2088 964806 : case coupling::temperature:
2089 : {
2090 964806 : _n_mapped_temp_elems++;
2091 964806 : break;
2092 : }
2093 11608 : case coupling::density:
2094 : {
2095 11608 : _n_mapped_density_elems++;
2096 11608 : break;
2097 : }
2098 568842 : case coupling::none:
2099 : {
2100 568842 : _uncoupled_volume += element_volume;
2101 568842 : _n_mapped_none_elems++;
2102 568842 : break;
2103 : }
2104 0 : default:
2105 0 : mooseError("Unhandled CouplingFields enum!");
2106 : }
2107 :
2108 2172310 : auto cell_index = _particle.coord(level).cell();
2109 2172310 : auto cell_instance = cell_instance_at_level(_particle, level);
2110 :
2111 : cellInfo cell_info = {cell_index, cell_instance};
2112 :
2113 2172310 : if (openmc::model::cells[cell_index]->type_ != openmc::Fill::MATERIAL)
2114 163584 : _material_cells_only = false;
2115 :
2116 : // store the map of cells to elements that will be coupled via feedback or a tally
2117 2172310 : if (requires_mapping)
2118 2082892 : _cell_to_elem[cell_info].push_back(local_elem);
2119 : }
2120 :
2121 3514 : _communicator.sum(_n_mapped_temp_elems);
2122 3514 : _communicator.sum(_n_mapped_temp_density_elems);
2123 3514 : _communicator.sum(_n_mapped_density_elems);
2124 3514 : _communicator.sum(_n_mapped_none_elems);
2125 3514 : _communicator.sum(_uncoupled_volume);
2126 :
2127 : // if ANY rank finds a non-material cell, they will hold 0 (false)
2128 3514 : _communicator.min(_material_cells_only);
2129 :
2130 : // store the local mapping of cells to elements for convenience
2131 : _local_cell_to_elem = _cell_to_elem;
2132 :
2133 : // flatten the cell IDs and instances
2134 15773 : for (const auto & c : _cell_to_elem)
2135 : {
2136 12259 : auto cell_info = c.first;
2137 12259 : _flattened_ids.push_back(cell_info.first);
2138 12259 : _flattened_instances.push_back(cell_info.second);
2139 : }
2140 :
2141 3514 : _communicator.allgather(_flattened_ids);
2142 3514 : _communicator.allgather(_flattened_instances);
2143 :
2144 : // collect the _cell_to_elem onto all ranks
2145 3514 : std::vector<unsigned int> n_elems;
2146 3514 : std::vector<unsigned int> elems;
2147 15773 : for (const auto & c : _cell_to_elem)
2148 : {
2149 12259 : n_elems.push_back(c.second.size());
2150 2093615 : for (const auto & e : c.second)
2151 2081356 : elems.push_back(_local_to_global_elem[e]);
2152 : }
2153 :
2154 3514 : gatherCellVector(elems, n_elems, _cell_to_elem);
2155 :
2156 : // fill out the elem_to_cell structure
2157 : // TODO: figure out how to shrink this so we only store the mapping for active
2158 : // elements as opposed to the entire element hierarchy.
2159 3514 : _elem_to_cell.resize(getMooseMesh().nElem(), {UNMAPPED, UNMAPPED});
2160 19127 : for (const auto & c : _cell_to_elem)
2161 4015477 : for (const auto & e : c.second)
2162 3999864 : _elem_to_cell[e] = c.first;
2163 3514 : }
2164 :
2165 : void
2166 3047 : OpenMCCellAverageProblem::getPointInCell()
2167 : {
2168 : std::vector<Real> x;
2169 : std::vector<Real> y;
2170 : std::vector<Real> z;
2171 14672 : for (const auto & c : _local_cell_to_elem)
2172 : {
2173 : // we are only dealing with local elements here, no need to check for nullptr
2174 11625 : const Elem * elem = getMooseMesh().queryElemPtr(globalElemID(c.second[0]));
2175 11625 : const Point & p = elem->vertex_average();
2176 :
2177 11625 : x.push_back(p(0));
2178 11625 : y.push_back(p(1));
2179 11625 : z.push_back(p(2));
2180 : }
2181 :
2182 3047 : _communicator.allgather(x);
2183 3047 : _communicator.allgather(y);
2184 3047 : _communicator.allgather(z);
2185 :
2186 : // this will get a point from the lowest rank in each cell
2187 : _cell_to_point.clear();
2188 30690 : for (unsigned int i = 0; i < _flattened_ids.size(); ++i)
2189 : {
2190 : cellInfo cell_info = {_flattened_ids[i], _flattened_instances[i]};
2191 : if (!_cell_to_point.count(cell_info))
2192 14819 : _cell_to_point[cell_info] = Point(x[i], y[i], z[i]);
2193 : }
2194 3047 : }
2195 :
2196 : void
2197 921 : OpenMCCellAverageProblem::resetTallies()
2198 : {
2199 921 : if (_local_tallies.size() == 0)
2200 : return;
2201 :
2202 : // We initialize [Problem/Tallies] by forward iterating this vector. We need to delete them in
2203 : // reverse.
2204 1185 : for (int i = _local_tallies.size() - 1; i >= 0; --i)
2205 628 : _local_tallies[i]->resetTally();
2206 : }
2207 :
2208 : void
2209 3015 : OpenMCCellAverageProblem::initializeTallies()
2210 : {
2211 : // add trigger information for k, if present
2212 3015 : openmc::settings::keff_trigger.metric = triggerMetric(_k_trigger);
2213 :
2214 3015 : if (_local_tallies.size() == 0)
2215 : return;
2216 :
2217 : // Initialize all of the [Problem/Tallies].
2218 5478 : for (auto & local_tally : _local_tallies)
2219 3280 : local_tally->initializeTally();
2220 :
2221 : // Ensure that any tally editors don't apply to mapped tallies
2222 2198 : checkTallyEditorIDs();
2223 : }
2224 :
2225 : void
2226 2 : OpenMCCellAverageProblem::latticeOuterError(const Point & c, int level) const
2227 : {
2228 2 : const auto & cell = openmc::model::cells[_particle.coord(level).cell()];
2229 2 : std::stringstream msg;
2230 2 : msg << "The point " << c << " mapped to cell " << cell->id_
2231 : << " in the OpenMC model is inside a universe "
2232 : "used as the 'outer' universe of a lattice. "
2233 : "All cells used for mapping in lattices must be explicitly set "
2234 : "on the 'universes' attribute of lattice objects. "
2235 : << "If you want to obtain feedback or cell tallies here, you "
2236 : "will need to widen your lattice to have universes covering all of the space you "
2237 : "want feedback or cell tallies.\n\nIn other words, re-build your OpenMC model but replace "
2238 : "lattice.outer by simply creating extra rings/rows in your lattice to cover all the space "
2239 : "needed. For more information, see: "
2240 2 : "https://github.com/openmc-dev/openmc/issues/551.";
2241 2 : mooseError(msg.str());
2242 0 : }
2243 :
2244 : void
2245 2082894 : OpenMCCellAverageProblem::latticeOuterCheck(const Point & c, int level) const
2246 : {
2247 4718240 : for (int i = 0; i <= level; ++i)
2248 : {
2249 : const auto & coord = _particle.coord(i);
2250 :
2251 : // if there is no lattice at this level, move on
2252 2635348 : if (coord.lattice() == openmc::C_NONE)
2253 2169134 : continue;
2254 :
2255 466214 : const auto & lat = openmc::model::lattices[coord.lattice()];
2256 :
2257 : // if the lattice's outer universe isn't set, move on
2258 466214 : if (lat->outer_ == openmc::NO_OUTER_UNIVERSE)
2259 1280 : continue;
2260 :
2261 464934 : if (coord.universe() != lat->outer_)
2262 464932 : continue;
2263 :
2264 : // move on if the lattice indices are valid (position is in the set of explicitly defined
2265 : // universes)
2266 2 : if (lat->are_valid_indices(coord.lattice_index()))
2267 0 : continue;
2268 :
2269 : // if we get here, the mapping is occurring in a universe that is not explicitly defined in the
2270 : // lattice
2271 2 : latticeOuterError(c, level);
2272 : }
2273 2082892 : }
2274 :
2275 : bool
2276 2198826 : OpenMCCellAverageProblem::findCell(const Point & point)
2277 : {
2278 2198826 : _particle.clear();
2279 : // Use a random direction to minimize "lost" virtual particles.
2280 2198826 : _particle.u() = {0.6339976, -0.538536, 0.555026};
2281 2198826 : _particle.u() /= _particle.u().norm();
2282 :
2283 2198826 : Point pt = transformPointToOpenMC(point);
2284 :
2285 2198826 : _particle.r() = {pt(0), pt(1), pt(2)};
2286 2198826 : return !openmc::exhaustive_find_cell(_particle);
2287 : }
2288 :
2289 : void
2290 2255 : OpenMCCellAverageProblem::addExternalVariables()
2291 : {
2292 : // We need to validate tallies here to we can add scores that may be missing.
2293 2255 : validateLocalTallies();
2294 :
2295 : // Add all of the auxvariables in which the [Tallies] block will store results.
2296 : unsigned int previous_valid_name_index = 0;
2297 5008 : for (unsigned int i = 0; i < _local_tallies.size(); ++i)
2298 : {
2299 2769 : _tally_var_ids.emplace_back();
2300 :
2301 : // Convert the subdomain ID map into a std::vector for addExternalVariable(...).
2302 : std::vector<SubdomainName> block_name_vec;
2303 6697 : for (const auto b : _local_tallies[i]->getBlocks())
2304 7856 : block_name_vec.emplace_back(mesh().getSubdomainName(b) != "" ? mesh().getSubdomainName(b)
2305 : : std::to_string(b));
2306 :
2307 : // We use this to check if a sequence of added tallies corresponds to a single translated mesh.
2308 : // If the number of names reported in getAuxVarNames is zero, the tally must store it's results
2309 : // in the variables added by the first mesh tally in the sequence.
2310 : bool is_instanced = _local_tallies[i]->getAuxVarNames().size() == 0;
2311 2769 : previous_valid_name_index = !is_instanced ? i : previous_valid_name_index;
2312 :
2313 2769 : const auto & names = _local_tallies[previous_valid_name_index]->getAuxVarNames();
2314 :
2315 2769 : _tally_ext_var_ids.emplace_back();
2316 2769 : if (_local_tallies[i]->hasOutputs())
2317 265 : _tally_ext_var_ids[i].resize(_local_tallies[i]->getOutputs().size());
2318 :
2319 7328 : for (unsigned int j = 0; j < names.size(); ++j)
2320 : {
2321 4563 : if (is_instanced)
2322 544 : _tally_var_ids[i].push_back(
2323 : _tally_var_ids[previous_valid_name_index][j]); // Use variables from first in sequence.
2324 : else
2325 8034 : _tally_var_ids[i].push_back(addExternalVariable(names[j], "Tally", &block_name_vec));
2326 :
2327 4559 : if (_local_tallies[i]->hasOutputs())
2328 : {
2329 : const auto & outs = _local_tallies[i]->getOutputs();
2330 872 : for (std::size_t k = 0; k < outs.size(); ++k)
2331 : {
2332 513 : std::string n = names[j] + "_" + outs[k];
2333 513 : if (is_instanced)
2334 16 : _tally_ext_var_ids[i][k].push_back(
2335 : _tally_ext_var_ids[previous_valid_name_index][k]
2336 : [j]); // Use variables from first in sequence.
2337 : else
2338 994 : _tally_ext_var_ids[i][k].push_back(addExternalVariable(n, "Tally", &block_name_vec));
2339 : }
2340 : }
2341 : }
2342 2765 : }
2343 :
2344 : // create the variable(s) that will be used to receive density
2345 : _subdomain_to_density_vars.clear();
2346 2836 : for (const auto & v : _density_vars_to_blocks)
2347 : {
2348 597 : auto number = addExternalVariable(v.first, "density feedback", &v.second);
2349 :
2350 597 : auto ids = getMooseMesh().getSubdomainIDs(v.second);
2351 1289 : for (const auto & s : ids)
2352 1384 : _subdomain_to_density_vars[s] = {number, v.first};
2353 597 : }
2354 :
2355 : // create the variable(s) that will be used to receive temperature
2356 : _subdomain_to_temp_vars.clear();
2357 3762 : for (const auto & v : _temp_vars_to_blocks)
2358 : {
2359 1523 : auto number = addExternalVariable(v.first, "temperature feedback", &v.second);
2360 :
2361 1523 : auto ids = getMooseMesh().getSubdomainIDs(v.second);
2362 3888 : for (const auto & s : ids)
2363 4730 : _subdomain_to_temp_vars[s] = {number, v.first};
2364 1523 : }
2365 :
2366 2239 : if (_output_cell_mapping && _needs_to_map_cells)
2367 : {
2368 1847 : std::string auxk_type = "CellIDAux";
2369 1847 : InputParameters params = _factory.getValidParams(auxk_type);
2370 3694 : addExternalVariable("cell_id", "cell mapping");
2371 3694 : params.set<AuxVariableName>("variable") = "cell_id";
2372 3694 : addAuxKernel(auxk_type, "cell_id", params);
2373 :
2374 : auxk_type = "CellInstanceAux";
2375 1847 : params = _factory.getValidParams(auxk_type);
2376 3694 : addExternalVariable("cell_instance", "cell mapping");
2377 3694 : params.set<AuxVariableName>("variable") = "cell_instance";
2378 1847 : addAuxKernel(auxk_type, "cell_instance", params);
2379 3694 : }
2380 : else
2381 : _console << "Skipping output of 'cell_id' and 'cell_instance' because 'temperature_blocks', "
2382 392 : "'density_blocks', and 'tally_blocks' are all empty"
2383 392 : << std::endl;
2384 2239 : }
2385 :
2386 : void
2387 2759 : OpenMCCellAverageProblem::externalSolve()
2388 : {
2389 : // if using Dufek-Gudowski acceleration and this is not the first iteration, update
2390 : // the number of particles; we put this here so that changing the number of particles
2391 : // doesn't intrude with any other postprocessing routines that happen outside this class's purview
2392 2759 : if (_relaxation == relaxation::dufek_gudowski && !firstSolve())
2393 32 : dufekGudowskiParticleUpdate();
2394 : else
2395 : {
2396 5454 : if (isParamValid("particles"))
2397 : {
2398 251 : if (*_particles <= 0.0)
2399 2 : mooseError(
2400 : "'particles' must be a positive integer. Try `execute_on = 'timestep_begin'` in "
2401 : "your postprocessor and check that the postprocessor value itself is not less than "
2402 : "or equal to zero.");
2403 249 : int64_t n = std::llround(*_particles);
2404 249 : openmc::settings::n_particles = n;
2405 : }
2406 : }
2407 :
2408 2757 : OpenMCProblemBase::externalSolve();
2409 2747 : }
2410 :
2411 : std::map<OpenMCCellAverageProblem::cellInfo, Real>
2412 1505 : OpenMCCellAverageProblem::computeVolumeWeightedCellInput(
2413 : const std::map<SubdomainID, std::pair<unsigned int, std::string>> & var_num,
2414 : const std::vector<coupling::CouplingFields> * phase,
2415 : const std::map<SubdomainID, Real> * scaling) const
2416 : {
2417 1505 : const auto & sys_number = _aux->number();
2418 :
2419 : // collect the volume-weighted product across local ranks
2420 : std::vector<Real> volume_product;
2421 12045 : for (const auto & c : _local_cell_to_elem)
2422 : {
2423 : // if a specific phase is passed in, only evaluate for those elements in the phase;
2424 : // in order to have the correct array sizes for gatherCellSum, we set zero values
2425 : // for any cells that aren't in the correct phase, and leave it up to the send...ToOpenMC()
2426 : // routines to properly shield against incorrect phases
2427 10540 : if (phase)
2428 : {
2429 10540 : if (std::find(phase->begin(), phase->end(), cellFeedback(c.first)) == phase->end())
2430 : {
2431 1692 : volume_product.push_back(0.0 /* dummy value */);
2432 1692 : continue;
2433 : }
2434 : }
2435 :
2436 8848 : Real product = 0.0;
2437 966593 : for (const auto & e : c.second)
2438 : {
2439 : // we are only accessing local elements here, so no need to check for nullptr
2440 957745 : const auto * elem = getMooseMesh().queryElemPtr(globalElemID(e));
2441 957745 : auto v = var_num.at(elem->subdomain_id()).first;
2442 957745 : auto dof_idx = elem->dof_number(sys_number, v, 0);
2443 957745 : const auto scale_val = scaling ? scaling->at(elem->subdomain_id()) : 1.0;
2444 957745 : product += _serialized_solution(dof_idx) * elem->volume() / scale_val;
2445 : }
2446 :
2447 8848 : volume_product.push_back(product);
2448 : }
2449 :
2450 : std::map<cellInfo, Real> global_volume_product;
2451 1505 : gatherCellSum(volume_product, global_volume_product);
2452 :
2453 1505 : return global_volume_product;
2454 1505 : }
2455 :
2456 : void
2457 2695 : OpenMCCellAverageProblem::sendTemperatureToOpenMC() const
2458 : {
2459 2695 : if (!_specified_temperature_feedback)
2460 1594 : return;
2461 :
2462 1101 : _console << "Sending temperature to OpenMC cells... " << std::endl;
2463 :
2464 1101 : double maximum = std::numeric_limits<double>::min();
2465 1101 : double minimum = std::numeric_limits<double>::max();
2466 :
2467 : // collect the volume-temperature product across local ranks
2468 : std::vector<coupling::CouplingFields> phase = {coupling::temperature,
2469 1101 : coupling::density_and_temperature};
2470 : std::map<cellInfo, Real> cell_vol_temp =
2471 1101 : computeVolumeWeightedCellInput(_subdomain_to_temp_vars, &phase);
2472 :
2473 10257 : for (const auto & c : _cell_to_elem)
2474 : {
2475 9160 : auto cell_info = c.first;
2476 9160 : if (!hasTemperatureFeedback(cell_info))
2477 32 : continue;
2478 :
2479 9128 : Real average_temp = cell_vol_temp.at(cell_info) / _cell_to_elem_volume.at(cell_info);
2480 :
2481 9128 : minimum = std::min(minimum, average_temp);
2482 9128 : maximum = std::max(maximum, average_temp);
2483 :
2484 9128 : if (_verbose)
2485 12824 : _console << "Setting cell " << printCell(cell_info) << " ["
2486 6412 : << _cell_to_n_contained.at(cell_info)
2487 6412 : << " contained cells] to temperature (K): " << std::setw(4) << average_temp
2488 6412 : << std::endl;
2489 :
2490 9128 : const bool identical_fill = cellHasIdenticalFill(cell_info);
2491 9128 : const auto & unshifted_contained_cells = unshiftedContainedCells(cell_info);
2492 192164 : for (auto & [cc_idx, cc_instances] : unshifted_contained_cells)
2493 : {
2494 11576908 : for (unsigned int cc_instance_idx = 0; cc_instance_idx < cc_instances.size();
2495 : ++cc_instance_idx)
2496 : {
2497 : // Shift the cell instances in-place if required for the identical cell fill optimization.
2498 : auto cc_instance = identical_fill
2499 11393872 : ? containedCellInstanceShift(cell_info, cc_idx, cc_instance_idx)
2500 11393872 : : cc_instances[cc_instance_idx];
2501 11393872 : setCellTemperature(cc_idx, cc_instance, average_temp, cell_info);
2502 : }
2503 : }
2504 : }
2505 :
2506 1097 : if (!_verbose)
2507 122 : _console << " Sent cell-averaged min/max (K): " << minimum << ", " << maximum << std::endl;
2508 1097 : }
2509 :
2510 : OpenMCCellAverageProblem::cellInfo
2511 6718976 : OpenMCCellAverageProblem::firstContainedMaterialCell(const cellInfo & cell_info) const
2512 : {
2513 : // this function is only used for displaying temperature and density in auxkernels; to avoid
2514 : // confusing the user, we return the first cell which is non-void fill - even if we set
2515 : // the density on a void cell via Cardinal, OpenMC will always be storing a density of zero
2516 : // because the density multiplier is zero. This could be confusing when reporting the density
2517 : // in CellDensityAux of a cell containing multiple nested cells and if void happens to be the
2518 : // first of those contained cells. So, we screen it out here.
2519 :
2520 6718976 : const bool identical_fill = cellHasIdenticalFill(cell_info);
2521 6718976 : const auto & unshifted_contained_cells = unshiftedContainedCells(cell_info);
2522 6740304 : for (auto & [cc_idx, cc_instances] : unshifted_contained_cells)
2523 : {
2524 6729216 : const auto & cell = openmc::model::cells[cc_idx];
2525 6842704 : for (unsigned int cc_instance_idx = 0; cc_instance_idx < cc_instances.size(); ++cc_instance_idx)
2526 : {
2527 : // Shift the cell instances in-place if required for the identical cell fill optimization.
2528 : auto cc_instance = identical_fill
2529 6821376 : ? containedCellInstanceShift(cell_info, cc_idx, cc_instance_idx)
2530 : : cc_instances[cc_instance_idx];
2531 :
2532 : const auto mat_index = cell->material(cc_instance);
2533 6821376 : if (mat_index != openmc::MATERIAL_VOID)
2534 : {
2535 : cellInfo first_cell = {cc_idx, cc_instance};
2536 6707888 : return first_cell;
2537 : }
2538 : }
2539 : }
2540 :
2541 : // if the cell only contains void, then we'll return that
2542 11088 : auto cc_idx = unshifted_contained_cells.begin()->first;
2543 : const auto & cc_instances = unshifted_contained_cells.begin()->second;
2544 : // Shift the cell instance in-place if required for the identical cell fill optimization.
2545 : auto cc_instance =
2546 11088 : identical_fill ? containedCellInstanceShift(cell_info, cc_idx, 0) : cc_instances[0];
2547 : cellInfo first_cell = {cc_idx, cc_instance};
2548 11088 : return first_cell;
2549 : }
2550 :
2551 : void
2552 2691 : OpenMCCellAverageProblem::sendDensityToOpenMC() const
2553 : {
2554 2691 : if (!_specified_density_feedback)
2555 2287 : return;
2556 :
2557 404 : _console << "Sending density to OpenMC cells... " << std::endl;
2558 :
2559 404 : double maximum = std::numeric_limits<double>::min();
2560 404 : double minimum = std::numeric_limits<double>::max();
2561 :
2562 : // collect the volume-density product across local ranks
2563 : std::vector<coupling::CouplingFields> phase = {coupling::density,
2564 404 : coupling::density_and_temperature};
2565 404 : const auto scaling = openmc::settings::run_CE ? nullptr : &_subdomain_to_ref_density;
2566 : std::map<cellInfo, Real> cell_vol_density =
2567 404 : computeVolumeWeightedCellInput(_subdomain_to_density_vars, &phase, scaling);
2568 :
2569 4726 : for (const auto & c : _cell_to_elem)
2570 : {
2571 4324 : auto cell_info = c.first;
2572 4324 : if (!hasDensityFeedback(cell_info))
2573 2096 : continue;
2574 :
2575 2228 : Real average_density = cell_vol_density.at(cell_info) / _cell_to_elem_volume.at(cell_info);
2576 :
2577 2228 : minimum = std::min(minimum, average_density);
2578 2228 : maximum = std::max(maximum, average_density);
2579 :
2580 2228 : if (_verbose)
2581 : {
2582 1650 : if (openmc::settings::run_CE)
2583 3236 : _console << "Setting cell " << printCell(cell_info)
2584 1618 : << " to density (kg/m3): " << std::setw(4) << average_density << std::endl;
2585 : else
2586 64 : _console << "Setting cell " << printCell(cell_info)
2587 32 : << " to MGXS density (-): " << std::setw(4) << average_density << std::endl;
2588 : }
2589 :
2590 2228 : auto & contained_cells = _cell_to_contained_material_cells.at(cell_info);
2591 :
2592 4552 : for (const auto & contained : contained_cells)
2593 4974 : for (const auto & instance : contained.second)
2594 2650 : setCellDensity(contained.first, instance, average_density, cell_info);
2595 : }
2596 :
2597 402 : if (!_verbose)
2598 : {
2599 50 : if (openmc::settings::run_CE)
2600 34 : _console << " Sent cell-averaged min/max (kg/m3): " << minimum << ", " << maximum
2601 34 : << std::endl;
2602 : else
2603 16 : _console << " Sent cell-averaged min/max (-): " << minimum << ", " << maximum << std::endl;
2604 : }
2605 402 : }
2606 :
2607 : Real
2608 479564 : OpenMCCellAverageProblem::tallyMultiplier(const std::string & score_name,
2609 : const Real & local_mean_tally) const
2610 : {
2611 479564 : if (!isHeatingScore(score_name))
2612 : {
2613 : // we need to get an effective source rate (particles / second) in order to
2614 : // normalize the tally
2615 192588 : Real source = local_mean_tally;
2616 192588 : if (_run_mode == openmc::RunMode::EIGENVALUE)
2617 178736 : source *= *_power / EV_TO_JOULE / _source_rate_norm_tally->getMean(_source_rate_score);
2618 : else
2619 13852 : source *= *_source_strength;
2620 :
2621 : // - Reaction rate scores have units of reactions/src (OpenMC) or reactions/s (Cardinal).
2622 : // - 'inverse-velocity' has units of particles*s/src (OpenMC) or particles (Cardinal).
2623 : // This score is flux-weighted, and must be divided by the flux to recover the true
2624 : // inverse velocity, which has units of s/cm.
2625 : // - 'decay-rate' has units of reactions/src/s (OpenMC) or reactions/s^2 (Cardinal).
2626 : // This score is weighted by the delayed fission rate, and must be divided by
2627 : // `delayed-nu-fission` to obtain the true decay rate, which has units of 1/s.
2628 : // - 'damage-energy' has units of eV/src (OpenMC) or eV/s (Cardinal). While the units of
2629 : // damage-energy are the same as a heating tally, we don't normalize it like one as it's
2630 : // used as an intermediate to compute DPA.
2631 369504 : if (isReactionRateScore(score_name) || score_name == "inverse-velocity" ||
2632 368392 : score_name == "decay-rate" || score_name == "damage-energy")
2633 : return source;
2634 :
2635 175268 : if (score_name == "flux")
2636 175268 : return source / _scaling;
2637 : else
2638 0 : mooseError("Unhandled tally score enum!");
2639 : }
2640 : else
2641 : {
2642 : // Heating tallies have units of eV / source particle
2643 286976 : if (_run_mode == openmc::RunMode::EIGENVALUE)
2644 286896 : return *_power;
2645 : else
2646 80 : return *_source_strength * EV_TO_JOULE * local_mean_tally;
2647 : }
2648 : }
2649 :
2650 : const Real
2651 1333120 : OpenMCCellAverageProblem::getReferenceDensity(const Elem * elem) const
2652 : {
2653 : // The element should never be null entering this function.
2654 : assert(elem != nullptr);
2655 :
2656 1333120 : return openmc::settings::run_CE ? 1.0 : _subdomain_to_ref_density.at(elem->subdomain_id());
2657 : }
2658 :
2659 : void
2660 32 : OpenMCCellAverageProblem::dufekGudowskiParticleUpdate()
2661 : {
2662 32 : int64_t n = (_n_particles_1 + std::sqrt(_n_particles_1 * _n_particles_1 +
2663 32 : 4.0 * _n_particles_1 * _total_n_particles)) /
2664 32 : 2.0;
2665 32 : openmc::settings::n_particles = n;
2666 32 : }
2667 :
2668 : void
2669 5526 : OpenMCCellAverageProblem::syncSolutions(ExternalProblem::Direction direction)
2670 : {
2671 5526 : OpenMCProblemBase::syncSolutions(direction);
2672 :
2673 : // We can skip syncronizing the solution when running with adaptivity
2674 : // and the mesh hasn't changed. This only applies to steady-state calculations
2675 : // as the mesh is adapted once per timestep in a transient calculation.
2676 5526 : if (_has_adaptivity && !_run_on_adaptivity_cycle)
2677 : return;
2678 :
2679 5490 : _aux->serializeSolution();
2680 :
2681 5490 : switch (direction)
2682 : {
2683 2759 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
2684 : {
2685 : // update the [Mesh] internally, so that if we have the skinner we then propagate those
2686 : // changes to the OpenMC geometry
2687 2759 : if (_use_displaced)
2688 : {
2689 99 : _console << "Updating the displaced mesh..." << std::endl;
2690 99 : _displaced_problem->updateMesh();
2691 : }
2692 :
2693 : // Reinitialize the MOOSE -> OpenMC coupling.
2694 2759 : reinitCouplingAndApplyFeedback();
2695 :
2696 2743 : break;
2697 : }
2698 2731 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
2699 : {
2700 2731 : _console << "Extracting OpenMC tallies..." << std::endl;
2701 :
2702 2731 : if (_local_tallies.size() == 0)
2703 : break;
2704 :
2705 : // Loop over all of the tallies and calculate their sums and averages.
2706 5809 : for (auto & local_tally : _local_tallies)
2707 3534 : local_tally->computeSumAndMean();
2708 :
2709 : // Recompute sums and means for tallies that are linked to other tallies.
2710 : // This is used to perform local normalization for translated copies of mesh tallies.
2711 : // These loops must be separate due to data dependencies.
2712 5809 : for (auto & local_tally : _local_tallies)
2713 3534 : local_tally->gatherLinkedSum();
2714 5809 : for (auto & local_tally : _local_tallies)
2715 3534 : local_tally->renormalizeLinkedTallies();
2716 :
2717 : // Loop over the tallies to relax and normalize their results score by score. Then, store the
2718 : // results.
2719 5803 : for (unsigned int i = 0; i < _local_tallies.size(); ++i)
2720 : {
2721 3532 : _local_tallies[i]->relaxAndNormalizeTally();
2722 :
2723 7572 : for (unsigned int score = 0; score < _local_tallies[i]->getScores().size(); ++score)
2724 : {
2725 : // Store the tally results.
2726 8088 : _local_tallies[i]->storeResults(_tally_var_ids[i], score, "relaxed");
2727 :
2728 : // Store additional tally outputs.
2729 4044 : if (_local_tallies[i]->hasOutputs())
2730 : {
2731 : const auto & outs = _local_tallies[i]->getOutputs();
2732 1172 : for (unsigned int j = 0; j < outs.size(); ++j)
2733 658 : _local_tallies[i]->storeResults(_tally_ext_var_ids[i][j], score, outs[j]);
2734 : }
2735 : }
2736 : }
2737 :
2738 : break;
2739 : }
2740 0 : default:
2741 0 : mooseError("Unhandled Direction enum in OpenMCCellAverageProblem!");
2742 : }
2743 :
2744 5470 : _first_transfer = false;
2745 5470 : _aux->solution().close();
2746 5470 : _aux->system().update();
2747 : }
2748 :
2749 : void
2750 3571 : OpenMCCellAverageProblem::reinitCouplingAndApplyFeedback()
2751 : {
2752 : #ifdef ENABLE_DAGMC
2753 1826 : if (_skinner)
2754 : {
2755 : // Update the OpenMC geometry to take into account skinning. This also calls
2756 : // _skinner->update().
2757 53 : updateOpenMCGeometry();
2758 :
2759 : // regenerate the DAGMC geometry
2760 49 : reloadDAGMC();
2761 : }
2762 : #endif
2763 :
2764 3567 : if (_need_to_reinit_coupling)
2765 : {
2766 921 : if (_volume_calc)
2767 7 : _volume_calc->resetVolumeCalculation();
2768 :
2769 921 : resetTallies();
2770 921 : setupProblem();
2771 : }
2772 :
2773 : // Change nuclide composition of material; we put this here so that we can still then change
2774 : // the _overall_ density (like due to thermal expansion, which does not change the relative
2775 : // amounts of the different nuclides)
2776 3567 : sendNuclideDensitiesToOpenMC();
2777 :
2778 3563 : if (_first_transfer && (_specified_temperature_feedback || _specified_density_feedback))
2779 : {
2780 : std::string incoming_transfer =
2781 2387 : _specified_density_feedback ? "temperature and density" : "temperature";
2782 :
2783 1457 : switch (_initial_condition)
2784 : {
2785 2 : case coupling::hdf5:
2786 : {
2787 : // if we're reading temperature and density from an existing HDF5 file,
2788 : // we don't need to send anything in to OpenMC, so we can leave.
2789 2 : importProperties();
2790 0 : _console << "Skipping " << incoming_transfer
2791 0 : << " transfer into OpenMC because 'initial_properties = hdf5'" << std::endl;
2792 0 : return;
2793 : }
2794 : case coupling::moose:
2795 : {
2796 : // transfer will happen from MOOSE - proceed normally
2797 : break;
2798 : }
2799 866 : case coupling::xml:
2800 : {
2801 : // if we're just using whatever temperature and density are already in the XML
2802 : // files, we don't need to send anything in to OpenMC, so we can leave.
2803 866 : _console << "Skipping " << incoming_transfer
2804 866 : << " transfer into OpenMC because 'initial_properties = xml'" << std::endl;
2805 866 : return;
2806 : }
2807 0 : default:
2808 0 : mooseError("Unhandled OpenMCInitialConditionEnum!");
2809 : }
2810 : }
2811 :
2812 : // Because we require at least one of fluid_blocks and solid_blocks, we are guaranteed
2813 : // to be setting the temperature of all of the cells in cell_to_elem - only for the density
2814 : // transfer do we need to filter for the fluid cells
2815 2695 : sendTemperatureToOpenMC();
2816 :
2817 2691 : sendDensityToOpenMC();
2818 :
2819 2689 : if (_export_properties)
2820 0 : openmc_properties_export("properties.h5");
2821 :
2822 : // After setting cell temperatures, we need to re-initialize MGXS data as temperature
2823 : // interpolation is performed on initialization. Verbosity is temporarily modified here
2824 : // as the user has seen the MGXS initialization info previously.
2825 2689 : if (!openmc::settings::run_CE)
2826 : {
2827 159 : auto initial_verbosity = openmc::settings::verbosity;
2828 159 : openmc::settings::verbosity = 1;
2829 : // Clear the MGXS manager.
2830 159 : openmc::data::mg = {};
2831 : // Reload the MGXS data.
2832 159 : openmc::data::mg.read_header(openmc::settings::path_cross_sections);
2833 159 : openmc::put_mgxs_header_data_to_globals();
2834 159 : openmc::finalize_cross_sections();
2835 159 : openmc::settings::verbosity = initial_verbosity;
2836 : }
2837 : }
2838 :
2839 : void
2840 812 : OpenMCCellAverageProblem::critSearchStep()
2841 : {
2842 812 : _aux->serializeSolution();
2843 :
2844 : // Reinitialize the OpenMC coupling prior to the execution of
2845 : // a criticality search step.
2846 812 : reinitCouplingAndApplyFeedback();
2847 :
2848 812 : _aux->solution().close();
2849 812 : _aux->system().update();
2850 812 : }
2851 :
2852 : void
2853 2239 : OpenMCCellAverageProblem::createQRules(QuadratureType type,
2854 : Order order,
2855 : Order volume_order,
2856 : Order face_order,
2857 : SubdomainID block,
2858 : const bool allow_negative_qweights)
2859 : {
2860 : // start copy: Copied from base class's createQRules in order to retain the same default behavior
2861 2239 : if (order == INVALID_ORDER)
2862 : {
2863 2239 : order = getNonlinearSystemBase(0).getMinQuadratureOrder();
2864 2239 : if (order < getAuxiliarySystem().getMinQuadratureOrder())
2865 2018 : order = getAuxiliarySystem().getMinQuadratureOrder();
2866 : }
2867 :
2868 2239 : if (volume_order == INVALID_ORDER)
2869 2239 : volume_order = order;
2870 :
2871 2239 : if (face_order == INVALID_ORDER)
2872 : face_order = order;
2873 : // end copy
2874 :
2875 : // The approximations made in elem->volume() are only valid for Gauss and Monomial quadratures
2876 : // if they are second order or above
2877 4478 : if (type == Moose::stringToEnum<QuadratureType>("GAUSS"))
2878 4478 : setMinimumVolumeQRules(volume_order, "GAUSS");
2879 4478 : if (type == Moose::stringToEnum<QuadratureType>("MONOMIAL"))
2880 0 : setMinimumVolumeQRules(volume_order, "MONOMIAL");
2881 4478 : if (type == Moose::stringToEnum<QuadratureType>("GAUSS_LOBATTO"))
2882 0 : setMinimumVolumeQRules(volume_order, "GAUSS_LOBATTO");
2883 :
2884 : // Some quadrature rules don't ever seem to give a matching elem->volume() with the MOOSE
2885 : // volume integrations
2886 6717 : if (type == Moose::stringToEnum<QuadratureType>("GRID") ||
2887 6717 : type == Moose::stringToEnum<QuadratureType>("TRAP"))
2888 0 : mooseError(
2889 : "The ",
2890 0 : std::to_string(type),
2891 : " quadrature set will never match the '_current_elem_volume' used to compute\n"
2892 : "integrals in MOOSE. This means that the tally computed by OpenMC is normalized by\n"
2893 : "a different volume than used for MOOSE volume integrations, such that the specified "
2894 : "'power' or 'source_strength'\n"
2895 : "would not be respected. Please switch to a different quadrature set.");
2896 :
2897 2239 : FEProblemBase::createQRules(
2898 : type, order, volume_order, face_order, block, allow_negative_qweights);
2899 2239 : }
2900 :
2901 : void
2902 2239 : OpenMCCellAverageProblem::setMinimumVolumeQRules(Order & volume_order,
2903 : const std::string & /* type */)
2904 : {
2905 4478 : if (volume_order < Moose::stringToEnum<Order>("SECOND"))
2906 2236 : volume_order = SECOND;
2907 2239 : }
2908 :
2909 : double
2910 188144 : OpenMCCellAverageProblem::cellMappedVolume(const cellInfo & cell_info) const
2911 : {
2912 188144 : return _cell_to_elem_volume.at(cell_info);
2913 : }
2914 :
2915 : double
2916 1333120 : OpenMCCellAverageProblem::cellDensity(const cellInfo & cell_info, const Elem * elem) const
2917 : {
2918 1333120 : auto material_cell = firstContainedMaterialCell(cell_info);
2919 :
2920 : double density;
2921 1333120 : int err = openmc_cell_get_density(material_cell.first, &material_cell.second, &density);
2922 1333120 : catchOpenMCError(err, "get density of cell " + printCell(cell_info));
2923 :
2924 : // Rescale by the reference density, if required.
2925 1333120 : const auto ref_den = getReferenceDensity(elem);
2926 :
2927 1333120 : return ref_den * density / densityConversionFactor();
2928 : }
2929 :
2930 : double
2931 4086016 : OpenMCCellAverageProblem::cellTemperature(const cellInfo & cell_info) const
2932 : {
2933 4086016 : auto material_cell = firstContainedMaterialCell(cell_info);
2934 :
2935 : double T;
2936 4086016 : int err = openmc_cell_get_temperature(material_cell.first, &material_cell.second, &T);
2937 4086016 : catchOpenMCError(err, "get temperature of cell " + printCell(cell_info));
2938 4086016 : return T;
2939 : }
2940 :
2941 : void
2942 49 : OpenMCCellAverageProblem::reloadDAGMC()
2943 : {
2944 : #ifdef ENABLE_DAGMC
2945 98 : _dagmc.reset(new moab::DagMC(_skinner->moabPtr(),
2946 : 0.0 /* overlap tolerance, default */,
2947 : 0.001 /* numerical precision, default */,
2948 147 : 0 /* verbosity */));
2949 :
2950 : // Set up geometry in DagMC from already-loaded mesh
2951 49 : _dagmc->load_existing_contents();
2952 :
2953 : // Initialize acceleration data structures
2954 49 : _dagmc->init_OBBTree();
2955 :
2956 : // Get an iterator to the DAGMC universe unique ptr
2957 : auto univ_it =
2958 49 : openmc::model::universes.begin() + openmc::model::universe_map.at(_dagmc_universe_id);
2959 :
2960 : // Remove the old universe
2961 : openmc::model::universes.erase(univ_it);
2962 :
2963 : // Create new DAGMC universe
2964 49 : openmc::model::universes.emplace_back(std::make_unique<openmc::DAGUniverse>(_dagmc, "", true));
2965 49 : _dagmc_universe_id = openmc::model::universes.back()->id_;
2966 :
2967 : openmc::model::universe_map.clear();
2968 103 : for (int32_t i = 0; i < openmc::model::universes.size(); ++i)
2969 54 : openmc::model::universe_map[openmc::model::universes[i]->id_] = i;
2970 :
2971 49 : if (!_dagmc_root_universe)
2972 5 : openmc::model::cells[openmc::model::cell_map.at(_cell_using_dagmc_universe_id)]->fill_ =
2973 5 : _dagmc_universe_id;
2974 :
2975 49 : _console << "Re-generating OpenMC model with " << openmc::model::cells.size() << " cells... "
2976 49 : << std::endl;
2977 :
2978 : // Clear cells on all surviving universes.
2979 103 : for (auto & universe : openmc::model::universes)
2980 54 : universe->cells_.clear();
2981 :
2982 : // Add cells to universes
2983 49 : openmc::populate_universes();
2984 :
2985 : // Set the root universe
2986 49 : openmc::model::root_universe = openmc::find_root_universe();
2987 49 : openmc::check_dagmc_root_univ();
2988 :
2989 : // Final geometry setup
2990 49 : openmc::finalize_geometry();
2991 :
2992 : // Finalize cross sections; we manually change the verbosity here because if skinning is
2993 : // enabled, we don't want to overwhelm the user with excess console output showing info
2994 : // which ultimately is no different from that shown on initialization
2995 49 : auto initial_verbosity = openmc::settings::verbosity;
2996 49 : openmc::settings::verbosity = 1;
2997 49 : openmc::finalize_cross_sections();
2998 :
2999 : // Finalize DAGMC cell densities after setting up the new geometry. CSG cells (and
3000 : // eventually non-skinned DAGMC cells) already have their densities finalized.
3001 511 : for (auto & c : openmc::model::cells)
3002 : {
3003 462 : if (c->geom_type() == openmc::GeometryType::CSG)
3004 5 : continue;
3005 :
3006 457 : c->density_mult_ = {1.0};
3007 : }
3008 :
3009 : // Needed to obtain correct cell instances
3010 49 : openmc::prepare_distribcell();
3011 49 : openmc::settings::verbosity = initial_verbosity;
3012 : #endif
3013 49 : }
3014 :
3015 : void
3016 632 : OpenMCCellAverageProblem::addFilter(const std::string & type,
3017 : const std::string & name,
3018 : InputParameters & moose_object_pars)
3019 : {
3020 1218 : auto filter = addObject<FilterBase>(type, name, moose_object_pars, false)[0];
3021 586 : _filters[name] = filter;
3022 586 : }
3023 :
3024 : std::shared_ptr<TallyBase>
3025 2847 : OpenMCCellAverageProblem::addTally(const std::string & type,
3026 : const std::string & name,
3027 : InputParameters & moose_object_pars)
3028 : {
3029 5643 : auto tally = addObject<TallyBase>(type, name, moose_object_pars, false)[0];
3030 2796 : _local_tallies.push_back(tally);
3031 :
3032 : // Set the relaxation scheme.
3033 8387 : tally->setRelaxation(_relaxation, getParam<Real>("relaxation_factor"));
3034 :
3035 : const auto & tally_scores = tally->getScores();
3036 6068 : for (unsigned int i = 0; i < tally_scores.size(); ++i)
3037 : {
3038 : // Populate a map which counts the number of times a score is referenced by local tallies.
3039 : // Used for error checking.
3040 : if (_score_count.count(tally_scores[i]) == 0)
3041 2805 : _score_count[tally_scores[i]] = 1;
3042 : else
3043 468 : _score_count[tally_scores[i]]++;
3044 :
3045 : // Add the local tally's score to the list of scores if we don't have it yet.
3046 3273 : if (std::find(_all_tally_scores.begin(), _all_tally_scores.end(), tally_scores[i]) ==
3047 : _all_tally_scores.end())
3048 2805 : _all_tally_scores.push_back(tally_scores[i]);
3049 : }
3050 :
3051 2795 : return tally;
3052 : }
3053 :
3054 : void
3055 2255 : OpenMCCellAverageProblem::validateLocalTallies()
3056 : {
3057 : // We can skip this check if we don't have tallies.
3058 2255 : if (_local_tallies.size() == 0)
3059 : return;
3060 :
3061 : // Make sure we can assume that tallies can be separate.
3062 1752 : if (_assume_separate_tallies)
3063 : {
3064 24 : for (const auto & tally : _local_tallies)
3065 14 : if (tally->addingGlobalTally())
3066 2 : paramError("assume_separate_tallies",
3067 : "Cannot assume separate tallies when either of 'check_tally_sum' or"
3068 : "'normalize_by_global_tally' is true!");
3069 :
3070 10 : if (_local_tallies.size() > 1)
3071 2 : paramError("assume_separate_tallies",
3072 : "Cannot assume separate tallies when there are multiple tallies added in the "
3073 : "[Tallies] block!");
3074 : }
3075 :
3076 : // need some special treatment for non-heating scores, in eigenvalue mode
3077 : bool has_non_heating_score = false;
3078 4539 : for (const auto & t : _all_tally_scores)
3079 2791 : if (!isHeatingScore(t))
3080 : has_non_heating_score = true;
3081 :
3082 1748 : if (has_non_heating_score && _run_mode == openmc::RunMode::EIGENVALUE)
3083 : {
3084 : std::string non_heating_scores;
3085 1947 : for (const auto & e : _all_tally_scores)
3086 : {
3087 1416 : if (!isHeatingScore(e))
3088 : {
3089 895 : std::string l = e;
3090 : std::replace(l.begin(), l.end(), '-', '_');
3091 1790 : non_heating_scores += "" + l + ", ";
3092 : }
3093 : }
3094 :
3095 531 : if (non_heating_scores.length() > 0)
3096 531 : non_heating_scores.erase(non_heating_scores.length() - 2);
3097 :
3098 1060 : checkRequiredParam(_pars,
3099 : "source_rate_normalization",
3100 531 : "using a non-heating tally (" + non_heating_scores + ") in eigenvalue mode");
3101 529 : const auto & norm = getParam<MooseEnum>("source_rate_normalization");
3102 1058 : std::string n = enumToTallyScore(norm);
3103 :
3104 529 : if (_local_tallies.size() > 1)
3105 : {
3106 : if (_score_count.count(n) == 0)
3107 2 : mooseError("The local tallies added in the [Tallies] block do not contain the requested "
3108 2 : "heating score " +
3109 0 : n +
3110 : ". You must either add this score in one of the tallies or choose a different "
3111 : "heating score.");
3112 :
3113 443 : if (_score_count.at(n) > 1)
3114 : {
3115 : // Edge case: multiple scores from linked MeshTally objects.
3116 36 : unsigned int linked = 0;
3117 : unsigned int num_with_score = 0;
3118 156 : for (auto tally : _local_tallies)
3119 : {
3120 120 : if (tally->hasScore(n))
3121 : {
3122 80 : linked = std::max(linked, static_cast<unsigned int>(tally->linkedTallies().size()) + 1);
3123 80 : num_with_score++;
3124 : }
3125 : }
3126 :
3127 : // Can only allow auto-detection of the normalization tally if there is a single linkage
3128 : // of every mesh tally with the normalization score.
3129 36 : if (_score_count.at(n) != linked || _score_count.at(n) != num_with_score)
3130 : {
3131 : // If there are more then one value of 'source_rate_normalization', the user needs
3132 : // to tell us which tally to use.
3133 54 : checkRequiredParam(
3134 : _pars,
3135 : "normalization_tally",
3136 28 : "using a non-heating tally (" + non_heating_scores +
3137 : ") in eigenvalue mode and adding more then one tally in the [Tallies] block");
3138 78 : const auto norm_tally_name = getParam<std::string>("normalization_tally");
3139 :
3140 : // Check to make sure the user provided a tally name for eigenvalue normalization
3141 : // that's been added.
3142 94 : for (auto tally : _local_tallies)
3143 68 : if (norm_tally_name == tally->name())
3144 : _source_rate_norm_tally = tally;
3145 :
3146 26 : if (!_source_rate_norm_tally)
3147 0 : paramError("normalization_tally",
3148 0 : "The tally " + norm_tally_name +
3149 : " does not exist in the problem! Please specify a tally added in the "
3150 : "[Tallies] block!");
3151 : }
3152 : else
3153 : {
3154 56 : for (auto tally : _local_tallies)
3155 48 : if (tally->hasScore(n))
3156 : _source_rate_norm_tally = tally;
3157 : }
3158 : }
3159 : else
3160 : {
3161 : // Otherwise, we can check the tallies added and find the one scoring the requested
3162 : // value of 'source_rate_normalization'.
3163 1485 : for (auto tally : _local_tallies)
3164 1078 : if (tally->hasScore(n))
3165 : _source_rate_norm_tally = tally;
3166 : }
3167 : }
3168 : else
3169 : _source_rate_norm_tally = _local_tallies[0];
3170 :
3171 : // If it's not in the specified source rate tally, we can add it for the user.
3172 525 : if (!_source_rate_norm_tally->hasScore(n))
3173 : {
3174 18 : if (_source_rate_norm_tally->renamesTallyVars())
3175 2 : mooseError("When specifying 'name', the score indicated in "
3176 : "'source_rate_normalization' must be\n"
3177 : "listed in 'score' so that we know what you want to name that score (",
3178 : norm,
3179 : ")");
3180 :
3181 : // We can add the requested normalization score if and only if a single tally was added by
3182 : // [Tallies].
3183 16 : _all_tally_scores.push_back(n);
3184 16 : _source_rate_norm_tally->addScore(n);
3185 16 : _source_rate_score = _source_rate_norm_tally->scoreIndex(n);
3186 : }
3187 : else
3188 507 : _source_rate_score = _source_rate_norm_tally->scoreIndex(n);
3189 523 : }
3190 2434 : else if (isParamValid("source_rate_normalization"))
3191 24 : mooseWarning(
3192 : "When either running in fixed-source mode, or all tallies have units of eV/src, the "
3193 : "'source_rate_normalization' parameter is unused!");
3194 : }
3195 :
3196 : void
3197 53 : OpenMCCellAverageProblem::updateOpenMCGeometry()
3198 : {
3199 : #ifdef ENABLE_DAGMC
3200 : // Need to swap array indices back to ids as OpenMC swapped these when preparing geometry.
3201 390 : for (const auto & cell : openmc::model::cells)
3202 : {
3203 337 : if (cell->type_ == openmc::Fill::MATERIAL)
3204 : {
3205 : std::vector<int32_t> mat_ids;
3206 664 : for (const auto & mat_index : cell->material_)
3207 332 : mat_ids.push_back(mat_index == openmc::MATERIAL_VOID
3208 : ? openmc::MATERIAL_VOID
3209 258 : : openmc::model::materials[mat_index]->id_);
3210 332 : cell->material_ = mat_ids;
3211 332 : }
3212 337 : if (cell->type_ == openmc::Fill::UNIVERSE && cell->fill_ != openmc::C_NONE)
3213 5 : cell->fill_ = openmc::model::universes[cell->fill_]->id_;
3214 337 : if (cell->type_ == openmc::Fill::LATTICE && cell->fill_ != openmc::C_NONE)
3215 0 : cell->fill_ = openmc::model::lattices[cell->fill_]->id_;
3216 :
3217 337 : cell->universe_ = openmc::model::universes[cell->universe_]->id_;
3218 : }
3219 :
3220 53 : for (const auto & lattice : openmc::model::lattices)
3221 : {
3222 0 : for (openmc::LatticeIter it = lattice->begin(); it != lattice->end(); ++it)
3223 : {
3224 0 : int u_index = *it;
3225 0 : *it = openmc::model::universes[u_index]->id_;
3226 : }
3227 :
3228 0 : if (lattice->outer_ != openmc::NO_OUTER_UNIVERSE)
3229 0 : lattice->outer_ = openmc::model::universes[lattice->outer_]->id_;
3230 : }
3231 :
3232 : // skin the mesh geometry according to contours in temperature, density, and subdomain
3233 53 : _skinner->update();
3234 :
3235 : openmc::model::universe_level_counts.clear();
3236 :
3237 : // Clear nuclides and elements, these will get reset in read_ce_cross_sections
3238 : // Horrible circular logic means that clearing nuclides clears nuclide_map, but
3239 : // which is needed before nuclides gets reset (similar for elements)
3240 : std::unordered_map<std::string, int> nuclide_map_copy = openmc::data::nuclide_map;
3241 49 : openmc::data::nuclides.clear();
3242 : openmc::data::nuclide_map = nuclide_map_copy;
3243 :
3244 : std::unordered_map<std::string, int> element_map_copy = openmc::data::element_map;
3245 49 : openmc::data::elements.clear();
3246 : openmc::data::element_map = element_map_copy;
3247 :
3248 : // Clear existing DAGMC cell data. Cells cannot be deleted in-place as that invalidates
3249 : // all pointers and iterators, so we loop over the cell map to store a list of DAGMC cells.
3250 : // Afterwards, the cells contained in the list can be deleted.
3251 : std::vector<int32_t> cells_to_delete;
3252 374 : for (auto [id, index] : openmc::model::cell_map)
3253 325 : if (openmc::model::cells[index]->geom_type() == openmc::GeometryType::DAG)
3254 320 : cells_to_delete.push_back(openmc::model::cells[index]->id_);
3255 :
3256 369 : for (auto cell : cells_to_delete)
3257 : {
3258 1850 : for (int32_t i = 0; i < openmc::model::cells.size(); ++i)
3259 : {
3260 1850 : if (openmc::model::cells[i]->id_ == cell)
3261 : {
3262 : openmc::model::cells.erase(openmc::model::cells.begin() + i);
3263 320 : break;
3264 : }
3265 : }
3266 : }
3267 49 : cells_to_delete.clear();
3268 :
3269 : // Clear existing surface data. Similar to cells, deletion of the DAGMC surfaces must be
3270 : // deferred.
3271 : std::vector<int> surfaces_to_delete;
3272 901 : for (auto [id, index] : openmc::model::surface_map)
3273 852 : if (openmc::model::surfaces[index]->geom_type() == openmc::GeometryType::DAG)
3274 822 : surfaces_to_delete.push_back(openmc::model::surfaces[index]->id_);
3275 :
3276 871 : for (auto surface : surfaces_to_delete)
3277 : {
3278 10289 : for (int i = 0; i < openmc::model::surfaces.size(); ++i)
3279 : {
3280 10289 : if (openmc::model::surfaces[i]->id_ == surface)
3281 : {
3282 : openmc::model::surface_map.erase(surface);
3283 : openmc::model::surfaces.erase(openmc::model::surfaces.begin() + i);
3284 822 : break;
3285 : }
3286 : }
3287 : }
3288 49 : surfaces_to_delete.clear();
3289 :
3290 : // Need to rebuild the cell_map and surface_map since the indices have changed.
3291 : openmc::model::cell_map.clear();
3292 54 : for (int32_t i = 0; i < openmc::model::cells.size(); ++i)
3293 5 : openmc::model::cell_map[openmc::model::cells[i]->id_] = i;
3294 :
3295 : // Horrible hack since we can't undo the surface id -> index swap that happens in
3296 : // CSGCell.region_.expression_, and so the 'surface_map' cannot be rebuilt. Intead, 'surfaces' is
3297 : // resized to the original length and the positions of each surface are shuffled such that they
3298 : // correspond to their indices in the original 'surface_map'. This results in the addition of N
3299 : // extra null 'DAGSurface' objects in 'surfaces', where N is the number of DAGMC surfaces in the
3300 : // geometry. These null surfaces aren't linked to a DAGMC universe and so they do not participate
3301 : // in particle transport, they just take up memory. CSGCell::region_ and Region::expression_ need
3302 : // to be made public in OpenMC to avoid this, or an appropriate series of C-API functions / member
3303 : // functions need to be added to OpenMC.
3304 49 : if (openmc::model::surfaces.size() > 0)
3305 : {
3306 38 : for (int i = openmc::model::surfaces.size(); i < _initial_num_openmc_surfaces; ++i)
3307 33 : openmc::model::surfaces.push_back(
3308 66 : std::move(std::make_unique<openmc::DAGSurface>(nullptr, 0)));
3309 35 : for (const auto & [id, index] : openmc::model::surface_map)
3310 : {
3311 : // If the surface at the index exists and the id is the same, do nothing.
3312 30 : if (openmc::model::surfaces[index]->id_ == id)
3313 30 : continue;
3314 : else
3315 : {
3316 : // Otherwise we need to find the filter and swap it with the filter at the current location.
3317 0 : for (int i = 0; i < openmc::model::surfaces.size(); ++i)
3318 : {
3319 0 : if (openmc::model::surfaces[i]->id_ == id)
3320 : {
3321 : auto temp = std::move(openmc::model::surfaces[index]);
3322 : openmc::model::surfaces[index] = std::move(openmc::model::surfaces[i]);
3323 : openmc::model::surfaces[i] = std::move(temp);
3324 : break;
3325 0 : }
3326 : }
3327 : }
3328 : }
3329 :
3330 : // Sanity check by looping over the surface_map to make sure the indices correspond to the
3331 : // surface ids.
3332 35 : for (const auto & [id, index] : openmc::model::surface_map)
3333 30 : if (openmc::model::surfaces[index]->id_ != id)
3334 0 : mooseError("Internal error: mismatch between surfaces[surface_map[id]]->id_ and id.");
3335 : }
3336 : #endif
3337 98 : }
3338 :
3339 : bool
3340 63612 : OpenMCCellAverageProblem::cellMapsToSubdomain(const cellInfo & cell_info,
3341 : const std::unordered_set<SubdomainID> & id) const
3342 : {
3343 63612 : auto s = _cell_to_elem_subdomain.at(cell_info);
3344 63636 : for (const auto & i : id)
3345 63612 : if (s.find(i) != s.end())
3346 : return true;
3347 :
3348 : return false;
3349 : }
3350 :
3351 : bool
3352 14801236 : OpenMCCellAverageProblem::cellHasIdenticalFill(const cellInfo & cell_info) const
3353 : {
3354 : // material cells are discounted as identical fill
3355 14801236 : const auto & cell = openmc::model::cells[cell_info.first];
3356 14801236 : if (!_has_identical_cell_fills || cell->type_ == openmc::Fill::MATERIAL)
3357 : return false;
3358 :
3359 63612 : return cellMapsToSubdomain(cell_info, _identical_cell_fill_blocks);
3360 : }
3361 :
3362 : int
3363 22905612 : OpenMCCellAverageProblem::containedCellInstanceShift(const cellInfo & cell_info,
3364 : int32_t cc_idx,
3365 : int32_t cc_instance_idx_to_shift) const
3366 : {
3367 22905612 : if (!_has_identical_cell_fills)
3368 0 : mooseError("Internal error: should not call containedCellInstanceShift!");
3369 :
3370 22905612 : auto offset = _n_offset.at(cell_info);
3371 : // All material filled cells in "cell_info".
3372 22905612 : const auto & first_cell_cc = _cell_to_contained_material_cells.at(_first_identical_cell);
3373 : // The cell instance we're shifting.
3374 22905612 : auto instance_to_shift = first_cell_cc.at(cc_idx)[cc_instance_idx_to_shift];
3375 : // The shift to apply.
3376 22905612 : auto shift = _instance_offsets.at(cc_idx)[cc_instance_idx_to_shift];
3377 :
3378 : // Compute the new instance.
3379 22905612 : return instance_to_shift + offset * shift;
3380 : }
3381 :
3382 : const OpenMCCellAverageProblem::containedCells &
3383 8042453 : OpenMCCellAverageProblem::unshiftedContainedCells(const cellInfo & cell_info) const
3384 : {
3385 8042453 : const bool identical_fill = cellHasIdenticalFill(cell_info);
3386 31536 : return identical_fill ? _cell_to_contained_material_cells.at(_first_identical_cell)
3387 8042453 : : _cell_to_contained_material_cells.at(cell_info);
3388 : }
3389 :
3390 : std::vector<int32_t>
3391 7589 : OpenMCCellAverageProblem::materialsInCells(const containedCells & contained_cells) const
3392 : {
3393 : std::vector<int32_t> mats;
3394 22070 : for (const auto & contained : contained_cells)
3395 : {
3396 14481 : if (_cell_material_modifiers.find(contained.first) != _cell_material_modifiers.end())
3397 : {
3398 : // find the iterator corresponding to the ModelModifier corresponding to the contained
3399 : // cell's index
3400 : std::vector<int32_t> modifier_mats =
3401 0 : _cell_material_modifiers.at(contained.first)->getMaterialIndices();
3402 : // insert exactly as many entries from the _material_indices vector as there are cell
3403 : // instances into the current mats vector, which in most cases is the same number
3404 :
3405 : // NOTE: for TRISO problems, when contained corresponds to a TRISO cell it is possible
3406 : // that modifier_mats has many more entries than the actual number of intances of the
3407 : // contained cell here.
3408 :
3409 : // This occurs when using the material_ids_file to do zoning, which assigns a material
3410 : // as many materials as there are instances of this cell in the containing cell.
3411 : // This case is why the below does not go to modifier_mats.end(), i.e. modifier_mats.end()
3412 : // does not always have to equal modifier_mats.begin() + contained.second.size().
3413 0 : mats.insert(
3414 : mats.end(), modifier_mats.begin(), modifier_mats.begin() + contained.second.size());
3415 0 : }
3416 : else
3417 : {
3418 447690 : for (const auto & instance : contained.second)
3419 : {
3420 : // we know this is a material cell, so we don't need to check that the fill is material
3421 : int32_t material_index;
3422 : cellInfo cell_info = {contained.first, instance};
3423 433209 : materialFill(cell_info, material_index);
3424 433209 : mats.push_back(material_index);
3425 : }
3426 : }
3427 : }
3428 :
3429 7589 : return mats;
3430 0 : }
3431 :
3432 : Point
3433 2214617 : OpenMCCellAverageProblem::transformPointToOpenMC(const Point & pt) const
3434 : {
3435 2214617 : Point pnt_out = transformPoint(pt);
3436 :
3437 : // scale point to OpenMC domain
3438 2214617 : pnt_out *= _scaling;
3439 :
3440 2214617 : return pnt_out;
3441 : }
3442 :
3443 : int
3444 1314349 : OpenMCCellAverageProblem::numContainedMaterialCells(const cellInfo & cell_info) const
3445 : {
3446 : int n_contained = 0;
3447 1314349 : const auto & contained_cells = unshiftedContainedCells(cell_info);
3448 2802712 : for (const auto & cell : contained_cells)
3449 1488363 : n_contained += cell.second.size();
3450 :
3451 1314349 : return n_contained;
3452 : }
3453 : #endif
|
