Add method on Material for computing photon contact dose rate (#3700)

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

docs/source/pythonapi/data.rst

@@ -71,6 +71,8 @@ Core Functions
     isotopes
     kalbach_slope
     linearize
+    mass_attenuation_coefficient
+    mass_energy_absorption_coefficient
     thin
     water_density
     zam

openmc/data/__init__.py

@@ -35,4 +35,6 @@ from .grid import *
 from .function import *
 from .vectfit import *
 
-from .effective_dose.dose import dose_coefficients
+from .dose.dose import dose_coefficients
+from .dose.mass_attenuation import \
+    mass_energy_absorption_coefficient, mass_attenuation_coefficient

openmc/data/dose/mass_attenuation.py

@@ -0,0 +1,153 @@
+from pathlib import Path
+
+import numpy as np
+import h5py
+
+import openmc.checkvalue as cv
+from openmc.data import EV_PER_MEV
+from ..data import ATOMIC_NUMBER
+from ..function import Tabulated1D
+
+# Embedded NIST-126 data
+# Air (Dry Near Sea Level) — NIST Standard Reference Database 126 Table 4 (doi: 10.18434/T4D01F)
+# Columns: Energy (MeV), μ_en/ρ (cm^2/g)
+_NIST126_AIR = np.array([
+    [1.00000e-03, 3.599e03],
+    [1.50000e-03, 1.188e03],
+    [2.00000e-03, 5.262e02],
+    [3.00000e-03, 1.614e02],
+    [3.20290e-03, 1.330e02],
+    [3.20290e-03, 1.460e02],
+    [4.00000e-03, 7.636e01],
+    [5.00000e-03, 3.931e01],
+    [6.00000e-03, 2.270e01],
+    [8.00000e-03, 9.446e00],
+    [1.00000e-02, 4.742e00],
+    [1.50000e-02, 1.334e00],
+    [2.00000e-02, 5.389e-01],
+    [3.00000e-02, 1.537e-01],
+    [4.00000e-02, 6.833e-02],
+    [5.00000e-02, 4.098e-02],
+    [6.00000e-02, 3.041e-02],
+    [8.00000e-02, 2.407e-02],
+    [1.00000e-01, 2.325e-02],
+    [1.50000e-01, 2.496e-02],
+    [2.00000e-01, 2.672e-02],
+    [3.00000e-01, 2.872e-02],
+    [4.00000e-01, 2.949e-02],
+    [5.00000e-01, 2.966e-02],
+    [6.00000e-01, 2.953e-02],
+    [8.00000e-01, 2.882e-02],
+    [1.00000e00, 2.789e-02],
+    [1.25000e00, 2.666e-02],
+    [1.50000e00, 2.547e-02],
+    [2.00000e00, 2.345e-02],
+    [3.00000e00, 2.057e-02],
+    [4.00000e00, 1.870e-02],
+    [5.00000e00, 1.740e-02],
+    [6.00000e00, 1.647e-02],
+    [8.00000e00, 1.525e-02],
+    [1.00000e01, 1.450e-02],
+    [1.50000e01, 1.353e-02],
+    [2.00000e01, 1.311e-02],
+])
+
+# Registry of embedded tables: (data_source, material) -> ndarray
+# Table shape: (N, 2) with columns [Energy (MeV), μen/ρ (cm^2/g)]
+_MUEN_TABLES = {
+    ("nist126", "air"): _NIST126_AIR,
+}
+
+
+def mass_energy_absorption_coefficient(
+    material: str, data_source: str = "nist126"
+) -> Tabulated1D:
+    """Return the mass energy-absorption coefficient as a function of energy.
+
+    The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is
+    defined as the fraction of incident photon energy absorbed in a material per
+    unit mass less the energy carried away by scattered photons. It is obtained
+    from `NIST Standard Reference Database 126
+    <https://doi.org/10.18434/T4D01F>`_: X-Ray Mass Attenuation Coefficients.
+
+    Parameters
+    ----------
+    material : {'air'}
+        Material compound for which to load coefficients.
+    data_source : {'nist126'}
+        Source library.
+
+    Returns
+    -------
+    Tabulated1D
+        Mass energy-absorption coefficient [cm^2/g] as a function of photon
+        energy [eV], using log-log interpolation.
+
+    """
+    cv.check_value("material", material, {"air"})
+    cv.check_value("data_source", data_source, {"nist126"})
+
+    key = (data_source, material)
+    if key not in _MUEN_TABLES:
+        available = sorted({m for (ds, m) in _MUEN_TABLES.keys() if ds == data_source})
+        raise ValueError(
+            f"No mass energy-absorption data for '{material}' in data source "
+            f"'{data_source}'. Available materials: {available}"
+        )
+
+    data = _MUEN_TABLES[key]
+    energy = data[:, 0].copy() * EV_PER_MEV  # MeV -> eV
+    mu_en_coeffs = data[:, 1].copy()
+    return Tabulated1D(energy, mu_en_coeffs,
+                       breakpoints=[len(energy)], interpolation=[5])
+
+
+# Used in mass_attenuation_coefficient function as a cache.
+# Maps atomic number Z (int) -> Tabulated1D of (mu/rho) [cm^2/g] vs E [eV]
+_MASS_ATTENUATION: dict[int, object] = {}
+
+
+def mass_attenuation_coefficient(element):
+    """Return the photon mass attenuation coefficient as a function of energy.
+
+    The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is
+    defined as the fraction of incident photon energy absorbed in a material per
+    unit mass. Values for each element are obtained from `NIST Standard
+    Reference Database 8 <https://doi.org/10.18434/T48G6X>`_: XCOM Photon Cross
+    Sections Database.
+
+    Parameters
+    ----------
+    element : str or int
+        Element symbol (e.g., 'Fe') or atomic number (e.g., 26).
+
+    Returns
+    -------
+    Tabulated1D
+        Mass attenuation coefficient [cm^2/g] as a function of photon energy
+        [eV], using log-log interpolation.
+
+    """
+    if not _MASS_ATTENUATION:
+        data_file = Path(__file__).with_name('mass_attenuation.h5')
+        with h5py.File(data_file, 'r') as f:
+            for key, dataset in f.items():
+                energies, mu_rho = dataset[()]  # shape (2, N)
+                _MASS_ATTENUATION[int(key)] = Tabulated1D(
+                    energies, mu_rho,
+                    breakpoints=[len(energies)],
+                    interpolation=[5]  # log-log
+                )
+
+    # Resolve element argument to atomic number
+    if isinstance(element, str):
+        if element not in ATOMIC_NUMBER:
+            raise ValueError(f"'{element}' is not a recognized element symbol")
+        Z = ATOMIC_NUMBER[element]
+    else:
+        Z = int(element)
+
+    if Z not in _MASS_ATTENUATION:
+        raise ValueError(f"No mass attenuation data available for Z={Z}")
+
+    return _MASS_ATTENUATION[Z]

openmc/material.py

@@ -2,6 +2,7 @@ from __future__ import annotations
 from collections import defaultdict, namedtuple, Counter
 from collections.abc import Iterable
 from copy import deepcopy
+from functools import reduce
 from numbers import Real
 from pathlib import Path
 import re
@@ -22,8 +23,10 @@ from .mixin import IDManagerMixin
 from .utility_funcs import input_path
 from . import waste
 from openmc.checkvalue import PathLike
-from openmc.stats import Univariate, Discrete, Mixture
-from openmc.data.data import _get_element_symbol
+from openmc.stats import Univariate, Discrete, Mixture, Tabular
+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
 
 
 # Units for density supported by OpenMC
@@ -409,6 +412,190 @@ class Material(IDManagerMixin):
 
         return combined
 
+    def get_photon_contact_dose_rate(
+        self,
+        dose_quantity: str = "absorbed-air",
+        build_up: float = 2.0,
+        by_nuclide: bool = False
+    ) -> float | dict[str, float]:
+        """Compute the photon contact dose rate (CDR) produced by radioactive decay
+        of the material.
+
+        The contact dose rate is calculated from decay photon energy spectra for
+        each nuclide in the material, combined with photon mass attenuation data
+        for the material and the appropriate response function for the dose quantity.
+        A slab-geometry approximation and a photon build-up factor are used.
+
+        Absorbed-air dose:
+            The approach follows the FISPACT-II manual (UKAEA-CCFE-RE(21)02 - May 2021).
+            Appendix C.7.1.
+            This method integrates over the photon energy:
+
+                (B/2) * (mu_en_air(E) / mu_material(E)) * E * S(E)
+
+        Effective dose:
+            The approach uses ICRP-116 effective dose coefficients to convert the photon
+            fluence due to decay photons to effective dose.
+            This method integrates over the photon energy:
+
+                (B/2) * (h_e(E) / mu_material(E)) * S(E)
+
+        where:
+            - mu_en_air(E) is the air mass energy-absorption coefficient,
+            - mu_material(E) is the photon mass attenuation coefficient of the material,
+            - S(E) is the photon emission spectrum per atom,
+            - h_e(E) is the ICRP-116 effective dose coefficient,
+            - B is the build-up factor,
+            - E is the photon energy.
+
+        Parameters
+        ----------
+        dose_quantity : {'absorbed-air', 'effective'}, optional
+            Specifies the dose quantity to be calculated.
+            The only supported options are 'absorbed-air' which implements the methodology
+            from FISPACT-II, and 'effective' which uses ICRP-116 effective dose coefficients.
+        build_up : float, optional. The default value is 2.0 as suggested in the FISPACT-II
+            manual.
+        by_nuclide : bool, optional
+            Specifies if the cdr should be returned for the material as a
+            whole or per nuclide. Default is False.
+
+        Limitations
+        ----------
+        This method does not implement correction from Bremsstrahlung particles which can be
+        relevant at close distances.
+        In addition, it computes the gamma contact dose rate only for the unstable nuclides
+        for which the radiation source specification is present in the chain file.
+
+        Returns
+        -------
+        cdr : float or dict[str, float]
+            Contact Dose Rate due to decay photons.
+            'absorbed-air': returns the absorbed dose in air [Gy/hr].
+            'effective': returns the effective dose [Sv/hr].
+        """
+
+        cv.check_type("by_nuclide", by_nuclide, bool)
+        cv.check_type("dose_quantity", dose_quantity, str)
+        cv.check_value("dose_quantity", dose_quantity, {'absorbed-air', 'effective'})
+        cv.check_type("build_up", build_up, Real)
+        cv.check_greater_than("build_up", build_up, 0.0)
+
+        nuc_densities = self.get_nuclide_atom_densities()
+        if not nuc_densities:
+            raise ValueError("Material has no nuclides; cannot compute mass attenuation")
+
+        # Collect partial mass densities ρ_i [g/cm³] and elemental mass
+        # attenuation coefficients µ_i/ρ_i [cm²/g] per nuclide
+        nuc_attenuation = []
+        for nuc, atom_density_bcm in nuc_densities.items():
+            Z = openmc.data.zam(nuc)[0]
+            mu_over_rho = openmc.data.mass_attenuation_coefficient(Z)
+            rho_i = (
+                atom_density_bcm * 1.0e24
+                * openmc.data.atomic_mass(nuc) / openmc.data.AVOGADRO
+            )
+            nuc_attenuation.append((rho_i, mu_over_rho))
+
+        # Build union energy grid across all nuclides
+        mu_e_vals = reduce(np.union1d, [t.x for _, t in nuc_attenuation])
+
+        # Build the material linear attenuation coefficient µ_material(E) [cm⁻¹]
+        # as the sum of ρ_i * (µ_i/ρ_i)(E) over all nuclides
+        mu_material_vals = np.zeros(len(mu_e_vals))
+        for rho_i, mu_over_rho in nuc_attenuation:
+            mu_material_vals += rho_i * mu_over_rho(mu_e_vals)
+        mu_material = Tabulated1D(
+            mu_e_vals, mu_material_vals, breakpoints=[len(mu_e_vals)], interpolation=[5])
+
+        # CDR computation
+        cdr = {}
+
+        geometry_factor_slab = 0.5
+
+        # ancillary conversion factors for clarity
+        seconds_per_hour = 3600.0
+        grams_per_kg = 1000.0
+        sv_per_psv = 1e-12
+
+        if dose_quantity == 'absorbed-air':
+            # mu_en/rho for air [cm²/g] as a function of energy [eV]
+            response_f = mass_energy_absorption_coefficient("air", data_source="nist126")
+
+            # Factor to convert [eV cm²/(b g s)] to [Gy/h]
+            multiplier = (build_up * geometry_factor_slab * seconds_per_hour
+                          * grams_per_kg * 1e24 * JOULE_PER_EV)
+
+        elif dose_quantity == 'effective':
+            # effective dose as a function of photon fluence [pSv cm²]
+            response_f_x, response_f_y = dose_coefficients(
+                "photon", geometry='AP', data_source='icrp116')
+            response_f = Tabulated1D(response_f_x, response_f_y, breakpoints=[
+                len(response_f_x)], interpolation=[5])
+
+            # Convert [pSv cm²/(b-s)] to [Sv/h]
+            multiplier = (build_up * geometry_factor_slab * seconds_per_hour
+                          * 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)
+
+            # nuclides with no contribution
+            if photon_source_per_atom is None or nuc_atoms_per_bcm <= 0.0:
+                cdr[nuc] = 0.0
+                continue
+
+            if not isinstance(photon_source_per_atom, (Discrete, Tabular)):
+                raise ValueError(
+                    f"Unknown decay photon energy data type for nuclide {nuc}"
+                    f"value returned: {type(photon_source_per_atom)}"
+                )
+
+            e_vals = photon_source_per_atom.x
+            p_vals = photon_source_per_atom.p
+
+            # Construct list of energies from (photon source, response function,
+            # mu_en_air) for clipping to common energy range
+            e_lists = [e_vals, response_f.x, mu_e_vals]
+
+            # clip distributions for values outside the tabulated values
+            left_bound = max(a.min() for a in e_lists)
+            right_bound = min(a.max() for a in e_lists)
+
+            mask = (e_vals >= left_bound) & (e_vals <= right_bound)
+            e_vals = e_vals[mask]
+            p_vals = p_vals[mask]
+
+            if isinstance(photon_source_per_atom, Tabular):
+                # limit the computation to the tabulated mu_en_air range
+                e_union = reduce(np.union1d, e_lists)
+                e_union = e_union[(e_union >= left_bound) & (e_union <= right_bound)]
+                if len(e_union) < 2:
+                    raise ValueError("Not enough overlapping energy points to compute CDR")
+
+                # Histogram interpolation: each new point inherits the value of
+                # the nearest original point to its left
+                p_vals = p_vals[np.searchsorted(e_vals, e_union, side='right') - 1]
+                e_vals = e_union
+
+            mu_vals = mu_material(e_vals)
+            if dose_quantity == 'absorbed-air':
+                # Compute (µ_en_air(E) / µ_material(E)) * E * S(E)
+                integrand = (response_f(e_vals) / mu_vals) * p_vals * e_vals
+            elif dose_quantity == 'effective':
+                # Compute (h_e(E) / µ_material(E)) * S(E)
+                integrand = (response_f(e_vals) / mu_vals) * p_vals
+
+            if isinstance(photon_source_per_atom, Discrete):
+                cdr_nuc = np.sum(integrand)
+            elif isinstance(photon_source_per_atom, Tabular):
+                cdr_nuc = np.trapezoid(integrand, e_vals)
+
+            # Compute air-absorbed dose [Gy/h] or effective dose [Sv/h]
+            cdr[nuc] = float(cdr_nuc * nuc_atoms_per_bcm * multiplier)
+
+        return cdr if by_nuclide else sum(cdr.values())
+
     @classmethod
     def from_hdf5(cls, group: h5py.Group) -> Material:
         """Create material from HDF5 group

pyproject.toml

@@ -69,7 +69,7 @@ include = ['openmc*']
 exclude = ['tests*']
 
 [tool.setuptools.package-data]
-"openmc.data.effective_dose" = ["**/*.txt"]
+"openmc.data.dose" = ["**/*.txt"]
 "openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"]
 "openmc.lib" = ["libopenmc.dylib", "libopenmc.so"]
 

tests/unit_tests/test_data_mass_attenuation.py

@@ -0,0 +1,53 @@
+from pytest import approx, raises
+
+from openmc.data import mass_energy_absorption_coefficient, mass_attenuation_coefficient
+from openmc.data.function import Tabulated1D
+
+
+def test_mass_attenuation_type():
+    mu = mass_attenuation_coefficient(26)  # Fe
+    assert isinstance(mu, Tabulated1D)
+
+
+def test_mass_attenuation_spot_values():
+    # Spot checks for Fe (Z=26) against NIST data: first/last tabulated points
+    # and a mid-range value at 1 MeV
+    mu = mass_attenuation_coefficient(26)
+    assert mu(1e3) == approx(9085.0)
+    assert mu(1e6) == approx(0.05995)
+    assert mu(2e7) == approx(0.03224)
+
+
+def test_mass_attenuation_caching():
+    # Repeated calls with the same Z should return the identical object
+    mu1 = mass_attenuation_coefficient(26)
+    mu2 = mass_attenuation_coefficient('Fe')
+    assert mu1 is mu2
+
+
+def test_mass_attenuation_invalid_z():
+    with raises(ValueError, match="Z=0"):
+        mass_attenuation_coefficient(0)
+    with raises(ValueError, match="Z=200"):
+        mass_attenuation_coefficient(200)
+
+
+def test_mass_energy_absorption_type():
+    # Spot checks on values from NIST tables
+    mu_en = mass_energy_absorption_coefficient("air")
+    assert isinstance(mu_en, Tabulated1D)
+
+
+def test_mass_energy_absorption_spot_values():
+    mu_en = mass_energy_absorption_coefficient("air")
+    assert mu_en(1e3) == approx(3.599e3)
+    assert mu_en(10.e3) == approx(4.742)
+    assert mu_en(2e7) == approx(1.311e-2)
+
+
+def test_mass_energy_absorption_invalid():
+    # Invalid material/data_source should raise an exception
+    with raises(ValueError):
+        mass_energy_absorption_coefficient("pasta")
+    with raises(ValueError):
+        mass_energy_absorption_coefficient("air", data_source="nist000")

tests/unit_tests/test_material.py

@@ -826,3 +826,58 @@ def test_material_from_constructor():
     assert mat2.density == 1e-7
     assert mat2.density_units == "g/cm3"
     assert mat2.nuclides == []
+
+
+def test_get_photon_contact_dose_rate():
+    # Set chain file for testing
+    openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    # A purely stable material (Fe) should give zero dose
+    m_stable = openmc.Material()
+    m_stable.add_element('Fe', 1.0)
+    m_stable.set_density('g/cm3', 7.87)
+    assert m_stable.get_photon_contact_dose_rate('absorbed-air') == 0.0
+    assert m_stable.get_photon_contact_dose_rate('effective') == 0.0
+
+    # I135 has a Discrete photon source (lines)
+    m_i135 = openmc.Material()
+    m_i135.add_nuclide('I135', 1.0)
+    m_i135.set_density('atom/b-cm', 1.0)
+
+    cdr_abs = m_i135.get_photon_contact_dose_rate('absorbed-air')
+    cdr_eff = m_i135.get_photon_contact_dose_rate('effective')
+    assert cdr_abs == pytest.approx(6.091547e10, rel=1e-4)   # [Gy/h]
+    assert cdr_eff == pytest.approx(6.102167e10, rel=1e-4)   # [Sv/h]
+
+    # Xe135 has a Tabular photon source (continuous distribution)
+    m_xe135 = openmc.Material()
+    m_xe135.add_nuclide('Xe135', 1.0)
+    m_xe135.set_density('atom/b-cm', 1.0)
+
+    cdr_xe_abs = m_xe135.get_photon_contact_dose_rate('absorbed-air')
+    cdr_xe_eff = m_xe135.get_photon_contact_dose_rate('effective')
+    assert cdr_xe_abs == pytest.approx(7.886077e8, rel=1e-4)  # [Gy/h]
+    assert cdr_xe_eff == pytest.approx(9.488298e8, rel=1e-4)  # [Sv/h]
+
+    # by_nuclide=True should return a dict whose values sum to the total
+    cdr_by_nuc = m_i135.get_photon_contact_dose_rate('absorbed-air', by_nuclide=True)
+    assert isinstance(cdr_by_nuc, dict)
+    assert 'I135' in cdr_by_nuc
+    assert sum(cdr_by_nuc.values()) == pytest.approx(cdr_abs)
+
+    # For a mixed material the sum over nuclides must equal the total
+    m_mix = openmc.Material()
+    m_mix.add_nuclide('I135', 0.5)
+    m_mix.add_nuclide('Xe135', 0.5)
+    m_mix.set_density('atom/b-cm', 1.0)
+    cdr_mix_total = m_mix.get_photon_contact_dose_rate('absorbed-air')
+    cdr_mix_nuc = m_mix.get_photon_contact_dose_rate('absorbed-air', by_nuclide=True)
+    assert sum(cdr_mix_nuc.values()) == pytest.approx(cdr_mix_total)
+
+    # Input validation
+    with pytest.raises(ValueError):
+        m_i135.get_photon_contact_dose_rate('invalid-quantity')
+    with pytest.raises(TypeError):
+        m_i135.get_photon_contact_dose_rate('absorbed-air', build_up='two')
+    with pytest.raises(ValueError):
+        m_i135.get_photon_contact_dose_rate('absorbed-air', build_up=-1.0)

tests/unit_tests/test_mesh.py

@@ -612,6 +612,7 @@ def test_mesh_get_homogenized_materials():
 
 @pytest.fixture
 def sphere_model():
+    openmc.reset_auto_ids()
     # Model with three materials separated by planes x=0 and z=0
     mats = []
     for i in range(3):