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