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