Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Decay photon source functionality #2234

Merged
merged 12 commits into from
Sep 30, 2022
1 change: 1 addition & 0 deletions docs/source/pythonapi/data.rst
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ Core Functions
atomic_weight
combine_distributions
decay_constant
decay_photon_energy
dose_coefficients
gnd_name
half_life
Expand Down
7 changes: 7 additions & 0 deletions openmc/config.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
import warnings

from openmc.data import DataLibrary
from openmc.data.decay import _DECAY_PHOTON_ENERGY

__all__ = ["config"]

Expand All @@ -22,6 +23,10 @@ def __delitem__(self, key):
del os.environ['OPENMC_CROSS_SECTIONS']
elif key == 'mg_cross_sections':
del os.environ['OPENMC_MG_CROSS_SECTIONS']
elif key == 'chain_file':
del os.environ['OPENMC_CHAIN_FILE']
# Reset photon source data since it relies on chain file
_DECAY_PHOTON_ENERGY.clear()

def __setitem__(self, key, value):
if key == 'cross_sections':
Expand All @@ -34,6 +39,8 @@ def __setitem__(self, key, value):
elif key == 'chain_file':
self._set_path(key, value)
os.environ['OPENMC_CHAIN_FILE'] = str(value)
# Reset photon source data since it relies on chain file
_DECAY_PHOTON_ENERGY.clear()
else:
raise KeyError(f'Unrecognized config key: {key}')

Expand Down
51 changes: 50 additions & 1 deletion openmc/data/decay.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,19 @@
from collections.abc import Iterable
from io import StringIO
from math import log
from pathlib import Path
import re
from warnings import warn
from xml.etree import ElementTree as ET

import numpy as np
from uncertainties import ufloat, UFloat

import openmc
import openmc.checkvalue as cv
from openmc.exceptions import DataError
from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular, combine_distributions
from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
from .function import INTERPOLATION_SCHEME
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
Expand Down Expand Up @@ -570,3 +574,48 @@ def sources(self):

self._sources = merged_sources
return self._sources


_DECAY_PHOTON_ENERGY = {}


def decay_photon_energy(nuclide: str) -> Univariate:
"""Get photon 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'].

.. versionadded:: 0.14.0

Parameters
----------
nuclide : str
Name of nuclide, e.g., 'Co58'

Returns
-------
openmc.stats.Univariate
Distribution of energies in [eV] of photons emitted from decay. Note
that the probabilities represent intensities, given as [decay/sec].
"""
if not _DECAY_PHOTON_ENERGY:
chain_file = openmc.config.get('chain_file')
shimwell marked this conversation as resolved.
Show resolved Hide resolved
if chain_file is None:
raise DataError(
"A depletion chain file must be specified with "
"openmc.config['chain_file'] in order to load decay data."
)
pshriwise marked this conversation as resolved.
Show resolved Hide resolved

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']

# 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.")

return _DECAY_PHOTON_ENERGY.get(nuclide)
2 changes: 2 additions & 0 deletions openmc/deplete/chain.py
Original file line number Diff line number Diff line change
Expand Up @@ -420,6 +420,8 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
# Append decay mode
nuclide.add_decay_mode(type_, target, br)

nuclide.sources = data.sources

fissionable = False
if parent in reactions:
reactions_available = set(reactions[parent].keys())
Expand Down
26 changes: 26 additions & 0 deletions openmc/deplete/nuclide.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
import numpy as np

from openmc.checkvalue import check_type
from openmc.stats import Univariate

__all__ = [
"DecayTuple", "ReactionTuple", "Nuclide", "FissionYield",
Expand Down Expand Up @@ -101,6 +102,9 @@ class Nuclide:
reactions : list of openmc.deplete.ReactionTuple
Reaction information. Each element of the list is a named tuple with
attribute 'type', 'target', 'Q', and 'branching_ratio'.
sources : dict
Dictionary mapping particle type as string to energy distribution of
decay source represented as :class:`openmc.stats.Univariate`
yield_data : FissionYieldDistribution or None
Fission product yields at tabulated energies for this nuclide. Can be
treated as a nested dictionary ``{energy: {product: yield}}``
Expand All @@ -120,6 +124,9 @@ def __init__(self, name=None):
# Reaction paths
self.reactions = []

# Decay sources
self.sources = {}

# Neutron fission yields, if present
self._yield_data = None

Expand Down Expand Up @@ -237,6 +244,12 @@ def from_xml(cls, element, root=None, fission_q=None):
branching_ratio = float(decay_elem.get('branching_ratio'))
nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))

# Check for sources
for src_elem in element.iter('source'):
particle = src_elem.get('particle')
distribution = Univariate.from_xml_element(src_elem)
nuc.sources[particle] = distribution

# Check for reaction paths
for reaction_elem in element.iter('reaction'):
r_type = reaction_elem.get('type')
Expand Down Expand Up @@ -302,6 +315,19 @@ def to_xml_element(self):
mode_elem.set('target', daughter)
mode_elem.set('branching_ratio', str(br))

# Write decay sources
if self.sources:
for particle, source in self.sources.items():
# TODO: Ugly hack to deal with the fact that
# 'source.to_xml_element' will return an xml.etree object
# whereas here lxml is being used preferentially. We should just
# switch to use lxml everywhere
import xml.etree.ElementTree as etree
src_elem_xmletree = source.to_xml_element('source')
src_elem = ET.fromstring(etree.tostring(src_elem_xmletree))
src_elem.set('particle', particle)
elem.append(src_elem)

elem.set('reactions', str(len(self.reactions)))
for rx, daughter, Q, br in self.reactions:
rx_elem = ET.SubElement(elem, 'reaction')
Expand Down
2 changes: 1 addition & 1 deletion openmc/deplete/openmc_operator.py
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@ def _set_number_from_results(self, mat, prev_res):
mat_id = str(mat.id)

# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind
depl_nuc = prev_res[-1].index_nuc
geom_nuc_densities = mat.get_nuclide_atom_densities()

# Merge lists of nuclides, with the same order for every calculation
Expand Down
4 changes: 2 additions & 2 deletions openmc/deplete/results.py
Original file line number Diff line number Diff line change
Expand Up @@ -435,7 +435,7 @@ def export_to_materials(
# results, and save them to the new depleted xml file.
for mat in mat_file:
mat_id = str(mat.id)
if mat_id in result.mat_to_ind:
if mat_id in result.index_mat:
mat.volume = result.volume[mat_id]

# Change density of all nuclides in material to atom/b-cm
Expand All @@ -447,7 +447,7 @@ def export_to_materials(

# For nuclides in chain that have cross sections, replace
# density in original material with new density from results
for nuc in result.nuc_to_ind:
for nuc in result.index_nuc:
if nuc not in available_cross_sections:
continue
atoms = result[0, mat_id, nuc]
Expand Down
Loading