Add distributed cell densities

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* rho)
+
+   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* rho: density of the cell in g/cc
+   :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
@@ -119,6 +130,20 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_cell_set_density(index index, double rho, const int32_t* instance, bool set_contained)
+
+   Set the density of a cell.
+
+   :param int32_t index: Index in the cells array
+   :param double rho: Density of the cell in g/cc
+   :param instance: Which instance of the cell. To set the density multiplier for all
+                    instances, pass a null pointer.
+   :param 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
+
 .. c:function:: int openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n)
 
    Return the bounding energies for an energy filter

docs/source/io_formats/properties.rst

@@ -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

@@ -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,19 @@ 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 multiplier 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 +239,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] rho 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 rho, int32_t instance = -1, bool set_contained = false);
+
   int32_t n_instances() const;
 
   //! Set the rotation matrix of a cell instance
@@ -334,6 +357,9 @@ class Cell {
   //! T. The units are sqrt(eV).
   vector<double> sqrtkT_;
 
+  //! \brief Unitless density multiplier(s) within this cell.
+  vector<double> rho_mult_ = {1.0};
+
   //! \brief Neighboring cells in the same universe.
   NeighborList neighbors_;
 
@@ -351,6 +377,10 @@ class Cell {
 
   // Right now, either CSG or DAGMC cells are used.
   virtual GeometryType geom_type() const = 0;
+
+  //! \brief A flag to indicate if this cell has it's density set in the
+  //! xml file.
+  bool xml_set_density_ = false;
 };
 
 struct CellInstanceItem {

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& rho_mult() { return rho_mult_; }
+  const double& rho_mult() const { return rho_mult_; }
+  double& rho_mult_last() { return rho_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 rho_mult_ {1.0};      //!< density multiplier
+  double rho_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,8 @@ 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 +264,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 +556,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 +683,13 @@ 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 +751,16 @@ 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:
+            if len(density) > 1:
+                c.density = density
+            else:
+                c.density = 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:

openmc/lib/cell.py

@@ -34,6 +34,10 @@
     c_int32, POINTER(c_int32), POINTER(c_double)]
 _dll.openmc_cell_get_temperature.restype = c_int
 _dll.openmc_cell_get_temperature.errcheck = _error_handler
+_dll.openmc_cell_get_density.argtypes = [
+    c_int32, POINTER(c_int32), POINTER(c_double)]
+_dll.openmc_cell_get_density.restype = c_int
+_dll.openmc_cell_get_density.errcheck = _error_handler
 _dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)]
 _dll.openmc_cell_get_name.restype = c_int
 _dll.openmc_cell_get_name.errcheck = _error_handler
@@ -58,6 +62,10 @@
     c_int32, c_double, POINTER(c_int32), c_bool]
 _dll.openmc_cell_set_temperature.restype = c_int
 _dll.openmc_cell_set_temperature.errcheck = _error_handler
+_dll.openmc_cell_set_density.argtypes = [
+    c_int32, c_double, POINTER(c_int32), c_bool]
+_dll.openmc_cell_set_density.restype = c_int
+_dll.openmc_cell_set_density.errcheck = _error_handler
 _dll.openmc_cell_set_translation.argtypes = [c_int32, POINTER(c_double)]
 _dll.openmc_cell_set_translation.restype = c_int
 _dll.openmc_cell_set_translation.errcheck = _error_handler
@@ -236,6 +244,43 @@ def set_temperature(self, T, instance=None, set_contained=False):
 
         _dll.openmc_cell_set_temperature(self._index, T, instance, set_contained)
 
+    def get_density(self, instance=None):
+        """Get the density of a cell (g/cc)
+
+        Parameters
+        ----------
+        instance: int or None
+            Which instance of the cell
+
+        """
+
+        if instance is not None:
+            instance = c_int32(instance)
+
+        rho = c_double()
+        _dll.openmc_cell_get_density(self._index, instance, rho)
+        return rho.value
+
+    def set_density(self, rho, instance=None, set_contained=False):
+        """Set the density of a cell
+
+        Parameters
+        ----------
+        rho : float
+            Density of the cell (g/cc)
+        instance : int or None
+            Which instance of the cell
+        set_contained: bool
+            If cell is not filled by a material, whether to set the density
+            of all filled cells
+
+        """
+
+        if instance is not None:
+            instance = c_int32(instance)
+
+        _dll.openmc_cell_set_density(self._index, rho, instance, set_contained)
+
     @property
     def translation(self):
         translation = np.zeros(3)

openmc/model/model.py

@@ -633,7 +633,7 @@ def import_properties(self, filename: PathLike):
                 raise ValueError("Number of cells in properties file doesn't "
                                  "match current model.")
 
-            # Update temperatures for cells filled with materials
+            # Update temperatures and densities for cells filled with materials
             for name, group in cells_group.items():
                 cell_id = int(name.split()[1])
                 cell = cells[cell_id]
@@ -648,6 +648,20 @@ def import_properties(self, filename: PathLike):
                         else:
                             lib_cell.set_temperature(temperature[0])
 
+                    if group['density']:
+                      density = group['density'][()]
+                      if density.size > 1:
+                          cell.density = [rho for rho in density]
+                      else:
+                          cell.density = density
+                      if self.is_initialized:
+                          lib_cell = openmc.lib.cells[cell_id]
+                          if density.size > 1:
+                              for i, rho in enumerate(density):
+                                  lib_cell.set_density(rho, i)
+                          else:
+                              lib_cell.set_density(density[0])
+
             # Make sure number of materials matches
             mats_group = fh['materials']
             n_cells = mats_group.attrs['n_materials']