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

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

Author: EdenRochmanSharabi

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 @@ def get_activity(
         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 @@ 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 +155,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 +171,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

@@ -8,7 +8,7 @@
 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.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 @@ 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: 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 @@ 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 +1457,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