Add chain parameter to Material.get_activity for half-life data (#3957)

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

openmc/data/data.py

@@ -1,14 +1,23 @@
+from __future__ import annotations
+
 import itertools
 import json
 import os
 import re
 from pathlib import Path
 from math import sqrt, log
+from typing import TYPE_CHECKING, Literal
 from warnings import warn
 
 from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL,
                        EV_PER_MEV, K_BOLTZMANN, gnds_name, zam)
 
+import openmc
+from openmc.checkvalue import PathLike
+
+if TYPE_CHECKING:
+    from openmc.deplete import Chain
+
 gnds_name.__module__ = __name__
 zam.__module__ = __name__
 
@@ -296,25 +305,47 @@ def atomic_weight(element):
         raise ValueError(f"No naturally-occurring isotopes for element '{element}'.")
 
 
-def half_life(isotope):
+def half_life(
+    isotope: str,
+    chain_file: Literal[False] | None | PathLike | Chain = False
+) -> float | None:
     """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 +355,10 @@ def half_life(isotope):
     return _HALF_LIFE.get(isotope.lower())
 
 
-def decay_constant(isotope):
+def decay_constant(
+    isotope: str,
+    chain_file: Literal[False] | None | PathLike | Chain = False
+) -> float:
     """Return decay constant of isotope in [s^-1]
 
     Decay constants are based on half-life values from the
@@ -333,10 +367,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 +392,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 +540,3 @@ def isotopes(element: str) -> list[tuple[str, float]]:
             result.append(kv)
 
     return result
-
-

openmc/deplete/results.py

@@ -1,12 +1,17 @@
+from __future__ import annotations
+
 import numbers
 import bisect
 import math
 from collections.abc import Iterable
+from typing import Literal
 from warnings import warn
 
 import h5py
 import numpy as np
 
+import openmc
+from .chain import Chain, _get_chain
 from .stepresult import StepResult, VERSION_RESULTS
 import openmc.checkvalue as cv
 from openmc.data import atomic_mass, AVOGADRO
@@ -103,7 +108,8 @@ class Results(list):
         mat: Material | str,
         units: str = "Bq/cm3",
         by_nuclide: bool = False,
-        volume: float | None = None
+        volume: float | None = None,
+        chain_file: Literal[False] | None | PathLike | Chain = None
     ) -> tuple[np.ndarray, np.ndarray | list[dict]]:
         """Get activity of material over time.
 
@@ -115,22 +121,31 @@ class Results(list):
             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 +155,13 @@ class Results(list):
         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 +171,8 @@ class Results(list):
         # 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

@@ -8,7 +8,7 @@ from pathlib import Path
 import re
 import sys
 import tempfile
-from typing import Sequence, Dict
+from typing import TYPE_CHECKING, Literal, Sequence, Dict
 import warnings
 
 import lxml.etree as ET
@@ -28,6 +28,9 @@ from openmc.data.data import _get_element_symbol, JOULE_PER_EV
 from openmc.data.function import Tabulated1D
 from openmc.data import mass_energy_absorption_coefficient, dose_coefficients
 
+if TYPE_CHECKING:
+    from openmc.deplete import Chain
+
 
 # Units for density supported by OpenMC
 DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
@@ -1385,8 +1388,13 @@ class Material(IDManagerMixin):
         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: Literal[False] | None | PathLike | Chain = None
+    ) -> dict[str, float] | float:
         """Return the activity of the material or each nuclide within.
 
         .. versionadded:: 0.13.1
@@ -1405,13 +1413,22 @@ class Material(IDManagerMixin):
             :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 +1457,8 @@ class Material(IDManagerMixin):
 
         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"
+    # 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 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)
+        # 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 @@ from pathlib import Path
 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 @@ from pathlib import Path
 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")

tests/unit_tests/test_material.py

@@ -7,7 +7,7 @@ import numpy as np
 
 import openmc
 from openmc.data import decay_photon_energy
-from openmc.deplete import Chain
+from openmc.deplete import Chain, Nuclide
 import openmc.examples
 import openmc.model
 import openmc.stats
@@ -614,6 +614,35 @@ def test_get_activity():
     assert m4.get_activity(units='Ci/m3') == pytest.approx(ci/m3)
 
 
+def test_get_activity_chain_file(tmp_path):
+    m = openmc.Material()
+    m.add_nuclide("H3", 1.0)
+    m.set_density('g/cm3', 1.0)
+
+    chain = Chain()
+    h3 = Nuclide("H3")
+    h3.half_life = 1.0
+    chain.add_nuclide(h3)
+
+    atoms_per_bcm = m.get_nuclide_atom_densities()["H3"]
+    expected = np.log(2.0) * 1e24 * atoms_per_bcm
+
+    assert m.get_activity(chain_file=chain) == pytest.approx(expected)
+
+    chain_path = tmp_path / "chain.xml"
+    chain.export_to_xml(chain_path)
+    assert m.get_activity(chain_file=chain_path) == pytest.approx(expected)
+
+    endf_activity = m.get_activity(chain_file=False)
+    with openmc.config.patch('chain_file', chain_path):
+        assert m.get_activity() == pytest.approx(expected)
+        assert m.get_activity(chain_file=False) == pytest.approx(endf_activity)
+
+    stable_chain = Chain()
+    stable_chain.add_nuclide(Nuclide("H3"))
+    assert m.get_activity(chain_file=stable_chain) == 0.0
+
+
 def test_get_decay_heat():
     # Set chain file for testing
     openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'