allow overlay material

docs/source/io_formats/tallies.rst

@@ -50,6 +50,22 @@ The ``<tally>`` element accepts the following sub-elements:
 
         <nuclides>U235 Pu239 total</nuclides>
 
+    An entry of the form ``material:<id>`` is a *virtual overlay material*
+    rather than a nuclide. Such a bin scores the macroscopic response of the
+    given material (the sum over its nuclides of atom density times microscopic
+    cross section) everywhere in the geometry, including in void and in regions
+    filled with other materials. This makes it possible to map, say, the
+    response of a silicon detector without placing any silicon in the model:
+
+    .. code-block:: xml
+
+        <nuclides>material:7</nuclides>
+
+    Overlay materials require ``multiply_density`` to be false. They are not
+    supported for the ``analog`` estimator, for surface or pulse-height
+    tallies, or for scores that do not depend on the nuclide (``flux``,
+    ``inverse-velocity``, ``events`` and the IFP scores).
+
     *Default*: total
 
   :estimator:
@@ -71,7 +87,9 @@ The ``<tally>`` element accepts the following sub-elements:
 
   :multiply_density:
     A boolean that indicates whether reaction rate scores should be computed by
-    multiplying by the atom density of a nuclide present in a material.
+    multiplying by the atom density of a nuclide present in a material. Setting
+    it to false decouples the tally from whatever fills the geometry, which is
+    what makes per-nuclide and ``material:<id>`` overlay bins meaningful.
 
     *Default*: true
 

docs/source/usersguide/tallies.rst

@@ -93,6 +93,52 @@ U238 (separately), we'd set::
 You can also list 'all' as a nuclide which will give you a separate reaction
 rate for every nuclide in the model.
 
+.. _usersguide_virtual_overlay:
+
+-------------------------
+Virtual Overlay Materials
+-------------------------
+
+By default, reaction rate scores are multiplied by the atom density of the
+nuclide in whatever material the particle is travelling through, so a tally
+only responds where that nuclide is actually present. Setting
+:attr:`Tally.multiply_density` to ``False`` removes that factor, which makes it
+possible to score a response for a nuclide that appears nowhere in the model::
+
+  tally = openmc.Tally()
+  tally.scores = ['heating']
+  tally.nuclides = ['Si28']
+  tally.multiply_density = False
+
+Each bin then holds a microscopic quantity for a fictitious density of one
+atom per barn-cm. To get the response of a whole *material* rather than a
+single nuclide, pass an :class:`openmc.Material` in the nuclide list::
+
+  silicon = openmc.Material()
+  silicon.add_element('Si', 1.0)
+  silicon.set_density('g/cm3', 2.33)
+
+  tally = openmc.Tally()
+  tally.scores = ['heating']
+  tally.nuclides = [silicon]
+  tally.multiply_density = False
+
+This produces a single bin holding the material's macroscopic response, i.e.
+the sum over the material's nuclides of atom density times microscopic cross
+section, evaluated everywhere in the geometry including in void. With a
+'heating' score the result is in eV/cm\ :sup:`3` per source particle at the
+overlay material's own density, so the absorbed dose in Gy per source particle
+is the tally result multiplied by :math:`1.602\times 10^{-19}` J/eV and divided
+by the material's mass density in kg/cm\ :sup:`3`. A typical use is mapping the
+dose in a silicon detector or in tissue without perturbing the transport
+problem by adding those materials to the geometry.
+
+The overlay material does not need to be assigned to any cell, and it is
+written to the model automatically. A few restrictions apply: overlay materials
+require ``multiply_density`` to be ``False``, they cannot be used with the
+'analog' estimator or on surface and pulse-height tallies, and they cannot be
+combined with scores that do not depend on the nuclide such as 'flux'.
+
 The following tables show all valid scores:
 
 .. table:: **Flux scores: units are particle-cm per source particle.**

include/openmc/tallies/tally.h

@@ -16,6 +16,31 @@
 
 namespace openmc {
 
+//==============================================================================
+// Nuclide bin encoding
+//==============================================================================
+
+//! Entry in Tally::nuclides_ indicating the total rate over the local material
+constexpr int NUCLIDE_BIN_TOTAL {-1};
+
+//! Encode a material index as an entry in Tally::nuclides_.
+//
+//! Such a bin makes the tally respond to a virtual overlay material rather
+//! than to a single nuclide. Entries less than NUCLIDE_BIN_TOTAL are material
+//! bins; entries of zero or greater are indices into data::nuclides.
+inline int nuclide_bin_from_material(int i_material)
+{
+  return NUCLIDE_BIN_TOTAL - 1 - i_material;
+}
+
+//! Decode an entry in Tally::nuclides_ into an index into model::materials,
+//! or C_NONE if the entry is not a virtual overlay material bin.
+inline int material_from_nuclide_bin(int nuclide_bin)
+{
+  return nuclide_bin < NUCLIDE_BIN_TOTAL ? NUCLIDE_BIN_TOTAL - 1 - nuclide_bin
+                                         : C_NONE;
+}
+
 //==============================================================================
 //! A user-specified flux-weighted (or current) measurement.
 //==============================================================================
@@ -54,6 +79,9 @@ class Tally {
 
   void set_nuclides(const vector<std::string>& nuclides);
 
+  //! Whether any nuclide bin refers to a virtual overlay material
+  bool has_material_bins() const;
+
   const tensor::Tensor<double>& results() const { return results_; }
 
   //! returns vector of indices corresponding to the tally this is called on
@@ -152,8 +180,10 @@ class Tally {
 
   vector<int> scores_; //!< Filter integrands (e.g. flux, fission)
 
-  //! Index of each nuclide to be tallied.  -1 indicates total material.
-  vector<int> nuclides_ {-1};
+  //! Index of each nuclide to be tallied.  NUCLIDE_BIN_TOTAL indicates the
+  //! total rate over the local material and values below it are virtual
+  //! overlay materials, see material_from_nuclide_bin.
+  vector<int> nuclides_ {NUCLIDE_BIN_TOTAL};
 
   //! Results for each bin -- the first dimension of the array is for the
   //! combination of filters (e.g. specific cell, specific energy group, etc.)

openmc/lib/tally.py

@@ -12,6 +12,7 @@
 from .core import _FortranObjectWithID
 from .error import _error_handler
 from .filter import _get_filter
+from .material import Material
 
 
 __all__ = ['Tally', 'tallies', 'global_tallies', 'num_realizations']
@@ -320,8 +321,16 @@ def nuclides(self):
         nucs = POINTER(c_int)()
         n = c_int()
         _dll.openmc_tally_get_nuclides(self._index, nucs, n)
-        return [Nuclide(nucs[i]).name if nucs[i] >= 0 else 'total'
-                for i in range(n.value)]
+
+        def name(bin):
+            if bin >= 0:
+                return Nuclide(bin).name
+            if bin == -1:
+                return 'total'
+            # Virtual overlay material bin, see material_from_nuclide_bin
+            return f'material:{Material(index=-2 - bin).id}'
+
+        return [name(nucs[i]) for i in range(n.value)]
 
     @nuclides.setter
     def nuclides(self, nuclides):

openmc/model/model.py

@@ -245,6 +245,43 @@ def _get_all_materials(self) -> dict[int, openmc.Material]:
 
         return materials
 
+    def _materials_for_export(self) -> openmc.Materials:
+        """Collect the materials that need to be written to XML.
+
+        This is the materials collection if one was specified, otherwise all
+        materials found in the geometry. Virtual overlay materials used as
+        tally nuclide bins are appended since they need to be defined for the
+        transport solver even though nothing in the geometry is filled by them.
+
+        Returns
+        -------
+        openmc.Materials
+            Materials to export
+
+        """
+        if self.materials:
+            materials = self.materials
+        else:
+            materials = openmc.Materials(
+                self.geometry.get_all_materials().values())
+
+        existing_ids = {mat.id for mat in materials}
+        overlay = [mat for tally in self.tallies
+                   for mat in tally.overlay_materials
+                   if mat.id not in existing_ids]
+        if not overlay:
+            return materials
+
+        # Don't mutate a user-supplied collection, but keep its settings
+        combined = openmc.Materials(materials)
+        combined.cross_sections = materials.cross_sections
+        for mat in overlay:
+            if mat.id not in existing_ids:
+                combined.append(mat)
+                existing_ids.add(mat.id)
+
+        return combined
+
     def add_kinetics_parameters_tallies(self, num_groups: int | None = None):
         """Add tallies for calculating kinetics parameters using the IFP method.
 
@@ -635,12 +672,8 @@ def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False,
         # If a materials collection was specified, export it. Otherwise, look
         # for all materials in the geometry and use that to automatically build
         # a collection.
-        if self.materials:
-            self.materials.export_to_xml(d, nuclides_to_ignore=nuclides_to_ignore)
-        else:
-            materials = openmc.Materials(self.geometry.get_all_materials()
-                                         .values())
-            materials.export_to_xml(d, nuclides_to_ignore=nuclides_to_ignore)
+        materials = self._materials_for_export()
+        materials.export_to_xml(d, nuclides_to_ignore=nuclides_to_ignore)
 
         if self.tallies:
             self.tallies.export_to_xml(d)
@@ -701,11 +734,7 @@ def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool =
         # If a materials collection was specified, export it. Otherwise, look
         # for all materials in the geometry and use that to automatically build
         # a collection.
-        if self.materials:
-            materials = self.materials
-        else:
-            materials = openmc.Materials(self.geometry.get_all_materials()
-                                         .values())
+        materials = self._materials_for_export()
 
         with open(xml_path, 'w', encoding='utf-8', errors='xmlcharrefreplace') as fh:
             # write the XML header

openmc/tallies.py

@@ -47,9 +47,10 @@
 _PRODUCT_TYPES = ['tensor', 'entrywise']
 
 # The following indicate acceptable types when setting Tally.scores,
-# Tally.nuclides, and Tally.filters
+# Tally.nuclides, and Tally.filters. A Material is only accepted as a nuclide
+# when multiply_density is False, where it acts as a virtual overlay material.
 _SCORE_CLASSES = (str, CrossScore, AggregateScore)
-_NUCLIDE_CLASSES = (str, CrossNuclide, AggregateNuclide)
+_NUCLIDE_CLASSES = (str, CrossNuclide, AggregateNuclide, openmc.Material)
 _FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
 
 # Valid types of estimators
@@ -71,8 +72,10 @@ class Tally(IDManagerMixin):
         List of scores, e.g. ['flux', 'fission']
     filters : list of openmc.Filter, optional
         List of filters for the tally
-    nuclides : list of str, optional
-        List of nuclides to score results for
+    nuclides : list of str or openmc.Material, optional
+        List of nuclides to score results for. When :attr:`multiply_density` is
+        False, an :class:`openmc.Material` may be given in place of a nuclide
+        name to score the macroscopic response of a virtual overlay material.
     estimator : {'analog', 'tracklength', 'collision'}, optional
         Type of estimator for the tally
     triggers : list of openmc.Trigger, optional
@@ -92,8 +95,17 @@ class Tally(IDManagerMixin):
         .. versionadded:: 0.14.0
     filters : list of openmc.Filter
         List of specified filters for the tally
-    nuclides : list of str
-        List of nuclides to score results for
+    nuclides : list of str or openmc.Material
+        List of nuclides to score results for. A :class:`openmc.Material`
+        entry is a virtual overlay material: the bin scores that material's
+        macroscopic response (sum over its nuclides of atom density times
+        microscopic cross section) everywhere in the geometry, including in
+        void and in regions filled with other materials. Overlay materials
+        require :attr:`multiply_density` to be False and are not supported for
+        the 'analog' estimator or for scores that do not depend on the
+        nuclide, such as 'flux'.
+
+        .. versionadded:: 0.15.4
     scores : list of str
         List of defined scores, e.g. 'flux', 'fission', etc.
     estimator : {'analog', 'tracklength', 'collision'}
@@ -237,7 +249,10 @@ def __repr__(self):
             parts.append('{: <15}=\t{}'.format('Derivative ID', self.derivative.id))
         filters = ', '.join(type(f).__name__ for f in self.filters)
         parts.append('{: <15}=\t{}'.format('Filters', filters))
-        nuclides = ' '.join(str(nuclide) for nuclide in self.nuclides)
+        nuclides = ' '.join(
+            f'material:{nuclide.id}'
+            if isinstance(nuclide, openmc.Material) else str(nuclide)
+            for nuclide in self.nuclides)
         parts.append('{: <15}=\t{}'.format('Nuclides', nuclides))
         parts.append('{: <15}=\t{}'.format('Scores', self.scores))
         parts.append('{: <15}=\t{}'.format('Estimator', self.estimator))
@@ -332,6 +347,12 @@ def nuclides(self, nuclides):
         self._nuclides = cv.CheckedList(_NUCLIDE_CLASSES, 'tally nuclides',
                                         nuclides)
 
+    @property
+    def overlay_materials(self):
+        """list of openmc.Material : Virtual overlay materials among the
+        tally's nuclide bins"""
+        return [n for n in self._nuclides if isinstance(n, openmc.Material)]
+
     @property
     def num_nuclides(self):
         return max(len(self._nuclides), 1)
@@ -1446,8 +1467,14 @@ def to_xml_element(self):
 
         # Optional Nuclides
         if self.nuclides:
+            if self.overlay_materials and self.multiply_density:
+                raise ValueError(
+                    f'Tally ID="{self.id}" specifies a material as a nuclide '
+                    'bin, which requires multiply_density to be False.')
             subelement = ET.SubElement(element, "nuclides")
-            subelement.text = ' '.join(str(n) for n in self.nuclides)
+            subelement.text = ' '.join(
+                f'material:{n.id}' if isinstance(n, openmc.Material) else str(n)
+                for n in self.nuclides)
 
         # Scores
         if len(self.scores) == 0:
@@ -1633,8 +1660,11 @@ def get_nuclide_index(self, nuclide):
 
         Parameters
         ----------
-        nuclide : str
-            The name of the Nuclide (e.g., 'H1', 'U238')
+        nuclide : str or openmc.Material
+            The name of the Nuclide (e.g., 'H1', 'U238'), or a virtual overlay
+            material. An overlay material may also be given by the string
+            ``'material:<id>'``, which is how it is labelled in statepoint
+            files.
 
         Returns
         -------
@@ -1648,10 +1678,20 @@ def get_nuclide_index(self, nuclide):
             in the Tally.
 
         """
+        # An overlay material is labelled 'material:<id>' once results have been
+        # read back, so match a Material against either representation
+        label = (f'material:{nuclide.id}'
+                 if isinstance(nuclide, openmc.Material) else nuclide)
+
         # Look for the user-requested nuclide in all of the Tally's nuclides
         for i, test_nuclide in enumerate(self.nuclides):
             if test_nuclide == nuclide:
                 return i
+            test_label = (f'material:{test_nuclide.id}'
+                          if isinstance(test_nuclide, openmc.Material)
+                          else test_nuclide)
+            if test_label == label:
+                return i
 
         msg = (f'Unable to get the nuclide index for Tally since "{nuclide}" '
                'is not one of the nuclides')
@@ -1760,8 +1800,8 @@ def get_nuclide_indices(self, nuclides):
 
         Parameters
         ----------
-        nuclides : list of str
-            A list of nuclide name strings
+        nuclides : list of str or openmc.Material
+            A list of nuclide name strings or virtual overlay materials
             (e.g., ['U235', 'U238']; default is [])
 
         Returns
@@ -1771,14 +1811,14 @@ def get_nuclide_indices(self, nuclides):
 
         """
 
-        cv.check_iterable_type('nuclides', nuclides, str)
+        cv.check_iterable_type('nuclides', nuclides, (str, openmc.Material))
 
         # If user did not specify any specific Nuclides, use them all
         if not nuclides:
             return np.arange(self.num_nuclides)
 
         # Determine the score indices from any of the requested scores
-        nuclide_indices = np.zeros_like(nuclides, dtype=int)
+        nuclide_indices = np.zeros(len(nuclides), dtype=int)
         for i, nuclide in enumerate(nuclides):
             nuclide_indices[i] = self.get_nuclide_index(nuclide)
         return nuclide_indices
@@ -1987,6 +2027,8 @@ def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True,
                 if isinstance(nuclide, openmc.AggregateNuclide):
                     nuclides.append(nuclide.name)
                     column_name = f'{nuclide.aggregate_op}(nuclide)'
+                elif isinstance(nuclide, openmc.Material):
+                    nuclides.append(f'material:{nuclide.id}')
                 else:
                     nuclides.append(nuclide)