Add chain parameter to Material.get_activity for half-life data

openmc/data/data.py

@@ -9,6 +9,8 @@
 from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL,
                        EV_PER_MEV, K_BOLTZMANN, gnds_name, zam)
 
+import openmc
+
 gnds_name.__module__ = __name__
 zam.__module__ = __name__
 
@@ -296,25 +298,44 @@ def atomic_weight(element):
         raise ValueError(f"No naturally-occurring isotopes for element '{element}'.")
 
 
-def half_life(isotope):
+def half_life(isotope, chain_file=False):
     """Return half-life of isotope in seconds or None if isotope is stable
 
-    Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
-    <https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
+    By default, half-life values are from the `ENDF/B-VIII.0 decay sublibrary
+    <https://www.nndc.bnl.gov/endf-b8.0/download.html>`_. A depletion chain can
+    also be used as the source of half-life values.
 
     .. versionadded:: 0.13.1
 
+    .. versionchanged:: 0.15.4
+        Added the ``chain_file`` argument.
+
     Parameters
     ----------
     isotope : str
         Name of isotope, e.g., 'Pu239'
+    chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
+        Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
+        used. If ``None``, the chain specified by
+        ``openmc.config['chain_file']`` is used when available. If a path or
+        :class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
+        or an explicit chain, nuclides absent from the chain fall back to
+        ENDF/B-VIII.0 data.
 
     Returns
     -------
-    float
-        Half-life of isotope in [s]
+    float or None
+        Half-life of isotope in [s], or None if the isotope is stable
 
     """
+    if chain_file is not False:
+        if chain_file is not None or openmc.config.get('chain_file') is not None:
+            # Local import avoids a circular dependency
+            from openmc.deplete.chain import _get_chain
+            chain = _get_chain(chain_file)
+            if isotope in chain:
+                return chain[isotope].half_life
+
     global _HALF_LIFE
     if not _HALF_LIFE:
         # Load ENDF/B-VIII.0 data from JSON file
@@ -324,7 +345,7 @@ def half_life(isotope):
     return _HALF_LIFE.get(isotope.lower())
 
 
-def decay_constant(isotope):
+def decay_constant(isotope, chain_file=False):
     """Return decay constant of isotope in [s^-1]
 
     Decay constants are based on half-life values from the
@@ -333,10 +354,20 @@ def decay_constant(isotope):
 
     .. versionadded:: 0.13.1
 
+    .. versionchanged:: 0.15.4
+        Added the ``chain_file`` argument.
+
     Parameters
     ----------
     isotope : str
         Name of isotope, e.g., 'Pu239'
+    chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
+        Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
+        used. If ``None``, the chain specified by
+        ``openmc.config['chain_file']`` is used when available. If a path or
+        :class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
+        or an explicit chain, nuclides absent from the chain fall back to
+        ENDF/B-VIII.0 data.
 
     Returns
     -------
@@ -348,7 +379,7 @@ def decay_constant(isotope):
     openmc.data.half_life
 
     """
-    t = half_life(isotope)
+    t = half_life(isotope, chain_file)
     return _LOG_TWO / t if t else 0.0
 
 
@@ -496,5 +527,3 @@ def isotopes(element: str) -> list[tuple[str, float]]:
             result.append(kv)
 
     return result
-
-

openmc/deplete/results.py

@@ -7,6 +7,8 @@
 import h5py
 import numpy as np
 
+import openmc
+from .chain import _get_chain
 from .stepresult import StepResult, VERSION_RESULTS
 import openmc.checkvalue as cv
 from openmc.data import atomic_mass, AVOGADRO
@@ -103,7 +105,8 @@ def get_activity(
         mat: Material | str,
         units: str = "Bq/cm3",
         by_nuclide: bool = False,
-        volume: float | None = None
+        volume: float | None = None,
+        chain_file=None
     ) -> tuple[np.ndarray, np.ndarray | list[dict]]:
         """Get activity of material over time.
 
@@ -115,22 +118,31 @@ def get_activity(
             Material object or material id to evaluate
         units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
             Specifies the type of activity to return, options include total
-            activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity [Bq/cm3].
+            activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity
+            [Bq/cm3].
         by_nuclide : bool
             Specifies if the activity should be returned for the material as a
             whole or per nuclide. Default is False.
         volume : float, optional
             Volume of the material. If not passed, defaults to using the
             :attr:`Material.volume` attribute.
+        chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
+            Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
+            used. If ``None``, the chain specified by
+            ``openmc.config['chain_file']`` is used when available. If a path or
+            :class:`openmc.deplete.Chain` is given, that chain is used. For
+            ``None`` or an explicit chain, nuclides absent from the chain fall
+            back to ENDF/B-VIII.0 data.
+
+            .. versionadded:: 0.15.4
 
         Returns
         -------
         times : numpy.ndarray
             Array of times in [s]
         activities : numpy.ndarray or List[dict]
-            Array of total activities if by_nuclide = False (default)
-            or list of dictionaries of activities by nuclide if
-            by_nuclide = True.
+            Array of total activities if by_nuclide = False (default) or list of
+            dictionaries of activities by nuclide if by_nuclide = True.
 
         """
         if isinstance(mat, Material):
@@ -140,6 +152,13 @@ def get_activity(
         else:
             raise TypeError('mat should be of type openmc.Material or str')
 
+        if chain_file is not False:
+            if chain_file is None:
+                if openmc.config.get('chain_file') is not None:
+                    chain_file = _get_chain(None)
+            else:
+                chain_file = _get_chain(chain_file)
+
         times = np.empty_like(self, dtype=float)
         if by_nuclide:
             activities = [None] * len(self)
@@ -149,7 +168,8 @@ def get_activity(
         # Evaluate activity for each depletion time
         for i, result in enumerate(self):
             times[i] = result.time[0]
-            activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume)
+            activities[i] = result.get_material(mat_id).get_activity(
+                units, by_nuclide, volume, chain_file=chain_file)
 
         return times, activities
 

openmc/material.py

@@ -1385,8 +1385,13 @@ def get_element_atom_densities(self, element: str | None = None) -> dict[str, fl
         return densities
 
 
-    def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
-                     volume: float | None = None) -> dict[str, float] | float:
+    def get_activity(
+        self,
+        units: str = 'Bq/cm3',
+        by_nuclide: bool = False,
+        volume: float | None = None,
+        chain_file=None
+    ) -> dict[str, float] | float:
         """Return the activity of the material or each nuclide within.
 
         .. versionadded:: 0.13.1
@@ -1405,13 +1410,22 @@ def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
             :attr:`Material.volume` attribute.
 
             .. versionadded:: 0.13.3
+        chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
+            Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
+            used. If ``None``, the chain specified by
+            ``openmc.config['chain_file']`` is used when available. If a path or
+            :class:`openmc.deplete.Chain` is given, that chain is used. For
+            ``None`` or an explicit chain, nuclides absent from the chain fall
+            back to ENDF/B-VIII.0 data.
+
+            .. versionadded:: 0.15.4
 
         Returns
         -------
         Union[dict, float]
-            If by_nuclide is True then a dictionary whose keys are nuclide
-            names and values are activity is returned. Otherwise the activity
-            of the material is returned as a float.
+            If by_nuclide is True then a dictionary whose keys are nuclide names
+            and values are activity is returned. Otherwise the activity of the
+            material is returned as a float.
         """
 
         cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'})
@@ -1440,7 +1454,8 @@ def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
 
         activity = {}
         for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items():
-            inv_seconds = openmc.data.decay_constant(nuclide)
+            inv_seconds = openmc.data.decay_constant(
+                nuclide, chain_file=chain_file)
             activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier
 
         return activity if by_nuclide else sum(activity.values())

tests/unit_tests/test_data_decay.py

@@ -155,16 +155,16 @@ def test_decay_photon_energy():
     with pytest.raises(DataError):
         openmc.data.decay_photon_energy('I135')
 
-    # 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')
-    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')
-    assert isinstance(src, openmc.stats.Tabular)
-    assert src.integral() == pytest.approx(2.076506258964966e-05)
+    # Temporarily Set chain file to simple chain
+    with openmc.config.patch('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')
+        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')
+        assert isinstance(src, openmc.stats.Tabular)
+        assert src.integral() == pytest.approx(2.076506258964966e-05)

tests/unit_tests/test_data_misc.py

@@ -7,6 +7,7 @@
 import numpy as np
 import pytest
 import openmc.data
+from openmc.deplete import Chain, Nuclide
 
 
 def test_data_library(tmpdir):
@@ -134,7 +135,8 @@ def test_zam():
     with pytest.raises(ValueError):
         openmc.data.zam('Am242-m1')
 
-def test_half_life():
+
+def test_half_life(tmp_path):
     assert openmc.data.half_life('H2') is None
     assert openmc.data.half_life('U235') == pytest.approx(2.22102e16)
     assert openmc.data.half_life('Am242') == pytest.approx(57672.0)
@@ -143,3 +145,32 @@ def test_half_life():
     assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16)
     assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0)
     assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0)
+
+    # Create minimal chain with H3 and Am242 to test half-life and decay
+    # constant retrieval from chain file
+    chain = Chain()
+    h3 = Nuclide("H3")
+    h3.half_life = 1.0
+    chain.add_nuclide(h3)
+    am242 = Nuclide("Am242")
+    chain.add_nuclide(am242)
+
+    assert openmc.data.half_life('H3', chain_file=chain) == 1.0
+    assert openmc.data.decay_constant('H3', chain_file=chain) == pytest.approx(log(2.0))
+
+    # Nuclides that are present but stable in the chain should not fall back to
+    # ENDF/B-VIII.0 data.
+    assert openmc.data.half_life('Am242', chain_file=chain) is None
+    assert openmc.data.decay_constant('Am242', chain_file=chain) == 0.0
+
+    # Nuclides missing from the chain fall back to ENDF/B-VIII.0 data.
+    assert openmc.data.half_life('U235', chain_file=chain) == pytest.approx(2.22102e16)
+
+    chain_path = tmp_path / "chain.xml"
+    chain.export_to_xml(chain_path)
+    assert openmc.data.half_life('H3', chain_file=chain_path) == 1.0
+
+    endf_h3 = openmc.data.half_life('H3')
+    with openmc.config.patch('chain_file', chain_path):
+        assert openmc.data.half_life('H3', chain_file=None) == 1.0
+        assert openmc.data.half_life('H3', chain_file=False) == endf_h3

tests/unit_tests/test_deplete_resultslist.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pytest
 
+import openmc
 import openmc.deplete
 
 
@@ -38,6 +39,36 @@ def test_get_activity(res):
     np.testing.assert_allclose(a_xe135, a_xe135_ref)
 
 
+def test_get_activity_chain_file(res, tmp_path):
+    """Tests evaluating activity with chain half-life data"""
+    _, a_endf = res.get_activity("1", by_nuclide=True, chain_file=False)
+    xe135_endf = np.array([a["Xe135"] for a in a_endf])
+
+    chain = openmc.deplete.Chain()
+    xe135 = openmc.deplete.Nuclide("Xe135")
+    xe135.half_life = openmc.data.half_life("Xe135") / 2.0
+    chain.add_nuclide(xe135)
+
+    t_chain, a_chain = res.get_activity("1", by_nuclide=True, chain_file=chain)
+    xe135_chain = np.array([a["Xe135"] for a in a_chain])
+
+    t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
+    np.testing.assert_allclose(t_chain, t_ref)
+    np.testing.assert_allclose(xe135_chain, 2.0 * xe135_endf)
+
+    chain_path = tmp_path / "chain.xml"
+    chain.export_to_xml(chain_path)
+    with openmc.config.patch('chain_file', chain_path):
+        _, a_config = res.get_activity("1", by_nuclide=True)
+    xe135_config = np.array([a["Xe135"] for a in a_config])
+    np.testing.assert_allclose(xe135_config, xe135_chain)
+
+    stable_chain = openmc.deplete.Chain()
+    stable_chain.add_nuclide(openmc.deplete.Nuclide("Xe135"))
+    _, a_stable = res.get_activity("1", by_nuclide=True, chain_file=stable_chain)
+    assert all(a["Xe135"] == 0.0 for a in a_stable)
+
+
 def test_get_atoms(res):
     """Tests evaluating single nuclide concentration."""
     t, n = res.get_atoms("1", "Xe135")