Add distributed cell densities (#3546)

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: 3 people

docs/source/capi/index.rst

@@ -84,6 +84,17 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)
+
+   Get the density of a cell
+
+   :param int32_t index: Index in the cells array
+   :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the density
+                             multiplier of the first instance is returned.
+   :param double* density: Density of the cell in [g/cm3]
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices)
 
    Set the fill for a cell
@@ -113,8 +124,22 @@ Functions
    :param double T: Temperature in Kelvin
    :param instance: Which instance of the cell. To set the temperature for all
                     instances, pass a null pointer.
-   :param set_contained: If the cell is not filled by a material, whether to set the temperatures
-                         of all filled cells
+   :param bool set_contained: If the cell is not filled by a material, whether
+                              to set the temperatures of all filled cells
+   :type instance: const int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_cell_set_density(index index, double density, const int32_t* instance, bool set_contained)
+
+   Set the density of a cell.
+
+   :param int32_t index: Index in the cells array
+   :param double density: Density of the cell in [g/cm3]
+   :param instance: Which instance of the cell. To set the density multiplier for all
+                    instances, pass a null pointer.
+   :param bool set_contained: If the cell is not filled by a material, whether
+                              to set the density multiplier of all filled cells
    :type instance: const int32_t*
    :return: Return status (negative if an error occurred)
    :rtype: int

docs/source/io_formats/properties.rst

@@ -4,7 +4,7 @@
 Properties File Format
 ======================
 
-The current version of the properties file format is 1.0.
+The current version of the properties file format is 1.1.
 
 **/**
 
@@ -25,6 +25,7 @@ The current version of the properties file format is 1.0.
 **/geometry/cells/cell <uid>/**
 
 :Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K].
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
 
 **/materials/**
 

docs/source/io_formats/summary.rst

@@ -4,7 +4,7 @@
 Summary File Format
 ===================
 
-The current version of the summary file format is 6.0.
+The current version of the summary file format is 6.1.
 
 **/**
 
@@ -38,6 +38,7 @@ The current version of the summary file format is 6.0.
              is an array if the cell uses distributed materials, otherwise it is
              a scalar.
            - **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
            - **translation** (*double[3]*) -- Translation applied to the fill
              universe. This dataset is present only if fill_type is set to
              'universe'.

include/openmc/capi.h

@@ -17,6 +17,8 @@ int openmc_cell_get_fill(
 int openmc_cell_get_id(int32_t index, int32_t* id);
 int openmc_cell_get_temperature(
   int32_t index, const int32_t* instance, double* T);
+int openmc_cell_get_density(
+  int32_t index, const int32_t* instance, double* rho);
 int openmc_cell_get_translation(int32_t index, double xyz[]);
 int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n);
 int openmc_cell_get_name(int32_t index, const char** name);
@@ -27,6 +29,8 @@ int openmc_cell_set_fill(
 int openmc_cell_set_id(int32_t index, int32_t id);
 int openmc_cell_set_temperature(
   int32_t index, double T, const int32_t* instance, bool set_contained = false);
+int openmc_cell_set_density(int32_t index, double rho, const int32_t* instance,
+  bool set_contained = false);
 int openmc_cell_set_translation(int32_t index, const double xyz[]);
 int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len);
 int openmc_dagmc_universe_get_cell_ids(

include/openmc/cell.h

@@ -216,6 +216,18 @@ class Cell {
   //! \return Temperature in [K]
   double temperature(int32_t instance = -1) const;
 
+  //! Get the density multiplier of a cell instance
+  //! \param[in] instance Instance index. If -1 is given, the density multiplier
+  //! for the first instance is returned.
+  //! \return Density multiplier
+  double density_mult(int32_t instance = -1) const;
+
+  //! Get the density of a cell instance in g/cm3
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //! for the first instance is returned.
+  //! \return Density in [g/cm3]
+  double density(int32_t instance = -1) const;
+
   //! Set the temperature of a cell instance
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
@@ -226,6 +238,16 @@ class Cell {
   void set_temperature(
     double T, int32_t instance = -1, bool set_contained = false);
 
+  //! Set the density of a cell instance
+  //! \param[in] density Density [g/cm3]
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //!   for all instances is set.
+  //! \param[in] set_contained If this cell is not filled with a material,
+  //!   collect all contained cells with material fills and set their
+  //!   densities.
+  void set_density(
+    double density, int32_t instance = -1, bool set_contained = false);
+
   int32_t n_instances() const;
 
   //! Set the rotation matrix of a cell instance
@@ -334,6 +356,9 @@ class Cell {
   //! T. The units are sqrt(eV).
   vector<double> sqrtkT_;
 
+  //! \brief Unitless density multiplier(s) within this cell.
+  vector<double> density_mult_;
+
   //! \brief Neighboring cells in the same universe.
   NeighborList neighbors_;
 

include/openmc/constants.h

@@ -28,11 +28,11 @@ constexpr int HDF5_VERSION[] {3, 0};
 constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
 constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr array<int, 2> VERSION_TRACK {3, 0};
-constexpr array<int, 2> VERSION_SUMMARY {6, 0};
+constexpr array<int, 2> VERSION_SUMMARY {6, 1};
 constexpr array<int, 2> VERSION_VOLUME {1, 0};
 constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
-constexpr array<int, 2> VERSION_PROPERTIES {1, 0};
+constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
 
 // ============================================================================

include/openmc/geometry_aux.h

@@ -37,6 +37,12 @@ void adjust_indices();
 
 void assign_temperatures();
 
+//==============================================================================
+//! Finalize densities (compute density multipliers).
+//==============================================================================
+
+void finalize_cell_densities();
+
 //==============================================================================
 //! \brief Obtain a list of temperatures that each nuclide/thermal scattering
 //! table appears at in the model. Later, this list is used to determine the

include/openmc/material.h

@@ -99,6 +99,13 @@ class Material {
   //----------------------------------------------------------------------------
   // Accessors
 
+  //! Get the atom density in [atom/b-cm]
+  //! \return Density in [atom/b-cm]
+  double atom_density(int32_t i, double rho_multiplier = 1.0) const
+  {
+    return atom_density_(i) * rho_multiplier;
+  }
+
   //! Get density in [atom/b-cm]
   //! \return Density in [atom/b-cm]
   double density() const { return density_; }

include/openmc/particle_data.h

@@ -389,6 +389,11 @@ class GeometryState {
   const double& sqrtkT() const { return sqrtkT_; }
   double& sqrtkT_last() { return sqrtkT_last_; }
 
+  // density multiplier of the current and last cell
+  double& density_mult() { return density_mult_; }
+  const double& density_mult() const { return density_mult_; }
+  double& density_mult_last() { return density_mult_last_; }
+
 private:
   int64_t id_ {-1}; //!< Unique ID
 
@@ -417,6 +422,9 @@ class GeometryState {
   double sqrtkT_ {-1.0};     //!< sqrt(k_Boltzmann * temperature) in eV
   double sqrtkT_last_ {0.0}; //!< last temperature
 
+  double density_mult_ {1.0};      //!< density multiplier
+  double density_mult_last_ {1.0}; //!< last density multiplier
+
 #ifdef OPENMC_DAGMC_ENABLED
   moab::DagMC::RayHistory history_;
   Direction last_dir_;

openmc/cell.py

@@ -72,6 +72,10 @@ class Cell(IDManagerMixin):
     temperature : float or iterable of float
         Temperature of the cell in Kelvin.  Multiple temperatures can be given
         to give each distributed cell instance a unique temperature.
+    density : float or iterable of float
+        Density of the cell in [g/cm3]. Multiple densities can be given to give
+        each distributed cell instance a unique density. Densities set here will
+        override the density set on materials used to fill the cell.
     translation : Iterable of float
         If the cell is filled with a universe, this array specifies a vector
         that is used to translate (shift) the universe.
@@ -109,6 +113,7 @@ def __init__(self, cell_id=None, name='', fill=None, region=None):
         self._rotation = None
         self._rotation_matrix = None
         self._temperature = None
+        self._density = None
         self._translation = None
         self._paths = None
         self._num_instances = None
@@ -141,6 +146,7 @@ def __repr__(self):
         if self.fill_type == 'material':
             string += '\t{0: <15}=\t{1}\n'.format('Temperature',
                                                   self.temperature)
+            string += '\t{0: <15}=\t{1}\n'.format('Density', self.density)
         string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
         string += '{: <16}=\t{}\n'.format('\tVolume', self.volume)
 
@@ -257,6 +263,30 @@ def temperature(self, temperature):
         else:
             self._temperature = temperature
 
+    @property
+    def density(self):
+        return self._density
+
+    @density.setter
+    def density(self, density):
+        # Make sure densities are greater than zero
+        cv.check_type('cell density', density, (Iterable, Real), none_ok=True)
+        if isinstance(density, Iterable):
+            cv.check_type('cell density', density, Iterable, Real)
+            for rho in density:
+                cv.check_greater_than('cell density', rho, 0.0, True)
+        elif isinstance(density, Real):
+            cv.check_greater_than('cell density', density, 0.0, True)
+
+        # If this cell is filled with a universe or lattice, propagate
+        # densities to all cells contained. Otherwise, simply assign it.
+        if self.fill_type in ('universe', 'lattice'):
+            for c in self.get_all_cells().values():
+                if c.fill_type == 'material':
+                    c._density = density
+        else:
+            self._density = density
+
     @property
     def translation(self):
         return self._translation
@@ -525,6 +555,8 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None):
             clone.volume = self.volume
             if self.temperature is not None:
                 clone.temperature = self.temperature
+            if self.density is not None:
+                clone.density = self.density
             if self.translation is not None:
                 clone.translation = self.translation
             if self.rotation is not None:
@@ -650,6 +682,12 @@ def create_surface_elements(node, element, memo=None):
             else:
                 element.set("temperature", str(self.temperature))
 
+        if self.density is not None:
+            if isinstance(self.density, Iterable):
+                element.set("density", ' '.join(str(t) for t in self.density))
+            else:
+                element.set("density", str(self.density))
+
         if self.translation is not None:
             element.set("translation", ' '.join(map(str, self.translation)))
 
@@ -711,10 +749,13 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe):
                 c.temperature = temperature
             else:
                 c.temperature = temperature[0]
+        density = get_elem_list(elem, 'density', float)
+        if density is not None:
+            c.density = density if len(density) > 1 else density[0]
         v = get_text(elem, 'volume')
         if v is not None:
             c.volume = float(v)
-        for key in ('temperature', 'rotation', 'translation'):
+        for key in ('temperature', 'density', 'rotation', 'translation'):
             values = get_elem_list(elem, key, float)
             if values is not None:
                 if key == 'rotation' and len(values) == 9: