diff --git a/openmc/data/data.py b/openmc/data/data.py index c22e54e7dc8..84e1669ae81 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -9,6 +9,8 @@ from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL, EV_PER_MEV, K_BOLTZMANN, gnds_name, zam) +import openmc + gnds_name.__module__ = __name__ zam.__module__ = __name__ @@ -296,25 +298,44 @@ def atomic_weight(element): raise ValueError(f"No naturally-occurring isotopes for element '{element}'.") -def half_life(isotope): +def half_life(isotope, chain_file=False): """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 - `_. + By default, half-life values are from the `ENDF/B-VIII.0 decay sublibrary + `_. 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 +345,7 @@ def half_life(isotope): return _HALF_LIFE.get(isotope.lower()) -def decay_constant(isotope): +def decay_constant(isotope, chain_file=False): """Return decay constant of isotope in [s^-1] Decay constants are based on half-life values from the @@ -333,10 +354,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 +379,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 +527,3 @@ def isotopes(element: str) -> list[tuple[str, float]]: result.append(kv) return result - - diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index adb0d3dbc04..7099b960959 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -7,6 +7,8 @@ import h5py import numpy as np +import openmc +from .chain import _get_chain from .stepresult import StepResult, VERSION_RESULTS import openmc.checkvalue as cv from openmc.data import atomic_mass, AVOGADRO @@ -103,7 +105,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=None ) -> tuple[np.ndarray, np.ndarray | list[dict]]: """Get activity of material over time. @@ -115,22 +118,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 +152,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 +168,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 diff --git a/openmc/material.py b/openmc/material.py index eded9c3faa6..f6313423a16 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1385,8 +1385,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=None + ) -> dict[str, float] | float: """Return the activity of the material or each nuclide within. .. versionadded:: 0.13.1 @@ -1405,13 +1410,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 +1454,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()) diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index a48c553255e..f81f857e0cd 100644 --- a/tests/unit_tests/test_data_decay.py +++ b/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) diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index 14db689130e..320d09eda5a 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/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 diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 39c532c549a..bc1a078b05a 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -6,6 +6,7 @@ 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") diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 89dfc03ddd8..12d3633767a 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -7,7 +7,7 @@ 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'