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