LCOV - code coverage report
Current view: top level - src/base - OpenMCCellAverageProblem.C (source / functions) Hit Total Coverage
Test: neams-th-coe/cardinal: f3518d Lines: 1351 1437 94.0 %
Date: 2025-10-01 10:06:53 Functions: 75 75 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14