diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 814de12fc4..a23ba2bc1c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -10,7 +10,7 @@ import math import re from collections import OrderedDict, defaultdict from collections.abc import Mapping, Iterable -from numbers import Real +from numbers import Real, Integral from warnings import warn from openmc.checkvalue import check_type, check_greater_than @@ -500,7 +500,7 @@ class Chain(object): # Gain for _, target, branching_ratio in nuc.decay_modes: # Allow for total annihilation for debug purposes - if target != 'Nothing': + if target is not None: branch_val = branching_ratio * decay_constant if branch_val != 0.0: @@ -525,17 +525,16 @@ class Chain(object): matrix[i, i] -= path_rate # Gain term; allow for total annihilation for debug purposes - if target != 'Nothing': - if r_type != 'fission': - if path_rate != 0.0: - k = self.nuclide_dict[target] - matrix[k, i] += path_rate * br - else: - for product, y in fission_yields[nuc.name].items(): - yield_val = y * path_rate - if yield_val != 0.0: - k = self.nuclide_dict[product] - matrix[k, i] += yield_val + 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 + else: + for product, y in fission_yields[nuc.name].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() @@ -821,3 +820,168 @@ class Chain(object): return stat valid = valid and stat return valid + + 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) + + nuclides = [] + nuclide_dict = {} + reactions = set() + + 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 = new_nuclide.decay_energy + + new_decay = [] + for mode in previous.decay_modes: + if mode.target in all_isotopes: + new_decay.append(mode) + else: + new_decay.append(DecayTuple( + mode.type, None, mode.branching_ratio)) + new_nuclide.decay_modes = new_decay + + new_reactions = [] + for rxn in previous.reactions: + if rxn.target in all_isotopes: + new_reactions.append(rxn) + reactions.add(rxn.type) + elif rxn.type == "fission": + new_yields = new_nuclide.yield_data = ( + previous.yield_data.restrict_products(name_sort)) + if new_yields is not None: + new_reactions.append(rxn) + reactions.add("fission") + # Maintain total destruction rates but set no target + else: + new_reactions.append(ReactionTuple( + rxn.type, None, rxn.Q, rxn.branching_ratio)) + reactions.add(rxn.type) + + new_nuclide.reactions = new_reactions + + nuclides.append(new_nuclide) + nuclide_dict[iso] = idx + + new_chain = type(self)() + new_chain.nuclides = nuclides + new_chain.nuclide_dict = nuclide_dict + + # Doesn't appear that the ordering matters for the reactions, + # just the contents + new_chain.reactions = sorted(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" or rxn.target is None: + continue + # Skip if we've already come across this isotope + elif (rxn.target in next_iso + or rxn.target in found or rxn.target in isotopes): + continue + next_iso.add(rxn.target) + + 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 diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 0aae5adb72..7ba72b6d56 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -13,7 +13,7 @@ try: except ImportError: import xml.etree.ElementTree as ET -from numpy import empty +from numpy import empty, searchsorted from openmc.checkvalue import check_type @@ -30,8 +30,10 @@ Parameters ---------- type : str Type of the decay mode, e.g., 'beta-' -target : str - Nuclide resulting from decay +target : str or None + Nuclide resulting from decay. A value of ``None`` implies the + target does not exist in the currently configured depletion + chain branching_ratio : float Branching ratio of the decay mode @@ -53,8 +55,11 @@ Parameters ---------- type : str Type of the reaction, e.g., 'fission' -target : str - nuclide resulting from reaction +target : str or None + Nuclide resulting from reaction. A value of ``None`` + implies either no single target, e.g. from fission, + or that the target nuclide is not considered + in the current depletion chain Q : float Q value of the reaction in [eV] branching_ratio : float @@ -179,6 +184,8 @@ class Nuclide(object): for decay_elem in element.iter('decay'): d_type = decay_elem.get('type') target = decay_elem.get('target') + if target is not None and target.lower() == "nothing": + target = None branching_ratio = float(decay_elem.get('branching_ratio')) nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio)) @@ -192,6 +199,8 @@ class Nuclide(object): # just set null values if r_type != 'fission': target = reaction_elem.get('target') + if target is not None and target.lower() == "nothing": + target = None else: target = None if fission_q is not None: @@ -226,7 +235,7 @@ class Nuclide(object): for mode, daughter, br in self.decay_modes: mode_elem = ET.SubElement(elem, 'decay') mode_elem.set('type', mode) - mode_elem.set('target', daughter) + mode_elem.set('target', daughter or "Nothing") mode_elem.set('branching_ratio', str(br)) elem.set('reactions', str(len(self.reactions))) @@ -234,7 +243,7 @@ class Nuclide(object): rx_elem = ET.SubElement(elem, 'reaction') rx_elem.set('type', rx) rx_elem.set('Q', str(Q)) - if rx != 'fission': + if rx != 'fission' or daughter is not None: rx_elem.set('target', daughter) if br != 1.0: rx_elem.set('branching_ratio', str(br)) @@ -397,9 +406,7 @@ class FissionYieldDistribution(Mapping): for g_index, energy in enumerate(energies): prod_map = fission_yields[energy] for prod_ix, product in enumerate(ordered_prod): - yield_val = prod_map.get(product) - yield_matrix[g_index, prod_ix] = ( - 0.0 if yield_val is None else yield_val) + yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0) self.energies = tuple(energies) self.products = tuple(ordered_prod) self.yield_matrix = yield_matrix @@ -435,7 +442,7 @@ class FissionYieldDistribution(Mapping): FissionYieldDistribution """ all_yields = {} - for elem_index, yield_elem in enumerate(element.iter("fission_yields")): + for yield_elem in element.iter("fission_yields"): energy = float(yield_elem.get("energy")) products = yield_elem.find("products").text.split() yields = map(float, yield_elem.find("data").text.split()) @@ -460,6 +467,37 @@ class FissionYieldDistribution(Mapping): data_elem = ET.SubElement(yield_element, "data") data_elem.text = " ".join(map(str, yield_obj.yields)) + def restrict_products(self, possible_products): + """Return a new distribution with select products + + Parameters + ---------- + possible_products : iterable of str + Candidate pool of fission products. Existing products + not contained here will not exist in the new instance + + Returns + ------- + FissionYieldDistribution or None + A value of None indicates no values in + ``possible_products`` exist in :attr:`products` + + """ + + overlap = set(self.products).intersection(possible_products) + if not overlap: + return None + + products = sorted(overlap) + indices = searchsorted(self.products, products) + + # coerce back to dictionary to pass back to __init__ + new_yields = {} + for ene, yields in zip(self.energies, self.yield_matrix.copy()): + new_yields[ene] = dict(zip(products, yields[indices])) + + return type(self)(new_yields) + class FissionYield(Mapping): """Mapping for fission yields of a parent at a specific energy diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 8eb033783f..ce7e709d9c 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem. import sys import copy from collections import OrderedDict -from itertools import chain import os -import time import xml.etree.ElementTree as ET from warnings import warn -import h5py import numpy as np from uncertainties import ufloat @@ -113,6 +110,13 @@ class Operator(TransportOperator): ``fission_yield_mode``. Will be passed directly on to the helper. Passing a value of None will use the defaults for the associated helper. + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. Default is False. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. Attributes ---------- @@ -158,7 +162,8 @@ class Operator(TransportOperator): def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, energy_mode="fission-q", fission_q=None, dilute_initial=1.0e3, - fission_yield_mode="constant", fission_yield_opts=None): + fission_yield_mode="constant", fission_yield_opts=None, + reduce_chain=False, reduce_chain_level=None): if fission_yield_mode not in self._fission_helpers: raise KeyError( "fission_yield_mode must be one of {}, not {}".format( @@ -179,6 +184,16 @@ class Operator(TransportOperator): self.geometry = geometry self.diff_burnable_mats = diff_burnable_mats + # Reduce the chain before we create more materials + if reduce_chain: + all_isotopes = set() + for material in geometry.get_all_materials().values(): + if not material.depletable: + continue + for name, _dens_percent, _dens_type in material.nuclides: + all_isotopes.add(name) + self.chain = self.chain.reduce(all_isotopes, reduce_chain_level) + # Differentiate burnable materials with multiple instances if self.diff_burnable_mats: self._differentiate_burnable_mats() diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 08c60d3ed1..ba3be4e37a 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -6,7 +6,6 @@ from pathlib import Path from itertools import product import numpy as np -from openmc.data import zam, ATOMIC_SYMBOL from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram import pytest @@ -405,7 +404,8 @@ def test_fission_yield_attribute(simple_chain): dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1) with pytest.raises( ValueError, match="fission yield.*not equal.*compositions"): - cram.deplete(empty_chain, dummy_conc, None, 0.5) + cram.deplete(empty_chain, dummy_conc, None, 0.5) + def test_validate(simple_chain): """Test the validate method""" @@ -457,3 +457,95 @@ def test_validate_inputs(): with pytest.raises(ValueError, match="tolerance"): c.validate(tolerance=-1) + + +@pytest.fixture +def gnd_simple_chain(): + chainfile = Path(__file__).parents[1] / "chain_simple.xml" + return Chain.from_xml(chainfile) + + +def test_reduce(gnd_simple_chain): + ref_U5 = gnd_simple_chain["U235"] + ref_iodine = gnd_simple_chain["I135"] + ref_U5_yields = ref_U5.yield_data + + no_depth = gnd_simple_chain.reduce(["U235", "I135"], 0) + # We should get a chain just containing U235 and I135 + assert len(no_depth) == 2 + assert set(no_depth.reactions) == set(gnd_simple_chain.reactions) + + u5_round0 = no_depth["U235"] + assert u5_round0.n_decay_modes == ref_U5.n_decay_modes + for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes): + assert newmode.target is None + assert newmode.type == refmode.type + assert newmode.branching_ratio == refmode.branching_ratio + + assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths + for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions): + assert newrxn.target is None + assert newrxn.type == refrxn.type + assert newrxn.Q == refrxn.Q + assert newrxn.branching_ratio == refrxn.branching_ratio + + assert u5_round0.yield_data is not None + assert u5_round0.yield_data.products == ("I135",) + assert u5_round0.yield_data.yield_matrix == ( + ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")] + ) + + bareI5 = no_depth["I135"] + assert bareI5.n_decay_modes == ref_iodine.n_decay_modes + for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes): + assert newmode.target is None + assert newmode.type == refmode.type + assert newmode.branching_ratio == refmode.branching_ratio + + assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths + for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions): + assert newrxn.target is None + assert newrxn.type == refrxn.type + assert newrxn.Q == refrxn.Q + assert newrxn.branching_ratio == refrxn.branching_ratio + + follow_u5 = gnd_simple_chain.reduce(["U235"], 1) + u5_round1 = follow_u5["U235"] + assert u5_round1.decay_modes == ref_U5.decay_modes + assert u5_round1.reactions == ref_U5.reactions + assert u5_round1.yield_data is not None + assert ( + u5_round1.yield_data.yield_matrix == ref_U5_yields.yield_matrix + ).all() + + # Per the chain_simple.xml + # I135 -> Xe135 -> Cs135 + # I135 -> Xe136 + # No limit on depth + iodine_chain = gnd_simple_chain.reduce(["I135"]) + truncated_iodine = gnd_simple_chain.reduce(["I135"], 1) + assert len(iodine_chain) == 4 + assert len(truncated_iodine) == 3 + assert set(iodine_chain.nuclide_dict) == { + "I135", "Xe135", "Xe136", "Cs135"} + assert set(truncated_iodine.nuclide_dict) == {"I135", "Xe135", "Xe136"} + assert iodine_chain.reactions == ["(n,gamma)"] + assert iodine_chain["I135"].decay_modes == ref_iodine.decay_modes + assert iodine_chain["I135"].reactions == ref_iodine.reactions + for mode in truncated_iodine["Xe135"].decay_modes: + assert mode.target is None + + # Test that no FissionYieldDistribution is made if there are no + # fission products + u5_noyields = gnd_simple_chain.reduce(["U235"], 0)["U235"] + assert u5_noyields.yield_data is None + + # Check early termination if the eventual full chain + # is specified by using the iodine isotopes + new_iodine = gnd_simple_chain.reduce(set(iodine_chain.nuclide_dict)) + assert set(iodine_chain.nuclide_dict) == set(new_iodine.nuclide_dict) + + # Failure if some requested isotopes not in chain + + with pytest.raises(IndexError, match=".*not found.*Xx999"): + gnd_simple_chain.reduce(["U235", "Xx999"]) diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index b8a258df00..126206cc70 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -176,6 +176,21 @@ def test_fission_yield_distribution(): with pytest.raises(TypeError): orig_yields *= similar + # Test restriction of fission products + strict_restrict = yield_dist.restrict_products(["Xe135", "Sm149"]) + with_extras = yield_dist.restrict_products( + ["Xe135", "Sm149", "H1", "U235"]) + + assert strict_restrict.products == ("Sm149", "Xe135") + assert strict_restrict.energies == yield_dist.energies + assert with_extras.products == ("Sm149", "Xe135") + assert with_extras.energies == yield_dist.energies + for ene, new_yields in strict_restrict.items(): + for product in strict_restrict.products: + assert new_yields[product] == yield_dist[ene][product] + assert with_extras[ene][product] == yield_dist[ene][product] + + assert yield_dist.restrict_products(["U235"]) is None def test_validate():