Fix CDR calculation with trace elements with Z>100

openmc/material.py

@@ -41,6 +41,10 @@
 
 _BECQUEREL_PER_CURIE = 3.7e10
 
+# Minimum mass fraction of nuclides without photon attenuation data that
+# results in a warning from Material.get_photon_contact_dose_rate()
+_MIN_ATTENUATION_MASS_FRACTION = 1e-6
+
 NuclideTuple = namedtuple('NuclideTuple', ['name', 'percent', 'percent_type'])
 
 
@@ -473,6 +477,9 @@ def get_photon_contact_dose_rate(
         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.
+        Photon attenuation data is only tabulated up to Z=100; nuclides with a
+        higher atomic number are neglected when building the material
+        attenuation coefficient.
 
         Returns
         -------
@@ -493,17 +500,48 @@ def get_photon_contact_dose_rate(
             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
+        # attenuation coefficients µ_i/ρ_i [cm²/g] per nuclide. Attenuation
+        # data is only tabulated up to Z=100, so nuclides beyond that -- which
+        # show up in trace quantities after depletion -- are left out of
+        # µ_material(E) instead of aborting the calculation.
         nuc_attenuation = []
+        missing_data = {}
+        total_rho = 0.0
         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
             )
+            total_rho += rho_i
+
+            Z = openmc.data.zam(nuc)[0]
+            try:
+                mu_over_rho = openmc.data.mass_attenuation_coefficient(Z)
+            except ValueError:
+                missing_data[nuc] = rho_i
+                continue
+
             nuc_attenuation.append((rho_i, mu_over_rho))
 
+        if not nuc_attenuation:
+            raise ValueError(
+                "No photon attenuation data is available for any nuclide in "
+                f"material ID={self.id}; cannot compute the contact dose rate."
+            )
+
+        # Only warn about neglected nuclides if they are more than a trace
+        if missing_data and total_rho > 0.0:
+            missing_frac = sum(missing_data.values()) / total_rho
+            if missing_frac > _MIN_ATTENUATION_MASS_FRACTION:
+                warnings.warn(
+                    'No photon attenuation data available for '
+                    f'{", ".join(sorted(missing_data))} in material '
+                    f'ID={self.id}. These nuclides make up a mass fraction of '
+                    f'{missing_frac:.3e} and are neglected in the material '
+                    'attenuation coefficient.',
+                    stacklevel=2,
+                )
+
         # Build union energy grid across all nuclides
         mu_e_vals = reduce(np.union1d, [t.x for _, t in nuc_attenuation])
 

tests/unit_tests/test_material.py

@@ -1,5 +1,6 @@
 from collections import defaultdict
 from pathlib import Path
+import warnings
 
 import pytest
 
@@ -922,3 +923,44 @@ def test_get_photon_contact_dose_rate():
         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)
+
+
+def test_get_photon_contact_dose_rate_missing_attenuation():
+    # Set chain file for testing
+    openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    m_i135 = openmc.Material()
+    m_i135.add_nuclide('I135', 1.0)
+    m_i135.set_density('atom/b-cm', 1.0)
+    reference = m_i135.get_photon_contact_dose_rate('absorbed-air')
+
+    # Depletion produces trace amounts of nuclides above Z=100, for which no
+    # photon attenuation data exists. They should be skipped rather than
+    # aborting the calculation, and are far too dilute to change the result.
+    m_trace = openmc.Material()
+    m_trace.add_nuclide('I135', 1.0)
+    for nuclide, percent in [('Rf265', 2.6e-36), ('Sg269', 1.0e-36),
+                             ('Hs273', 1.0e-36)]:
+        m_trace.add_nuclide(nuclide, percent)
+    m_trace.set_density('atom/b-cm', 1.0)
+
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter('always')
+        cdr = m_trace.get_photon_contact_dose_rate('absorbed-air')
+    assert not caught
+    assert cdr == pytest.approx(reference, rel=1e-10)
+
+    # A non-trace amount of such nuclides is worth telling the user about
+    m_bulk = openmc.Material()
+    m_bulk.add_nuclide('I135', 0.5)
+    m_bulk.add_nuclide('Rf265', 0.5)
+    m_bulk.set_density('atom/b-cm', 1.0)
+    with pytest.warns(UserWarning, match='Rf265'):
+        m_bulk.get_photon_contact_dose_rate('absorbed-air')
+
+    # If nothing in the material has attenuation data, give a clear error
+    m_none = openmc.Material()
+    m_none.add_nuclide('Rf265', 1.0)
+    m_none.set_density('atom/b-cm', 1.0e-10)
+    with pytest.raises(ValueError, match='No photon attenuation data'):
+        m_none.get_photon_contact_dose_rate('absorbed-air')