Source code for openmc.deplete.chain

"""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 OrderedDict, defaultdict, namedtuple
from import Mapping, Iterable
from numbers import Real, Integral
from warnings import warn

from openmc.checkvalue import check_type, check_greater_than
from import gnd_name, zam, DataLibrary
from openmc.exceptions import DataError
from .nuclide import FissionYieldDistribution

# Try to use lxml if it is available. It preserves the order of attributes and
# provides a pretty-printer by default. If not available,
# use OpenMC function to pretty print.
    import lxml.etree as ET
    _have_lxml = True
except ImportError:
    import xml.etree.ElementTree as ET
    _have_lxml = False
import scipy.sparse as sp

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

    '(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.

    product : str
        Name of product in GND format, e.g. 'Y86_m1'.
    decay_data : dict
        Dictionary of decay data

    product : str
        Replacement for missing product in GND format.

    # Determine atomic number, mass number, and metastable state
    Z, A, state =
    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(
            if data.half_life.nominal_value > half_life:
                mass_longest_lived = int(
                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
            # Otherwise assume a beta- or beta+
            if beta_minus:
                Z += 1
                Z -= 1
        product = '{}{}'.format([Z], A)

    return product

def replace_missing_fpy(actinide, fpy_data, decay_data):
    """Replace missing fission product yields

    actinide : str
        Name of actinide missing FPY data
    fpy_data : dict
        Dictionary of FPY data
    decay_data : dict
        Dictionary of decay data

        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 = gnd_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 = gnd_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 = gnd_name(Z, A, 0)
        if isotone in fpy_data:
            return isotone

    # If all else fails, use U235 yields
    return 'U235'

def _find_chain_file(cross_sections=None):
    # First check deprecated OPENMC_DEPLETE_CHAIN environment variable
    chain_file = os.environ.get("OPENMC_DEPLETE_CHAIN")
    if chain_file is not None:
        warn("Use of OPENMC_DEPLETE_CHAIN is deprecated in favor of adding "
             "depletion_chain to OPENMC_CROSS_SECTIONS", FutureWarning)
        return chain_file

    # Check for depletion chain in cross_sections.xml
    data = DataLibrary.from_xml(cross_sections)
    for lib in reversed(data.libraries):
        if lib['type'] == 'depletion_chain':
            return lib['path']

    raise DataError("No depletion chain specified and could not find depletion "
                    f"chain in {cross_sections}")

[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 the ``depletion_chain`` item in the :envvar:`OPENMC_CROSS_SECTIONS` environment variable. 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 : OrderedDict 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 = OrderedDict() 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[] = 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:`` instances can be included in these arguments. Parameters ---------- decay_files : list of str or List of ENDF decay sub-library files fpy_files : list of str or List of ENDF neutron-induced fission product yield sub-library files neutron_files : list of str or 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 = name = evaluation.gnd_name reactions[name] = {} for mf, mt, nc, mod in evaluation.reaction_list: if mf == 3: file_obj = StringIO(evaluation.section[3, mt]) q_value =[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 = # 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 = 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, 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) 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([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 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) if _have_lxml: tree.write(str(filename), encoding='utf-8', pretty_print=True) else: clean_indentation(root_elem) tree.write(str(filename), encoding='utf-8')
[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[] = dict(yield_obj) return out
[docs] def form_matrix(self, rates, fission_yields=None): """Forms depletion matrix. Parameters ---------- rates : numpy.ndarray 2D array indexed by (nuclide, reaction) fission_yields : dict, optional Option to use a custom set of fission yields. Expected to be of the form ``{parent : {product : f_yield}}`` with string nuclide names for ``parent`` and ``product``, and ``f_yield`` as the respective fission yield Returns ------- scipy.sparse.csr_matrix Sparse matrix representing depletion. See Also -------- :meth:`get_default_fission_yields` """ matrix = defaultdict(float) reactions = set() if fission_yields is None: fission_yields = self.get_default_fission_yields() for i, nuc in enumerate(self.nuclides): # Loss from radioactive decay if nuc.half_life is not None: decay_constant = math.log(2) / nuc.half_life if decay_constant != 0.0: matrix[i, i] -= decay_constant # Gain from radioactive decay if nuc.n_decay_modes != 0: for _, target, branching_ratio in nuc.decay_modes: # Allow for total annihilation for debug purposes if target is not None: branch_val = branching_ratio * decay_constant if branch_val != 0.0: k = self.nuclide_dict[target] matrix[k, i] += branch_val if in rates.index_nuc: # Extract all reactions for this nuclide in this cell nuc_ind = rates.index_nuc[] nuc_rates = rates[nuc_ind, :] for r_type, target, _, br in nuc.reactions: # Extract reaction index, and then final reaction rate r_id = rates.index_rx[r_type] path_rate = nuc_rates[r_id] # Loss term -- make sure we only count loss once for # reactions with branching ratios if r_type not in reactions: reactions.add(r_type) if path_rate != 0.0: matrix[i, i] -= path_rate # Gain term; allow for total annihilation for debug purposes if r_type != 'fission': if target is not None and path_rate != 0.0: k = self.nuclide_dict[target] matrix[k, i] += path_rate * br # Determine light nuclide production, e.g., (n,d) should # produce H2 light_nucs = REACTIONS[r_type].secondaries for light_nuc in light_nucs: k = self.nuclide_dict.get(light_nuc) if k is not None: matrix[k, i] += path_rate * br else: for product, y in fission_yields[].items(): yield_val = y * path_rate if yield_val != 0.0: k = self.nuclide_dict[product] matrix[k, i] += yield_val # Clear set of reactions reactions.clear() # Use DOK matrix as intermediate representation, then convert to CSR and return n = len(self) matrix_dok = sp.dok_matrix((n, n)) dict.update(matrix_dok, matrix) return matrix_dok.tocsr()
[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.branching_ratio if len(nuc_capt) > 0: capt[] = 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 not in secondary: indexes.append(ix) if "_m" not in grounds[parent] = 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 = gnd_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, previous = self[iso] new_nuclide = Nuclide( new_nuclide.half_life = previous.half_life new_nuclide.decay_energy = previous.decay_energy if hasattr(previous, '_fpy'): new_nuclide._fpy = previous._fpy for mode in previous.decay_modes: if 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 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([], 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