implementation of the generalized methods to fetch particle sources from chain file

Author: marco-de-pietri

openmc/config.py

@@ -20,7 +20,7 @@
 from typing import Any, Dict, Iterator
 
 from openmc.data import DataLibrary
-from openmc.data.decay import _DECAY_ENERGY, _DECAY_PHOTON_ENERGY
+from openmc.data.decay import _DECAY_ENERGY, _DECAY_PARTICLE_ENERGY
 
 __all__ = ["config"]
 
@@ -77,7 +77,7 @@ def __delitem__(self, key: str):
             if env_var in os.environ:
                 del os.environ[env_var]
         if key == 'chain_file':
-            _DECAY_PHOTON_ENERGY.clear()
+            _DECAY_PARTICLE_ENERGY.clear()
             _DECAY_ENERGY.clear()
 
     def __setitem__(self, key: str, value: Any):
@@ -103,7 +103,7 @@ def __setitem__(self, key: str, value: Any):
             os.environ[self._PATH_KEYS[key]] = str(stored_path)
 
             if key == 'chain_file':
-                _DECAY_PHOTON_ENERGY.clear()
+                _DECAY_PARTICLE_ENERGY.clear()
                 _DECAY_ENERGY.clear()
 
             if not stored_path.exists():

openmc/data/decay.py

@@ -561,31 +561,32 @@ def sources(self):
         return merged_sources
 
 
-_DECAY_PHOTON_ENERGY = {}
+_DECAY_PARTICLE_ENERGY = {}
 
 
-def decay_photon_energy(nuclide: str) -> Univariate | None:
-    """Get photon energy distribution resulting from the decay of a nuclide
+def decay_particle_energy(nuclide: str, particle: str) -> Univariate | None:
+    """Get a decay particle energy distribution resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
     for the first time, you need to ensure that a depletion chain has been
     specified in openmc.config['chain_file'].
 
-    .. versionadded:: 0.13.2
-
     Parameters
     ----------
     nuclide : str
         Name of nuclide, e.g., 'Co58'
+    particle : str
+        Name of the decay particle, e.g., 'photon', 'neutron'
 
     Returns
     -------
     openmc.stats.Univariate or None
-        Distribution of energies in [eV] of photons emitted from decay, or None
-        if no photon source exists. Note that the probabilities represent
+        Distribution of energies in [eV] of decay particles emitted from decay, or None
+        if no decay source of that particle exists. Note that the probabilities represent
         intensities, given as [Bq/atom] (in other words, decay constants).
     """
-    if not _DECAY_PHOTON_ENERGY:
+
+    if not _DECAY_PARTICLE_ENERGY:
         chain_file = openmc.config.get('chain_file')
         if chain_file is None:
             raise DataError(
@@ -596,15 +597,15 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
         from openmc.deplete import Chain
         chain = Chain.from_xml(chain_file)
         for nuc in chain.nuclides:
-            if 'photon' in nuc.sources:
-                _DECAY_PHOTON_ENERGY[nuc.name] = nuc.sources['photon']
+            for part, dist in nuc.sources.items():
+                _DECAY_PARTICLE_ENERGY[(nuc.name, part)] = dist
 
         # If the chain file contained no sources at all, warn the user
-        if not _DECAY_PHOTON_ENERGY:
-            warn(f"Chain file '{chain_file}' does not have any decay photon "
-                 "sources listed.")
+        if not _DECAY_PARTICLE_ENERGY:
+            warn(f"Chain file '{chain_file}' does not have any decay particle "
+                 "source listed.")
 
-    return _DECAY_PHOTON_ENERGY.get(nuclide)
+    return _DECAY_PARTICLE_ENERGY.get((nuclide, particle))
 
 
 _DECAY_ENERGY = {}

openmc/material.py

@@ -332,20 +332,21 @@ def decay_photon_energy(self) -> Univariate | None:
             "version.", FutureWarning)
         return self.get_decay_photon_energy(0.0)
 
-    def get_decay_photon_energy(
+    def get_decay_particle_energy(
         self,
+        particle: str = 'photon',
         clip_tolerance: float = 1e-6,
         units: str = 'Bq',
         volume: float | None = None,
         exclude_nuclides: list[str] | None = None,
         include_nuclides: list[str] | None = None
     ) -> Univariate | None:
-        r"""Return energy distribution of decay photons from unstable nuclides.
-
-        .. versionadded:: 0.14.0
+        r"""Return energy distribution of decay particles from unstable nuclides.
 
         Parameters
         ----------
+        particle : string
+            Particle type to return the energy distribution for.
         clip_tolerance : float
             Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions
             that will be discarded.
@@ -355,19 +356,21 @@ def get_decay_photon_energy(
             Volume of the material. If not passed, defaults to using the
             :attr:`Material.volume` attribute.
         exclude_nuclides : list of str, optional
-            Nuclides to exclude from the photon source calculation.
+            Nuclides to exclude from the particle source calculation.
         include_nuclides : list of str, optional
-            Nuclides to include in the photon source calculation. If specified,
+            Nuclides to include in the particle source calculation. If specified,
             only these nuclides are used.
 
         Returns
         -------
         Univariate or None
-            Decay photon energy distribution. The integral of this distribution is
-            the total intensity of the photon source in the requested units.
+            Decay particle energy distribution. The integral of this distribution is
+            the total intensity of the particle source in the requested units.
 
         """
         cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'})
+        cv.check_type('particle', particle, str)
+        cv.check_value('particle', particle, {'alpha', 'fragment', 'neutron','photon', 'electron', 'positron','proton'})
 
         if exclude_nuclides is not None and include_nuclides is not None:
             raise ValueError("Cannot specify both exclude_nuclides and include_nuclides")
@@ -393,12 +396,12 @@ def get_decay_photon_energy(
             if include_nuclides is not None and nuc not in include_nuclides:
                 continue
 
-            source_per_atom = openmc.data.decay_photon_energy(nuc)
+            source_per_atom = openmc.data.decay_particle_energy(nuc,particle)
             if source_per_atom is not None and atoms_per_bcm > 0.0:
                 dists.append(source_per_atom)
                 probs.append(1e24 * atoms_per_bcm * multiplier)
 
-        # If no photon sources, exit early
+        # If no particle sources, exit early
         if not dists:
             return None
 
@@ -414,6 +417,54 @@ def get_decay_photon_energy(
 
         return combined
 
+    def get_decay_photon_energy(
+        self,
+        clip_tolerance: float = 1e-6,
+        units: str = 'Bq',
+        volume: float | None = None,
+        exclude_nuclides: list[str] | None = None,
+        include_nuclides: list[str] | None = None
+    ) -> Univariate | None:
+        r"""Return energy distribution of decay photons from unstable nuclides.
+
+        This is a convenience wrapper around :meth:`Material.get_decay_particle_energy`
+        with ``particle='photon'``.
+
+        .. versionadded:: 0.14.0
+
+        Parameters
+        ----------
+        clip_tolerance : float
+            Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions
+            that will be discarded.
+        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
+            Specifies the units on the integral of the distribution.
+        volume : float, optional
+            Volume of the material. If not passed, defaults to using the
+            :attr:`Material.volume` attribute.
+        exclude_nuclides : list of str, optional
+            Nuclides to exclude from the photon source calculation.
+        include_nuclides : list of str, optional
+            Nuclides to include in the photon source calculation. If specified,
+            only these nuclides are used.
+
+        Returns
+        -------
+        Univariate or None
+            Decay photon energy distribution. The integral of this distribution is
+            the total intensity of the photon source in the requested units.
+
+        """
+
+        return self.get_decay_particle_energy(
+            particle='photon',
+            clip_tolerance=clip_tolerance,
+            units=units,
+            volume=volume,
+            exclude_nuclides=exclude_nuclides,
+            include_nuclides=include_nuclides,
+        )
+
     def get_photon_contact_dose_rate(
         self,
         dose_quantity: str = "absorbed-air",
@@ -540,7 +591,7 @@ def get_photon_contact_dose_rate(
                           * sv_per_psv * 1e24)
 
         for nuc, nuc_atoms_per_bcm in self.get_nuclide_atom_densities().items():
-            photon_source_per_atom = openmc.data.decay_photon_energy(nuc)
+            photon_source_per_atom = openmc.data.decay_particle_energy(nuc,'photon')
 
             # nuclides with no contribution
             if photon_source_per_atom is None or nuc_atoms_per_bcm <= 0.0:

tests/unit_tests/test_config.py

@@ -99,6 +99,6 @@ def test_config_chain_side_effect(tmp_path):
     """Test that modifying chain_file clears decay data caches."""
     chain_file = tmp_path / "chain.xml"; chain_file.touch()
     decay._DECAY_ENERGY['U235'] = (1.0, 2.0)
-    decay._DECAY_PHOTON_ENERGY['PU239'] = {}
+    decay._DECAY_PARTICLE_ENERGY['PU239','photon'] = {}
     openmc.config['chain_file'] = chain_file
-    assert not decay._DECAY_ENERGY and not decay._DECAY_PHOTON_ENERGY
+    assert not decay._DECAY_ENERGY and not decay._DECAY_PARTICLE_ENERGY

tests/unit_tests/test_data_decay.py

@@ -131,18 +131,38 @@ def test_decay_photon_energy():
     if 'chain_file' in openmc.config:
         del openmc.config['chain_file']
     with pytest.raises(DataError):
-        openmc.data.decay_photon_energy('I135')
+        openmc.data.decay_particle_energy('I135','photon')
 
     # Set chain file to simple chain
     openmc.config['chain_file'] = Path(__file__).parents[1] / "chain_simple.xml"
 
     # Check strength of I135 source and presence of specific spectral line
-    src = openmc.data.decay_photon_energy('I135')
+    src = openmc.data.decay_particle_energy('I135','photon')
     assert isinstance(src, openmc.stats.Discrete)
     assert src.integral() == pytest.approx(3.920996223799345e-05)
     assert 1260409. in src.x
 
     # Check Xe135 source, which should be tabular
-    src = openmc.data.decay_photon_energy('Xe135')
+    src = openmc.data.decay_particle_energy('Xe135','photon')
     assert isinstance(src, openmc.stats.Tabular)
     assert src.integral() == pytest.approx(2.076506258964966e-05)
+
+
+def test_decay_particle_energy():
+    # If chain file is not set, we should get a data error
+    if 'chain_file' in openmc.config:
+        del openmc.config['chain_file']
+    with pytest.raises(DataError):
+        openmc.data.decay_particle_energy('Fe55', 'electron')
+
+    # Set chain file to Ni chain and check electron/positron sources
+    with openmc.config.patch('chain_file', Path(__file__).parents[1] / "chain_ni.xml"):
+        src = openmc.data.decay_particle_energy('Fe55', 'electron')
+        assert isinstance(src, openmc.stats.Discrete)
+        assert src.integral() == pytest.approx(4.225437849490689e-08)
+        assert 6377.571 in src.x
+
+        src = openmc.data.decay_particle_energy('Fe55', 'positron')
+        assert isinstance(src, openmc.stats.Discrete)
+        assert src.integral() == pytest.approx(8.004558990612365e-09)
+        assert 231210.0 in src.x

tests/unit_tests/test_material.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 import openmc
-from openmc.data import decay_photon_energy
+from openmc.data import decay_particle_energy
 from openmc.deplete import Chain
 import openmc.examples
 import openmc.model
@@ -699,13 +699,13 @@ def intensity(src):
         return src.integral() if src is not None else 0.0
 
     assert src.integral() == pytest.approx(sum(
-        intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides()
+        intensity(decay_particle_energy(nuc,"photon")) for nuc in m.get_nuclides()
     ), rel=1e-3)
 
     # When the clipping threshold is zero, the intensities should match exactly
     src = m.get_decay_photon_energy(0.0)
     assert src.integral() == pytest.approx(sum(
-        intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides()
+        intensity(decay_particle_energy(nuc,"photon")) for nuc in m.get_nuclides()
     ))
 
     # A material with no unstable nuclides should have no decay photon source
@@ -715,6 +715,54 @@ def intensity(src):
     assert stable.get_decay_photon_energy() is None
 
 
+def test_decay_particle_energy():
+    with openmc.config.patch('chain_file', Path(__file__).parents[1] / 'chain_ni.xml'):
+        # Material representing single atom of Fe55 and Co58
+        m = openmc.Material()
+        m.add_nuclide('Fe55', 1.0e-24)
+        m.add_nuclide('Co58', 1.0e-24)
+        m.volume = 1.0
+
+        # Get decay positron source and make sure it's the right type
+        src = m.get_decay_particle_energy('positron')
+        assert isinstance(src, openmc.stats.Discrete)
+
+        # Make sure units/volume work as expected
+        src_v2 = m.get_decay_particle_energy('positron', volume=2.0)
+        assert src.p * 2.0 == pytest.approx(src_v2.p)
+        src_per_cm3 = m.get_decay_particle_energy('positron', units='Bq/cm3',
+                                                  volume=100.0)
+        assert (src.p == src_per_cm3.p).all()
+        src_per_bqg = m.get_decay_particle_energy('positron', units='Bq/g')
+        src_per_bqkg = m.get_decay_particle_energy('positron', units='Bq/kg')
+        assert pytest.approx(src_per_bqg.integral()) == src_per_bqkg.integral() / 1000.
+        src_per_bqm3 = m.get_decay_particle_energy('positron', units='Bq/m3')
+        assert pytest.approx(src_per_bqm3.integral()) == src_per_cm3.integral() * 1e6
+
+        # With a single atom of each, the intensity of the positron source
+        # should be equal to the sum of the intensities for each nuclide
+        def intensity(src):
+            return src.integral() if src is not None else 0.0
+
+        assert src.integral() == pytest.approx(sum(
+            intensity(decay_particle_energy(nuc, 'positron'))
+            for nuc in m.get_nuclides()
+        ), rel=1e-3)
+
+        # When the clipping threshold is zero, the intensities should match exactly
+        src = m.get_decay_particle_energy('positron', 0.0)
+        assert src.integral() == pytest.approx(sum(
+            intensity(decay_particle_energy(nuc, 'positron'))
+            for nuc in m.get_nuclides()
+        ))
+
+        # A material with no unstable nuclides should have no decay positron source
+        stable = openmc.Material()
+        stable.add_nuclide('Gd156', 1.0)
+        stable.volume = 1.0
+        assert stable.get_decay_particle_energy('positron') is None
+
+
 def test_avoid_subnormal(run_in_tmpdir):
     # Write a materials.xml with a material that has a nuclide density that is
     # represented as a subnormal floating point value