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