"""chain module.
This module contains information about a depletion chain. A depletion chain is
loaded from an .xml file and all the nuclides are linked together.
"""
from io import StringIO
from itertools import chain
import math
import os
import re
from collections import defaultdict, namedtuple
from collections.abc import Mapping, Iterable
from numbers import Real, Integral
from warnings import warn
from openmc.checkvalue import check_type, check_greater_than
from openmc.data import gnds_name, zam, DataLibrary
from openmc.exceptions import DataError
from .nuclide import FissionYieldDistribution
import lxml.etree as ET
import scipy.sparse as sp
import openmc.data
from openmc._xml import clean_indentation
from .nuclide import Nuclide, DecayTuple, ReactionTuple
# tuple of (possible MT values, (dA, dZ), secondaries) where dA is the change in
# the mass number and dZ is the change in the atomic number
ReactionInfo = namedtuple('ReactionInfo', ('mts', 'dadz', 'secondaries'))
REACTIONS = {
'(n,2nd)': ReactionInfo({11}, (-3, -1), ('H2',)),
'(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), (-1, 0), ()),
'(n,3n)': ReactionInfo({17}, (-2, 0), ()),
'(n,na)': ReactionInfo({22}, (-4, -2), ('He4',)),
'(n,n3a)': ReactionInfo({23}, (-12, -6), ('He4', 'He4', 'He4')),
'(n,2na)': ReactionInfo({24}, (-5, -2), ('He4',)),
'(n,3na)': ReactionInfo({25}, (-6, -2), ('He4',)),
'(n,np)': ReactionInfo({28}, (-1, -1), ('H1',)),
'(n,n2a)': ReactionInfo({29}, (-8, -4), ('He4', 'He4')),
'(n,2n2a)': ReactionInfo({30}, (-9, -4), ('He4', 'He4')),
'(n,nd)': ReactionInfo({32}, (-2, -1), ('H2',)),
'(n,nt)': ReactionInfo({33}, (-3, -1), ('H3',)),
'(n,n3He)': ReactionInfo({34}, (-3, -2), ('He3',)),
'(n,nd2a)': ReactionInfo({35}, (-10, -5), ('H2', 'He4', 'He4')),
'(n,nt2a)': ReactionInfo({36}, (-11, -5), ('H3', 'He4', 'He4')),
'(n,4n)': ReactionInfo({37}, (-3, 0), ()),
'(n,2np)': ReactionInfo({41}, (-2, -1), ('H1',)),
'(n,3np)': ReactionInfo({42}, (-3, -1), ('H1',)),
'(n,n2p)': ReactionInfo({44}, (-2, -2), ('H1', 'H1')),
'(n,npa)': ReactionInfo({45}, (-5, -3), ('H1', 'He4')),
'(n,gamma)': ReactionInfo({102}, (1, 0), ()),
'(n,p)': ReactionInfo(set(chain([103], range(600, 650))), (0, -1), ('H1',)),
'(n,d)': ReactionInfo(set(chain([104], range(650, 700))), (-1, -1), ('H2',)),
'(n,t)': ReactionInfo(set(chain([105], range(700, 750))), (-2, -1), ('H3',)),
'(n,3He)': ReactionInfo(set(chain([106], range(750, 800))), (-2, -2), ('He3',)),
'(n,a)': ReactionInfo(set(chain([107], range(800, 850))), (-3, -2), ('He4',)),
'(n,2a)': ReactionInfo({108}, (-7, -4), ('He4', 'He4')),
'(n,3a)': ReactionInfo({109}, (-11, -6), ('He4', 'He4', 'He4')),
'(n,2p)': ReactionInfo({111}, (-1, -2), ('H1', 'H1')),
'(n,pa)': ReactionInfo({112}, (-4, -3), ('H1', 'He4')),
'(n,t2a)': ReactionInfo({113}, (-10, -5), ('H3', 'He4', 'He4')),
'(n,d2a)': ReactionInfo({114}, (-9, -5), ('H2', 'He4', 'He4')),
'(n,pd)': ReactionInfo({115}, (-2, -2), ('H1', 'H2')),
'(n,pt)': ReactionInfo({116}, (-3, -2), ('H1', 'H3')),
'(n,da)': ReactionInfo({117}, (-5, -3), ('H2', 'He4')),
'(n,5n)': ReactionInfo({152}, (-4, 0), ()),
'(n,6n)': ReactionInfo({153}, (-5, 0), ()),
'(n,2nt)': ReactionInfo({154}, (-4, -1), ('H3',)),
'(n,ta)': ReactionInfo({155}, (-6, -3), ('H3', 'He4')),
'(n,4np)': ReactionInfo({156}, (-4, -1), ('H1',)),
'(n,3nd)': ReactionInfo({157}, (-4, -1), ('H2',)),
'(n,nda)': ReactionInfo({158}, (-6, -3), ('H2', 'He4')),
'(n,2npa)': ReactionInfo({159}, (-6, -3), ('H1', 'He4')),
'(n,7n)': ReactionInfo({160}, (-6, 0), ()),
'(n,8n)': ReactionInfo({161}, (-7, 0), ()),
'(n,5np)': ReactionInfo({162}, (-5, -1), ('H1',)),
'(n,6np)': ReactionInfo({163}, (-6, -1), ('H1',)),
'(n,7np)': ReactionInfo({164}, (-7, -1), ('H1',)),
'(n,4na)': ReactionInfo({165}, (-7, -2), ('He4',)),
'(n,5na)': ReactionInfo({166}, (-8, -2), ('He4',)),
'(n,6na)': ReactionInfo({167}, (-9, -2), ('He4',)),
'(n,7na)': ReactionInfo({168}, (-10, -2), ('He4',)),
'(n,4nd)': ReactionInfo({169}, (-5, -1), ('H2',)),
'(n,5nd)': ReactionInfo({170}, (-6, -1), ('H2',)),
'(n,6nd)': ReactionInfo({171}, (-7, -1), ('H2',)),
'(n,3nt)': ReactionInfo({172}, (-5, -1), ('H3',)),
'(n,4nt)': ReactionInfo({173}, (-6, -1), ('H3',)),
'(n,5nt)': ReactionInfo({174}, (-7, -1), ('H3',)),
'(n,6nt)': ReactionInfo({175}, (-8, -1), ('H3',)),
'(n,2n3He)': ReactionInfo({176}, (-4, -2), ('He3',)),
'(n,3n3He)': ReactionInfo({177}, (-5, -2), ('He3',)),
'(n,4n3He)': ReactionInfo({178}, (-6, -2), ('He3',)),
'(n,3n2p)': ReactionInfo({179}, (-4, -2), ('H1', 'H1')),
'(n,3n2a)': ReactionInfo({180}, (-10, -4), ('He4', 'He4')),
'(n,3npa)': ReactionInfo({181}, (-7, -3), ('H1', 'He4')),
'(n,dt)': ReactionInfo({182}, (-4, -2), ('H2', 'H3')),
'(n,npd)': ReactionInfo({183}, (-3, -2), ('H1', 'H2')),
'(n,npt)': ReactionInfo({184}, (-4, -2), ('H1', 'H3')),
'(n,ndt)': ReactionInfo({185}, (-5, -2), ('H2', 'H3')),
'(n,np3He)': ReactionInfo({186}, (-4, -3), ('H1', 'He3')),
'(n,nd3He)': ReactionInfo({187}, (-5, -3), ('H2', 'He3')),
'(n,nt3He)': ReactionInfo({188}, (-6, -3), ('H3', 'He3')),
'(n,nta)': ReactionInfo({189}, (-7, -3), ('H3', 'He4')),
'(n,2n2p)': ReactionInfo({190}, (-3, -2), ('H1', 'H1')),
'(n,p3He)': ReactionInfo({191}, (-4, -3), ('H1', 'He3')),
'(n,d3He)': ReactionInfo({192}, (-5, -3), ('H2', 'He3')),
'(n,3Hea)': ReactionInfo({193}, (-6, -4), ('He3', 'He4')),
'(n,4n2p)': ReactionInfo({194}, (-5, -2), ('H1', 'H1')),
'(n,4n2a)': ReactionInfo({195}, (-11, -4), ('He4', 'He4')),
'(n,4npa)': ReactionInfo({196}, (-8, -3), ('H1', 'He4')),
'(n,3p)': ReactionInfo({197}, (-2, -3), ('H1', 'H1', 'H1')),
'(n,n3p)': ReactionInfo({198}, (-3, -3), ('H1', 'H1', 'H1')),
'(n,3n2pa)': ReactionInfo({199}, (-8, -4), ('H1', 'H1', 'He4')),
'(n,5n2p)': ReactionInfo({200}, (-6, -2), ('H1', 'H1')),
}
__all__ = ["Chain", "REACTIONS"]
def replace_missing(product, decay_data):
"""Replace missing product with suitable decay daughter.
Parameters
----------
product : str
Name of product in GNDS format, e.g. 'Y86_m1'.
decay_data : dict
Dictionary of decay data
Returns
-------
product : str
Replacement for missing product in GNDS format.
"""
# Determine atomic number, mass number, and metastable state
Z, A, state = openmc.data.zam(product)
symbol = openmc.data.ATOMIC_SYMBOL[Z]
# Replace neutron with nothing
if Z == 0:
return None
# First check if ground state is available
if state:
product = '{}{}'.format(symbol, A)
# Find isotope with longest half-life
half_life = 0.0
for nuclide, data in decay_data.items():
m = re.match(r'{}(\d+)(?:_m\d+)?'.format(symbol), nuclide)
if m:
# If we find a stable nuclide, stop search
if data.nuclide['stable']:
mass_longest_lived = int(m.group(1))
break
if data.half_life.nominal_value > half_life:
mass_longest_lived = int(m.group(1))
half_life = data.half_life.nominal_value
# If mass number of longest-lived isotope is less than that of missing
# product, assume it undergoes beta-. Otherwise assume beta+.
beta_minus = (mass_longest_lived < A)
# Iterate until we find an existing nuclide
while product not in decay_data:
if Z > 98:
# Assume alpha decay occurs for Z=99 and above
Z -= 2
A -= 4
else:
# Otherwise assume a beta- or beta+
if beta_minus:
Z += 1
else:
Z -= 1
product = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
return product
def replace_missing_fpy(actinide, fpy_data, decay_data):
"""Replace missing fission product yields
Parameters
----------
actinide : str
Name of actinide missing FPY data
fpy_data : dict
Dictionary of FPY data
decay_data : dict
Dictionary of decay data
Returns
-------
str
Actinide that can be used as replacement for FPY purposes
"""
# Check if metastable state has data (e.g., Am242m)
Z, A, m = zam(actinide)
if m == 0:
metastable = gnds_name(Z, A, 1)
if metastable in fpy_data:
return metastable
# Try increasing Z, holding N constant
isotone = actinide
while isotone in decay_data:
Z += 1
A += 1
isotone = gnds_name(Z, A, 0)
if isotone in fpy_data:
return isotone
# Try decreasing Z, holding N constant
isotone = actinide
while isotone in decay_data:
Z -= 1
A -= 1
isotone = gnds_name(Z, A, 0)
if isotone in fpy_data:
return isotone
# If all else fails, use U235 yields
return 'U235'
[docs]class Chain:
"""Full representation of a depletion chain.
A depletion chain can be created by using the :meth:`from_endf` method which
requires a list of ENDF incident neutron, decay, and neutron fission product
yield sublibrary files. The depletion chain used during a depletion
simulation is indicated by either an argument to
:class:`openmc.deplete.CoupledOperator` or
:class:`openmc.deplete.IndependentOperator`, or through
openmc.config['chain_file'].
Attributes
----------
nuclides : list of openmc.deplete.Nuclide
Nuclides present in the chain.
reactions : list of str
Reactions that are tracked in the depletion chain
nuclide_dict : dict of str to int
Maps a nuclide name to an index in nuclides.
fission_yields : None or iterable of dict
List of effective fission yields for materials. Each dictionary
should be of the form ``{parent: {product: yield}}`` with
types ``{str: {str: float}}``, where ``yield`` is the fission product
yield for isotope ``parent`` producing isotope ``product``.
A single entry indicates yields are constant across all materials.
Otherwise, an entry can be added for each material to be burned.
Ordering should be identical to how the operator orders reaction
rates for burnable materials.
"""
def __init__(self):
self.nuclides = []
self.reactions = []
self.nuclide_dict = {}
self._fission_yields = None
def __contains__(self, nuclide):
return nuclide in self.nuclide_dict
def __getitem__(self, name):
"""Get a Nuclide by name."""
return self.nuclides[self.nuclide_dict[name]]
def __len__(self):
"""Number of nuclides in chain."""
return len(self.nuclides)
[docs] def add_nuclide(self, nuclide):
"""Add a nuclide to the depletion chain
Parameters
----------
nuclide : openmc.deplete.Nuclide
Nuclide to add
"""
self.nuclide_dict[nuclide.name] = len(self.nuclides)
self.nuclides.append(nuclide)
# Check for reaction paths
for rx in nuclide.reactions:
if rx.type not in self.reactions:
self.reactions.append(rx.type)
[docs] @classmethod
def from_endf(cls, decay_files, fpy_files, neutron_files,
reactions=('(n,2n)', '(n,3n)', '(n,4n)', '(n,gamma)', '(n,p)', '(n,a)'),
progress=True
):
"""Create a depletion chain from ENDF files.
String arguments in ``decay_files``, ``fpy_files``, and
``neutron_files`` will be treated as file names to be read.
Alternatively, :class:`openmc.data.endf.Evaluation` instances
can be included in these arguments.
Parameters
----------
decay_files : list of str or openmc.data.endf.Evaluation
List of ENDF decay sub-library files
fpy_files : list of str or openmc.data.endf.Evaluation
List of ENDF neutron-induced fission product yield sub-library files
neutron_files : list of str or openmc.data.endf.Evaluation
List of ENDF neutron reaction sub-library files
reactions : iterable of str, optional
Transmutation reactions to include in the depletion chain, e.g.,
`["(n,2n)", "(n,gamma)"]`. Note that fission is always included if
it is present. A complete listing of transmutation reactions can be
found in :data:`openmc.deplete.chain.REACTIONS`.
.. versionadded:: 0.12.1
progress : bool, optional
Flag to print status messages during processing. Does not
effect warning messages
Returns
-------
Chain
Notes
-----
When an actinide is missing fission product yield (FPY) data, yields will
copied from a parent isotope, found according to:
1. If the nuclide is in a ground state and a metastable state exists with
fission yields, copy the yields from the metastable
2. Find an isotone (same number of neutrons) and copy those yields
3. Copy the yields of U235 if the previous two checks fail
"""
transmutation_reactions = reactions
# Create dictionary mapping target to filename
if progress:
print('Processing neutron sub-library files...')
reactions = {}
for f in neutron_files:
evaluation = openmc.data.endf.Evaluation(f)
name = evaluation.gnds_name
reactions[name] = {}
for mf, mt, nc, mod in evaluation.reaction_list:
if mf == 3:
file_obj = StringIO(evaluation.section[3, mt])
openmc.data.endf.get_head_record(file_obj)
q_value = openmc.data.endf.get_cont_record(file_obj)[1]
reactions[name][mt] = q_value
# Determine what decay and FPY nuclides are available
if progress:
print('Processing decay sub-library files...')
decay_data = {}
for f in decay_files:
data = openmc.data.Decay(f)
# Skip decay data for neutron itself
if data.nuclide['atomic_number'] == 0:
continue
decay_data[data.nuclide['name']] = data
if progress:
print('Processing fission product yield sub-library files...')
fpy_data = {}
for f in fpy_files:
data = openmc.data.FissionProductYields(f)
fpy_data[data.nuclide['name']] = data
if progress:
print('Creating depletion_chain...')
missing_daughter = []
missing_rx_product = []
missing_fpy = []
missing_fp = []
chain = cls()
for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)):
data = decay_data[parent]
nuclide = Nuclide(parent)
if not data.nuclide['stable'] and data.half_life.nominal_value != 0.0:
nuclide.half_life = data.half_life.nominal_value
nuclide.decay_energy = data.decay_energy.nominal_value
sum_br = 0.0
for i, mode in enumerate(data.modes):
type_ = ','.join(mode.modes)
if mode.daughter in decay_data:
target = mode.daughter
else:
print('missing {} {} {}'.format(
parent, ','.join(mode.modes), mode.daughter))
target = replace_missing(mode.daughter, decay_data)
# Write branching ratio, taking care to ensure sum is unity
br = mode.branching_ratio.nominal_value
sum_br += br
if i == len(data.modes) - 1 and sum_br != 1.0:
br = 1.0 - sum(m.branching_ratio.nominal_value
for m in data.modes[:-1])
# 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())
for name in transmutation_reactions:
mts, changes, _ = REACTIONS[name]
if mts & reactions_available:
delta_A, delta_Z = changes
A = data.nuclide['mass_number'] + delta_A
Z = data.nuclide['atomic_number'] + delta_Z
daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
if daughter not in decay_data:
daughter = replace_missing(daughter, decay_data)
if daughter is None:
missing_rx_product.append((parent, name, daughter))
# Store Q value
for mt in sorted(mts):
if mt in reactions[parent]:
q_value = reactions[parent][mt]
break
else:
q_value = 0.0
nuclide.add_reaction(name, daughter, q_value, 1.0)
if any(mt in reactions_available for mt in openmc.data.FISSION_MTS):
q_value = reactions[parent][18]
nuclide.add_reaction('fission', None, q_value, 1.0)
fissionable = True
if fissionable:
if parent in fpy_data:
fpy = fpy_data[parent]
if fpy.energies is not None:
yield_energies = fpy.energies
else:
yield_energies = [0.0]
yield_data = {}
for E, yield_table in zip(yield_energies, fpy.independent):
yield_replace = 0.0
yields = defaultdict(float)
for product, y in yield_table.items():
# Handle fission products that have no decay data
if product not in decay_data:
daughter = replace_missing(product, decay_data)
product = daughter
yield_replace += y.nominal_value
yields[product] += y.nominal_value
if yield_replace > 0.0:
missing_fp.append((parent, E, yield_replace))
yield_data[E] = yields
nuclide.yield_data = FissionYieldDistribution(yield_data)
else:
nuclide._fpy = replace_missing_fpy(parent, fpy_data, decay_data)
missing_fpy.append((parent, nuclide._fpy))
# Add nuclide to chain
chain.add_nuclide(nuclide)
# Replace missing FPY data
for nuclide in chain.nuclides:
if hasattr(nuclide, '_fpy'):
nuclide.yield_data = chain[nuclide._fpy].yield_data
# Display warnings
if missing_daughter:
print('The following decay modes have daughters with no decay data:')
for mode in missing_daughter:
print(' {}'.format(mode))
print('')
if missing_rx_product:
print('The following reaction products have no decay data:')
for vals in missing_rx_product:
print('{} {} -> {}'.format(*vals))
print('')
if missing_fpy:
print('The following fissionable nuclides have no fission product yields:')
for parent, replacement in missing_fpy:
print(' {}, replaced with {}'.format(parent, replacement))
print('')
if missing_fp:
print('The following nuclides have fission products with no decay data:')
for vals in missing_fp:
print(' {}, E={} eV (total yield={})'.format(*vals))
return chain
[docs] @classmethod
def from_xml(cls, filename, fission_q=None):
"""Reads a depletion chain XML file.
Parameters
----------
filename : str
The path to the depletion chain XML file.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
If not given, values will be pulled from ``filename``
"""
chain = cls()
if fission_q is not None:
check_type("fission_q", fission_q, Mapping)
else:
fission_q = {}
# Load XML tree
root = ET.parse(str(filename))
for i, nuclide_elem in enumerate(root.findall('nuclide')):
this_q = fission_q.get(nuclide_elem.get("name"))
nuc = Nuclide.from_xml(nuclide_elem, root, this_q)
chain.add_nuclide(nuc)
return chain
[docs] def export_to_xml(self, filename):
"""Writes a depletion chain XML file.
Parameters
----------
filename : str
The path to the depletion chain XML file.
"""
root_elem = ET.Element('depletion_chain')
for nuclide in self.nuclides:
root_elem.append(nuclide.to_xml_element())
tree = ET.ElementTree(root_elem)
tree.write(str(filename), encoding='utf-8', pretty_print=True)
[docs] def get_default_fission_yields(self):
"""Return fission yields at lowest incident neutron energy
Used as the default set of fission yields for :meth:`form_matrix`
if ``fission_yields`` are not provided
Returns
-------
fission_yields : dict
Dictionary of ``{parent: {product: f_yield}}``
where ``parent`` and ``product`` are both string
names of nuclides with yield data and ``f_yield``
is a float for the fission yield.
"""
out = defaultdict(dict)
for nuc in self.nuclides:
if nuc.yield_data is None:
continue
yield_obj = nuc.yield_data[min(nuc.yield_energies)]
out[nuc.name] = dict(yield_obj)
return out
[docs] def get_branch_ratios(self, reaction="(n,gamma)"):
"""Return a dictionary with reaction branching ratios
Parameters
----------
reaction : str, optional
Reaction name like ``"(n,gamma)"`` [default], or
``"(n,alpha)"``.
Returns
-------
branches : dict
nested dict of parent nuclide keys with reaction targets and
branching ratios. Consider the capture, ``"(n,gamma)"``,
reaction for Am241::
{"Am241": {"Am242": 0.91, "Am242_m1": 0.09}}
See Also
--------
:meth:`set_branch_ratios`
"""
capt = {}
for nuclide in self.nuclides:
nuc_capt = {}
for rx in nuclide.reactions:
if rx.type == reaction and rx.branching_ratio != 1.0:
nuc_capt[rx.target] = rx.branching_ratio
if len(nuc_capt) > 0:
capt[nuclide.name] = nuc_capt
return capt
[docs] def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)",
strict=True, tolerance=1e-5):
"""Set the branching ratios for a given reactions
Parameters
----------
branch_ratios : dict of {str: {str: float}}
Capture branching ratios to be inserted.
First layer keys are names of parent nuclides, e.g.
``"Am241"``. The branching ratios for these
parents will be modified. Corresponding values are
dictionaries of ``{target: branching_ratio}``
reaction : str, optional
Reaction name like ``"(n,gamma)"`` [default], or
``"(n, alpha)"``.
strict : bool, optional
Error control. If this evalutes to ``True``, then errors will
be raised if inconsistencies are found. Otherwise, warnings
will be raised for most issues.
tolerance : float, optional
Tolerance on the sum of all branching ratios for a
single parent. Will be checked with::
1 - tol < sum_br < 1 + tol
Raises
------
IndexError
If no isotopes were found on the chain that have the requested
reaction
KeyError
If ``strict`` evaluates to ``False`` and a parent isotope in
``branch_ratios`` does not exist on the chain
AttributeError
If ``strict`` evaluates to ``False`` and a parent isotope in
``branch_ratios`` does not have the requested reaction
ValueError
If ``strict`` evalutes to ``False`` and the sum of one parents
branch ratios is outside 1 +/- ``tolerance``
See Also
--------
:meth:`get_branch_ratios`
"""
# Store some useful information through the validation stage
sums = {}
rxn_ix_map = {}
grounds = {}
tolerance = abs(tolerance)
missing_parents = set()
missing_products = {}
missing_reaction = set()
bad_sums = {}
# Secondary products, like alpha particles, should not be modified
secondary = REACTIONS[reaction].secondaries
# Check for validity before manipulation
for parent, sub in branch_ratios.items():
if parent not in self:
if strict:
raise KeyError(parent)
missing_parents.add(parent)
continue
# Make sure all products are present in the chain
prod_flag = False
for product in sub:
if product not in self:
if strict:
raise KeyError(product)
missing_products[parent] = product
prod_flag = True
break
if prod_flag:
continue
# Make sure this nuclide has the reaction
indexes = []
for ix, rx in enumerate(self[parent].reactions):
if rx.type == reaction and rx.target not in secondary:
indexes.append(ix)
if "_m" not in rx.target:
grounds[parent] = rx.target
if len(indexes) == 0:
if strict:
raise AttributeError(
"Nuclide {} does not have {} reactions".format(
parent, reaction))
missing_reaction.add(parent)
continue
this_sum = sum(sub.values())
# sum of branching ratios can be lower than 1 if no ground
# target is given, but never greater
if (this_sum >= 1 + tolerance or (grounds[parent] in sub
and this_sum <= 1 - tolerance)):
if strict:
msg = ("Sum of {} branching ratios for {} "
"({:7.3f}) outside tolerance of 1 +/- "
"{:5.3e}".format(
reaction, parent, this_sum, tolerance))
raise ValueError(msg)
bad_sums[parent] = this_sum
else:
rxn_ix_map[parent] = indexes
sums[parent] = this_sum
if len(rxn_ix_map) == 0:
raise IndexError(
"No {} reactions found in this {}".format(
reaction, self.__class__.__name__))
if len(missing_parents) > 0:
warn("The following nuclides were not found in {}: {}".format(
self.__class__.__name__, ", ".join(sorted(missing_parents))))
if len(missing_reaction) > 0:
warn("The following nuclides did not have {} reactions: "
"{}".format(reaction, ", ".join(sorted(missing_reaction))))
if len(missing_products) > 0:
tail = ("{} -> {}".format(k, v)
for k, v in sorted(missing_products.items()))
warn("The following products were not found in the {} and "
"parents were unmodified: \n{}".format(
self.__class__.__name__, ", ".join(tail)))
if len(bad_sums) > 0:
tail = ("{}: {:5.3f}".format(k, s)
for k, s in sorted(bad_sums.items()))
warn("The following parent nuclides were given {} branch ratios "
"with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format(
reaction, tolerance, "\n".join(tail)))
# Insert new ReactionTuples with updated branch ratios
for parent_name, rxn_index in rxn_ix_map.items():
parent = self[parent_name]
new_ratios = branch_ratios[parent_name]
rxn_index = rxn_ix_map[parent_name]
# Assume Q value is independent of target state
rxn_Q = parent.reactions[rxn_index[0]].Q
# Remove existing reactions
for ix in reversed(rxn_index):
parent.reactions.pop(ix)
# Add new reactions
all_meta = True
for target, br in new_ratios.items():
all_meta = all_meta and ("_m" in target)
parent.add_reaction(reaction, target, rxn_Q, br)
# If branching ratios don't add to unity, add reaction to ground
# with remainder of branching ratio
if all_meta and sums[parent_name] != 1.0:
ground_br = 1.0 - sums[parent_name]
ground_target = grounds.get(parent_name)
if ground_target is None:
pz, pa, pm = zam(parent_name)
ground_target = gnds_name(pz, pa + 1, 0)
new_ratios[ground_target] = ground_br
parent.add_reaction(reaction, ground_target, rxn_Q, ground_br)
@property
def fission_yields(self):
if self._fission_yields is None:
self._fission_yields = [self.get_default_fission_yields()]
return self._fission_yields
@fission_yields.setter
def fission_yields(self, yields):
if yields is not None:
if isinstance(yields, Mapping):
yields = [yields]
check_type("fission_yields", yields, Iterable, Mapping)
self._fission_yields = yields
[docs] def validate(self, strict=True, quiet=False, tolerance=1e-4):
"""Search for possible inconsistencies
The following checks are performed for all nuclides present:
1) For all non-fission reactions, does the sum of branching
ratios equal about one?
2) For fission reactions, does the sum of fission yield
fractions equal about two?
Parameters
----------
strict : bool, optional
Raise exceptions at the first inconsistency if true.
Otherwise mark a warning
quiet : bool, optional
Flag to suppress warnings and return immediately at
the first inconsistency. Used only if
``strict`` does not evaluate to ``True``.
tolerance : float, optional
Absolute tolerance for comparisons. Used to compare computed
value ``x`` to intended value ``y`` as::
valid = (y - tolerance <= x <= y + tolerance)
Returns
-------
valid : bool
True if no inconsistencies were found
Raises
------
ValueError
If ``strict`` evaluates to ``True`` and an inconistency was
found
See Also
--------
openmc.deplete.Nuclide.validate
"""
check_type("tolerance", tolerance, Real)
check_greater_than("tolerance", tolerance, 0.0, True)
valid = True
# Sort through nuclides by name
for name in sorted(self.nuclide_dict):
stat = self[name].validate(strict, quiet, tolerance)
if quiet and not stat:
return stat
valid = valid and stat
return valid
[docs] def reduce(self, initial_isotopes, level=None):
"""Reduce the size of the chain by following transmutation paths
As an example, consider a simple chain with the following
isotopes and transmutation paths::
U235 (n,gamma) U236
(n,fission) (Xe135, I135, Cs135)
I135 (beta decay) Xe135 (beta decay) Cs135
Xe135 (n,gamma) Xe136
Calling ``chain.reduce(["I135"])`` will produce a depletion
chain that contains only isotopes that would originate from
I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not
be included, but multiple isotopes can be used to start
the search.
The ``level`` value controls the depth of the search.
``chain.reduce(["U235"], level=1)`` would return a chain
with all isotopes except Xe136, since it is two transmutations
removed from U235 in this case.
While targets will not be included in the new chain, the
total destruction rate and decay rate of included isotopes
will be preserved.
Parameters
----------
initial_isotopes : iterable of str
Start the search based on the contents of these isotopes
level : int, optional
Depth of transmuation path to follow. Must be greater than
or equal to zero. A value of zero returns a chain with
``initial_isotopes``. The default value of None implies
that all isotopes that appear in the transmutation paths
of the initial isotopes and their progeny should be
explored
Returns
-------
Chain
Depletion chain containing isotopes that would appear
after following up to ``level`` reactions and decay paths
"""
check_type("initial_isotopes", initial_isotopes, Iterable, str)
if level is None:
level = math.inf
else:
check_type("level", level, Integral)
check_greater_than("level", level, 0, equality=True)
all_isotopes = self._follow(set(initial_isotopes), level)
# Avoid re-sorting for fission yields
name_sort = sorted(all_isotopes)
new_chain = type(self)()
for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
previous = self[iso]
new_nuclide = Nuclide(previous.name)
new_nuclide.half_life = previous.half_life
new_nuclide.decay_energy = previous.decay_energy
new_nuclide.sources = previous.sources.copy()
if hasattr(previous, '_fpy'):
new_nuclide._fpy = previous._fpy
for mode in previous.decay_modes:
if mode.target in all_isotopes:
new_nuclide.add_decay_mode(*mode)
else:
new_nuclide.add_decay_mode(mode.type, None, mode.branching_ratio)
for rx in previous.reactions:
if rx.target in all_isotopes:
new_nuclide.add_reaction(*rx)
elif rx.type == "fission":
new_yields = new_nuclide.yield_data = (
previous.yield_data.restrict_products(name_sort))
if new_yields is not None:
new_nuclide.add_reaction(*rx)
# Maintain total destruction rates but set no target
else:
new_nuclide.add_reaction(rx.type, None, rx.Q, rx.branching_ratio)
new_chain.add_nuclide(new_nuclide)
# Doesn't appear that the ordering matters for the reactions,
# just the contents
new_chain.reactions = sorted(new_chain.reactions)
return new_chain
def _follow(self, isotopes, level):
"""Return all isotopes present up to depth level"""
found = isotopes.copy()
remaining = set(self.nuclide_dict)
if not found.issubset(remaining):
raise IndexError(
"The following isotopes were not found in the chain: "
"{}".format(", ".join(found - remaining)))
if level == 0:
return found
remaining -= found
depth = 0
next_iso = set()
while depth < level and remaining:
# Exhaust all isotopes at this level
while isotopes:
iso = isotopes.pop()
found.add(iso)
nuclide = self[iso]
# Follow all transmutation paths for this nuclide
for rxn in nuclide.reactions + nuclide.decay_modes:
if rxn.type == "fission":
continue
# Figure out if this reaction produces light nuclides
if rxn.type in REACTIONS:
secondaries = REACTIONS[rxn.type].secondaries
else:
secondaries = []
# Only include secondaries if they are present in original chain
secondaries = [x for x in secondaries if x in self]
for product in chain([rxn.target], secondaries):
if product is None:
continue
# Skip if we've already come across this isotope
elif (product in next_iso or product in found
or product in isotopes):
continue
next_iso.add(product)
if nuclide.yield_data is not None:
for product in nuclide.yield_data.products:
if (product in next_iso
or product in found or product in isotopes):
continue
next_iso.add(product)
if not next_iso:
# No additional isotopes to process, nor to update the
# current set of discovered isotopes
return found
# Prepare for next dig
depth += 1
isotopes |= next_iso
remaining -= next_iso
next_iso.clear()
# Process isotope that would have started next depth
found.update(isotopes)
return found