diff --git a/openmc/material.py b/openmc/material.py index ece4cba2b1b..0e18552ce55 100644 --- a/openmc/material.py +++ b/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]) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 9cb8405f1e8..d8e7f6f9003 100644 --- a/tests/unit_tests/test_material.py +++ b/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')