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