Add skinning to OpenMC wrapping.

Makefile

@@ -285,6 +285,9 @@ endif
 
 ifeq ($(ENABLE_DAGMC), yes)
   libmesh_CXXFLAGS    += -DENABLE_DAGMC
+
+  # this flag is used in OpenMC
+  libmesh_CXXFLAGS    += -DDAGMC
 endif
 
 # ======================================================================================

doc/content/source/problems/OpenMCCellAverageProblem.md

@@ -432,7 +432,76 @@ element in the `[Mesh]`.
 
 While this class mainly facilitates data transfers to and from OpenMC, a number of
 other features are implemented in order to enable more convenient input file
-setup and achieve better iterative performance. These are described in this section.
+setup, achieve better iterative performance, and iteratively update the geometry
+(for CAD tracking only). These are described in this section.
+
+#### CAD Geometry Skinning
+  id=skinning
+
+For all cell-based OpenMC models, OpenMC is currently limited to settings a single
+constant temperature, and a single constant density, in each cell. For CSG geometries,
+the user needs to manually set up sub-divisions in the geometry in order to capture
+spatial variation in temperature and density (computed by some other MOOSE application).
+For DAGMC geometries, you can instead optionally re-generate the OpenMC cells after
+each Picard iteration according to contours in temperature and/or density - in other
+words, automatically setting up the OpenMC model so that it can receive spatially-varying
+temperature and density without the user needing to manually subdivide regions of space
+*a priori*.
+
+!alert note
+This skinning feature is only available for DAGMC geometries for which the `[Mesh]`
+obeys the same material boundaries in the starting `.h5m` file.
+
+To use this feature, provide the `skinner` user
+object parameter, of type [MoabSkinner](/userobjects/MoabSkinner.md). After each OpenMC
+run, this object will group elements in the `[Mesh]` into "bins" for temperature, density,
+and/or subdomain (mesh block). A new DAGMC cell will then be created by "skinning" the mesh
+elements according to the boundaries between the unique bin regions.
+For example, suppose the starting DAGMC model consists of two
+cells, each with a different material, as shown in [dagmc_model].
+
+!media dagmc_model.png
+  id=dagmc_model
+  style=width:50%;margin-left:auto;margin-right:auto
+  caption=Example illustration of a DAGMC model before a skinning operation is applied.
+
+Suppose this geometry is skinned with 4 temperature bins between 500 K and 700 K.
+The input syntax would look something like
+
+!listing test/tests/neutronics/dagmc/mesh_tallies/openmc.i
+  start=Problem
+  end=Postprocessors
+
+Elements will then be categorized as falling into one of these bins:
+
+- Bin 0, for $500\ K \leq T< 550\ K$.
+- Bin 1, for $550\ K \leq T< 600\ K$.
+- Bin 2, for $600\ K \leq T< 650\ K$.
+- Bin 3, for $650\ K \leq T\leq 700\ K$.
+
+The geometry will also be skinned according to the subdomain IDs, which align with the
+materials in the original DAGMC geometry. For this problem with two materials, elements
+will be categorized according to falling into one of two subdomain bins:
+
+- Bin 0, for elements in material A
+- Bin 1, for elements in material B
+
+Then, if the following temperature distribution is sent to OpenMC, the DAGMC geometry
+will be re-skinned into 8 unique cells, as shown in [skinning]. Depending on the
+temperature and material definitions, the problem would be skinned into more or less
+cells, depending on which temperatures are observed and whether those element "clumps"
+are contiguous or disjoint from one another.
+A "graveyard" volume is automatically built around the new DAGMC model between two bounding boxes. This
+region is used to apply vacuum boundary conditions by killing any particles that
+enter this region.
+
+!media skinning.png
+  id=skinning
+  caption=Example illustration of a DAGMC model after skinning
+
+Once the geometry has been re-created, the rest of the coupling proceeds as normal -
+the average temperature and density of each new DAGMC cell is set as a volume-average
+of the elements corresponding to each cell.
 
 #### Mesh Scaling
   id=scaling

include/base/OpenMCCellAverageProblem.h

@@ -25,6 +25,11 @@
 #include "openmc/tallies/filter_mesh.h"
 #include "openmc/mesh.h"
 
+#ifdef ENABLE_DAGMC
+#include "MoabSkinner.h"
+#include "DagMC.hpp"
+#endif
+
 /**
  * Mapping of OpenMC to a collection of MOOSE elements, with temperature feedback
  * on solid cells and both temperature and density feedback on fluid cells. The
@@ -307,10 +312,16 @@ class OpenMCCellAverageProblem : public OpenMCProblemBase
    */
   double cellMappedVolume(const cellInfo & cell_info);
 
+  /// Reconstruct the DAGMC geometry after skinning
+  void reloadDAGMC();
+
   /// Constant flag to indicate that a cell/element was unmapped
   static constexpr int32_t UNMAPPED{-1};
 
 protected:
+  /// Calculate the number of unique OpenMC cells (each with a unique ID & instance)
+  void calculateNumCells();
+
   /// Loop over the mapped cells, and build a map between subdomains to OpenMC materials
   void subdomainsToMaterials();
 
@@ -935,7 +946,7 @@ class OpenMCCellAverageProblem : public OpenMCProblemBase
   std::map<cellInfo, std::unordered_set<SubdomainID>> _cell_to_elem_subdomain;
 
   /// Mapping of elem subdomains to materials
-  std::map<SubdomainID, std::unordered_set<int32_t>> _subdomain_to_material;
+  std::map<SubdomainID, std::set<int32_t>> _subdomain_to_material;
 
   /**
    * A point inside the cell, taken simply as the centroid of the first global
@@ -1117,6 +1128,20 @@ class OpenMCCellAverageProblem : public OpenMCProblemBase
   /// Index in tally_score pointing to the score used for normalizing flux tallies in eigenvalue mode
   unsigned int _source_rate_index;
 
+#ifdef ENABLE_DAGMC
+  /// Optional skinner to re-generate the OpenMC geometry on-the-fly for DAGMC models
+  MoabSkinner * _skinner;
+
+  /// Pointer to DAGMC
+  std::shared_ptr<moab::DagMC> _dagmc;
+#endif
+
+  /// Total number of unique OpenMC cell IDs + instances combinations
+  long unsigned int _n_openmc_cells;
+
+  /// Index in the OpenMC universes corresponding to the DAGMC universe
+  int32_t _dagmc_universe_index;
+
   /// Conversion rate from eV to Joule
   static constexpr Real EV_TO_JOULE = 1.6022e-19;
 

include/base/OpenMCProblemBase.h

@@ -326,6 +326,12 @@ class OpenMCProblemBase : public CardinalProblem, public PostprocessorInterface
   /// Total number of unique OpenMC cell IDs + instances combinations
   long unsigned int _n_openmc_cells;
 
+  /**
+   * Whether the OpenMC model consists of a single coordinate level; this can
+   * in some cases be used for more verbose error messages
+   */
+  const bool _single_coord_level;
+
   /**
    * Fixed point iteration index used in relaxation; because we sometimes run OpenMC
    * in a pseudo-transient coupling with NekRS, we simply increment this by 1 each

include/userobjects/MoabSkinner.h

@@ -37,6 +37,8 @@ class MoabSkinner : public GeneralUserObject
    */
   virtual moab::ErrorCode check(const moab::ErrorCode input) const;
 
+  std::string materialName(const unsigned int & block, const unsigned int & density, const unsigned int & temp) const;
+
   /// Perform the skinning operation
   virtual void update();
 
@@ -118,6 +120,31 @@ class MoabSkinner : public GeneralUserObject
                               const unsigned int & density_bin,
                               const unsigned int & subdomain_bin) const;
 
+  /**
+   * Whether the skinner builds a graveyard
+   * @return whether a graveyard is built
+   */
+  virtual const bool & hasGraveyard() const { return _build_graveyard; }
+
+  /**
+   * Set the graveyard setting
+   * @param[in] build whether to build a graveyard
+   */
+  void setGraveyard(bool build) { _build_graveyard = build; }
+
+  /**
+   * Number of density bins; if greater than 1, this means we must be re-generating
+   * OpenMC materials during the course of the simulation.
+   * @return number of density bins
+   */
+  virtual unsigned int nDensityBins() const { return _n_density_bins; }
+
+  /**
+   * Get pointer to underlying moab interface
+   * @return pointer to moab interface
+   */
+  const std::shared_ptr<moab::Interface> & moabPtr() const { return _moab; }
+
 protected:
   std::unique_ptr<NumericVector<Number>> _serialized_solution;
 
@@ -127,12 +154,6 @@ class MoabSkinner : public GeneralUserObject
   /// Whether to print diagnostic information
   bool _verbose;
 
-  /**
-   * Whether to build a graveyard as two additional cube surfaces surrounding the mesh.
-   * This is only needed if the skinned geometry is fed into a Monte Carlo code.
-   */
-  const bool & _build_graveyard;
-
   /// Name of the temperature variable
   const std::string & _temperature_name;
 
@@ -173,6 +194,12 @@ class MoabSkinner : public GeneralUserObject
   /// Whether to output the MOAB mesh to a .h5m file
   const bool & _output_full;
 
+  /**
+   * Whether to build a graveyard as two additional cube surfaces surrounding the mesh.
+   * This is only needed if the skinned geometry is fed into a Monte Carlo code.
+   */
+  bool _build_graveyard;
+
   /// Length multiplier to get from [Mesh] units into OpenMC's centimeters
   Real _scaling;
 
@@ -309,14 +336,6 @@ class MoabSkinner : public GeneralUserObject
   /// NB elems in param is a copy, localElems is a reference
   void groupLocalElems(std::set<dof_id_type> elems, std::vector<moab::Range> & localElems);
 
-  /**
-   * Get a unique index for an OpenMC material (which we can simply take as the density bin,
-   * because in OpenMC you can set a unique temperature on different cells filled with the
-   * same material, but this is not the case for density).
-   * @return OpenMC material index
-   */
-  unsigned int getMatBin(const unsigned int & iDenBin) const;
-
   /// Clear MOAB entity sets
   bool resetMOAB();