diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index 74413e7b610..e00a895ddf8 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -60,8 +60,8 @@ attributes: ```` Element -------------------- -The ```` element represents photon and electron sources associated with -the decay of a nuclide and contains information to construct an +The ```` element represents decay particle sources associated with the +decay of a nuclide and contains information to construct an :class:`openmc.stats.Univariate` object that represents this emission as an energy distribution. This element has the following attributes: @@ -69,7 +69,8 @@ energy distribution. This element has the following attributes: The type of :class:`openmc.stats.Univariate` source term. :particle: - The type of particle emitted, e.g., 'photon' or 'electron' + The type of particle emitted, e.g., ``photon``, ``electron``, + ``positron``, ``alpha``, ``neutron``, ``proton``, or ``fragment`` :parameters: The parameters of the source term, e.g., for a diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 9d47430f75f..170d4a6fd7b 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -64,6 +64,7 @@ Core Functions combine_distributions decay_constant decay_energy + decay_particle_energy decay_photon_energy dose_coefficients gnds_name diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index 398680e7464..cdd504fafb5 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -39,9 +39,13 @@ Generally, it involves the following steps: with the timestep index returns a :class:`~openmc.deplete.StepResult` object, which itself has a :meth:`~openmc.deplete.StepResult.get_material` method returning an activated material. -4. The :meth:`openmc.Material.get_decay_photon_energy` method is used to obtain - the energy spectrum of the decay photon source. The integral of the spectrum - also indicates the intensity of the source in units of [Bq]. +4. The :meth:`openmc.Material.get_decay_particle_energy` method can be used to + obtain the energy spectrum of decay particles from the activated material. + For R2S, the photon spectrum is needed and can be obtained either with + :meth:`openmc.Material.get_decay_particle_energy` using + ``particle='photon'`` or with the wrapper + :meth:`openmc.Material.get_decay_photon_energy`. The integral of the + spectrum indicates the intensity of the source in units of [Bq]. 5. A new photon source is defined using one of OpenMC's source classes with the energy distribution set equal to the object returned by the :meth:`openmc.Material.get_decay_photon_energy` method. The source is then @@ -76,9 +80,10 @@ Altogether, the workflow looks as follows:: model.settings.source = photon_source model.run() -Note that by default, the :meth:`~openmc.Material.get_decay_photon_energy` -method will eliminate spectral lines with very low intensity, but this behavior -can be configured with the ``clip_tolerance`` argument. +Note that :meth:`~openmc.Material.get_decay_particle_energy` and +:meth:`~openmc.Material.get_decay_photon_energy` will, by default, eliminate +spectral lines with very low intensity, but this behavior can be configured +with the ``clip_tolerance`` argument. Cell-based R2S -------------- @@ -236,4 +241,3 @@ relevant tallies. This can be done with the aid of the # Apply time correction factors tally = d1s.apply_time_correction(dose_tally, factors, time_index) - diff --git a/openmc/config.py b/openmc/config.py index 23d8e23a7b9..be64a285019 100644 --- a/openmc/config.py +++ b/openmc/config.py @@ -20,7 +20,7 @@ from typing import Any, Dict, Iterator from openmc.data import DataLibrary -from openmc.data.decay import _DECAY_ENERGY, _DECAY_PHOTON_ENERGY +from openmc.data.decay import _DECAY_ENERGY, _DECAY_PARTICLE_ENERGY __all__ = ["config"] @@ -77,7 +77,7 @@ def __delitem__(self, key: str): if env_var in os.environ: del os.environ[env_var] if key == 'chain_file': - _DECAY_PHOTON_ENERGY.clear() + _DECAY_PARTICLE_ENERGY.clear() _DECAY_ENERGY.clear() def __setitem__(self, key: str, value: Any): @@ -103,7 +103,7 @@ def __setitem__(self, key: str, value: Any): os.environ[self._PATH_KEYS[key]] = str(stored_path) if key == 'chain_file': - _DECAY_PHOTON_ENERGY.clear() + _DECAY_PARTICLE_ENERGY.clear() _DECAY_ENERGY.clear() if not stored_path.exists(): diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 7cd4bf43d41..40ce8027193 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -561,7 +561,7 @@ def sources(self): return merged_sources -_DECAY_PHOTON_ENERGY = {} +_DECAY_PARTICLE_ENERGY = {} def decay_photon_energy(nuclide: str) -> Univariate | None: @@ -571,7 +571,8 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: for the first time, you need to ensure that a depletion chain has been specified in openmc.config['chain_file']. - .. versionadded:: 0.13.2 + This is a backward-compatible convenience wrapper around + :func:`decay_particle_energy` with ``particle='photon'``. Parameters ---------- @@ -585,7 +586,34 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: if no photon source exists. Note that the probabilities represent intensities, given as [Bq/atom] (in other words, decay constants). """ - if not _DECAY_PHOTON_ENERGY: + return decay_particle_energy(nuclide, 'photon') + + +def decay_particle_energy(nuclide: str, particle: str) -> Univariate | None: + """Get a decay particle energy distribution resulting from the decay of a nuclide + + This function relies on data stored in a depletion chain. Before calling it + for the first time, you need to ensure that a depletion chain has been + specified in openmc.config['chain_file']. + + Parameters + ---------- + nuclide : str + Name of nuclide, e.g., 'Co58' + particle : str + Name of the decay particle, e.g., ``photon``, ``electron``, + ``positron``, ``alpha``, ``neutron``, ``proton``, or ``fragment``. + + Returns + ------- + openmc.stats.Univariate or None + Distribution of energies in [eV] of decay particles emitted from decay, + or None if no decay source of that particle exists. Note that the + probabilities represent intensities, given as [Bq/atom] (in other + words, decay constants). + """ + + if not _DECAY_PARTICLE_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: raise DataError( @@ -596,15 +624,15 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: - if 'photon' in nuc.sources: - _DECAY_PHOTON_ENERGY[nuc.name] = nuc.sources['photon'] + for part, dist in nuc.sources.items(): + _DECAY_PARTICLE_ENERGY[(nuc.name, part)] = dist # If the chain file contained no sources at all, warn the user - if not _DECAY_PHOTON_ENERGY: - warn(f"Chain file '{chain_file}' does not have any decay photon " - "sources listed.") + if not _DECAY_PARTICLE_ENERGY: + warn(f"Chain file '{chain_file}' does not have any decay particle " + "source listed.") - return _DECAY_PHOTON_ENERGY.get(nuclide) + return _DECAY_PARTICLE_ENERGY.get((nuclide, particle)) _DECAY_ENERGY = {} @@ -649,5 +677,3 @@ def decay_energy(nuclide: str): warn(f"Chain file '{chain_file}' does not have any decay energy.") return _DECAY_ENERGY.get(nuclide, 0.0) - - diff --git a/openmc/material.py b/openmc/material.py index 2dbe691f525..9db9dea2902 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -332,20 +332,23 @@ def decay_photon_energy(self) -> Univariate | None: "version.", FutureWarning) return self.get_decay_photon_energy(0.0) - def get_decay_photon_energy( + def get_decay_particle_energy( self, + particle: str = 'photon', clip_tolerance: float = 1e-6, units: str = 'Bq', volume: float | None = None, exclude_nuclides: list[str] | None = None, include_nuclides: list[str] | None = None ) -> Univariate | None: - r"""Return energy distribution of decay photons from unstable nuclides. - - .. versionadded:: 0.14.0 + r"""Return energy distribution of decay particles from unstable nuclides. Parameters ---------- + particle : string + Particle type to return the energy distribution for. Supported + values are ``'photon'``, ``'electron'``, ``'positron'``, + ``'alpha'``, ``'neutron'``, ``'proton'``, and ``'fragment'``. clip_tolerance : float Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions that will be discarded. @@ -355,19 +358,21 @@ def get_decay_photon_energy( Volume of the material. If not passed, defaults to using the :attr:`Material.volume` attribute. exclude_nuclides : list of str, optional - Nuclides to exclude from the photon source calculation. + Nuclides to exclude from the particle source calculation. include_nuclides : list of str, optional - Nuclides to include in the photon source calculation. If specified, + Nuclides to include in the particle source calculation. If specified, only these nuclides are used. Returns ------- Univariate or None - Decay photon energy distribution. The integral of this distribution is - the total intensity of the photon source in the requested units. + Decay particle energy distribution. The integral of this distribution is + the total intensity of the particle source in the requested units. """ cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}) + cv.check_type('particle', particle, str) + cv.check_value('particle', particle, {'alpha', 'fragment', 'neutron','photon', 'electron', 'positron','proton'}) if exclude_nuclides is not None and include_nuclides is not None: raise ValueError("Cannot specify both exclude_nuclides and include_nuclides") @@ -393,12 +398,12 @@ def get_decay_photon_energy( if include_nuclides is not None and nuc not in include_nuclides: continue - source_per_atom = openmc.data.decay_photon_energy(nuc) + source_per_atom = openmc.data.decay_particle_energy(nuc,particle) if source_per_atom is not None and atoms_per_bcm > 0.0: dists.append(source_per_atom) probs.append(1e24 * atoms_per_bcm * multiplier) - # If no photon sources, exit early + # If no particle sources, exit early if not dists: return None @@ -414,6 +419,53 @@ def get_decay_photon_energy( return combined + def get_decay_photon_energy( + self, + clip_tolerance: float = 1e-6, + units: str = 'Bq', + volume: float | None = None, + exclude_nuclides: list[str] | None = None, + include_nuclides: list[str] | None = None + ) -> Univariate | None: + r"""Return energy distribution of decay photons from unstable nuclides. + + This is a backward-compatible convenience wrapper around + :meth:`Material.get_decay_particle_energy` with + ``particle='photon'``. + + Parameters + ---------- + clip_tolerance : float + Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions + that will be discarded. + units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'} + Specifies the units on the integral of the distribution. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + exclude_nuclides : list of str, optional + Nuclides to exclude from the photon source calculation. + include_nuclides : list of str, optional + Nuclides to include in the photon source calculation. If specified, + only these nuclides are used. + + Returns + ------- + Univariate or None + Decay photon energy distribution. The integral of this distribution is + the total intensity of the photon source in the requested units. + + """ + + return self.get_decay_particle_energy( + particle='photon', + clip_tolerance=clip_tolerance, + units=units, + volume=volume, + exclude_nuclides=exclude_nuclides, + include_nuclides=include_nuclides, + ) + def get_photon_contact_dose_rate( self, dose_quantity: str = "absorbed-air", @@ -540,7 +592,7 @@ def get_photon_contact_dose_rate( * 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) + photon_source_per_atom = openmc.data.decay_particle_energy(nuc,'photon') # nuclides with no contribution if photon_source_per_atom is None or nuc_atoms_per_bcm <= 0.0: diff --git a/tests/unit_tests/test_config.py b/tests/unit_tests/test_config.py index 45e6a7e5a1b..fe9112cff83 100644 --- a/tests/unit_tests/test_config.py +++ b/tests/unit_tests/test_config.py @@ -99,6 +99,6 @@ def test_config_chain_side_effect(tmp_path): """Test that modifying chain_file clears decay data caches.""" chain_file = tmp_path / "chain.xml"; chain_file.touch() decay._DECAY_ENERGY['U235'] = (1.0, 2.0) - decay._DECAY_PHOTON_ENERGY['PU239'] = {} + decay._DECAY_PARTICLE_ENERGY['PU239','photon'] = {} openmc.config['chain_file'] = chain_file - assert not decay._DECAY_ENERGY and not decay._DECAY_PHOTON_ENERGY + assert not decay._DECAY_ENERGY and not decay._DECAY_PARTICLE_ENERGY diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index de8d90a434d..1b3ac35cb8f 100644 --- a/tests/unit_tests/test_data_decay.py +++ b/tests/unit_tests/test_data_decay.py @@ -146,3 +146,23 @@ def test_decay_photon_energy(): src = openmc.data.decay_photon_energy('Xe135') assert isinstance(src, openmc.stats.Tabular) assert src.integral() == pytest.approx(2.076506258964966e-05) + + +def test_decay_particle_energy(): + # If chain file is not set, we should get a data error + if 'chain_file' in openmc.config: + del openmc.config['chain_file'] + with pytest.raises(DataError): + openmc.data.decay_particle_energy('Fe55', 'electron') + + # Set chain file to Ni chain and check electron/positron sources + with openmc.config.patch('chain_file', Path(__file__).parents[1] / "chain_ni.xml"): + src = openmc.data.decay_particle_energy('Fe55', 'electron') + assert isinstance(src, openmc.stats.Discrete) + assert src.integral() == pytest.approx(4.225437849490689e-08) + assert 6377.571 in src.x + + src = openmc.data.decay_particle_energy('Fe55', 'positron') + assert isinstance(src, openmc.stats.Discrete) + assert src.integral() == pytest.approx(8.004558990612365e-09) + assert 231210.0 in src.x diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 911b8867f96..86d19e4b110 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -6,7 +6,7 @@ import numpy as np import openmc -from openmc.data import decay_photon_energy +from openmc.data import decay_particle_energy from openmc.deplete import Chain import openmc.examples import openmc.model @@ -699,13 +699,13 @@ def intensity(src): return src.integral() if src is not None else 0.0 assert src.integral() == pytest.approx(sum( - intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides() + intensity(decay_particle_energy(nuc,"photon")) for nuc in m.get_nuclides() ), rel=1e-3) # When the clipping threshold is zero, the intensities should match exactly src = m.get_decay_photon_energy(0.0) assert src.integral() == pytest.approx(sum( - intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides() + intensity(decay_particle_energy(nuc,"photon")) for nuc in m.get_nuclides() )) # A material with no unstable nuclides should have no decay photon source @@ -715,6 +715,54 @@ def intensity(src): assert stable.get_decay_photon_energy() is None +def test_decay_particle_energy(): + with openmc.config.patch('chain_file', Path(__file__).parents[1] / 'chain_ni.xml'): + # Material representing single atom of Fe55 and Co58 + m = openmc.Material() + m.add_nuclide('Fe55', 1.0e-24) + m.add_nuclide('Co58', 1.0e-24) + m.volume = 1.0 + + # Get decay positron source and make sure it's the right type + src = m.get_decay_particle_energy('positron') + assert isinstance(src, openmc.stats.Discrete) + + # Make sure units/volume work as expected + src_v2 = m.get_decay_particle_energy('positron', volume=2.0) + assert src.p * 2.0 == pytest.approx(src_v2.p) + src_per_cm3 = m.get_decay_particle_energy('positron', units='Bq/cm3', + volume=100.0) + assert (src.p == src_per_cm3.p).all() + src_per_bqg = m.get_decay_particle_energy('positron', units='Bq/g') + src_per_bqkg = m.get_decay_particle_energy('positron', units='Bq/kg') + assert pytest.approx(src_per_bqg.integral()) == src_per_bqkg.integral() / 1000. + src_per_bqm3 = m.get_decay_particle_energy('positron', units='Bq/m3') + assert pytest.approx(src_per_bqm3.integral()) == src_per_cm3.integral() * 1e6 + + # With a single atom of each, the intensity of the positron source + # should be equal to the sum of the intensities for each nuclide + def intensity(src): + return src.integral() if src is not None else 0.0 + + assert src.integral() == pytest.approx(sum( + intensity(decay_particle_energy(nuc, 'positron')) + for nuc in m.get_nuclides() + ), rel=1e-3) + + # When the clipping threshold is zero, the intensities should match exactly + src = m.get_decay_particle_energy('positron', 0.0) + assert src.integral() == pytest.approx(sum( + intensity(decay_particle_energy(nuc, 'positron')) + for nuc in m.get_nuclides() + )) + + # A material with no unstable nuclides should have no decay positron source + stable = openmc.Material() + stable.add_nuclide('Gd156', 1.0) + stable.volume = 1.0 + assert stable.get_decay_particle_energy('positron') is None + + def test_avoid_subnormal(run_in_tmpdir): # Write a materials.xml with a material that has a nuclide density that is # represented as a subnormal floating point value