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