From f3e7c7c623d2aba1a524b98a89467d1002b17fc1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 24 Jul 2019 11:58:54 -0500 Subject: [PATCH 01/45] Support for passing fission yields to Chain.form_matrix In support of using energy dependent fission yields, this commit allows a dictionary of fission yields to be passed into Chain.form_matrix. The dictionary is expected to be of the form {source_name: {target_name: yield_fraction}} where source_name and target_name are string GND names for fissionable nuclide and fission yield products, respectively. Currently, the fission yield dictionary does not have to be passed, defaulting to using the thermal yield values. This is done to make testing easier, and for back compatability. The goal of this feature, however, is to pass material [and thus spectrum] specific yields into this method. The test test_deplete_chain::test_form_matrix has been modified to pass the exact fission yield dictionary into Chain.form_matrix. The resulting matrix is compared against the matrix when no yields are passed. --- openmc/deplete/chain.py | 25 ++++++++++++++++++++++--- tests/unit_tests/test_deplete_chain.py | 7 +++++++ 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 1f16d9cafa..865f533ceb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -129,6 +129,7 @@ class Chain(object): self.nuclides = [] self.reactions = [] self.nuclide_dict = OrderedDict() + self._default_fsn_yields = None def __contains__(self, nuclide): return nuclide in self.nuclide_dict @@ -374,13 +375,27 @@ class Chain(object): clean_indentation(root_elem) tree.write(str(filename), encoding='utf-8') - def form_matrix(self, rates): + def _build_default_yields(self): + """Return dictionary {str: {str: float}}""" + # Take lowest energy for back compatability + # Should be removed by end of this feature + out = {} + for nuc in self.nuclides: + if len(nuc.yield_data) == 0: + continue + _energy, yield_data = sorted(nuc.yield_data.items())[0] + out[nuc.name] = {prod: frac for prod, frac in yield_data} + return out + + def form_matrix(self, rates, fission_yields=None): """Forms depletion matrix. Parameters ---------- rates : numpy.ndarray 2D array indexed by (nuclide, reaction) + fission_yields : dictionary of tuple to float, optional + Option to use a custom set of fission yields Returns ------- @@ -391,6 +406,11 @@ class Chain(object): matrix = defaultdict(float) reactions = set() + if fission_yields is None: + if self._default_fsn_yields is None: + self._default_fsn_yields = self._build_default_yields() + fission_yields = self._default_fsn_yields + for i, nuc in enumerate(self.nuclides): if nuc.n_decay_modes != 0: @@ -438,8 +458,7 @@ class Chain(object): # Assume that we should always use thermal fission # yields. At some point it would be nice to account # for the energy-dependence.. - energy, data = sorted(nuc.yield_data.items())[0] - for product, y in data: + for product, y in fission_yields[nuc.name].items(): yield_val = y * path_rate if yield_val != 0.0: k = self.nuclide_dict[product] diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 33b13b596d..df82eab1b0 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -3,6 +3,7 @@ from collections.abc import Mapping import os from pathlib import Path +from itertools import product import numpy as np from openmc.data import zam, ATOMIC_SYMBOL @@ -220,6 +221,12 @@ def test_form_matrix(simple_chain): assert mat[1, 2] == mat12 assert mat[2, 2] == mat22 + # Pass equivalent fission yields directly + # Ensure identical matrix is formed + f_yields = {"C": {"A": 0.0292737, "B": 0.002566345}} + new_mat = chain.form_matrix(react[0], f_yields) + for r, c in product(range(3), range(3)): + assert new_mat[r, c] == mat[r, c] def test_getitem(): """ Test nuc_by_ind converter function. """ From 72bf884d9bdfe62af065642207fe7fc6b79df9c5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 2 Aug 2019 17:00:38 -0500 Subject: [PATCH 02/45] Pass fission yields to depletion matrix_func Takes a single set of fission yields and passes them as an additional argument to matrix_func: >>> A = matrix_func(chain, rates, fission_yields) Applied to cf4, epc_rk4, celi, and leqi functions. Assumes that fission yields will not change during a depletion event. This change is probably overshadowed by how much the reaction rates may change, but still worth pointing out. --- openmc/deplete/_matrix_funcs.py | 83 +++++++++++++++++---------------- openmc/deplete/chain.py | 11 ++--- openmc/deplete/cram.py | 21 +++++++-- tests/dummy_operator.py | 69 ++++++++++++++------------- 4 files changed, 99 insertions(+), 85 deletions(-) diff --git a/openmc/deplete/_matrix_funcs.py b/openmc/deplete/_matrix_funcs.py index 1cacbd2c06..a606d739f7 100644 --- a/openmc/deplete/_matrix_funcs.py +++ b/openmc/deplete/_matrix_funcs.py @@ -1,77 +1,78 @@ """Functions to form the special matrix for depletion""" -def celi_f1(chain, rates): - return (5 / 12 * chain.form_matrix(rates[0]) - + 1 / 12 * chain.form_matrix(rates[1])) +def celi_f1(chain, rates, fission_yields=None): + return (5 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 12 * chain.form_matrix(rates[1], fission_yields)) -def celi_f2(chain, rates): - return (1 / 12 * chain.form_matrix(rates[0]) - + 5 / 12 * chain.form_matrix(rates[1])) +def celi_f2(chain, rates, fission_yields=None): + return (1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 5 / 12 * chain.form_matrix(rates[1], fission_yields)) -def cf4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def cf4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def cf4_f2(chain, rates): - return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1]) +def cf4_f2(chain, rates, fission_yields=None): + return (-1 / 2 * chain.form_matrix(rates[0], fission_yields) + + chain.form_matrix(rates[1], fission_yields)) -def cf4_f3(chain, rates): - return (1 / 4 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - - 1 / 12 * chain.form_matrix(rates[3])) +def cf4_f3(chain, rates, fission_yields=None): + return (1 / 4 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + - 1 / 12 * chain.form_matrix(rates[3], fission_yields)) -def cf4_f4(chain, rates): - return (-1 / 12 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - + 1 / 4 * chain.form_matrix(rates[3])) +def cf4_f4(chain, rates, fission_yields=None): + return (-1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + + 1 / 4 * chain.form_matrix(rates[3], fission_yields)) -def rk4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def rk4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def rk4_f4(chain, rates): - return (1 / 6 * chain.form_matrix(rates[0]) - + 1 / 3 * chain.form_matrix(rates[1]) - + 1 / 3 * chain.form_matrix(rates[2]) - + 1 / 6 * chain.form_matrix(rates[3])) +def rk4_f4(chain, rates, fission_yields=None): + return (1 / 6 * chain.form_matrix(rates[0], fission_yields) + + 1 / 3 * chain.form_matrix(rates[1], fission_yields) + + 1 / 3 * chain.form_matrix(rates[2], fission_yields) + + 1 / 6 * chain.form_matrix(rates[3], fission_yields)) -def leqi_f1(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f1(chain, inputs, fission_yields): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f2(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f2(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f3(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f3(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3) -def leqi_f4(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f4(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 2 * dt * dt_l + dt_l ** 2) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 865f533ceb..62158255a7 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -129,7 +129,6 @@ class Chain(object): self.nuclides = [] self.reactions = [] self.nuclide_dict = OrderedDict() - self._default_fsn_yields = None def __contains__(self, nuclide): return nuclide in self.nuclide_dict @@ -375,7 +374,7 @@ class Chain(object): clean_indentation(root_elem) tree.write(str(filename), encoding='utf-8') - def _build_default_yields(self): + def get_thermal_fission_yields(self): """Return dictionary {str: {str: float}}""" # Take lowest energy for back compatability # Should be removed by end of this feature @@ -383,8 +382,8 @@ class Chain(object): for nuc in self.nuclides: if len(nuc.yield_data) == 0: continue - _energy, yield_data = sorted(nuc.yield_data.items())[0] - out[nuc.name] = {prod: frac for prod, frac in yield_data} + yield_obj = nuc.yield_data[min(nuc.yield_energies)] + out[nuc.name] = dict(yield_obj) return out def form_matrix(self, rates, fission_yields=None): @@ -407,9 +406,7 @@ class Chain(object): reactions = set() if fission_yields is None: - if self._default_fsn_yields is None: - self._default_fsn_yields = self._build_default_yields() - fission_yields = self._default_fsn_yields + fission_yields = self.get_thermal_fission_yields() for i, nuc in enumerate(self.nuclides): diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index 41a4e4c710..d0c868004c 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -30,7 +30,9 @@ def deplete(chain, x, rates, dt, matrix_func=None): dt : float Time in [s] to deplete for maxtrix_func : Callable, optional - Function of two variables: ``chain`` and ``rates``. + Function to form the depletion matrix after calling + ``matrix_func(chain, rates, fission_yields)``, where + ``fission_yields = {parent: {product: yield_frac}}`` Expected to return the depletion matrix required by :func:`CRAM48`. @@ -40,9 +42,15 @@ def deplete(chain, x, rates, dt, matrix_func=None): Updated atom number vectors for each material """ + if not hasattr(chain, "fission_yields"): + fission_yields = repeat(chain.get_thermal_fission_yields()) + else: + fission_yields = chain.fission_yields + # Use multiprocessing pool to distribute work with Pool() as pool: - iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func)) + iters = zip(repeat(chain), x, rates, repeat(dt), + fission_yields, repeat(matrix_func)) x_result = list(pool.starmap(_cram_wrapper, iters)) return x_result @@ -67,7 +75,7 @@ def timed_deplete(*args, **kwargs): return time.time() - start, results -def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): +def _cram_wrapper(chain, n0, rates, dt, fission_yields, matrix_func=None): """Wraps depletion matrix creation / CRAM solve for multiprocess execution Parameters @@ -82,6 +90,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): Time to integrate to. maxtrix_func : function, optional Function to form the depletion matrix + fission_yields : dict + Single-energy fission yields of the form + ``{parent: {product: fission_yield}}`` Returns ------- @@ -90,9 +101,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): """ if matrix_func is None: - A = chain.form_matrix(rates) + A = chain.form_matrix(rates, fission_yields) else: - A = matrix_func(chain, rates) + A = matrix_func(chain, rates, fission_yields) return CRAM48(A, n0, dt) diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index c23abef631..acfdc5b583 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -6,6 +6,42 @@ from openmc.deplete.reaction_rates import ReactionRates from openmc.deplete.abc import TransportOperator, OperatorResult +class TestChain(object): + + @staticmethod + def get_thermal_fission_yields(): + return None + + def form_matrix(self, rates, _fission_yields=None): + """Forms the f(y) matrix in y' = f(y)y. + + Nominally a depletion matrix, this is abstracted on the off chance + that the function f has nothing to do with depletion at all. + + Parameters + ---------- + rates : numpy.ndarray + Slice of reaction rates for a single material + _fission_yields : optional + Not used + + Returns + ------- + scipy.sparse.csr_matrix + Sparse matrix representing f(y). + """ + + y_1 = rates[0, 0] + y_2 = rates[1, 0] + + a11 = np.sin(y_2) + a12 = np.cos(y_1) + a21 = -np.cos(y_2) + a22 = np.sin(y_1) + + return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) + + class DummyOperator(TransportOperator): """This is a dummy operator class with no statistical uncertainty. @@ -21,6 +57,7 @@ class DummyOperator(TransportOperator): """ def __init__(self, previous_results=None): self.prev_res = previous_results + self.chain = TestChain() def __call__(self, vec, power, print_out=False): """Evaluates F(y) @@ -52,38 +89,6 @@ class DummyOperator(TransportOperator): # Create a fake rates object return OperatorResult(ufloat(0.0, 0.0), reaction_rates) - @property - def chain(self): - return self - - def form_matrix(self, rates): - """Forms the f(y) matrix in y' = f(y)y. - - Nominally a depletion matrix, this is abstracted on the off chance - that the function f has nothing to do with depletion at all. - - Parameters - ---------- - rates : numpy.ndarray - Slice of reaction rates for a single material - - Returns - ------- - scipy.sparse.csr_matrix - Sparse matrix representing f(y). - """ - - y_1 = rates[0, 0] - y_2 = rates[1, 0] - - mat = np.zeros((2, 2)) - a11 = np.sin(y_2) - a12 = np.cos(y_1) - a21 = -np.cos(y_2) - a22 = np.sin(y_1) - - return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) - @property def volume(self): """ From 8cd71be8641790d0d7a96e0359b6eecd165f8176 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 2 Aug 2019 17:48:06 -0500 Subject: [PATCH 03/45] Store depletion fission yields with common yield matrix Add two new classes for working with fission yield data on the Chain. The primary class is FissionYieldDistribution. This class stores distributions for one nuclide with potentially many energies, with the assumption that the products don't change __too__ much across energy. Looking at the CASL chain and one generated by ENDF data, this appears to be the case. Most distributions are "full" in the sense that energies produce the same products. There are some cases where one or two products may not exist for a given energy. The FissionYieldDistribution retains a dictionary-like behavior, e.g. d[0.0253]["Xe135"] is a valid command and returns the yield for Xe-135 at a "thermal" spectrum, due to yields provided at 0.0253 eV. This is done with a helper class, _FissionYield, implemented first to support this behavior, but also to allow vector-operations in combining fission yields. These _FissionYield objects store the same product vector and a view into the underlying yield_matrix for a single energy. To support simple weighting of yields, the __mull__ and __iadd__ methods are implemented. A copy method is provided to remove the chance of modifying the original yield data. test_deplete_chain and test_deplete_nuclide have been modified in order to utilize these classes, without ruining the validity of the tests. Tests for the view / copy methods and vector operations on _FissionYield objects are included. --- openmc/deplete/nuclide.py | 201 +++++++++++++++++++++-- tests/unit_tests/test_deplete_chain.py | 6 +- tests/unit_tests/test_deplete_nuclide.py | 48 +++++- 3 files changed, 232 insertions(+), 23 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 88d5ac6eae..80ec648802 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -3,12 +3,19 @@ Contains the per-nuclide components of a depletion chain. """ +from numbers import Real from collections import namedtuple +from collections.abc import Iterable, Mapping + try: import lxml.etree as ET except ImportError: import xml.etree.ElementTree as ET +from numpy import asarray, fromiter, empty + +from openmc.checkvalue import check_type, check_length + DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio') DecayTuple.__doc__ = """\ @@ -83,7 +90,7 @@ class Nuclide(object): Maps tabulated energy to list of (product, yield) for all neutron-induced fission products. yield_energies : list of float - Energies at which fission product yiels exist + Energies at which fission product yields exist """ @@ -165,12 +172,7 @@ class Nuclide(object): fpy_elem = element.find('neutron_fission_yields') if fpy_elem is not None: - for yields_elem in fpy_elem.iter('fission_yields'): - E = float(yields_elem.get('energy')) - products = yields_elem.find('products').text.split() - yields = [float(y) for y in - yields_elem.find('data').text.split()] - nuc.yield_data[E] = list(zip(products, yields)) + nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem) nuc.yield_energies = list(sorted(nuc.yield_data.keys())) return nuc @@ -211,14 +213,181 @@ class Nuclide(object): fpy_elem = ET.SubElement(elem, 'neutron_fission_yields') energy_elem = ET.SubElement(fpy_elem, 'energies') energy_elem.text = ' '.join(str(E) for E in self.yield_energies) - - for E in self.yield_energies: - yields_elem = ET.SubElement(fpy_elem, 'fission_yields') - yields_elem.set('energy', str(E)) - - products_elem = ET.SubElement(yields_elem, 'products') - products_elem.text = ' '.join(x[0] for x in self.yield_data[E]) - data_elem = ET.SubElement(yields_elem, 'data') - data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E]) + self.yield_data.to_xml_element(fpy_elem) return elem + + +class FissionYieldDistribution(Mapping): + """Class for storing energy-dependent fission yields for a single nuclide + + Parameters + ---------- + ordered_energies : iterable of real + Energies for which fission yield data exist + orderded_products : iterable of str + Fission products produced by this parent at all energies + group_fission_yields : numpy.ndarray or iterable of iterable of float + Array of shape ``(n_energy, n_products)`` where + ``group_fission_yields[g][j]`` is the yield of + ``ordered_products[j]`` due to a fission in energy region ``g``. + + Attributes + ---------- + energies : tuple + Energies for which fission yields exist. Converted for + indexing + products : tuple + Fission products produced at all energies. Converted + for indexing + yield_matrix : numpy.ndarray + Array ``(n_energy, n_products)`` where + ``yield_matrix[g, j]`` is the fission yield of product + ``j`` for energy group ``g``. + + See Also + -------- + :meth:`from_xml`, :meth:`from_dict` + """ + + def __init__(self, ordered_energies, ordered_products, group_fission_yields): + check_type("energies", ordered_energies, Iterable, Real) + self.energies = tuple(ordered_energies) + check_type("products", ordered_products, Iterable, str) + self.products = tuple(ordered_products) + yield_matrix = asarray(group_fission_yields, dtype=float) + if yield_matrix.shape != (len(self.energies), len(self.products)): + raise ValueError( + "Shape of yield matrix inconsistent. " + "Should be ({}, {}), is {}".format( + len(energy_map), len(ordered_products), yield_matrix.shape)) + self.yield_matrix = yield_matrix + + def __len__(self): + return len(self.energies) + + def __getitem__(self, energy): + if energy not in self.energies: + raise KeyError(energy) + return _FissionYield( + self.products, self.yield_matrix[self.energies.index(energy)]) + + def __iter__(self): + return iter(self.energies) + + @classmethod + def from_xml_element(cls, element): + """Construct a distribution from a depletion chain xml file + + Parameters + ---------- + element : xml.etree.ElementTree.Element + XML element to pull fission yield data from + + Returns + ------- + FissionYieldDistribution + """ + yields = {} + for elem_index, yield_elem in enumerate(element.iter("fission_yields")): + energy = float(yield_elem.get("energy")) + products = yield_elem.find("products").text.split() + yield_mapobj = map(float, yield_elem.find("data").text.split()) + # Get a map of products to their corresponding yield + yields[energy] = dict(zip(products, yield_mapobj)) + + return cls.from_dict(yields) + + @classmethod + def from_dict(cls, fission_yields): + """Construct a distribution from a dictionary of yields + + Parameters + ----------- + fission_yields : dict + Dictionary ``{energy: {product: yield}}`` + + Returns + ------- + FissionYieldDistribution + """ + # mapping {energy: {product: value}} + energies = tuple(sorted(fission_yields)) + + # Get a consistent set of products to produce a matrix of yields + shared_prod = set() + for prod_set in map(set, fission_yields.values()): + shared_prod |= prod_set + ordered_prod = tuple(sorted(shared_prod)) + + yield_matrix = empty((len(energies), len(shared_prod))) + + for g_index, energy in enumerate(sorted(energies)): + prod_map = fission_yields[energy] + for prod_ix, product in enumerate(ordered_prod): + yield_val = prod_map.get(product) + if yield_val is None: + yield_matrix[g_index, prod_ix] = 0.0 + else: + yield_matrix[g_index, prod_ix] = yield_val + + return cls(energies, ordered_prod, yield_matrix) + + def to_xml_element(self, root): + """Write fission yield data to an xml element + + Parameters + ---------- + root : xml.etree.ElementTree.Element + Element to write distribution data to + """ + for energy, yield_obj in self.items(): + yield_element = ET.SubElement(root, "fission_yields") + yield_element.set("energy", str(energy)) + product_elem = ET.SubElement(yield_element, "products") + product_elem.text = " ".join(map(str, yield_obj.products)) + data_elem = ET.SubElement(yield_element, "data") + data_elem.text = " ".join(map(str, yield_obj.yields)) + + +class _FissionYield(Mapping): + """Mapping for fission yields of a parent at a specific energy + + Abstracted to support nested dictionary-like behavior for + :class:`FissionYieldDistribution`, and allowing math operations + on a single vector of yields + + Parameters + ---------- + products : tuple of str + Products for this specific distribution + yields : numpy.ndarray + View into associated :attr:`FissionYieldDistribution.yield_matrix` + """ + + def __init__(self, products, yields): + self.products = products + self.yields = yields + + def __getitem__(self, product): + if product not in self.products: + raise KeyError(product) + return self.yields[self.products.index(product)] + + def __len__(self): + return len(self.products) + + def __iter__(self): + return iter(self.products) + + def __iadd__(self, other): + """Increment value from other fission yield""" + self.yields += other.yields + return self + + def __mul__(self, value): + return _FissionYield(self.products, self.yields * value) + + def copy(self): + """Return an identical yield object, with unique yields""" + return _FissionYield(self.products, self.yields.copy()) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index df82eab1b0..d31933c34a 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -131,7 +131,8 @@ def test_from_xml(simple_chain): # Yield tests assert nuc.yield_energies == [0.0253] assert list(nuc.yield_data) == [0.0253] - assert nuc.yield_data[0.0253] == [("A", 0.0292737), ("B", 0.002566345)] + assert nuc.yield_data[0.0253].products == ("A", "B") + assert (nuc.yield_data[0.0253].yields == [0.0292737, 0.002566345]).all() def test_export_to_xml(run_in_tmpdir): @@ -163,7 +164,8 @@ def test_export_to_xml(run_in_tmpdir): nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3) ] C.yield_energies = [0.0253] - C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]} + C.yield_data = nuclide.FissionYieldDistribution.from_dict({ + 0.0253: {"A": 0.0292737, "B": 0.002566345}}) chain = Chain() chain.nuclides = [A, B, C] diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index f2a101d2a9..07f2edb85c 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -2,6 +2,8 @@ import xml.etree.ElementTree as ET +import numpy + from openmc.deplete import nuclide @@ -69,11 +71,10 @@ def test_from_xml(): nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0), nuclide.ReactionTuple('fission', None, 193405400.0, 1.0), ] + expected_yield_data = nuclide.FissionYieldDistribution.from_dict({ + 0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}}) assert u235.yield_energies == [0.0253] - assert u235.yield_data == { - 0.0253: [('Te134', 0.062155), ('Zr100', 0.0497641), - ('Xe138', 0.0481413)] - } + assert u235.yield_data == expected_yield_data def test_to_xml_element(): @@ -91,7 +92,8 @@ def test_to_xml_element(): nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0) ] C.yield_energies = [0.0253] - C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]} + C.yield_data = nuclide.FissionYieldDistribution.from_dict( + {0.0253: {"A": 0.0292737, "B": 0.002566345}}) element = C.to_xml_element() assert element.get("half_life") == "0.123" @@ -114,3 +116,39 @@ def test_to_xml_element(): assert float(rx_elems[1].get("Q")) == 0.0 assert element.find('neutron_fission_yields') is not None + + +def test_fission_yield_distribution(): + """Test an energy-dependent yield distribution""" + yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8, "Sm149": 2.69e-8}, + 5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}, # drop Sm149 + } + yield_dist = nuclide.FissionYieldDistribution.from_dict(yield_dict) + assert len(yield_dist) == len(yield_dict) + assert yield_dist.energies == tuple(sorted(yield_dict.keys())) + for exp_ene, exp_dist in yield_dict.items(): + act_dist = yield_dict[exp_ene] + for exp_prod, exp_yield in exp_dist.items(): + assert act_dist[exp_prod] == exp_yield + exp_yield_matrix = numpy.array([ + [4.08e-12, 1.71e-12, 7.85e-4], + [1.32e-12, 0.0, 1.12e-3], + [5.83e-8, 2.69e-8, 4.54e-3]]) + assert numpy.array_equal(yield_dist.yield_matrix, exp_yield_matrix) + + # Test the operations / special methods for fission yield + orig_yield_obj = yield_dist[0.0253] + # __getitem__ return yields as a view into yield matrix + assert orig_yield_obj.yields.base is yield_dist.yield_matrix + copied_yield = orig_yield_obj.copy() + # copied yields own their own memory -> not a view + assert copied_yield.yields.base is None + + # Fission yield feature uses scaled and incremented + mod_yields = orig_yield_obj * 2 + assert numpy.array_equal(orig_yield_obj.yields * 2, mod_yields.yields) + mod_yields += orig_yield_obj + assert numpy.array_equal(orig_yield_obj.yields * 3, mod_yields.yields) + From 599d473888be6e443cac4465987078a58eec9c19 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 2 Aug 2019 17:51:33 -0500 Subject: [PATCH 04/45] Add FissionYieldHelper for energy-dependent fission yields Operator now has a new helper for tallying fission rates and computing fission yields using energy dependency. The FissionYieldHelper creates fission rate tallies across burnable materials with an energy filter. The energy filter corresponds to each energy group for which there are fission yields. During the operator's unpacking of tallies, this helper computes a relative fission contribution by normalizing fission rates, such that the sum of all group fission rates in each material and each nuclide is unity. A single fission yield for each parent nuclide is computed by weighting the energy-dependent fission yields by these relative contributions. A nested dictionary {parent: {product: yield}} is added to a library, and then attached to the depletion chain at the end of the unpacking. These libraries are used during the depletion process instead of the only-thermal yields from before. A test was added that creates a proxy helper that doesn't use the C-API, but unpacks the "fission tally" and computes the fission yield library in an identical manner as the actual class. Since the proxy helper is a subclass of the one used in real transport, the functionality is consistent. --- openmc/deplete/chain.py | 3 - openmc/deplete/helpers.py | 193 ++++++++++++++++++++++- openmc/deplete/operator.py | 22 ++- tests/unit_tests/test_deplete_helpers.py | 99 ++++++++++++ 4 files changed, 309 insertions(+), 8 deletions(-) create mode 100644 tests/unit_tests/test_deplete_helpers.py diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 62158255a7..31b074bddb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -452,9 +452,6 @@ class Chain(object): k = self.nuclide_dict[target] matrix[k, i] += path_rate * br else: - # Assume that we should always use thermal fission - # yields. At some point it would be nice to account - # for the energy-dependence.. for product, y in fission_yields[nuc.name].items(): yield_val = y * path_rate if yield_val != 0.0: diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 66722c00a2..07c8e20242 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,11 +1,13 @@ """ Class for normalizing fission energy deposition """ +from collections import defaultdict from itertools import product -from numpy import dot, zeros +from numpy import dot, zeros, newaxis, divide, asarray -from openmc.capi import Tally, MaterialFilter +from openmc.capi import ( + Tally, MaterialFilter, EnergyFilter) from .abc import ReactionRateHelper, EnergyHelper # ------------------------------------- @@ -142,3 +144,190 @@ class ChainFissionHelper(EnergyHelper): isotopes in all materials have the same Q value. """ self._energy += dot(fission_rates, self._fission_q_vector) + + +# ------------------------------------ +# Helper for collapsing fission yields +# ------------------------------------ + + +class FissionYieldHelper(object): + """Class for using energy-dependent fission yields in depletion chain + + Creates a tally across all burnable materials to score the fission + rate in nuclides with yield data. An energy filter is used to + compute this rates in a group structure corresponding to the + fission yield data. This tally data is used to compute the + relative number of fission events in each energy region, + which serve as the weights for each energy-dependent fission + yield distribution. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + n_bmats : int + Number of burnable materials tracked in the problem + + Attributes + ---------- + energy_bounds : tuple of float + Sorted energy bounds from the tally filter + results : numpy.ndarray + Array of tally results for this process with shape + ``(n_local_mat, n_energy, n_nucs)`` + libraries : list of dict + List of fission yield dictionaries of the form + ``{parent: {product: yield}}``. Populated in + :meth:`compute_yields` and reset during + :meth:`unpack`. + """ + + def __init__(self, chain_nuclides, n_bmats): + self.libraries = [] + self._chain_nuclides = {} + self._chain_set = set() + self._tally_map = {} + # TODO Support user-requested minimum energy? + yield_energies = {0.0} + + # Get nuclides with fission yield data, names + # and all energy points + # Names are provided from operator tally nuclides + for nuc in chain_nuclides: + if len(nuc.yield_data) == 0: + continue + self._chain_nuclides[nuc.name] = nuc + self._chain_set.add(nuc.name) + yield_energies.update(nuc.yield_energies) + + # Create energy grid + self._energy_bounds = tuple(sorted(yield_energies)) + self.n_bmats = n_bmats + + self._reaction_tally = None + self.results = None + self.local_indexes = None + + @property + def energy_bounds(self): + return self._energy_bounds + + def generate_tallies(self, materials, mat_indexes): + """Construct the fission rate tally + + Parameters + ---------- + materials : iterable of C-API materials + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indexes for materials in ``materials`` tracked on this + process + """ + # Tally group-wise fission reaction rates + self._reaction_tally = Tally() + self._reaction_tally.scores = ['fission'] + + # Tally energy-weighted group-wise fission reaction rate + # Used to evaluated linear interpolation between fission yield points + self._weighted_reaction_tally = Tally() + self._weighted_reaction_tally.scores = ['fission'] + + filters = [ + MaterialFilter(materials), EnergyFilter(self._energy_bounds)] + + self._reaction_tally.filters = filters + self.local_indexes = asarray(mat_indexes) + + def set_fissionable_nuclides(self, nuclides): + """List of string of nuclides with data to be tallied + + Parameters + ---------- + nuclides : iterable of str + nuclides with non-zero densities that are candidates + for the fission tally. Not necessary that all are nuclides + with fission yields, but at least one must be + + Returns + ------- + nuclides : tuple of str + Nuclides ordered as they appear in the tally and in + the nuclide column of :attr:`results` + + Raises + ------ + ValueError + If no nuclides in ``nuclides`` are tracked on this + object + """ + # Set of all nuclides with positive density + # and fission yield data + nuc_set = self._chain_set & set(nuclides) + if len(nuc_set) == 0: + raise ValueError( + "No overlap between chain nuclides with fission yields and " + "requested tally nuclides") + nuclides = tuple(sorted(nuc_set)) + self._tally_index = [self._chain_nuclides[n] for n in nuclides] + self._reaction_tally.nuclides = nuclides + return nuclides + + def unpack(self): + """Unpack fission rate tallies to produce :attr:`results` + + Resets :attr:`libraries` under the assumption this is called + during the :class:`openmc.deplete.Operator` unpackign process + """ + # clear old libraries + self.libraries = [] + + # get view into tally results + # new shape: [material, energy, parent nuclide] + result_view = self._reaction_tally.results[..., 1].reshape( + self.n_bmats, len(self._energy_bounds) - 1, + len(self._reaction_tally.nuclides)) + + # Get results specific to this process + self.results = result_view[self.local_indexes, ...] + + # scale fission yields proportional to total fission rate + fsn_rate = self.results.sum(axis=1) + # TODO Guard against divide by zero + self.results /= fsn_rate[:, newaxis, :] + + def compute_yields(self, local_mat_index): + """Compute single fission yields using :attr:`results` + + Produces a new library in :attr:`self.libraries` + + Parameters + ---------- + local_mat_index : int + Index for material tracked on this process that + exists in :attr:`local_mat_index` and fits within + the first axis in :attr:`results` + """ + tally_results = self.results[local_mat_index] + + # Dictionary {parent_nuclide : [product, yield_vector]} + initial_library = {} + for i_energy, energy in enumerate(self._energy_bounds[1:]): + for i_nuc, fiss_frac in enumerate(tally_results[i_energy]): + parent = self._tally_index[i_nuc] + yield_data = parent.yield_data.get(energy) + if yield_data is None: + continue + if parent not in initial_library: + initial_library[parent] = yield_data * fiss_frac + continue + initial_library[parent] += yield_data * fiss_frac + + # convert to dictionary that can be passed to Chain.form_matrix + # {parent: {product: yield}} + library = {} + for k, yield_obj in initial_library.items(): + library[k.name] = dict(zip(yield_obj.products, yield_obj.yields)) + + self.libraries.append(library) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index f8288841f7..0c8c72987b 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -25,7 +25,8 @@ from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results_list import ResultsList -from .helpers import DirectReactionRateHelper, ChainFissionHelper +from .helpers import ( + DirectReactionRateHelper, ChainFissionHelper, FissionYieldHelper) def _distribute(items): @@ -177,6 +178,8 @@ class Operator(TransportOperator): self._rate_helper = DirectReactionRateHelper( self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() + self._fsn_yield_helper = FissionYieldHelper( + self.chain.nuclides, len(self.burnable_mats)) def __call__(self, vec, power): """Runs a simulation. @@ -204,8 +207,10 @@ class Operator(TransportOperator): # Update material compositions and tally nuclides self._update_materials() - self._rate_helper.nuclides = self._get_tally_nuclides() - self._energy_helper.nuclides = self._rate_helper.nuclides + nuclides = self._get_tally_nuclides() + self._rate_helper.nuclides = nuclides + self._energy_helper.nuclides = nuclides + self._fsn_yield_helper.set_fissionable_nuclides(nuclides) # Run OpenMC openmc.capi.reset() @@ -409,6 +414,10 @@ class Operator(TransportOperator): self._rate_helper.generate_tallies(materials, self.chain.reactions) self._energy_helper.prepare( self.chain.nuclides, self.reaction_rates.index_nuc, materials) + # Tell fission yield helper what materials this process is + # responsible for + self._fsn_yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) @@ -554,6 +563,7 @@ class Operator(TransportOperator): # Keep track of energy produced from all reactions in eV per source # particle self._energy_helper.reset() + self._fsn_yield_helper.unpack() # Create arrays to store fission Q values, reaction rates, and nuclide # numbers, zeroed out in material iteration @@ -576,6 +586,9 @@ class Operator(TransportOperator): tally_rates = self._rate_helper.get_material_rates( mat_index, nuc_ind, react_ind) + # Compute fission yields for this material + self._fsn_yield_helper.compute_yields(i) + # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) @@ -589,6 +602,9 @@ class Operator(TransportOperator): # Scale reaction rates to obtain units of reactions/sec rates *= power / energy + # Store new fission yields on the chain + self.chain.fission_yields = self._fsn_yield_helper.libraries + return OperatorResult(k_combined, rates) def _get_nuclides_with_data(self): diff --git a/tests/unit_tests/test_deplete_helpers.py b/tests/unit_tests/test_deplete_helpers.py new file mode 100644 index 0000000000..5db979098b --- /dev/null +++ b/tests/unit_tests/test_deplete_helpers.py @@ -0,0 +1,99 @@ +"""Test the Operator helpers""" + +from unittest.mock import Mock + +import pytest +import numpy +from numpy.testing import assert_array_equal + +from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution +from openmc.deplete.helpers import FissionYieldHelper + + +class FissionYieldHelperProxy(FissionYieldHelper): + + def __init__(self, chain_nuclides, n_bmats): + super().__init__(chain_nuclides, n_bmats) + self._reaction_tally = Mock() + self.local_indexes = numpy.array([0]) + + def generate_tallies(self, *args, **kwargs): + # Avoid calls to the C-API + pass + + +def test_fission_yield_helper(): + """Test the collection of fission yield data using approximated tallies + """ + u5yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} + u235 = Nuclide() + u235.name = "U235" + u235.yield_energies = (0.0253, 1.40e7) + u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) + + u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} + u238 = Nuclide() + u238.name = "U238" + u238.yield_energies = (5.00e5, ) + u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) + + xe135 = Nuclide() + xe135.name = "Xe135" + + n_bmats = 2 + helper = FissionYieldHelperProxy([u235, u238, xe135], n_bmats) + + assert helper.energy_bounds == (0, 0.0253, 5.00e5, 1.40e7) + + # test that tally must be created with nuclides with yields + with pytest.raises(ValueError, match="No overlap"): + helper.set_fissionable_nuclides(["Xe135", ]) + + act_nucs = helper.set_fissionable_nuclides(["U235", "U238", "Xe135"]) + assert act_nucs == ("U235", "U238") + + # Emulate getting tally data from transport run + # Test as if this Helper is responsible for one of two materials + # Tally results ordered [n_mat * n_ene, n_fiss_nuc, 3] + + tally_res = numpy.zeros((3 * n_bmats, 2, 3)) + u5_fiss_rates = numpy.array([1.0, 1.5, 2.0]) + u8_fiss_rates = numpy.array([0.0, 1.0, 1.0]) + tally_res[:3, 0, 1] = u5_fiss_rates + tally_res[:3, 1, 1] = u8_fiss_rates + + helper._reaction_tally.results = tally_res + helper.unpack() # compute yield fractions + + # Compare fraction fission rate from helper.reset + exp_results = numpy.empty((1, 3, 2)) + # Fraction of fission events in each energy range + u5_vec = u5_fiss_rates / u5_fiss_rates.sum() + u8_vec = u8_fiss_rates / u8_fiss_rates.sum() + exp_results[..., 0] = u5_vec + exp_results[..., 1] = u8_vec + assert_array_equal(helper.results, exp_results) + + # Compute and compare the library of fission yields + exp_lib = { + "U235": { + "Xe135": (u5yield_dict[0.0253]["Xe135"] * u5_vec[0] + + u5yield_dict[1.40e7]["Xe135"] * u5_vec[2]), + "Gd155": (u5yield_dict[0.0253]["Gd155"] * u5_vec[0] + + u5yield_dict[1.40e7]["Gd155"] * u5_vec[2]), + "Sm149": u5yield_dict[0.0253]["Sm149"] * u5_vec[0], + }, + "U238": { + "Xe135": u8yield_dict[5.00e5]["Xe135"] * u8_vec[1], + "Gd155": u8yield_dict[5.00e5]["Gd155"] * u8_vec[1], + } + } + + assert len(helper.libraries) == 0 + helper.compute_yields(0) + assert len(helper.libraries) == 1 + act_library = helper.libraries[0] + for parent, sub_yields in exp_lib.items(): + assert act_library[parent] == pytest.approx(sub_yields) From 27f2c6b15fe951104d48fda726e4ddc8764b2abb Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 5 Aug 2019 09:28:47 -0500 Subject: [PATCH 05/45] Guard against divide by zero in fission yield helper Create an empty matrix to fit energy-dependent fission rates for [materials, energies, nuclides]. Use the nonzero method on total fission rate to find indicies along the first and last axis for non-zero total fission rate. Populate corresponding fractional fission rates by dividing group fission rates by total fission rate, using the indices for materials and nuclides with non-zero total fission rate. Use numpy.where to find where total fission rate is zero, and directly set these locations to zero in the final result matrix. Also clean up some lint in openmc.deplete.nuclide related to old variable names and unused imports. --- openmc/deplete/helpers.py | 20 +++++++++++++------- openmc/deplete/nuclide.py | 7 ++++--- tests/unit_tests/test_deplete_nuclide.py | 1 - 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 07c8e20242..43a835908c 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,10 +1,9 @@ """ Class for normalizing fission energy deposition """ -from collections import defaultdict from itertools import product -from numpy import dot, zeros, newaxis, divide, asarray +from numpy import dot, zeros, newaxis, asarray, empty_like, where from openmc.capi import ( Tally, MaterialFilter, EnergyFilter) @@ -290,12 +289,19 @@ class FissionYieldHelper(object): len(self._reaction_tally.nuclides)) # Get results specific to this process - self.results = result_view[self.local_indexes, ...] + fission_rates = result_view[self.local_indexes] + self.results = empty_like(fission_rates) - # scale fission yields proportional to total fission rate - fsn_rate = self.results.sum(axis=1) - # TODO Guard against divide by zero - self.results /= fsn_rate[:, newaxis, :] + # scale group fission rates proportional to total fission rate + fission_total = fission_rates.sum(axis=1) + nz_mat, nz_nuc = fission_total.nonzero() + self.results[nz_mat, :, nz_nuc] = ( + fission_rates[nz_mat, :, nz_nuc] + / fission_total[nz_mat, newaxis, nz_nuc]) + + # directly set values to zero where total fission rate is zero + z_mat, z_nuc = where(fission_total == 0.0) + self.results[z_mat, :, z_nuc] = 0.0 def compute_yields(self, local_mat_index): """Compute single fission yields using :attr:`results` diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 80ec648802..8a0a9dc642 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -12,9 +12,9 @@ try: except ImportError: import xml.etree.ElementTree as ET -from numpy import asarray, fromiter, empty +from numpy import asarray, empty -from openmc.checkvalue import check_type, check_length +from openmc.checkvalue import check_type DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio') @@ -260,7 +260,8 @@ class FissionYieldDistribution(Mapping): raise ValueError( "Shape of yield matrix inconsistent. " "Should be ({}, {}), is {}".format( - len(energy_map), len(ordered_products), yield_matrix.shape)) + len(ordered_energies), len(ordered_products), + yield_matrix.shape)) self.yield_matrix = yield_matrix def __len__(self): diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 07f2edb85c..2a3993921f 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -151,4 +151,3 @@ def test_fission_yield_distribution(): assert numpy.array_equal(orig_yield_obj.yields * 2, mod_yields.yields) mod_yields += orig_yield_obj assert numpy.array_equal(orig_yield_obj.yields * 3, mod_yields.yields) - From c732336536669e5232fef19b0261fb2d941e0ef7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 6 Aug 2019 13:38:13 -0500 Subject: [PATCH 06/45] Avert bad numpy slices unpacking fission yields with no materials For the case where there are more processes than burnable materials, e.g. single pin cell or assembly with one burnable material, then the local_indexes attribute on some FissionYieldHelpers will be an empty array. This causes errors trying to use this as an indexer into the tally results. Instead, simply return from the unpacking if this is the case. The compute_yields function is called inside the Operator's unpacking of it's local materials, so that __should__ not be in a position to call compute_yields for a process with no burnable materials. --- openmc/deplete/helpers.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 43a835908c..8f35fac433 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -277,8 +277,14 @@ class FissionYieldHelper(object): """Unpack fission rate tallies to produce :attr:`results` Resets :attr:`libraries` under the assumption this is called - during the :class:`openmc.deplete.Operator` unpackign process + during the :class:`openmc.deplete.Operator` unpacking process """ + # if this process is not responsible for depleting anything + # [more processes than burnable materials] + # don't do anything + if self.local_indexes.size == 0: + return + # clear old libraries self.libraries = [] From 5fb8ec2d7aeb0af1d408898b6f57b620132bf39e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 8 Aug 2019 08:30:23 -0500 Subject: [PATCH 07/45] Document and test the default Chain fission yields --- openmc/deplete/chain.py | 26 +++++++++++++++++++++----- tests/unit_tests/test_deplete_chain.py | 8 ++++++++ 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 31b074bddb..2354158562 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -375,9 +375,19 @@ class Chain(object): tree.write(str(filename), encoding='utf-8') def get_thermal_fission_yields(self): - """Return dictionary {str: {str: float}}""" - # Take lowest energy for back compatability - # Should be removed by end of this feature + """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 = {} for nuc in self.nuclides: if len(nuc.yield_data) == 0: @@ -393,14 +403,20 @@ class Chain(object): ---------- rates : numpy.ndarray 2D array indexed by (nuclide, reaction) - fission_yields : dictionary of tuple to float, optional - Option to use a custom set of fission yields + 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_thermal_fission_yields` """ matrix = defaultdict(float) reactions = set() diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index d31933c34a..546fa99c3e 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -230,6 +230,7 @@ def test_form_matrix(simple_chain): for r, c in product(range(3), range(3)): assert new_mat[r, c] == mat[r, c] + def test_getitem(): """ Test nuc_by_ind converter function. """ chain = Chain() @@ -338,3 +339,10 @@ def test_capture_branch_failures(simple_chain): br = {"C": {"A": 1.0, "B": 1.0}} with pytest.raises(ValueError, match="C ratios"): simple_chain.set_capture_branches(br) + + +def test_simple_fission_yields(simple_chain): + """Check the default fission yields that can be used to form the matrix + """ + fission_yields = simple_chain.get_thermal_fission_yields() + assert fission_yields == {"C": {"A": 0.0292737, "B": 0.002566345}} From 3f42fc8505a37af535cf5786e14c38891da3a034 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 8 Aug 2019 08:50:32 -0500 Subject: [PATCH 08/45] Remove FissionYieldHelper.libraries; Return from compute_yields Previously the yields were appended into the libraries attribute. Now, the yields are returned directly and not appended at all. It is up to the caller to place these yields in the correct location. The Operator maintains a list of of the libraries that is updated through the unpacking method, and then set onto the Chain after working through all local materials --- openmc/deplete/helpers.py | 26 ++++++++---------------- openmc/deplete/operator.py | 7 +++++-- tests/unit_tests/test_deplete_helpers.py | 5 +---- 3 files changed, 15 insertions(+), 23 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 8f35fac433..ac81d0ac7a 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -163,7 +163,7 @@ class FissionYieldHelper(object): Parameters ---------- - chain_nuclides : iterable of openmc.deplete.nuclide + chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. n_bmats : int @@ -176,15 +176,9 @@ class FissionYieldHelper(object): results : numpy.ndarray Array of tally results for this process with shape ``(n_local_mat, n_energy, n_nucs)`` - libraries : list of dict - List of fission yield dictionaries of the form - ``{parent: {product: yield}}``. Populated in - :meth:`compute_yields` and reset during - :meth:`unpack`. """ def __init__(self, chain_nuclides, n_bmats): - self.libraries = [] self._chain_nuclides = {} self._chain_set = set() self._tally_map = {} @@ -245,7 +239,7 @@ class FissionYieldHelper(object): Parameters ---------- nuclides : iterable of str - nuclides with non-zero densities that are candidates + Nuclides with non-zero densities that are candidates for the fission tally. Not necessary that all are nuclides with fission yields, but at least one must be @@ -275,9 +269,6 @@ class FissionYieldHelper(object): def unpack(self): """Unpack fission rate tallies to produce :attr:`results` - - Resets :attr:`libraries` under the assumption this is called - during the :class:`openmc.deplete.Operator` unpacking process """ # if this process is not responsible for depleting anything # [more processes than burnable materials] @@ -285,9 +276,6 @@ class FissionYieldHelper(object): if self.local_indexes.size == 0: return - # clear old libraries - self.libraries = [] - # get view into tally results # new shape: [material, energy, parent nuclide] result_view = self._reaction_tally.results[..., 1].reshape( @@ -312,7 +300,7 @@ class FissionYieldHelper(object): def compute_yields(self, local_mat_index): """Compute single fission yields using :attr:`results` - Produces a new library in :attr:`self.libraries` + Produces a new library in :attr:`libraries` Parameters ---------- @@ -320,6 +308,11 @@ class FissionYieldHelper(object): Index for material tracked on this process that exists in :attr:`local_mat_index` and fits within the first axis in :attr:`results` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` """ tally_results = self.results[local_mat_index] @@ -337,9 +330,8 @@ class FissionYieldHelper(object): initial_library[parent] += yield_data * fiss_frac # convert to dictionary that can be passed to Chain.form_matrix - # {parent: {product: yield}} library = {} for k, yield_obj in initial_library.items(): library[k.name] = dict(zip(yield_obj.products, yield_obj.yields)) - self.libraries.append(library) + return library diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 0c8c72987b..d99c96f185 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -565,6 +565,9 @@ class Operator(TransportOperator): self._energy_helper.reset() self._fsn_yield_helper.unpack() + # Store fission yield dictionaries + fission_yields = [] + # Create arrays to store fission Q values, reaction rates, and nuclide # numbers, zeroed out in material iteration number = np.empty(rates.n_nuc) @@ -587,7 +590,7 @@ class Operator(TransportOperator): mat_index, nuc_ind, react_ind) # Compute fission yields for this material - self._fsn_yield_helper.compute_yields(i) + fission_yields.append(self._fsn_yield_helper.compute_yields(i)) # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) @@ -603,7 +606,7 @@ class Operator(TransportOperator): rates *= power / energy # Store new fission yields on the chain - self.chain.fission_yields = self._fsn_yield_helper.libraries + self.chain.fission_yields = fission_yields return OperatorResult(k_combined, rates) diff --git a/tests/unit_tests/test_deplete_helpers.py b/tests/unit_tests/test_deplete_helpers.py index 5db979098b..15ae1550ed 100644 --- a/tests/unit_tests/test_deplete_helpers.py +++ b/tests/unit_tests/test_deplete_helpers.py @@ -91,9 +91,6 @@ def test_fission_yield_helper(): } } - assert len(helper.libraries) == 0 - helper.compute_yields(0) - assert len(helper.libraries) == 1 - act_library = helper.libraries[0] + act_library = helper.compute_yields(0) for parent, sub_yields in exp_lib.items(): assert act_library[parent] == pytest.approx(sub_yields) From d10b83e06454b573e19e77375e46cb2a2040b423 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 8 Aug 2019 08:58:37 -0500 Subject: [PATCH 09/45] Document FissionYieldDistribution Pretty internal, but people may be curious --- docs/source/pythonapi/deplete.rst | 2 ++ openmc/deplete/nuclide.py | 14 +++++++++++--- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index a72504ac92..243dcad430 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -65,6 +65,7 @@ for a depletion chain: DecayTuple Nuclide ReactionTuple + FissionYieldDistribution The following classes are used during a depletion simulation and store auxiliary data, such as number densities and reaction rates for each material. @@ -77,6 +78,7 @@ data, such as number densities and reaction rates for each material. AtomNumber ChainFissionHelper DirectReactionRateHelper + FissionYieldHelper OperatorResult ReactionRates Results diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 8a0a9dc642..23f14331e2 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -221,9 +221,17 @@ class Nuclide(object): class FissionYieldDistribution(Mapping): """Class for storing energy-dependent fission yields for a single nuclide + Can be used as a dictionary mapping energies and products to fission + yields:: + + >>> fydist = FissionYieldDistribution.from_dict({ + ... {0.0253: {"Xe135": 0.021}}) + >>> fydist[0.0253]["Xe135"] + 0.021 + Parameters ---------- - ordered_energies : iterable of real + ordered_energies : iterable of float Energies for which fission yield data exist orderded_products : iterable of str Fission products produced by this parent at all energies @@ -247,7 +255,7 @@ class FissionYieldDistribution(Mapping): See Also -------- - :meth:`from_xml`, :meth:`from_dict` + :meth:`from_xml_element`, :meth:`from_dict` """ def __init__(self, ordered_energies, ordered_products, group_fission_yields): @@ -354,7 +362,7 @@ class FissionYieldDistribution(Mapping): class _FissionYield(Mapping): """Mapping for fission yields of a parent at a specific energy - Abstracted to support nested dictionary-like behavior for + Separated to support nested dictionary-like behavior for :class:`FissionYieldDistribution`, and allowing math operations on a single vector of yields From 726811394781a2b3b88440a3e4b36172aa596139 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 9 Aug 2019 17:43:34 -0500 Subject: [PATCH 10/45] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/deplete/nuclide.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 23f14331e2..7353cb9023 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -331,7 +331,7 @@ class FissionYieldDistribution(Mapping): yield_matrix = empty((len(energies), len(shared_prod))) - for g_index, energy in enumerate(sorted(energies)): + 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) From fab859982f867dc64f73678d38832cb84d435bc9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 12 Aug 2019 15:15:26 -0500 Subject: [PATCH 11/45] Document openmc.deplete.nuclide.FissionYield --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/nuclide.py | 81 +++++++++++++++++++++---------- 2 files changed, 56 insertions(+), 26 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 243dcad430..d53db063e4 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -66,6 +66,7 @@ for a depletion chain: Nuclide ReactionTuple FissionYieldDistribution + FissionYield The following classes are used during a depletion simulation and store auxiliary data, such as number densities and reaction rates for each material. diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 7353cb9023..86a10d4294 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -231,13 +231,15 @@ class FissionYieldDistribution(Mapping): Parameters ---------- - ordered_energies : iterable of float - Energies for which fission yield data exist - orderded_products : iterable of str - Fission products produced by this parent at all energies - group_fission_yields : numpy.ndarray or iterable of iterable of float + energies : iterable of float + Energies for which fission yield data exist. Must be ordered + by increasing energy + products : iterable of str + Fission products produced by this parent at all energies. + Must be ordered alphabetically + yield_matrix : numpy.ndarray or iterable of iterable of float Array of shape ``(n_energy, n_products)`` where - ``group_fission_yields[g][j]`` is the yield of + ``yield_matrix[g][j]`` is the yield of ``ordered_products[j]`` due to a fission in energy region ``g``. Attributes @@ -255,20 +257,21 @@ class FissionYieldDistribution(Mapping): See Also -------- - :meth:`from_xml_element`, :meth:`from_dict` + * :meth:`from_xml_element`, :meth:`from_dict` - Construction methods + * :class:`FissionYield` - Class used for storing yields at a given energy """ - def __init__(self, ordered_energies, ordered_products, group_fission_yields): - check_type("energies", ordered_energies, Iterable, Real) - self.energies = tuple(ordered_energies) - check_type("products", ordered_products, Iterable, str) - self.products = tuple(ordered_products) - yield_matrix = asarray(group_fission_yields, dtype=float) + def __init__(self, energies, products, yield_matrix): + check_type("energies", energies, Iterable, Real) + self.energies = tuple(energies) + check_type("products", products, Iterable, str) + self.products = tuple(products) + yield_matrix = asarray(yield_matrix, dtype=float) if yield_matrix.shape != (len(self.energies), len(self.products)): raise ValueError( "Shape of yield matrix inconsistent. " "Should be ({}, {}), is {}".format( - len(ordered_energies), len(ordered_products), + len(energies), len(products), yield_matrix.shape)) self.yield_matrix = yield_matrix @@ -278,7 +281,7 @@ class FissionYieldDistribution(Mapping): def __getitem__(self, energy): if energy not in self.energies: raise KeyError(energy) - return _FissionYield( + return FissionYield( self.products, self.yield_matrix[self.energies.index(energy)]) def __iter__(self): @@ -321,13 +324,13 @@ class FissionYieldDistribution(Mapping): FissionYieldDistribution """ # mapping {energy: {product: value}} - energies = tuple(sorted(fission_yields)) + energies = sorted(fission_yields) # Get a consistent set of products to produce a matrix of yields shared_prod = set() for prod_set in map(set, fission_yields.values()): shared_prod |= prod_set - ordered_prod = tuple(sorted(shared_prod)) + ordered_prod = sorted(shared_prod) yield_matrix = empty((len(energies), len(shared_prod))) @@ -335,10 +338,8 @@ class FissionYieldDistribution(Mapping): prod_map = fission_yields[energy] for prod_ix, product in enumerate(ordered_prod): yield_val = prod_map.get(product) - if yield_val is None: - yield_matrix[g_index, prod_ix] = 0.0 - else: - yield_matrix[g_index, prod_ix] = yield_val + yield_matrix[g_index, prod_ix] = ( + 0.0 if yield_val is None else yield_val) return cls(energies, ordered_prod, yield_matrix) @@ -359,19 +360,47 @@ class FissionYieldDistribution(Mapping): data_elem.text = " ".join(map(str, yield_obj.yields)) -class _FissionYield(Mapping): +class FissionYield(Mapping): """Mapping for fission yields of a parent at a specific energy Separated to support nested dictionary-like behavior for :class:`FissionYieldDistribution`, and allowing math operations - on a single vector of yields + on a single vector of yields. Can in turn be used like a + dictionary to fetch fission yields. + + Does not support resizing / inserting new products that do + not exist. Parameters ---------- products : tuple of str Products for this specific distribution yields : numpy.ndarray - View into associated :attr:`FissionYieldDistribution.yield_matrix` + Fission product yields for each product in ``products`` + + Attributes + ---------- + products : tuple of str + Products for this specific distribution + yields : numpy.ndarray + Fission product yields for each product in ``products`` + + Examples + -------- + >>> import numpy + >>> fy_vector = FissionYield( + ... ("Xe135", "I129", "Sm149"), + ... numpy.array((0.002, 0.001, 0.0003))) + >>> fy_vector["Xe135"] + 0.002 + >>> new = fy_vector.copy() + >>> fy_vector *= 2 + >>> fy_vector["Xe135"] + 0.004 + >>> new["Xe135"] + 0.002 + >>> (new + fy_vector)["Sm149"] + 0.0009 """ def __init__(self, products, yields): @@ -395,8 +424,8 @@ class _FissionYield(Mapping): return self def __mul__(self, value): - return _FissionYield(self.products, self.yields * value) + return FissionYield(self.products, self.yields * value) def copy(self): """Return an identical yield object, with unique yields""" - return _FissionYield(self.products, self.yields.copy()) + return FissionYield(self.products, self.yields.copy()) From 6d4d25acb2a53df609df76eb0a79ebb5e77b4804 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 15:16:41 -0500 Subject: [PATCH 12/45] Implement add, imul, rmul for openmc.deplete.FissionYield Addition methods __add__ and __iadd__ expect the other argument to a FissionYield object and will add the yields from one to the other [in place for __iadd__]. Multiplication methods __imul__, __mul__, and __rmul__ expect the other argument to be a scalar and will scale the fission yields by this value [in place for __imul__] Also added a __repr__ method that returns a dictionary-like representation --- openmc/deplete/nuclide.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 86a10d4294..615761cde9 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -367,6 +367,8 @@ class FissionYield(Mapping): :class:`FissionYieldDistribution`, and allowing math operations on a single vector of yields. Can in turn be used like a dictionary to fetch fission yields. + Supports multiplication of a scalar to scale the fission + yields and addition of another set of yields. Does not support resizing / inserting new products that do not exist. @@ -401,6 +403,8 @@ class FissionYield(Mapping): 0.002 >>> (new + fy_vector)["Sm149"] 0.0009 + >>> dict(new) + {"Xe135": 0.002, "I129": 0.001, "Sm149": 0.0003} """ def __init__(self, products, yields): @@ -418,13 +422,32 @@ class FissionYield(Mapping): def __iter__(self): return iter(self.products) + def __add__(self, other): + new = self.copy() + new += other + return new + def __iadd__(self, other): """Increment value from other fission yield""" self.yields += other.yields return self - def __mul__(self, value): - return FissionYield(self.products, self.yields * value) + def __imul__(self, scalar): + self.yields *= scalar + return self + + def __mul__(self, scalar): + new = self.copy() + new *= scalar + return new + + def __rmul__(self, scalar): + new = self.copy() + new *= scalar + return new + + def __repr__(self): + return "<{}: {}>".format(self.__class__.__name__, repr(dict(self))) def copy(self): """Return an identical yield object, with unique yields""" From 2cfbb21858283e78bfbfeab2ac8fc8fb5b6ec91a Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 15:23:19 -0500 Subject: [PATCH 13/45] Provide abstract FissionYieldHelper class API used by the Operator: - generate_tallies - weighted_yields [abstract] - update_nuclides_from_operator - unpack generate_tallies and unpack are empty methods, provided so that a consistent API can be found on helpers with tallies. Sorts nuclides into two dictionaries: those with a single set of fission yields, and those with multiple sets. The former can is exposed through the constant_yields property, returning a copy of the fission yield dictionary {parent: {product: yield}} The Operator now looks for/uses the weighted_yields and update_nuclides_from_operator methods instead of the old compute_yields and set_fissionable_nuclides methods --- docs/source/pythonapi/deplete.rst | 8 +-- openmc/deplete/abc.py | 102 ++++++++++++++++++++++++++++++ openmc/deplete/operator.py | 4 +- 3 files changed, 108 insertions(+), 6 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index d53db063e4..260838439c 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -6,7 +6,7 @@ .. module:: openmc.deplete -Several functions are provided that implement different time-integration +Several classes are provided that implement different time-integration algorithms for depletion calculations, which are described in detail in Colin Josey's thesis, `Development and analysis of high order neutron transport-depletion coupling algorithms `_. @@ -25,7 +25,7 @@ transport-depletion coupling algorithms `_. SICELIIntegrator SILEQIIntegrator -Each of these functions expects a "transport operator" to be passed. An operator +Each of these classes expects a "transport operator" to be passed. An operator specific to OpenMC is available using the following class: .. autosummary:: @@ -79,7 +79,6 @@ data, such as number densities and reaction rates for each material. AtomNumber ChainFissionHelper DirectReactionRateHelper - FissionYieldHelper OperatorResult ReactionRates Results @@ -94,8 +93,9 @@ The following classes are abstract classes that can be used to extend the :nosignatures: :template: myclass.rst - ReactionRateHelper EnergyHelper + FissionYieldHelper + ReactionRateHelper TransportOperator Custom integrators can be developed by subclassing from the following abstract diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 0f1203aeb8..356e9186df 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -364,6 +364,108 @@ class EnergyHelper(ABC): self._nuclides = nuclides +class FissionYieldHelper(ABC): + """Abstract class for processing energy dependent fission yields + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + + Attributes + ---------- + n_bmats : int + Number of burnable materials tracked in the problem + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + """ + + def __init__(self, chain_nuclides): + self._chain_nuclides = {} + self._constant_yields = {} + + # Get all nuclides with fission yield data + for nuc in chain_nuclides: + if len(nuc.yield_data) == 1: + self._constant_yields[nuc.name] = ( + nuc.yield_data[nuc.yield_energies[0]]) + elif len(nuc.yield_data) > 1: + self._chain_nuclides[nuc.name] = nuc + self._chain_set = set(self._chain_nuclides) | set(self._constant_yields) + + @property + def constant_yields(self): + out = {} + for key, sub in self._constant_yields.items(): + out[key] = sub.copy() + return out + + @abstractmethod + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + Parameters + ---------- + local_mat_index : int + Index for material tracked on this process that + exists in :attr:`local_mat_index` and fits within + the first axis in :attr:`results` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + + @staticmethod + def unpack(): + """Unpack tally data prior to compute fission yields. + + Called after a :meth:`openmc.deplete.Operator.__call__` + routine during the normalization of reaction rates. + + Not necessary for all subclasses to implement, unless tallies + are used. + """ + + @staticmethod + def generate_tallies(materials, mat_indexes): + """Construct tallies necessary for computing fission yields + + Called during the operator set up phase prior to depleting. + Not necessary for subclasses to implement + + Parameters + ---------- + materials : iterable of C-API materials + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indexes for materials in ``materials`` tracked on this + process + """ + + def update_nuclides_from_operator(self, nuclides): + """Return nuclides with non-zero densities and yield data + + Parameters + ---------- + nuclides : iterable of str + Nuclides with non-zero densities from the + :class:`openmc.deplete.Operator` + + Returns + ------- + nuclides : tuple of str + Union of nuclides that the :class:`openmc.deplete.Operator` + says have non-zero densities at this stage and those that + have yield data. Sorted by nuclide name + + """ + return tuple(sorted(self._chain_set & set(nuclides))) + + class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index d99c96f185..15f85aa4b2 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -210,7 +210,7 @@ class Operator(TransportOperator): nuclides = self._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._energy_helper.nuclides = nuclides - self._fsn_yield_helper.set_fissionable_nuclides(nuclides) + self._fsn_yield_helper.update_nuclides_from_operator(nuclides) # Run OpenMC openmc.capi.reset() @@ -590,7 +590,7 @@ class Operator(TransportOperator): mat_index, nuc_ind, react_ind) # Compute fission yields for this material - fission_yields.append(self._fsn_yield_helper.compute_yields(i)) + fission_yields.append(self._fsn_yield_helper.weighted_yields(i)) # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) From d8b9f5db54d28c6dc16330cd056e47eed90bedec Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 15:39:13 -0500 Subject: [PATCH 14/45] Add ConstantFissionYieldHelper concrete class Given a requested energy, will take fission yields from that energy on all nuclides and hold them as constant throughout the simulation. If the requested energy is not found for a specific nuclide, then the closest energy is used. The default is to take the thermal 0.0253 eV yields. This is now the helper attached to the Operator. The user can select what energy of yields they would like to use. --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/helpers.py | 173 +++++------------- openmc/deplete/operator.py | 12 +- .../unit_tests/test_deplete_fission_yields.py | 42 +++++ tests/unit_tests/test_deplete_helpers.py | 96 ---------- 5 files changed, 99 insertions(+), 225 deletions(-) create mode 100644 tests/unit_tests/test_deplete_fission_yields.py delete mode 100644 tests/unit_tests/test_deplete_helpers.py diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 260838439c..8e8cf755b8 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -78,6 +78,7 @@ data, such as number densities and reaction rates for each material. AtomNumber ChainFissionHelper + ConstantFissionYieldHelper DirectReactionRateHelper OperatorResult ReactionRates diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index ac81d0ac7a..44ed9293d0 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -2,12 +2,19 @@ Class for normalizing fission energy deposition """ from itertools import product +from numbers import Real from numpy import dot, zeros, newaxis, asarray, empty_like, where +from openmc.checkvalue import check_type, check_greater_than from openmc.capi import ( Tally, MaterialFilter, EnergyFilter) -from .abc import ReactionRateHelper, EnergyHelper +from .abc import ( + ReactionRateHelper, EnergyHelper, FissionYieldHelper) + +__all__ = ( + "DirectReactionRateHelper", "ChainFissionHelper", + "ConstantFissionYieldHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -150,152 +157,68 @@ class ChainFissionHelper(EnergyHelper): # ------------------------------------ -class FissionYieldHelper(object): - """Class for using energy-dependent fission yields in depletion chain - - Creates a tally across all burnable materials to score the fission - rate in nuclides with yield data. An energy filter is used to - compute this rates in a group structure corresponding to the - fission yield data. This tally data is used to compute the - relative number of fission events in each energy region, - which serve as the weights for each energy-dependent fission - yield distribution. +class ConstantFissionYieldHelper(FissionYieldHelper): + """Class that uses a single set of fission yields on each isotope Parameters ---------- chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int - Number of burnable materials tracked in the problem + energy : float, optional + Key in :attr:`openmc.deplete.Nuclide.yield_data` corresponding + to the desired set of fission yield data. Typically one of + ``{0.0253, 500000, 14000000}`` corresponding to 0.0253 eV, + 500 keV, and 14 MeV yield libraries. If the specific key is not + found, will fall back to closest energy present. + Default: 0.0253 eV for thermal yields Attributes ---------- - energy_bounds : tuple of float - Sorted energy bounds from the tally filter - results : numpy.ndarray - Array of tally results for this process with shape - ``(n_local_mat, n_energy, n_nucs)`` + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + energy : float + Energy of fission yield libraries. """ - def __init__(self, chain_nuclides, n_bmats): - self._chain_nuclides = {} - self._chain_set = set() - self._tally_map = {} - # TODO Support user-requested minimum energy? - yield_energies = {0.0} - - # Get nuclides with fission yield data, names - # and all energy points - # Names are provided from operator tally nuclides - for nuc in chain_nuclides: - if len(nuc.yield_data) == 0: + def __init__(self, chain_nuclides, energy=0.0253): + check_type("energy", energy, Real) + check_greater_than("energy", energy, 0.0, equality=True) + self._energy = energy + super().__init__(chain_nuclides) + # Iterate over all nuclides with > 1 set of yields + for name, nuc in self._chain_nuclides.items(): + yield_data = nuc.yield_data.get(energy) + if yield_data is not None: + self._constant_yields[name] = yield_data continue - self._chain_nuclides[nuc.name] = nuc - self._chain_set.add(nuc.name) - yield_energies.update(nuc.yield_energies) - - # Create energy grid - self._energy_bounds = tuple(sorted(yield_energies)) - self.n_bmats = n_bmats - - self._reaction_tally = None - self.results = None - self.local_indexes = None + # Specific energy not found, use closest energy + distances = [abs(energy - ene) for ene in nuc.yield_energies] + min_index = min( + range(len(nuc.yield_energies)), key=distances.__getitem__) + self._constant_yields[name] = ( + nuc.yield_data[nuc.yield_energies[min_index]]) @property - def energy_bounds(self): - return self._energy_bounds + def energy(self): + return self._energy - def generate_tallies(self, materials, mat_indexes): - """Construct the fission rate tally + def weighted_yields(self, _local_mat_index=None): + """Return fission yields for all nuclides requested Parameters ---------- - materials : iterable of C-API materials - Materials to be used in :class:`openmc.capi.MaterialFilter` - mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process - """ - # Tally group-wise fission reaction rates - self._reaction_tally = Tally() - self._reaction_tally.scores = ['fission'] - - # Tally energy-weighted group-wise fission reaction rate - # Used to evaluated linear interpolation between fission yield points - self._weighted_reaction_tally = Tally() - self._weighted_reaction_tally.scores = ['fission'] - - filters = [ - MaterialFilter(materials), EnergyFilter(self._energy_bounds)] - - self._reaction_tally.filters = filters - self.local_indexes = asarray(mat_indexes) - - def set_fissionable_nuclides(self, nuclides): - """List of string of nuclides with data to be tallied - - Parameters - ---------- - nuclides : iterable of str - Nuclides with non-zero densities that are candidates - for the fission tally. Not necessary that all are nuclides - with fission yields, but at least one must be + _local_mat_index : int, optional + Current material index. Not used since all yields are + constant Returns ------- - nuclides : tuple of str - Nuclides ordered as they appear in the tally and in - the nuclide column of :attr:`results` - - Raises - ------ - ValueError - If no nuclides in ``nuclides`` are tracked on this - object + library : dict + Dictionary of ``{parent: {product: fyield}}`` """ - # Set of all nuclides with positive density - # and fission yield data - nuc_set = self._chain_set & set(nuclides) - if len(nuc_set) == 0: - raise ValueError( - "No overlap between chain nuclides with fission yields and " - "requested tally nuclides") - nuclides = tuple(sorted(nuc_set)) - self._tally_index = [self._chain_nuclides[n] for n in nuclides] - self._reaction_tally.nuclides = nuclides - return nuclides - - def unpack(self): - """Unpack fission rate tallies to produce :attr:`results` - """ - # if this process is not responsible for depleting anything - # [more processes than burnable materials] - # don't do anything - if self.local_indexes.size == 0: - return - - # get view into tally results - # new shape: [material, energy, parent nuclide] - result_view = self._reaction_tally.results[..., 1].reshape( - self.n_bmats, len(self._energy_bounds) - 1, - len(self._reaction_tally.nuclides)) - - # Get results specific to this process - fission_rates = result_view[self.local_indexes] - self.results = empty_like(fission_rates) - - # scale group fission rates proportional to total fission rate - fission_total = fission_rates.sum(axis=1) - nz_mat, nz_nuc = fission_total.nonzero() - self.results[nz_mat, :, nz_nuc] = ( - fission_rates[nz_mat, :, nz_nuc] - / fission_total[nz_mat, newaxis, nz_nuc]) - - # directly set values to zero where total fission rate is zero - z_mat, z_nuc = where(fission_total == 0.0) - self.results[z_mat, :, z_nuc] = 0.0 + return self.constant_yields def compute_yields(self, local_mat_index): """Compute single fission yields using :attr:`results` diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 15f85aa4b2..660539af42 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -26,7 +26,7 @@ from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results_list import ResultsList from .helpers import ( - DirectReactionRateHelper, ChainFissionHelper, FissionYieldHelper) + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper) def _distribute(items): @@ -85,6 +85,10 @@ class Operator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. + fission_yield_energy : float, optional + Energy [eV] to pull fission product yields from. Passed to the + :class:`openmc.deplete.ConstantFissionYieldHelper`. Default: + 0.0253 eV Attributes ---------- @@ -123,7 +127,7 @@ class Operator(TransportOperator): """ def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, fission_q=None, - dilute_initial=1.0e3): + dilute_initial=1.0e3, fission_yield_energy=0.0253): super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None @@ -178,8 +182,8 @@ class Operator(TransportOperator): self._rate_helper = DirectReactionRateHelper( self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() - self._fsn_yield_helper = FissionYieldHelper( - self.chain.nuclides, len(self.burnable_mats)) + self._fsn_yield_helper = ConstantFissionYieldHelper( + self.chain.nuclides, energy=fission_yield_energy) def __call__(self, vec, power): """Runs a simulation. diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py new file mode 100644 index 0000000000..b9eaf9d2fb --- /dev/null +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -0,0 +1,42 @@ +"""Test the FissionYieldHelpers""" + +from collections import namedtuple + +import pytest + +from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution +from openmc.deplete.helpers import ConstantFissionYieldHelper + + +@pytest.fixture(scope="module") +def nuclide_bundle(): + u5yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} + u235 = Nuclide() + u235.name = "U235" + u235.yield_energies = (0.0253, 5.0e5, 1.40e7) + u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) + + u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} + u238 = Nuclide() + u238.name = "U238" + u238.yield_energies = (5.00e5, ) + u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) + + xe135 = Nuclide() + xe135.name = "Xe135" + + NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") + return NuclideBundle(u235, u238, xe135) +@pytest.mark.parametrize("input_energy, u5_yield_energy", ( + (0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) +def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): + helper = ConstantFissionYieldHelper( + nuclide_bundle, energy=input_energy) + assert helper.energy == input_energy + assert helper.constant_yields == { + "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], + "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal + assert helper.constant_yields == helper.weighted_yields(1) diff --git a/tests/unit_tests/test_deplete_helpers.py b/tests/unit_tests/test_deplete_helpers.py deleted file mode 100644 index 15ae1550ed..0000000000 --- a/tests/unit_tests/test_deplete_helpers.py +++ /dev/null @@ -1,96 +0,0 @@ -"""Test the Operator helpers""" - -from unittest.mock import Mock - -import pytest -import numpy -from numpy.testing import assert_array_equal - -from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution -from openmc.deplete.helpers import FissionYieldHelper - - -class FissionYieldHelperProxy(FissionYieldHelper): - - def __init__(self, chain_nuclides, n_bmats): - super().__init__(chain_nuclides, n_bmats) - self._reaction_tally = Mock() - self.local_indexes = numpy.array([0]) - - def generate_tallies(self, *args, **kwargs): - # Avoid calls to the C-API - pass - - -def test_fission_yield_helper(): - """Test the collection of fission yield data using approximated tallies - """ - u5yield_dict = { - 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, - 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} - u235 = Nuclide() - u235.name = "U235" - u235.yield_energies = (0.0253, 1.40e7) - u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) - - u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} - u238 = Nuclide() - u238.name = "U238" - u238.yield_energies = (5.00e5, ) - u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) - - xe135 = Nuclide() - xe135.name = "Xe135" - - n_bmats = 2 - helper = FissionYieldHelperProxy([u235, u238, xe135], n_bmats) - - assert helper.energy_bounds == (0, 0.0253, 5.00e5, 1.40e7) - - # test that tally must be created with nuclides with yields - with pytest.raises(ValueError, match="No overlap"): - helper.set_fissionable_nuclides(["Xe135", ]) - - act_nucs = helper.set_fissionable_nuclides(["U235", "U238", "Xe135"]) - assert act_nucs == ("U235", "U238") - - # Emulate getting tally data from transport run - # Test as if this Helper is responsible for one of two materials - # Tally results ordered [n_mat * n_ene, n_fiss_nuc, 3] - - tally_res = numpy.zeros((3 * n_bmats, 2, 3)) - u5_fiss_rates = numpy.array([1.0, 1.5, 2.0]) - u8_fiss_rates = numpy.array([0.0, 1.0, 1.0]) - tally_res[:3, 0, 1] = u5_fiss_rates - tally_res[:3, 1, 1] = u8_fiss_rates - - helper._reaction_tally.results = tally_res - helper.unpack() # compute yield fractions - - # Compare fraction fission rate from helper.reset - exp_results = numpy.empty((1, 3, 2)) - # Fraction of fission events in each energy range - u5_vec = u5_fiss_rates / u5_fiss_rates.sum() - u8_vec = u8_fiss_rates / u8_fiss_rates.sum() - exp_results[..., 0] = u5_vec - exp_results[..., 1] = u8_vec - assert_array_equal(helper.results, exp_results) - - # Compute and compare the library of fission yields - exp_lib = { - "U235": { - "Xe135": (u5yield_dict[0.0253]["Xe135"] * u5_vec[0] - + u5yield_dict[1.40e7]["Xe135"] * u5_vec[2]), - "Gd155": (u5yield_dict[0.0253]["Gd155"] * u5_vec[0] - + u5yield_dict[1.40e7]["Gd155"] * u5_vec[2]), - "Sm149": u5yield_dict[0.0253]["Sm149"] * u5_vec[0], - }, - "U238": { - "Xe135": u8yield_dict[5.00e5]["Xe135"] * u8_vec[1], - "Gd155": u8yield_dict[5.00e5]["Gd155"] * u8_vec[1], - } - } - - act_library = helper.compute_yields(0) - for parent, sub_yields in exp_lib.items(): - assert act_library[parent] == pytest.approx(sub_yields) From ca435b0bda717d1e4764a2907c42369019034d3b Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 16:10:55 -0500 Subject: [PATCH 15/45] Add TalliedFissionYieldHelper abstract helper Designed to be subclassed by helpers that aim to compute effective fission yields with tallied data. Constructs a tally in all burnable materials, scoring the fission rate in all nuclides that have non-zero density (as reported by the Operator) and have multiple sets of yields. For nuclides with non-zero densities and a single set of yields, those yields are assumed to be constant across incident neutron energies. The unpack method is now abstract, as each subclass should have its own data layout, some with different energy structures perhaps?? Maybe some weighting functions? --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/abc.py | 92 ++++++++++++++++++++++++++++++- 2 files changed, 92 insertions(+), 1 deletion(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 8e8cf755b8..7efc8c6393 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -97,6 +97,7 @@ The following classes are abstract classes that can be used to extend the EnergyHelper FissionYieldHelper ReactionRateHelper + TalliedFissionYieldHelper TransportOperator Custom integrators can be developed by subclassing from the following abstract diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 356e9186df..3e1dbde599 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -13,10 +13,11 @@ from copy import deepcopy from warnings import warn from numbers import Real, Integral -from numpy import nonzero, empty +from numpy import nonzero, empty, asarray from uncertainties import ufloat from openmc.data import DataLibrary, JOULE_PER_EV +from openmc.capi import MaterialFilter, Tally from openmc.checkvalue import check_type, check_greater_than from .results import Results from .chain import Chain @@ -466,6 +467,95 @@ class FissionYieldHelper(ABC): return tuple(sorted(self._chain_set & set(nuclides))) +class TalliedFissionYieldHelper(FissionYieldHelper): + """Abstract class for computing fission yields with tallies + + Generates a basic fission rate tally in all burnable materials with + :meth:`generate_tallies`, and set nuclides to be tallied with + :meth:`update_nuclides_from_operator`. Subclasses will need to implement + :meth:`unpack` and :meth:`weighted_yields`. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + n_bmats : int + Number of burnable materials tracked in the problem. + + Attributes + ---------- + n_bmats : int + Number of burnable materials tracked in the problem + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + """ + + _upper_energy = 20.0e6 # upper energy for tallies + + def __init__(self, chain_nuclides, n_bmats): + super().__init__(chain_nuclides) + self.n_bmats = n_bmats + self._local_indexes = None + self._fission_rate_tally = None + self._tally_index = {} + + def generate_tallies(self, materials, mat_indexes): + """Construct the fission rate tally + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indexes for materials in ``materials`` tracked on this + process + """ + self._local_indexes = asarray(mat_indexes) + + # Tally group-wise fission reaction rates + self._fission_rate_tally = Tally() + self._fission_rate_tally.scores = ['fission'] + + self._fission_rate_tally.filters = [MaterialFilter(materials)] + + def update_nuclides_from_operator(self, nuclides): + """Tally nuclides with non-zero density and multiple yields + + Parameters + ---------- + nuclides : iterable of str + Nuclides with non-zero densities from the + :class:`openmc.deplete.Operator` + + Returns + ------- + nuclides : tuple of str + Union of nuclides that the :class:`openmc.deplete.Operator` + says have non-zero densities at this stage and those that + have multiple sets of yield data. Sorted by nuclide name + """ + overlap = set(self._chain_nuclides).intersection(set(nuclides)) + if len(overlap) == 0: + # tally no nuclides, but keep the Tally alive + self._fission_rate_tally.nuclides = None + self._tally_index = [] + return tuple() + nuclides = tuple(sorted(overlap)) + self._tally_index = [self._chain_nuclides[n] for n in nuclides] + self._fission_rate_tally.nuclides = nuclides + return nuclides + + @abstractmethod + def unpack(self): + """Unpack tallies after a transport run. + + Abstract because each subclass will need to arrange its + tally data. + """ + + class Integrator(ABC): """Abstract class for solving the time-integration for depletion From f82dc83473cd27c80820077beb5b013444187aaf Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 16:38:24 -0500 Subject: [PATCH 16/45] Add FissionYieldCutoffHelper for weighting yields by a cutoff energy This helper performs the following tasks to weight fission yields: 1) Determine a set of "fast" and "thermal" yields based on a user specified cutoff energy, 2) Tally fission rates above and below the cutoff, and 3) Compute the effective yield by weighting the fast and thermal yields by the fraction of fast and thermal fissions for all nuclides in all burnable materials. The user is allowed to specify the cutoff energy and energies preferred for the fast and thermal yields. By default, the cutoff is 112 eV, chosen as it is the logarithmic mean of 0.0253 eV and 500 keV. These two values are the default energies for thermal and fast yields. Tests are added to check for failure in constructing this helper (thermal energy > cutoff, non-real values passed, etc.), selection of yields given a range of user data, and a proxy class that emulates a variety of thermal and fast splits to check the eventual computed yield libraries. --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/helpers.py | 189 +++++++++++++++--- .../unit_tests/test_deplete_fission_yields.py | 143 +++++++++++++ 3 files changed, 306 insertions(+), 27 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 7efc8c6393..44f8b70e05 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -80,6 +80,7 @@ data, such as number densities and reaction rates for each material. ChainFissionHelper ConstantFissionYieldHelper DirectReactionRateHelper + FissionYieldCutoffHelper OperatorResult ReactionRates Results diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 44ed9293d0..37d7359e83 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -3,18 +3,20 @@ Class for normalizing fission energy deposition """ from itertools import product from numbers import Real +import operator -from numpy import dot, zeros, newaxis, asarray, empty_like, where +from numpy import dot, zeros, newaxis from openmc.checkvalue import check_type, check_greater_than from openmc.capi import ( Tally, MaterialFilter, EnergyFilter) from .abc import ( - ReactionRateHelper, EnergyHelper, FissionYieldHelper) + ReactionRateHelper, EnergyHelper, FissionYieldHelper, + TalliedFissionYieldHelper) __all__ = ( "DirectReactionRateHelper", "ChainFissionHelper", - "ConstantFissionYieldHelper") + "ConstantFissionYieldHelper", "FissionYieldCutoffHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -220,41 +222,174 @@ class ConstantFissionYieldHelper(FissionYieldHelper): """ return self.constant_yields - def compute_yields(self, local_mat_index): - """Compute single fission yields using :attr:`results` - Produces a new library in :attr:`libraries` +class FissionYieldCutoffHelper(TalliedFissionYieldHelper): + """Helper that computes fission yields based on a cutoff energy + + Tally fission rates above and below the cutoff energy. + Assume that all fissions below cutoff energy have use thermal fission + product yield distributions, while all fissions above use a faster + set of yield distributions. + + Uses a limit of 20 MeV for tallying fission. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + n_bmats : int + Number of burnable materials tracked in the problem + cutoff : float, optional + Cutoff energy in [eV] below which all fissions will be + use thermal yields. All other fissions will use a + faster set of yields. Default: 112 [eV] + thermal_energy : float, optional + Energy of yield data corresponding to thermal yields. + Default: 0.0253 [eV] + fast_energy : float, optional + Energy of yield data corresponding to fast yields. + Default: 500 [kev] + + Attributes + ---------- + n_bmats : int + Number of burnable materials tracked in the problem + thermal_yields : dict + Dictionary of the form ``{parent: {product: yield}}`` + with thermal yields + fast_yields : dict + Dictionary of the form ``{parent: {product: yield}}`` + with fast yields + results : numpy.ndarray + Array of fission rate fractions with shape + ``(n_mats, 2, n_nucs)``. ``results[:, 0]`` + corresponds to the fraction of all fissions + that occured below ``cutoff``. The number + of materials in the first axis corresponds + to the number of materials burned by the + :class:`openmc.deplete.Operator` + """ + + def __init__(self, chain_nuclides, n_bmats, cutoff=112.0, + thermal_energy=0.0253, fast_energy=500.0e3): + check_type("cutoff", cutoff, Real) + check_type("thermal_energy", thermal_energy, Real) + check_type("fast_energy", fast_energy, Real) + check_greater_than("thermal_energy", thermal_energy, 0.0, equality=True) + check_greater_than("cutoff", cutoff, thermal_energy, equality=False) + check_greater_than("fast_energy", fast_energy, cutoff, equality=False) + super().__init__(chain_nuclides, n_bmats) + self._cutoff = cutoff + self._thermal_yields = {} + self._fast_yields = {} + for name, nuc in self._chain_nuclides.items(): + yields = nuc.yield_data + energies = nuc.yield_energies + thermal = yields.get(thermal_energy) + if thermal is None: + # find first index >= cutoff + ix = self._find_fallback_energy( + name, energies, cutoff, True) + thermal = yields[energies[ix - 1]] + fast = yields.get(fast_energy) + if fast is None: + # find first index <= cutoff + rev_ix = self._find_fallback_energy( + name, list(reversed(energies)), cutoff, False) + fast = yields[energies[-rev_ix]] + self._thermal_yields[name] = thermal + self._fast_yields[name] = fast + + @staticmethod + def _find_fallback_energy(name, energies, cutoff, check_under): + cutoff_func = operator.ge if check_under else operator.le + found = False + for ix, ene in enumerate(energies): + if cutoff_func(ene, cutoff): + found = True + break + if found and ix != 0: + return ix + domain = "thermal" if check_under else "fast" + raise ValueError("Could not find {} yields for {} " + "with cutoff {} eV".format(domain, name, cutoff)) + + def generate_tallies(self, materials, mat_indexes): + """Use C API to produce a fission rate tally in burnable materials + + Include a :class:`openmc.capi.EnergyFilter` to tally fission rates + above and below cutoff energy. + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indexes for materials in ``materials`` tracked on this + process + """ + super().generate_tallies(materials, mat_indexes) + energy_filter = EnergyFilter() + energy_filter.bins = (0.0, self._cutoff, self._upper_energy) + self._fission_rate_tally.filters.append(energy_filter) + + def unpack(self): + """Obtain fast and thermal fission fractions from tally""" + fission_rates = self._fission_rate_tally.results[..., 1].reshape( + self.n_bmats, 2, len(self._tally_index)) + self.results = fission_rates[self._local_indexes] + total_fission = self.results.sum(axis=1) + nz_mat, nz_nuc = total_fission.nonzero() + self.results[nz_mat, :, nz_nuc] /= total_fission[nz_mat, newaxis, nz_nuc] + + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + For nuclides with both yield data above and below + the cutoff energy, the effective yield for nuclide ``A`` + will be a weighted sum of fast and thermal yields. The + weights will be the fraction of ``A``s fission events + in the above and below the cutoff energy. + + If ``A`` has fission product distribution ``F`` + for fast fissions and ``T`` for thermal fissions, and + 70% of ``A``'s fissions are considered thermal, then + the effective fission product yield distributions + for ``A`` is ``0.7 * T + 0.3 * F`` Parameters ---------- local_mat_index : int - Index for material tracked on this process that - exists in :attr:`local_mat_index` and fits within - the first axis in :attr:`results` + Index for specific burnable material. Effective + yields will be produced using + ``self.results[local_mat_index]`` Returns ------- library : dict Dictionary of ``{parent: {product: fyield}}`` """ - tally_results = self.results[local_mat_index] + rates = self.results[local_mat_index] + yields = self.constant_yields + # iterate over thermal then fast yields, prefer __mul__ to __rmul__ + for therm_frac, nuc in zip(rates[0], self._tally_index): + yields[nuc.name] = self._thermal_yields[nuc.name] * therm_frac - # Dictionary {parent_nuclide : [product, yield_vector]} - initial_library = {} - for i_energy, energy in enumerate(self._energy_bounds[1:]): - for i_nuc, fiss_frac in enumerate(tally_results[i_energy]): - parent = self._tally_index[i_nuc] - yield_data = parent.yield_data.get(energy) - if yield_data is None: - continue - if parent not in initial_library: - initial_library[parent] = yield_data * fiss_frac - continue - initial_library[parent] += yield_data * fiss_frac + for fast_frac, nuc in zip(rates[1], self._tally_index): + yields[nuc.name] += self._fast_yields[nuc.name] * fast_frac + return yields - # convert to dictionary that can be passed to Chain.form_matrix - library = {} - for k, yield_obj in initial_library.items(): - library[k.name] = dict(zip(yield_obj.products, yield_obj.yields)) + @property + def thermal_yields(self): + out = {} + for key, sub in self._thermal_yields.items(): + out[key] = sub.copy() + return out - return library + @property + def fast_yields(self): + out = {} + for key, sub in self._fast_yields.items(): + out[key] = sub.copy() + return out diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index b9eaf9d2fb..f3911d80d0 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -1,5 +1,148 @@ """Test the FissionYieldHelpers""" +from collections import namedtuple +from unittest.mock import Mock + +import pytest +import numpy +from numpy.testing import assert_array_equal + +from openmc.capi import Material +from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution +from openmc.deplete.helpers import ( + FissionYieldCutoffHelper, ConstantFissionYieldHelper) + + +@pytest.fixture(scope="module") +def nuclide_bundle(): + u5yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} + u235 = Nuclide() + u235.name = "U235" + u235.yield_energies = (0.0253, 5.0e5, 1.40e7) + u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) + + u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} + u238 = Nuclide() + u238.name = "U238" + u238.yield_energies = (5.00e5, ) + u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) + + xe135 = Nuclide() + xe135.name = "Xe135" + + NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") + return NuclideBundle(u235, u238, xe135) +@pytest.mark.parametrize("input_energy, u5_yield_energy", ( + (0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) +def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): + helper = ConstantFissionYieldHelper( + nuclide_bundle, energy=input_energy) + assert helper.energy == input_energy + assert helper.constant_yields == { + "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], + "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal + assert helper.constant_yields == helper.weighted_yields(1) + + +# --------------------------------- +# Test the FissionYieldCutoffHelper +# --------------------------------- + + +def test_cutoff_construction(nuclide_bundle): + # defaults + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1) + assert helper.constant_yields == { + "U238": nuclide_bundle.u238.yield_data[5.0e5]} + assert helper.thermal_yields == { + "U235": nuclide_bundle.u235.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": nuclide_bundle.u235.yield_data[5e5]} + # use 14 MeV yields + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, fast_energy=14e6) + assert helper.constant_yields == { + "U238": nuclide_bundle.u238.yield_data[5.0e5]} + assert helper.thermal_yields == { + "U235": nuclide_bundle.u235.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": nuclide_bundle.u235.yield_data[14e6]} + # specify missing thermal yields -> use 0.0253 + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, thermal_energy=1) + assert helper.thermal_yields == { + "U235": nuclide_bundle.u235.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": nuclide_bundle.u235.yield_data[5e5]} + # request missing fast yields -> use epithermal + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, fast_energy=1e4) + assert helper.thermal_yields == { + "U235": nuclide_bundle.u235.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": nuclide_bundle.u235.yield_data[5e5]} + # test failures in cutoff: super low, super high + with pytest.raises(ValueError, match="thermal yields"): + FissionYieldCutoffHelper( + nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002) + with pytest.raises(ValueError, match="fast yields"): + FissionYieldCutoffHelper( + nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6) + + +@pytest.mark.parametrize("key", ("cutoff", "thermal_energy", "fast_energy")) +def test_cutoff_failure(key): + with pytest.raises(TypeError, match=key): + FissionYieldCutoffHelper(None, None, **{key: None}) + with pytest.raises(ValueError, match=key): + FissionYieldCutoffHelper(None, None, **{key: -1}) + + +class CutoffProxy(FissionYieldCutoffHelper): + """Proxy that supplies a set of tallies""" + + def generate_tallies(self, materials, mat_indexes): + self._fission_rate_tally = Mock() + self._local_indexes = numpy.asarray(mat_indexes) + + +# emulate some split between fast and thermal U235 fissions +@pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) +def test_cutoff_helper(nuclide_bundle, therm_frac): + materials = ["1", "2"] # TODO Use real C API materials? + n_bmats = len(materials) + local_mats = [0] + proxy = CutoffProxy(nuclide_bundle, n_bmats) + proxy.generate_tallies(materials, local_mats) + non_zero_nucs = [n.name for n in nuclide_bundle] + tally_nucs = proxy.update_nuclides_from_operator(non_zero_nucs) + assert tally_nucs == ("U235", ) + # Emulate building tallies + # material x energy, tallied_nuclides, 3 + proxy_flux = 1e6 + tally_data = numpy.empty((n_bmats * 2, 1, 3)) + tally_data[0, 0, 1] = therm_frac * proxy_flux + tally_data[1, 0, 1] = (1 - therm_frac) * proxy_flux + proxy._fission_rate_tally.results = tally_data + + proxy.unpack() + # expected results of shape (n_mats, 2, n_tnucs) + expected_results = numpy.empty((1, 2, 1)) + expected_results[:, 0] = therm_frac + expected_results[:, 1] = (1 - therm_frac) + assert proxy.results == pytest.approx(expected_results) + + actual_yields = proxy.weighted_yields(0) + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + assert actual_yields["U235"] == ( + proxy.thermal_yields["U235"] * therm_frac + + proxy.fast_yields["U235"] * (1 - therm_frac)) +"""Test the FissionYieldHelpers""" + from collections import namedtuple import pytest From ebf2f25ab1d7e55e5cef3089380af4ef0a7f7b2a Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 16:57:27 -0500 Subject: [PATCH 17/45] Make number of burnable materials optional for TalliedFissionYieldHelper Must be supplied prior to generating tallies. This makes the Operator's job a little easier when picking the fission yield mode. --- openmc/deplete/abc.py | 22 +++++++++++++++++++--- openmc/deplete/helpers.py | 5 +++-- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 3e1dbde599..7f6de6d62e 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -480,13 +480,14 @@ class TalliedFissionYieldHelper(FissionYieldHelper): chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int + n_bmats : int, optional Number of burnable materials tracked in the problem. Attributes ---------- n_bmats : int - Number of burnable materials tracked in the problem + Number of burnable materials tracked in the problem. + Must be set prior to generating tallies constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` @@ -494,13 +495,26 @@ class TalliedFissionYieldHelper(FissionYieldHelper): _upper_energy = 20.0e6 # upper energy for tallies - def __init__(self, chain_nuclides, n_bmats): + def __init__(self, chain_nuclides, n_bmats=None): super().__init__(chain_nuclides) self.n_bmats = n_bmats self._local_indexes = None self._fission_rate_tally = None self._tally_index = {} + @property + def n_bmats(self): + return self._n_bmats + + @n_bmats.setter + def n_bmats(self, value): + if value is None: + self._n_bmats = None + return + check_type("n_bmats", value, Integral) + check_greater_than("n_bmats", value, 0) + sef._n_bmats = value + def generate_tallies(self, materials, mat_indexes): """Construct the fission rate tally @@ -512,6 +526,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): Indexes for materials in ``materials`` tracked on this process """ + if self._n_bmat is None: + raise AttributeError("Number of burnable materials is not set") self._local_indexes = asarray(mat_indexes) # Tally group-wise fission reaction rates diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 37d7359e83..2dae7c0362 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -238,7 +238,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int + n_bmats : int, optional Number of burnable materials tracked in the problem cutoff : float, optional Cutoff energy in [eV] below which all fissions will be @@ -254,7 +254,8 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): Attributes ---------- n_bmats : int - Number of burnable materials tracked in the problem + Number of burnable materials tracked in the problem. + Must be set prior to generating tallies thermal_yields : dict Dictionary of the form ``{parent: {product: yield}}`` with thermal yields From 1d7f79e341c96a8c563abe38702b3a7e66c3c7ac Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 17:12:31 -0500 Subject: [PATCH 18/45] Determine fission yield helper based on Operator arguments A single argument, fission_yield_mode, is used to determine what type of FissionYield to use on the Operator. The options are * "constant": ConstantFissionYieldHelper [default] * "cutoff": FissionYieldCutoffHelper An additional argument, fission_yield_opts, can be supplied as a dictionary of keyword arguments to pass to the new helper. This allows the user to have control over how the private helper is created. --- openmc/deplete/helpers.py | 4 ++-- openmc/deplete/operator.py | 48 +++++++++++++++++++++++++++++++------- 2 files changed, 42 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 2dae7c0362..d6051c6dd3 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -350,12 +350,12 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): For nuclides with both yield data above and below the cutoff energy, the effective yield for nuclide ``A`` will be a weighted sum of fast and thermal yields. The - weights will be the fraction of ``A``s fission events + weights will be the fraction of ``A`` fission events in the above and below the cutoff energy. If ``A`` has fission product distribution ``F`` for fast fissions and ``T`` for thermal fissions, and - 70% of ``A``'s fissions are considered thermal, then + 70% of ``A`` fissions are considered thermal, then the effective fission product yield distributions for ``A`` is ``0.7 * T + 0.3 * F`` diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 660539af42..11fded8e44 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -26,7 +26,8 @@ from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results_list import ResultsList from .helpers import ( - DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper) + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, + FissionYieldCutoffHelper) def _distribute(items): @@ -85,10 +86,23 @@ class Operator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. - fission_yield_energy : float, optional - Energy [eV] to pull fission product yields from. Passed to the - :class:`openmc.deplete.ConstantFissionYieldHelper`. Default: - 0.0253 eV + fission_yield_mode : str, {"constant", "cutoff"} + Key indicating what fission product yield scheme to use. The + key determines what fission energy helper is used:: + + * "constant": :class:`openmc.deplete.ConstantFissionYieldHelper` + * "cutoff": :class:`openmc.deplete.FissionYieldCutoffHelper` + + The documentation on these classes describe their methodology + and differences. ``"constant"`` will treat fission yields as + constant with respect to energy, while ``"cutoff"`` will + compute effective yields based on the number of fission events + above and below a cutoff. Default: ``"constant"`` + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the helper determined by + ``fission_yield_mode``. Will be passed directly on to the + helper. Passing a value of None will use the defaults for + the associated helper. Attributes ---------- @@ -125,9 +139,19 @@ class Operator(TransportOperator): diff_burnable_mats : bool Whether to differentiate burnable materials with multiple instances """ + _fission_helpers_ = { + "constant": ConstantFissionYieldHelper, + "cutoff": FissionYieldCutoffHelper, + } + def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, fission_q=None, - dilute_initial=1.0e3, fission_yield_energy=0.0253): + dilute_initial=1.0e3, fission_yield_mode="constant", + fission_yield_opts=None): + if fission_yield_mode not in self._fission_helpers_: + raise KeyError( + "fission_yield_mode must be one of {}, not {}".format( + ", ".join(self._fission_helpers_), fission_yield_mode)) super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None @@ -182,8 +206,16 @@ class Operator(TransportOperator): self._rate_helper = DirectReactionRateHelper( self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() - self._fsn_yield_helper = ConstantFissionYieldHelper( - self.chain.nuclides, energy=fission_yield_energy) + + # Select and create fission yield helper + fission_helper = self._fission_helpers_[fission_yield_mode] + if fission_yield_opts is None: + self._fsn_yield_helper = fission_helper(self.chain.nuclides) + else: + self._fsn_yield_helper = fission_yield_mode( + self.chain.nuclides, **fission_yield_opts) + if hasattr(self._fsn_yield_helper, "n_bmats"): + self._fsn_yield_helper.n_bmats = len(self.burnable_mats) def __call__(self, vec, power): """Runs a simulation. From 46ce147a81c963f23150c9bf4b1d77baf520eb67 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 17:47:34 -0500 Subject: [PATCH 19/45] Revert ebf2f25ab and add from_operator class method TalliedFissionYieldHelper instances must be passed the number of burnable materials at construction again. To make things easier for the Operator, a from_operator class method has been added to the FissionYieldHelper and TalliedFissionYieldHelper abstract classes, as well as more explicit versions on the ConstantFissionYieldHelper and FissionYieldCutoffHelper. All have the call signature cls.from_operator(operator, **kwargs) The concrete classes have explicitely declared what the allowed key word arguments are, while the abstract classes forward kwargs directly onto __init__. All explicit keyword arguments should have a counterpart in the __init__ method. --- openmc/deplete/abc.py | 51 ++++++++++++++++++++++++-------------- openmc/deplete/helpers.py | 48 +++++++++++++++++++++++++++++++++++ openmc/deplete/operator.py | 11 +++----- 3 files changed, 85 insertions(+), 25 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 7f6de6d62e..5efda87029 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -466,6 +466,20 @@ class FissionYieldHelper(ABC): """ return tuple(sorted(self._chain_set & set(nuclides))) + @classmethod + def from_operator(cls, operator, **kwargs): + """Create a new instance by pulling data from the operator + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator with a depletion chain + kwargs: optional + Optional keyword arguments to be passed to the underlying + ``__init__`` method + """ + return cls(operator.chain.nuclides, **kwargs) + class TalliedFissionYieldHelper(FissionYieldHelper): """Abstract class for computing fission yields with tallies @@ -480,14 +494,13 @@ class TalliedFissionYieldHelper(FissionYieldHelper): chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int, optional + n_bmats : int Number of burnable materials tracked in the problem. Attributes ---------- n_bmats : int Number of burnable materials tracked in the problem. - Must be set prior to generating tallies constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` @@ -495,26 +508,13 @@ class TalliedFissionYieldHelper(FissionYieldHelper): _upper_energy = 20.0e6 # upper energy for tallies - def __init__(self, chain_nuclides, n_bmats=None): + def __init__(self, chain_nuclides, n_bmats): super().__init__(chain_nuclides) self.n_bmats = n_bmats self._local_indexes = None self._fission_rate_tally = None self._tally_index = {} - @property - def n_bmats(self): - return self._n_bmats - - @n_bmats.setter - def n_bmats(self, value): - if value is None: - self._n_bmats = None - return - check_type("n_bmats", value, Integral) - check_greater_than("n_bmats", value, 0) - sef._n_bmats = value - def generate_tallies(self, materials, mat_indexes): """Construct the fission rate tally @@ -526,8 +526,6 @@ class TalliedFissionYieldHelper(FissionYieldHelper): Indexes for materials in ``materials`` tracked on this process """ - if self._n_bmat is None: - raise AttributeError("Number of burnable materials is not set") self._local_indexes = asarray(mat_indexes) # Tally group-wise fission reaction rates @@ -571,6 +569,23 @@ class TalliedFissionYieldHelper(FissionYieldHelper): tally data. """ + @classmethod + def from_operator(cls, operator, **kwargs): + """Create a new instance by pulling data from the operator + + Require the more concrete :class:`openmc.deplete.Operator` + because this needs knowledge of burnable materials + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator with a depletion chain + kwargs: optional + Optional keyword arguments to be passed to the underlying + ``__init__`` method + """ + return cls(operator.chain.nuclides, len(operator.burnable_mats), + **kwargs) class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index d6051c6dd3..b66c9ba5e7 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -202,6 +202,24 @@ class ConstantFissionYieldHelper(FissionYieldHelper): self._constant_yields[name] = ( nuc.yield_data[nuc.yield_energies[min_index]]) + @classmethod + def from_operator(cls, operator, energy=0.0253): + """Return a new ConstantFissionYieldHelper using operator data + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + operator with a depletion chain + energy : float, optional + Energy for default fission yield libraries for nuclides with + multiple sets of yield data + + Returns + ------- + ConstantFissionYieldHelper + """ + return cls(operator.chain.nuclides, energy=energy) + @property def energy(self): return self._energy @@ -302,6 +320,36 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): self._thermal_yields[name] = thermal self._fast_yields[name] = fast + @classmethod + def from_operator(cls, operator, cutoff=112.0, + thermal_energy=0.0253, fast_energy=500e3): + """Construct a helper from an operator + + All keyword arguments are identical to their counterpart + in the main ``__init__`` method + + Parameters + ---------- + operator : openmc.deplete.Operator + Operator with a chain and burnable materials + cutoff : float, optional + Cutoff energy for tallying fast and thermal fissions + thermal_energy : float, optional + Energy to use when pulling thermal fission yields from + nuclides with multiple sets of yields + fast_energy : float, optional + Energy to use when pulling fast fission yields from + nuclides with multiple sets of yields + + Returns + ------- + FissionYieldCutoffHelper + + """ + return cls(operator.chain.nuclides, len(operator.burnable_mats), + cutoff=cutoff, thermal_energy=thermal_energy, + fast_energy=fast_energy) + @staticmethod def _find_fallback_energy(name, energies, cutoff, check_under): cutoff_func = operator.ge if check_under else operator.le diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 11fded8e44..e967dc85d1 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -209,13 +209,10 @@ class Operator(TransportOperator): # Select and create fission yield helper fission_helper = self._fission_helpers_[fission_yield_mode] - if fission_yield_opts is None: - self._fsn_yield_helper = fission_helper(self.chain.nuclides) - else: - self._fsn_yield_helper = fission_yield_mode( - self.chain.nuclides, **fission_yield_opts) - if hasattr(self._fsn_yield_helper, "n_bmats"): - self._fsn_yield_helper.n_bmats = len(self.burnable_mats) + fission_yield_opts = ( + {} if fission_yield_opts is None else fission_yield_opts) + self._fsn_yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) def __call__(self, vec, power): """Runs a simulation. From 0944503464aec68046a9407dd19078286d937fdf Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 09:47:36 -0500 Subject: [PATCH 20/45] Document Chain.fission_yields attribute The user is allowed to pass a single dictionary of yields to be applied to all materials. Otherwise, each element in the iterable should correspond to a material ordered by the TransportOperator. Sparse checking is performed on the supplied input. Validation is possible, but could become a burden for potentially many burnable materials representing a lengthy list of nested dictionaries to dig through and validate. --- openmc/deplete/chain.py | 31 ++++++++++++++++++++++++-- tests/unit_tests/test_deplete_chain.py | 13 +++++++++++ 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 2354158562..a204e4821c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -9,7 +9,8 @@ from itertools import chain import math import re from collections import OrderedDict, defaultdict -from collections.abc import Mapping +from collections.abc import Mapping, Iterable +from numbers import Real from warnings import warn from openmc.checkvalue import check_type, check_less_than @@ -122,13 +123,22 @@ class Chain(object): 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: Real}}``, 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 @@ -642,3 +652,20 @@ class Chain(object): new_ratios[ground_tgt] = ground_br parent.reactions.append(ReactionTuple( "(n,gamma)", ground_tgt, capt_Q, ground_br)) + + @property + def fission_yields(self): + if self._fission_yields is None: + self._fission_yields = [self.get_thermal_fission_yields()] + return self._fission_yields + + @fission_yields.setter + def fission_yields(self, yields): + if yields is None: + self._fission_yields = None + return + if isinstance(yields, Mapping): + self._fission_yields = [yields] + return + check_type("fission_yields", yields, Iterable, Mapping) + self._fission_yields = yields diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 546fa99c3e..2e5ad09832 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -346,3 +346,16 @@ def test_simple_fission_yields(simple_chain): """ fission_yields = simple_chain.get_thermal_fission_yields() assert fission_yields == {"C": {"A": 0.0292737, "B": 0.002566345}} + + +def test_fission_yield_attribute(simple_chain): + """Test the fission_yields property""" + thermal_yields = simple_chain.get_thermal_fission_yields() + # generate default with property + assert simple_chain.fission_yields[0] == thermal_yields + empty_chain = Chain() + empty_chain.fission_yields = thermal_yields + assert empty_chain.fission_yields[0] == thermal_yields + empty_chain.fission_yields = [thermal_yields] * 2 + assert empty_chain.fission_yields[0] == thermal_yields + assert empty_chain.fission_yields[1] == thermal_yields From c1d66bc0226b721620e1aad78b684074acb0c080 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 10:49:46 -0500 Subject: [PATCH 21/45] Ensure # fission yields == # burnable materials in cram.deplete If the fission yields are a single entry, then itertools.repeat is used to apply the fission yields to all burnable materials. Otherwise, check that the number of fission yield libraries is equal to the number of materials to be burned. --- openmc/deplete/cram.py | 12 ++++++++---- openmc/deplete/nuclide.py | 21 +++++++++++---------- tests/unit_tests/test_deplete_chain.py | 9 ++++++++- 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index d0c868004c..184c411bf3 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -42,10 +42,14 @@ def deplete(chain, x, rates, dt, matrix_func=None): Updated atom number vectors for each material """ - if not hasattr(chain, "fission_yields"): - fission_yields = repeat(chain.get_thermal_fission_yields()) - else: - fission_yields = chain.fission_yields + fission_yields = chain.fission_yields + if len(fission_yields) == 1: + fission_yields = repeat(fission_yields[0]) + elif len(fission_yields) != len(x): + raise ValueError( + "Number of material fission yield distributions {} is not equal " + "to the number of compositions {}".format(len(fission_yields), + len(x))) # Use multiprocessing pool to distribute work with Pool() as pool: diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 615761cde9..5217831f49 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -86,9 +86,9 @@ class Nuclide(object): 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'. - yield_data : dict of float to list - Maps tabulated energy to list of (product, yield) for all - neutron-induced fission products. + yield_data : dict of float to :class:`openmc.deplete.FissionYieldDistribution` + Fission product yields at tabulated energies for this nuclide. Can be + treated as a nested dictionary ``{energy: {product: yield}}`` yield_energies : list of float Energies at which fission product yields exist @@ -300,15 +300,15 @@ class FissionYieldDistribution(Mapping): ------- FissionYieldDistribution """ - yields = {} + all_yields = {} for elem_index, yield_elem in enumerate(element.iter("fission_yields")): energy = float(yield_elem.get("energy")) products = yield_elem.find("products").text.split() - yield_mapobj = map(float, yield_elem.find("data").text.split()) + yields = map(float, yield_elem.find("data").text.split()) # Get a map of products to their corresponding yield - yields[energy] = dict(zip(products, yield_mapobj)) + all_yields[energy] = dict(zip(products, yields)) - return cls.from_dict(yields) + return cls.from_dict(all_yields) @classmethod def from_dict(cls, fission_yields): @@ -327,9 +327,7 @@ class FissionYieldDistribution(Mapping): energies = sorted(fission_yields) # Get a consistent set of products to produce a matrix of yields - shared_prod = set() - for prod_set in map(set, fission_yields.values()): - shared_prod |= prod_set + shared_prod = set.union(*(set(x) for x in fission_yields.values())) ordered_prod = sorted(shared_prod) yield_matrix = empty((len(energies), len(shared_prod))) @@ -422,6 +420,9 @@ class FissionYield(Mapping): def __iter__(self): return iter(self.products) + def items(self): + return zip(self.products, self.yields) + def __add__(self, other): new = self.copy() new += other diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 2e5ad09832..dcb5127cba 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -7,7 +7,7 @@ from itertools import product import numpy as np from openmc.data import zam, ATOMIC_SYMBOL -from openmc.deplete import comm, Chain, reaction_rates, nuclide +from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram import pytest from tests import cdtemp @@ -359,3 +359,10 @@ def test_fission_yield_attribute(simple_chain): empty_chain.fission_yields = [thermal_yields] * 2 assert empty_chain.fission_yields[0] == thermal_yields assert empty_chain.fission_yields[1] == thermal_yields + + # test failure with deplete function + # number fission yields != number of materials + 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) From 6556abc2a869edd970c4adec247f95521df6828f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 10:53:49 -0500 Subject: [PATCH 22/45] Add empty fission yields to tests.TestChain This object is used in unit testing the depletion schemes, where the depletion matrix is pre-defined and there is no need for fission yields. However, the depletion interface expects the chain to have a fission_yields attribute. --- tests/dummy_operator.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index acfdc5b583..cf76f86b59 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -7,6 +7,13 @@ from openmc.deplete.abc import TransportOperator, OperatorResult class TestChain(object): + """Empty chain to assist with unit testing depletion routines + + Only really provides the form_matrix function, but acts like + a real Chain + """ + + fission_yields = [None] @staticmethod def get_thermal_fission_yields(): From 0931b58269f1e0e324c81048503680764411f7d4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 15:15:47 -0500 Subject: [PATCH 23/45] Fix wonky layout in test_deplete_fission_yields --- .../unit_tests/test_deplete_fission_yields.py | 89 +++++-------------- 1 file changed, 20 insertions(+), 69 deletions(-) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index f3911d80d0..bf0bc6b8fc 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -5,9 +5,7 @@ from unittest.mock import Mock import pytest import numpy -from numpy.testing import assert_array_equal -from openmc.capi import Material from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( FissionYieldCutoffHelper, ConstantFissionYieldHelper) @@ -27,7 +25,7 @@ def nuclide_bundle(): u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} u238 = Nuclide() u238.name = "U238" - u238.yield_energies = (5.00e5, ) + u238.yield_energies = (5.00e5,) u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) xe135 = Nuclide() @@ -35,11 +33,13 @@ def nuclide_bundle(): NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") return NuclideBundle(u235, u238, xe135) -@pytest.mark.parametrize("input_energy, u5_yield_energy", ( - (0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) + + +@pytest.mark.parametrize( + "input_energy, u5_yield_energy", + ((0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): - helper = ConstantFissionYieldHelper( - nuclide_bundle, energy=input_energy) + helper = ConstantFissionYieldHelper(nuclide_bundle, energy=input_energy) assert helper.energy == input_energy assert helper.constant_yields == { "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], @@ -54,37 +54,29 @@ def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): def test_cutoff_construction(nuclide_bundle): # defaults - helper = FissionYieldCutoffHelper( - nuclide_bundle, 1) + helper = FissionYieldCutoffHelper(nuclide_bundle, 1) assert helper.constant_yields == { "U238": nuclide_bundle.u238.yield_data[5.0e5]} assert helper.thermal_yields == { "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": nuclide_bundle.u235.yield_data[5e5]} + assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} # use 14 MeV yields - helper = FissionYieldCutoffHelper( - nuclide_bundle, 1, fast_energy=14e6) + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6) assert helper.constant_yields == { "U238": nuclide_bundle.u238.yield_data[5.0e5]} assert helper.thermal_yields == { "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": nuclide_bundle.u235.yield_data[14e6]} + assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[14e6]} # specify missing thermal yields -> use 0.0253 - helper = FissionYieldCutoffHelper( - nuclide_bundle, 1, thermal_energy=1) + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1) assert helper.thermal_yields == { "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": nuclide_bundle.u235.yield_data[5e5]} + assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} # request missing fast yields -> use epithermal - helper = FissionYieldCutoffHelper( - nuclide_bundle, 1, fast_energy=1e4) + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4) assert helper.thermal_yields == { "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": nuclide_bundle.u235.yield_data[5e5]} + assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} # test failures in cutoff: super low, super high with pytest.raises(ValueError, match="thermal yields"): FissionYieldCutoffHelper( @@ -120,66 +112,25 @@ def test_cutoff_helper(nuclide_bundle, therm_frac): proxy.generate_tallies(materials, local_mats) non_zero_nucs = [n.name for n in nuclide_bundle] tally_nucs = proxy.update_nuclides_from_operator(non_zero_nucs) - assert tally_nucs == ("U235", ) + assert tally_nucs == ("U235",) # Emulate building tallies # material x energy, tallied_nuclides, 3 proxy_flux = 1e6 tally_data = numpy.empty((n_bmats * 2, 1, 3)) tally_data[0, 0, 1] = therm_frac * proxy_flux - tally_data[1, 0, 1] = (1 - therm_frac) * proxy_flux + tally_data[1, 0, 1] = (1 - therm_frac) * proxy_flux proxy._fission_rate_tally.results = tally_data proxy.unpack() # expected results of shape (n_mats, 2, n_tnucs) expected_results = numpy.empty((1, 2, 1)) expected_results[:, 0] = therm_frac - expected_results[:, 1] = (1 - therm_frac) + expected_results[:, 1] = 1 - therm_frac assert proxy.results == pytest.approx(expected_results) actual_yields = proxy.weighted_yields(0) assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] assert actual_yields["U235"] == ( proxy.thermal_yields["U235"] * therm_frac - + proxy.fast_yields["U235"] * (1 - therm_frac)) -"""Test the FissionYieldHelpers""" - -from collections import namedtuple - -import pytest - -from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution -from openmc.deplete.helpers import ConstantFissionYieldHelper - - -@pytest.fixture(scope="module") -def nuclide_bundle(): - u5yield_dict = { - 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, - 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, - 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} - u235 = Nuclide() - u235.name = "U235" - u235.yield_energies = (0.0253, 5.0e5, 1.40e7) - u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) - - u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} - u238 = Nuclide() - u238.name = "U238" - u238.yield_energies = (5.00e5, ) - u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) - - xe135 = Nuclide() - xe135.name = "Xe135" - - NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") - return NuclideBundle(u235, u238, xe135) -@pytest.mark.parametrize("input_energy, u5_yield_energy", ( - (0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) -def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): - helper = ConstantFissionYieldHelper( - nuclide_bundle, energy=input_energy) - assert helper.energy == input_energy - assert helper.constant_yields == { - "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], - "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal - assert helper.constant_yields == helper.weighted_yields(1) + + proxy.fast_yields["U235"] * (1 - therm_frac) + ) From f2fa86fd5264432c89907c40109bbbb429746baf Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 16:09:41 -0500 Subject: [PATCH 24/45] Delegate Nuclide yield properties to FissionYieldDistribution For consistency across the various fission yield helpers, the Nuclide object now sets the yield data to always be a FissionYieldDistribution or an empty dictionary. Nuclide.yield_energies is now a property that returns the energies from the FissionYieldDistribution. When setting the yield_data, the Nuclide enforces the object to be a Mapping instance. If the incoming yield_data is a FissionYieldDistribution, no conversions are made. Otherwise, a new distribution is created with the from_dict method. The Chain.from_xml method now creates the yield_data for nuclides using FissionYieldDistribution.from_dict Changes were made to some unit tests that attempted to set the yield_energies, which do not have a setter. --- openmc/deplete/chain.py | 14 +++++----- openmc/deplete/nuclide.py | 27 +++++++++++++++---- tests/unit_tests/test_deplete_chain.py | 3 +-- .../unit_tests/test_deplete_fission_yields.py | 2 -- tests/unit_tests/test_deplete_nuclide.py | 8 ++++-- 5 files changed, 37 insertions(+), 17 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index a204e4821c..dd13a338be 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -15,6 +15,7 @@ from warnings import warn from openmc.checkvalue import check_type, check_less_than from openmc.data import gnd_name, zam +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, @@ -276,11 +277,12 @@ class Chain(object): fpy = fpy_data[parent] if fpy.energies is not None: - nuclide.yield_energies = fpy.energies + yield_energies = fpy.energies else: - nuclide.yield_energies = [0.0] + yield_energies = [0.0] - for E, table in zip(nuclide.yield_energies, fpy.independent): + yield_data = {} + for E, table in zip(yield_energies, fpy.independent): yield_replace = 0.0 yields = defaultdict(float) for product, y in table.items(): @@ -294,10 +296,10 @@ class Chain(object): if yield_replace > 0.0: missing_fp.append((parent, E, yield_replace)) + yield_data[E] = yields - nuclide.yield_data[E] = [] - for k in sorted(yields, key=openmc.data.zam): - nuclide.yield_data[E].append((k, yields[k])) + nuclide.yield_data = ( + FissionYieldDistribution.from_dict(yield_data)) # Display warnings if missing_daughter: diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 5217831f49..649b6700f8 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -86,10 +86,10 @@ class Nuclide(object): 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'. - yield_data : dict of float to :class:`openmc.deplete.FissionYieldDistribution` + yield_data : FissionYieldDistribution or None Fission product yields at tabulated energies for this nuclide. Can be treated as a nested dictionary ``{energy: {product: yield}}`` - yield_energies : list of float + yield_energies : tuple of float or None Energies at which fission product yields exist """ @@ -107,8 +107,7 @@ class Nuclide(object): self.reactions = [] # Neutron fission yields, if present - self.yield_data = {} - self.yield_energies = [] + self._yield_data = None @property def n_decay_modes(self): @@ -118,6 +117,25 @@ class Nuclide(object): def n_reaction_paths(self): return len(self.reactions) + @property + def yield_data(self): + if self._yield_data is None: + return {} + return self._yield_data + + @yield_data.setter + def yield_data(self, fission_yields): + check_type("fission_yields", fission_yields, Mapping) + if isinstance(fission_yields, FissionYieldDistribution): + self._yield_data = fission_yields + self._yield_data = FissionYieldDistribution.from_dict(fission_yields) + + @property + def yield_energies(self): + if self._yield_data is None: + return None + return self.yield_data.energies + @classmethod def from_xml(cls, element, fission_q=None): """Read nuclide from an XML element. @@ -173,7 +191,6 @@ class Nuclide(object): fpy_elem = element.find('neutron_fission_yields') if fpy_elem is not None: nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem) - nuc.yield_energies = list(sorted(nuc.yield_data.keys())) return nuc diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index dcb5127cba..98ca023f00 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -129,7 +129,7 @@ def test_from_xml(simple_chain): assert [r.branching_ratio for r in nuc.reactions] == [1.0, 0.7, 0.3] # Yield tests - assert nuc.yield_energies == [0.0253] + assert nuc.yield_energies == (0.0253, ) assert list(nuc.yield_data) == [0.0253] assert nuc.yield_data[0.0253].products == ("A", "B") assert (nuc.yield_data[0.0253].yields == [0.0292737, 0.002566345]).all() @@ -163,7 +163,6 @@ def test_export_to_xml(run_in_tmpdir): nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7), nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3) ] - C.yield_energies = [0.0253] C.yield_data = nuclide.FissionYieldDistribution.from_dict({ 0.0253: {"A": 0.0292737, "B": 0.002566345}}) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index bf0bc6b8fc..2330742361 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -19,13 +19,11 @@ def nuclide_bundle(): 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} u235 = Nuclide() u235.name = "U235" - u235.yield_energies = (0.0253, 5.0e5, 1.40e7) u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} u238 = Nuclide() u238.name = "U238" - u238.yield_energies = (5.00e5,) u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) xe135 = Nuclide() diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 2a3993921f..ef5c3d1ca7 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -3,6 +3,7 @@ import xml.etree.ElementTree as ET import numpy +import pytest from openmc.deplete import nuclide @@ -73,8 +74,12 @@ def test_from_xml(): ] expected_yield_data = nuclide.FissionYieldDistribution.from_dict({ 0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}}) - assert u235.yield_energies == [0.0253] assert u235.yield_data == expected_yield_data + # test accessing the yield energies through the FissionYieldDistribution + assert u235.yield_energies == (0.0253, ) + assert u235.yield_energies is u235.yield_data.energies + with pytest.raises(AttributeError): # not settable + u235.yield_energies = [0.0253, 5e5] def test_to_xml_element(): @@ -91,7 +96,6 @@ def test_to_xml_element(): nuclide.ReactionTuple('fission', None, 2.0e8, 1.0), nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0) ] - C.yield_energies = [0.0253] C.yield_data = nuclide.FissionYieldDistribution.from_dict( {0.0253: {"A": 0.0292737, "B": 0.002566345}}) element = C.to_xml_element() From 454e6cc6b6739cd6f358951757354f654171ce80 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 17:39:57 -0500 Subject: [PATCH 25/45] Add AveragedFissionYieldHelper operator helper Computes the effective fission yields for nuclides with multiple sets of yields by 1) Tallying fission rate and energy-weighted fission rate 2) Determining average energy at which fission events occur 3) Interpolating [lin-lin] between adjacent fission yields based on this average energy Uses a second tally with an EnergyFunctionFilter to tally E * sigma_f * phi. Added unit tests with the proxy-style class that doesn't use the C API to generate tallies but tests all other aspects of the implementation. Tests cover three possible average energies (below minimum supplied yield energy, above max supplied yield energy, and between two) The Operator uses fission_yield_mode="average" to select this helper. The helper also has a from_operator class method. --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/abc.py | 10 +- openmc/deplete/helpers.py | 157 +++++++++++++++++- openmc/deplete/operator.py | 4 +- .../unit_tests/test_deplete_fission_yields.py | 69 ++++++-- 5 files changed, 217 insertions(+), 24 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 44f8b70e05..5699f3f2a8 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -77,6 +77,7 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber + AveragedFissionYieldHelper ChainFissionHelper ConstantFissionYieldHelper DirectReactionRateHelper diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 5efda87029..d1dc2cdc4d 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -504,6 +504,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` + results : None or numpy.ndarray + Tally results shaped in a manner useful to this helper. """ _upper_energy = 20.0e6 # upper energy for tallies @@ -513,7 +515,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): self.n_bmats = n_bmats self._local_indexes = None self._fission_rate_tally = None - self._tally_index = {} + self._tally_nucs = [] + self.results = None def generate_tallies(self, materials, mat_indexes): """Construct the fission rate tally @@ -554,10 +557,10 @@ class TalliedFissionYieldHelper(FissionYieldHelper): if len(overlap) == 0: # tally no nuclides, but keep the Tally alive self._fission_rate_tally.nuclides = None - self._tally_index = [] + self._tally_nucs = [] return tuple() nuclides = tuple(sorted(overlap)) - self._tally_index = [self._chain_nuclides[n] for n in nuclides] + self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] self._fission_rate_tally.nuclides = nuclides return nuclides @@ -587,6 +590,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): return cls(operator.chain.nuclides, len(operator.burnable_mats), **kwargs) + class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index b66c9ba5e7..19e51628d9 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -9,14 +9,15 @@ from numpy import dot, zeros, newaxis from openmc.checkvalue import check_type, check_greater_than from openmc.capi import ( - Tally, MaterialFilter, EnergyFilter) + Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) from .abc import ( ReactionRateHelper, EnergyHelper, FissionYieldHelper, TalliedFissionYieldHelper) __all__ = ( "DirectReactionRateHelper", "ChainFissionHelper", - "ConstantFissionYieldHelper", "FissionYieldCutoffHelper") + "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", + "AveragedFissionYieldHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -386,7 +387,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): def unpack(self): """Obtain fast and thermal fission fractions from tally""" fission_rates = self._fission_rate_tally.results[..., 1].reshape( - self.n_bmats, 2, len(self._tally_index)) + self.n_bmats, 2, len(self._tally_nucs)) self.results = fission_rates[self._local_indexes] total_fission = self.results.sum(axis=1) nz_mat, nz_nuc = total_fission.nonzero() @@ -422,10 +423,10 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): rates = self.results[local_mat_index] yields = self.constant_yields # iterate over thermal then fast yields, prefer __mul__ to __rmul__ - for therm_frac, nuc in zip(rates[0], self._tally_index): + for therm_frac, nuc in zip(rates[0], self._tally_nucs): yields[nuc.name] = self._thermal_yields[nuc.name] * therm_frac - for fast_frac, nuc in zip(rates[1], self._tally_index): + for fast_frac, nuc in zip(rates[1], self._tally_nucs): yields[nuc.name] += self._fast_yields[nuc.name] * fast_frac return yields @@ -442,3 +443,149 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): for key, sub in self._fast_yields.items(): out[key] = sub.copy() return out + + +class AveragedFissionYieldHelper(TalliedFissionYieldHelper): + r"""Class that computes fission yields based on average fission energy + + Computes average energy at which fission events occured + reactions for all nuclides with multiple sets of fission yields + by + + .. math:: + + \bar{E} = \frac{ + \int_0^\infty E\sigma_f(E)\phi(E)dE + }{ + \int_0^\infty\sigma_f(E)\phi(E)dE + } + + If the average energy for a nuclide is below the lowest energy + with yield data, that set of fission yields is taken. + Conversely, if the average energy is above the highest energy + with yield data, that set of fission yields is used. + For the case where the average energy is between two sets + of yields, a geometric mean of the yield distributions is + used. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + n_bmats : int + Number of burnable materials tracked in the problem. + + Attributes + ---------- + n_bmats : int + Number of burnable materials tracked in the problem. + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + results : None or numpy.ndarray + If tallies have been generated and unpacked, then the array will + have shape ``(n_mats, n_tnucs)``, where ``n_mats`` is the number + of materials where fission reactions were tallied and ``n_tnucs`` + is the number of nuclides with multiple sets of fission yields. + """ + + def __init__(self, chain_nuclides, n_bmats): + super().__init__(chain_nuclides, n_bmats) + self._weighted_tally = None + + def generate_tallies(self, materials, mat_indexes): + """Construct tallies to determine average energy of fissions + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indexes for materials in ``materials`` tracked on this + process + """ + super().generate_tallies(materials, mat_indexes) + fission_tally = self._fission_rate_tally + + weighted_tally = Tally() + weighted_tally.filters = fission_tally.filters.copy() + weighted_tally.nuclides = fission_tally.nuclides + weighted_tally.scores = ['fission'] + + ene_bin = EnergyFilter() + ene_bin.bins = (0, self._upper_energy) + fission_tally.filters.append(ene_bin) + + ene_filter = EnergyFunctionFilter() + ene_filter.set_data((0, self._upper_energy), (0, self._upper_energy)) + weighted_tally.filters.append(ene_filter) + self._weighted_tally = weighted_tally + + def unpack(self): + """Unpack tallies and populate :attr:`results` with average energies""" + fission_results = ( + self._fission_rate_tally.results[self._local_indexes, :, 1]) + self.results = ( + self._weighted_tally.results[self._local_indexes, :, 1]).copy() + nz_mat, nz_nuc = fission_results.nonzero() + self.results[nz_mat, nz_nuc] /= fission_results[nz_mat, nz_nuc] + + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + Use the computed average energy of fission + events to determine fission yields. If average + energy is between two sets of yields, linearly + interpolate bewteen the two. + Otherwise take the closet set of yields. + + Parameters + ---------- + local_mat_index : int + Index for specific burnable material. Effective + yields will be produced using + ``self.results[local_mat_index]`` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + mat_yields = {} + average_energies = self.results[local_mat_index] + for avg_e, nuc in zip(average_energies, self._tally_nucs): + nuc_energies = nuc.yield_energies + if avg_e <= nuc_energies[0]: + mat_yields[nuc.name] = nuc.yield_data[nuc_energies[0]] + continue + if avg_e >= nuc_energies[-1]: + mat_yields[nuc.name] = nuc.yield_data[nuc_energies[-1]] + continue + # in-between two energies + # linear search since there are usually ~3 energies + for ix, ene in enumerate(nuc_energies[:-1]): + if nuc_energies[ix + 1] > avg_e: + break + lower, upper = nuc_energies[ix:ix + 2] + fast_frac = (avg_e - lower) / (upper - lower) + mat_yields[nuc.name] = ( + nuc.yield_data[lower] * (1 - fast_frac) + + nuc.yield_data[upper] * fast_frac) + mat_yields.update(self.constant_yields) + return mat_yields + + @classmethod + def from_operator(cls, operator): + """Return a new helper with data from an operator + + Parameters + ---------- + operator : openmc.deplete.Operator + Operator with a depletion chain + + Returns + ------- + AveragedFissionYieldHelper + """ + return cls(operator.chain.nuclides, len(operator.burnable_mats)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index e967dc85d1..b7df9b8fc7 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -27,7 +27,7 @@ from .reaction_rates import ReactionRates from .results_list import ResultsList from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, - FissionYieldCutoffHelper) + FissionYieldCutoffHelper, AveragedFissionYieldHelper) def _distribute(items): @@ -92,6 +92,7 @@ class Operator(TransportOperator): * "constant": :class:`openmc.deplete.ConstantFissionYieldHelper` * "cutoff": :class:`openmc.deplete.FissionYieldCutoffHelper` + * "average": :class:`openmc.deplete.AveragedFissionYieldHelper` The documentation on these classes describe their methodology and differences. ``"constant"`` will treat fission yields as @@ -140,6 +141,7 @@ class Operator(TransportOperator): Whether to differentiate burnable materials with multiple instances """ _fission_helpers_ = { + "average": AveragedFissionYieldHelper, "constant": ConstantFissionYieldHelper, "cutoff": FissionYieldCutoffHelper, } diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 2330742361..0a87c4fb6c 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -8,7 +8,11 @@ import numpy from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( - FissionYieldCutoffHelper, ConstantFissionYieldHelper) + FissionYieldCutoffHelper, ConstantFissionYieldHelper, + AveragedFissionYieldHelper) + + +MATERIALS = ["1", "2"] @pytest.fixture(scope="module") @@ -45,11 +49,6 @@ def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): assert helper.constant_yields == helper.weighted_yields(1) -# --------------------------------- -# Test the FissionYieldCutoffHelper -# --------------------------------- - - def test_cutoff_construction(nuclide_bundle): # defaults helper = FissionYieldCutoffHelper(nuclide_bundle, 1) @@ -92,22 +91,23 @@ def test_cutoff_failure(key): FissionYieldCutoffHelper(None, None, **{key: -1}) -class CutoffProxy(FissionYieldCutoffHelper): - """Proxy that supplies a set of tallies""" - +class ProxyMixin(object): + """Mixing that overloads the tally generation""" def generate_tallies(self, materials, mat_indexes): self._fission_rate_tally = Mock() self._local_indexes = numpy.asarray(mat_indexes) +class CutoffProxy(ProxyMixin, FissionYieldCutoffHelper): + """Proxy that supplies a set of tallies""" + + # emulate some split between fast and thermal U235 fissions @pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) def test_cutoff_helper(nuclide_bundle, therm_frac): - materials = ["1", "2"] # TODO Use real C API materials? - n_bmats = len(materials) - local_mats = [0] + n_bmats = len(MATERIALS) proxy = CutoffProxy(nuclide_bundle, n_bmats) - proxy.generate_tallies(materials, local_mats) + proxy.generate_tallies(MATERIALS, [0]) non_zero_nucs = [n.name for n in nuclide_bundle] tally_nucs = proxy.update_nuclides_from_operator(non_zero_nucs) assert tally_nucs == ("U235",) @@ -130,5 +130,44 @@ def test_cutoff_helper(nuclide_bundle, therm_frac): assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] assert actual_yields["U235"] == ( proxy.thermal_yields["U235"] * therm_frac - + proxy.fast_yields["U235"] * (1 - therm_frac) - ) + + proxy.fast_yields["U235"] * (1 - therm_frac)) + + +class AverageProxy(ProxyMixin, AveragedFissionYieldHelper): + """Proxy for generating mock set of tallies""" + def generate_tallies(self, materials, mat_indexes): + super().generate_tallies(materials, mat_indexes) + self._weighted_tally = Mock() + + +@pytest.mark.parametrize("avg_energy", (0.01, 100, 15e6)) +def test_averaged_helper(nuclide_bundle, avg_energy): + proxy = AverageProxy(nuclide_bundle, len(MATERIALS)) + proxy.generate_tallies(MATERIALS, [0]) + tallied_nucs = proxy.update_nuclides_from_operator( + [n.name for n in nuclide_bundle]) + assert tallied_nucs == ("U235", ) + # enforce some average energy + proxy_flux = 1e16 + fission_results = numpy.ones((len(MATERIALS), 1, 3)) * proxy_flux + weighted_results = fission_results * avg_energy + proxy._fission_rate_tally.results = fission_results + proxy._weighted_tally.results = weighted_results + proxy.unpack() + expected_results = numpy.ones((1, 1)) * avg_energy + assert proxy.results == pytest.approx(expected_results) + + actual_yields = proxy.weighted_yields(0) + # constant U238 => no interpolation + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + # construct expected yields + if avg_energy < 0.0253: # take thermal U235 yields + exp_u235_yields = nuclide_bundle.u235.yield_data[0.0253] + elif avg_energy > 14e6: # take fastest U235 yields + exp_u235_yields = nuclide_bundle.u235.yield_data[14e6] + else: # reconstruct between thermal and epithermal + thermal = nuclide_bundle.u235.yield_data[0.0253] + epithermal = nuclide_bundle.u235.yield_data[5e5] + split = (avg_energy - 0.0253) / (5e5 - 0.0253) + exp_u235_yields = thermal * (1 - split) + epithermal * split + assert actual_yields["U235"] == exp_u235_yields From aa86aa5262690e9a3e5b3ac08bac799ad642209f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 18:02:27 -0500 Subject: [PATCH 26/45] Improve documentation for FissionYieldHelpers Add a separate section in the documentation describing the helpers. Reword purpose of mat_indexes in generate_tallies --- docs/source/pythonapi/deplete.rst | 19 ++++++++++++++----- openmc/deplete/abc.py | 14 ++++++++++---- openmc/deplete/helpers.py | 19 +++++++++++++------ openmc/deplete/nuclide.py | 2 +- openmc/deplete/operator.py | 15 ++++++--------- 5 files changed, 44 insertions(+), 25 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 5699f3f2a8..5b26eb38aa 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -77,16 +77,25 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber - AveragedFissionYieldHelper - ChainFissionHelper - ConstantFissionYieldHelper - DirectReactionRateHelper - FissionYieldCutoffHelper OperatorResult ReactionRates Results ResultsList +The following classes are used to help the :class:`openmc.deplete.Operator` +compute quantities like effective fission yields, reaction rates, and +total system energy. + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + helpers.AveragedFissionYieldHelper + helpers.ChainFissionHelper + helpers.ConstantFissionYieldHelper + helpers.DirectReactionRateHelper + helpers.FissionYieldCutoffHelper The following classes are abstract classes that can be used to extend the :mod:`openmc.deplete` capabilities: diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index d1dc2cdc4d..fd5278d894 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -443,8 +443,11 @@ class FissionYieldHelper(ABC): materials : iterable of C-API materials Materials to be used in :class:`openmc.capi.MaterialFilter` mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. """ def update_nuclides_from_operator(self, nuclides): @@ -526,8 +529,11 @@ class TalliedFissionYieldHelper(FissionYieldHelper): materials : iterable of :class:`openmc.capi.Material` Materials to be used in :class:`openmc.capi.MaterialFilter` mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. """ self._local_indexes = asarray(mat_indexes) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 19e51628d9..8dc9035b3d 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -376,8 +376,11 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): materials : iterable of :class:`openmc.capi.Material` Materials to be used in :class:`openmc.capi.MaterialFilter` mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. """ super().generate_tallies(materials, mat_indexes) energy_filter = EnergyFilter() @@ -465,8 +468,9 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): Conversely, if the average energy is above the highest energy with yield data, that set of fission yields is used. For the case where the average energy is between two sets - of yields, a geometric mean of the yield distributions is - used. + of yields, the effective fission yield computed by + linearly interpolating between yields provided at the + nearest energies above and below the average. Parameters ---------- @@ -502,8 +506,11 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): materials : iterable of :class:`openmc.capi.Material` Materials to be used in :class:`openmc.capi.MaterialFilter` mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. """ super().generate_tallies(materials, mat_indexes) fission_tally = self._fission_rate_tally diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 649b6700f8..d2b0acc169 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -236,7 +236,7 @@ class Nuclide(object): class FissionYieldDistribution(Mapping): - """Class for storing energy-dependent fission yields for a single nuclide + """Energy-dependent fission product yields for a single nuclide Can be used as a dictionary mapping energies and products to fission yields:: diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index b7df9b8fc7..5b9d2fca67 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -86,19 +86,16 @@ class Operator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. - fission_yield_mode : str, {"constant", "cutoff"} + fission_yield_mode : ("constant", "cutoff", "average") Key indicating what fission product yield scheme to use. The - key determines what fission energy helper is used:: + key determines what fission energy helper is used: - * "constant": :class:`openmc.deplete.ConstantFissionYieldHelper` - * "cutoff": :class:`openmc.deplete.FissionYieldCutoffHelper` - * "average": :class:`openmc.deplete.AveragedFissionYieldHelper` + * "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper` + * "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper` + * "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper` The documentation on these classes describe their methodology - and differences. ``"constant"`` will treat fission yields as - constant with respect to energy, while ``"cutoff"`` will - compute effective yields based on the number of fission events - above and below a cutoff. Default: ``"constant"`` + and differences. Default: ``"constant"`` fission_yield_opts : dict of str to option, optional Optional arguments to pass to the helper determined by ``fission_yield_mode``. Will be passed directly on to the From d8b5752b59da10eb4d0a902da9e902e0c83ce5d4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 15 Aug 2019 11:20:04 -0500 Subject: [PATCH 27/45] Make n_bmat required only for FissionYieldCutoffHelper Base class TalliedFissionYieldHelper and AveragedFissionYieldHelper instances do not need to know the number of burnable materials, and thus it is no longer a required input argument. Changes propagated into the from_operator methods and tests/unit_tests/test_deplete_fission_yields.py Removed explicit TalliedFissionYieldHelper.from_operator method as it is identical and now inherited from the base FissionYieldHelper.from_operator method --- openmc/deplete/abc.py | 25 +------------------ openmc/deplete/helpers.py | 15 +++++------ .../unit_tests/test_deplete_fission_yields.py | 2 +- 3 files changed, 8 insertions(+), 34 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index fd5278d894..56914367e2 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -497,13 +497,9 @@ class TalliedFissionYieldHelper(FissionYieldHelper): chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int - Number of burnable materials tracked in the problem. Attributes ---------- - n_bmats : int - Number of burnable materials tracked in the problem. constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` @@ -513,9 +509,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): _upper_energy = 20.0e6 # upper energy for tallies - def __init__(self, chain_nuclides, n_bmats): + def __init__(self, chain_nuclides): super().__init__(chain_nuclides) - self.n_bmats = n_bmats self._local_indexes = None self._fission_rate_tally = None self._tally_nucs = [] @@ -578,24 +573,6 @@ class TalliedFissionYieldHelper(FissionYieldHelper): tally data. """ - @classmethod - def from_operator(cls, operator, **kwargs): - """Create a new instance by pulling data from the operator - - Require the more concrete :class:`openmc.deplete.Operator` - because this needs knowledge of burnable materials - - Parameters - ---------- - operator : openmc.deplete.TransportOperator - Operator with a depletion chain - kwargs: optional - Optional keyword arguments to be passed to the underlying - ``__init__`` method - """ - return cls(operator.chain.nuclides, len(operator.burnable_mats), - **kwargs) - class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 8dc9035b3d..3303dd4e40 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -299,7 +299,8 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): check_greater_than("thermal_energy", thermal_energy, 0.0, equality=True) check_greater_than("cutoff", cutoff, thermal_energy, equality=False) check_greater_than("fast_energy", fast_energy, cutoff, equality=False) - super().__init__(chain_nuclides, n_bmats) + self.n_bmats = n_bmats + super().__init__(chain_nuclides) self._cutoff = cutoff self._thermal_yields = {} self._fast_yields = {} @@ -477,13 +478,9 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): chain_nuclides : iterable of openmc.deplete.Nuclide Nuclides tracked in the depletion chain. Not necessary that all have yield data. - n_bmats : int - Number of burnable materials tracked in the problem. Attributes ---------- - n_bmats : int - Number of burnable materials tracked in the problem. constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` @@ -494,8 +491,8 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): is the number of nuclides with multiple sets of fission yields. """ - def __init__(self, chain_nuclides, n_bmats): - super().__init__(chain_nuclides, n_bmats) + def __init__(self, chain_nuclides): + super().__init__(chain_nuclides) self._weighted_tally = None def generate_tallies(self, materials, mat_indexes): @@ -588,11 +585,11 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): Parameters ---------- - operator : openmc.deplete.Operator + operator : openmc.deplete.TransportOperator Operator with a depletion chain Returns ------- AveragedFissionYieldHelper """ - return cls(operator.chain.nuclides, len(operator.burnable_mats)) + return cls(operator.chain.nuclides) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 0a87c4fb6c..6285fe5b9b 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -142,7 +142,7 @@ class AverageProxy(ProxyMixin, AveragedFissionYieldHelper): @pytest.mark.parametrize("avg_energy", (0.01, 100, 15e6)) def test_averaged_helper(nuclide_bundle, avg_energy): - proxy = AverageProxy(nuclide_bundle, len(MATERIALS)) + proxy = AverageProxy(nuclide_bundle) proxy.generate_tallies(MATERIALS, [0]) tallied_nucs = proxy.update_nuclides_from_operator( [n.name for n in nuclide_bundle]) From 73a3ad3e5ee26386259dbeb66f17f2a9f6445e61 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 16 Aug 2019 16:13:33 -0500 Subject: [PATCH 28/45] Guard against using FPY tally results if not needed If there are no nuclides with multiple sets of yields, then the tallies needed for fission yields may or or may not exist. Regardless, this means that using constant_yields is absolutely appropriate and there is no need for unpacking on FissionYieldCutoffHelper nor AveragedFissionYieldHelper. The weighted_yields methods for these two classes simply return the set of constant yields if this is the case. --- openmc/deplete/abc.py | 5 ++++- openmc/deplete/helpers.py | 12 +++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 56914367e2..4406614eb4 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -556,8 +556,11 @@ class TalliedFissionYieldHelper(FissionYieldHelper): """ overlap = set(self._chain_nuclides).intersection(set(nuclides)) if len(overlap) == 0: + if self._fission_rate_tally is None: + # No tally to update + return tuple() # tally no nuclides, but keep the Tally alive - self._fission_rate_tally.nuclides = None + self._fission_rate_tally.nuclides = [] self._tally_nucs = [] return tuple() nuclides = tuple(sorted(overlap)) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 3303dd4e40..1336af7e69 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -390,6 +390,9 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): def unpack(self): """Obtain fast and thermal fission fractions from tally""" + if not self._tally_nucs: + self.results = None + return fission_rates = self._fission_rate_tally.results[..., 1].reshape( self.n_bmats, 2, len(self._tally_nucs)) self.results = fission_rates[self._local_indexes] @@ -424,8 +427,10 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): library : dict Dictionary of ``{parent: {product: fyield}}`` """ - rates = self.results[local_mat_index] yields = self.constant_yields + if not self._tally_nucs: + return yields + rates = self.results[local_mat_index] # iterate over thermal then fast yields, prefer __mul__ to __rmul__ for therm_frac, nuc in zip(rates[0], self._tally_nucs): yields[nuc.name] = self._thermal_yields[nuc.name] * therm_frac @@ -528,6 +533,9 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): def unpack(self): """Unpack tallies and populate :attr:`results` with average energies""" + if not self._tally_nucs: + self.results = None + return fission_results = ( self._fission_rate_tally.results[self._local_indexes, :, 1]) self.results = ( @@ -556,6 +564,8 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): library : dict Dictionary of ``{parent: {product: fyield}}`` """ + if not self._tally_nucs: + return self.constant_yields mat_yields = {} average_energies = self.results[local_mat_index] for avg_e, nuc in zip(average_energies, self._tally_nucs): From 0cf3ea3b16e65e928af1b66e9b0b3e38567a8106 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 19 Aug 2019 09:43:40 -0500 Subject: [PATCH 29/45] Improve updating of nuclides on tallied FPY helpers Renamed update_nuclides_from_operator -> update_tally_nuclides to be more general. The base method TalliedFissionYieldHelper.update_tally_nuclides will raise an AttributeError if this method is called before the tallies are generated. The inner logic for setting the tallied nuclides can be streamlined given the assumption that the tally exists. AveragedFissionYieldHelper.update_tally_nuclides has been improved to set the nuclides for the weighted tally now. This was not captured prior because the nuclides were copied over from the fission rate tally during tally generation, prior to when nuclides were set. --- openmc/deplete/abc.py | 33 ++++++++++--------- openmc/deplete/helpers.py | 27 ++++++++++++++- openmc/deplete/operator.py | 2 +- .../unit_tests/test_deplete_fission_yields.py | 4 +-- 4 files changed, 46 insertions(+), 20 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 4406614eb4..da9bb5244e 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -450,7 +450,7 @@ class FissionYieldHelper(ABC): in parallel mode. """ - def update_nuclides_from_operator(self, nuclides): + def update_tally_nuclides(self, nuclides): """Return nuclides with non-zero densities and yield data Parameters @@ -489,7 +489,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): Generates a basic fission rate tally in all burnable materials with :meth:`generate_tallies`, and set nuclides to be tallied with - :meth:`update_nuclides_from_operator`. Subclasses will need to implement + :meth:`update_tally_nuclides`. Subclasses will need to implement :meth:`unpack` and :meth:`weighted_yields`. Parameters @@ -538,31 +538,32 @@ class TalliedFissionYieldHelper(FissionYieldHelper): self._fission_rate_tally.filters = [MaterialFilter(materials)] - def update_nuclides_from_operator(self, nuclides): + def update_tally_nuclides(self, nuclides): """Tally nuclides with non-zero density and multiple yields + Must be run after :meth:`generate_tallies`. + Parameters ---------- nuclides : iterable of str - Nuclides with non-zero densities from the - :class:`openmc.deplete.Operator` + Potential nuclides to be tallied, such as those with + non-zero density at this stage. Returns ------- nuclides : tuple of str - Union of nuclides that the :class:`openmc.deplete.Operator` - says have non-zero densities at this stage and those that - have multiple sets of yield data. Sorted by nuclide name + Union of input nuclides and those that have multiple sets + of yield data. Sorted by nuclide name + + Raises + ------ + AttributeError + If tallies not generated """ + if self._fission_rate_tally is None: + raise AttributeError( + "Tallies not built. Run generate_tallies first") overlap = set(self._chain_nuclides).intersection(set(nuclides)) - if len(overlap) == 0: - if self._fission_rate_tally is None: - # No tally to update - return tuple() - # tally no nuclides, but keep the Tally alive - self._fission_rate_tally.nuclides = [] - self._tally_nucs = [] - return tuple() nuclides = tuple(sorted(overlap)) self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] self._fission_rate_tally.nuclides = nuclides diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 1336af7e69..aeafd412e8 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -519,7 +519,6 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): weighted_tally = Tally() weighted_tally.filters = fission_tally.filters.copy() - weighted_tally.nuclides = fission_tally.nuclides weighted_tally.scores = ['fission'] ene_bin = EnergyFilter() @@ -531,6 +530,32 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): weighted_tally.filters.append(ene_filter) self._weighted_tally = weighted_tally + def update_tally_nuclides(self, nuclides): + """Tally nuclides with non-zero density and multiple yields + + Must be run after :meth:`generate_tallies`. + + Parameters + ---------- + nuclides : iterable of str + Potential nuclides to be tallied, such as those with + non-zero density at this stage. + + Returns + ------- + nuclides : tuple of str + Union of input nuclides and those that have multiple sets + of yield data. Sorted by nuclide name + + Raises + ------ + AttributeError + If tallies not generated + """ + tally_nucs = super().update_tally_nuclides(nuclides) + self._weighted_tally.nuclides = tally_nucs + return tally_nucs + def unpack(self): """Unpack tallies and populate :attr:`results` with average energies""" if not self._tally_nucs: diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 5b9d2fca67..b7a4d2c8f0 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -242,7 +242,7 @@ class Operator(TransportOperator): nuclides = self._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._energy_helper.nuclides = nuclides - self._fsn_yield_helper.update_nuclides_from_operator(nuclides) + self._fsn_yield_helper.update_tally_nuclides(nuclides) # Run OpenMC openmc.capi.reset() diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 6285fe5b9b..3b29aac96e 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -109,7 +109,7 @@ def test_cutoff_helper(nuclide_bundle, therm_frac): proxy = CutoffProxy(nuclide_bundle, n_bmats) proxy.generate_tallies(MATERIALS, [0]) non_zero_nucs = [n.name for n in nuclide_bundle] - tally_nucs = proxy.update_nuclides_from_operator(non_zero_nucs) + tally_nucs = proxy.update_tally_nuclides(non_zero_nucs) assert tally_nucs == ("U235",) # Emulate building tallies # material x energy, tallied_nuclides, 3 @@ -144,7 +144,7 @@ class AverageProxy(ProxyMixin, AveragedFissionYieldHelper): def test_averaged_helper(nuclide_bundle, avg_energy): proxy = AverageProxy(nuclide_bundle) proxy.generate_tallies(MATERIALS, [0]) - tallied_nucs = proxy.update_nuclides_from_operator( + tallied_nucs = proxy.update_tally_nuclides( [n.name for n in nuclide_bundle]) assert tallied_nucs == ("U235", ) # enforce some average energy From 2aa1849edaf01ef763a3da27355a4523a7a79ef7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 08:46:55 -0500 Subject: [PATCH 30/45] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/deplete/chain.py | 6 +++--- openmc/deplete/helpers.py | 4 ++-- tests/unit_tests/test_deplete_chain.py | 2 +- tests/unit_tests/test_deplete_fission_yields.py | 2 +- tests/unit_tests/test_deplete_nuclide.py | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index dd13a338be..27b73134b3 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -127,7 +127,7 @@ class Chain(object): 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: Real}}``, where ``yield`` is the fission product yield + 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. @@ -395,14 +395,14 @@ class Chain(object): Returns ------- fission_yields : dict - Dictionary of ``{parent : {product : f_yield}}`` + 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 = {} for nuc in self.nuclides: - if len(nuc.yield_data) == 0: + if not nuc.yield_data: continue yield_obj = nuc.yield_data[min(nuc.yield_energies)] out[nuc.name] = dict(yield_obj) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index aeafd412e8..a07630df76 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -317,7 +317,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): if fast is None: # find first index <= cutoff rev_ix = self._find_fallback_energy( - name, list(reversed(energies)), cutoff, False) + name, reversed(energies), cutoff, False) fast = yields[energies[-rev_ix]] self._thermal_yields[name] = thermal self._fast_yields[name] = fast @@ -518,7 +518,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): fission_tally = self._fission_rate_tally weighted_tally = Tally() - weighted_tally.filters = fission_tally.filters.copy() + weighted_tally.filters = fission_tally.filters weighted_tally.scores = ['fission'] ene_bin = EnergyFilter() diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 98ca023f00..7a5a128450 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -129,7 +129,7 @@ def test_from_xml(simple_chain): assert [r.branching_ratio for r in nuc.reactions] == [1.0, 0.7, 0.3] # Yield tests - assert nuc.yield_energies == (0.0253, ) + assert nuc.yield_energies == (0.0253,) assert list(nuc.yield_data) == [0.0253] assert nuc.yield_data[0.0253].products == ("A", "B") assert (nuc.yield_data[0.0253].yields == [0.0292737, 0.002566345]).all() diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 3b29aac96e..547f34c227 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -91,7 +91,7 @@ def test_cutoff_failure(key): FissionYieldCutoffHelper(None, None, **{key: -1}) -class ProxyMixin(object): +class ProxyMixin: """Mixing that overloads the tally generation""" def generate_tallies(self, materials, mat_indexes): self._fission_rate_tally = Mock() diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index ef5c3d1ca7..2cfb98e6ca 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -76,7 +76,7 @@ def test_from_xml(): 0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}}) assert u235.yield_data == expected_yield_data # test accessing the yield energies through the FissionYieldDistribution - assert u235.yield_energies == (0.0253, ) + assert u235.yield_energies == (0.0253,) assert u235.yield_energies is u235.yield_data.energies with pytest.raises(AttributeError): # not settable u235.yield_energies = [0.0253, 5e5] From 1d046ab4f62eec2c699999fb921c893240b98257 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 09:34:44 -0500 Subject: [PATCH 31/45] Support kwargs for FissionYieldHelper.from_operator Applied to concrete classes AveragedFissionYieldHelper, FissionYieldCutoffHelper, ConstantFissionYieldHelper. --- openmc/deplete/abc.py | 17 +++++++------- openmc/deplete/helpers.py | 49 +++++++++++++++++++-------------------- 2 files changed, 32 insertions(+), 34 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index da9bb5244e..7cccad3695 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -371,13 +371,11 @@ class FissionYieldHelper(ABC): Parameters ---------- chain_nuclides : iterable of openmc.deplete.Nuclide - Nuclides tracked in the depletion chain. Not necessary - that all have yield data. + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. Attributes ---------- - n_bmats : int - Number of burnable materials tracked in the problem constant_yields : dict of str to :class:`openmc.deplete.FissionYield` Fission yields for all nuclides that only have one set of fission yield data. Can be accessed as ``{parent: {product: yield}}`` @@ -473,13 +471,15 @@ class FissionYieldHelper(ABC): def from_operator(cls, operator, **kwargs): """Create a new instance by pulling data from the operator + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + Parameters ---------- operator : openmc.deplete.TransportOperator Operator with a depletion chain kwargs: optional - Optional keyword arguments to be passed to the underlying - ``__init__`` method + Additional keyword arguments to be used in constuction """ return cls(operator.chain.nuclides, **kwargs) @@ -560,9 +560,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): AttributeError If tallies not generated """ - if self._fission_rate_tally is None: - raise AttributeError( - "Tallies not built. Run generate_tallies first") + assert self._fission_rate_tally is not None, ( + "Run generate_tallies first") overlap = set(self._chain_nuclides).intersection(set(nuclides)) nuclides = tuple(sorted(overlap)) self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index a07630df76..ced3a2523a 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -166,8 +166,8 @@ class ConstantFissionYieldHelper(FissionYieldHelper): Parameters ---------- chain_nuclides : iterable of openmc.deplete.Nuclide - Nuclides tracked in the depletion chain. Not necessary - that all have yield data. + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. energy : float, optional Key in :attr:`openmc.deplete.Nuclide.yield_data` corresponding to the desired set of fission yield data. Typically one of @@ -204,22 +204,24 @@ class ConstantFissionYieldHelper(FissionYieldHelper): nuc.yield_data[nuc.yield_energies[min_index]]) @classmethod - def from_operator(cls, operator, energy=0.0253): + def from_operator(cls, operator, **kwargs): """Return a new ConstantFissionYieldHelper using operator data + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + Parameters ---------- operator : openmc.deplete.TransportOperator operator with a depletion chain - energy : float, optional - Energy for default fission yield libraries for nuclides with - multiple sets of yield data + kwargs: + Additional keyword arguments to be used in construction Returns ------- ConstantFissionYieldHelper """ - return cls(operator.chain.nuclides, energy=energy) + return cls(operator.chain.nuclides, **kwargs) @property def energy(self): @@ -255,8 +257,8 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): Parameters ---------- chain_nuclides : iterable of openmc.deplete.Nuclide - Nuclides tracked in the depletion chain. Not necessary - that all have yield data. + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. n_bmats : int, optional Number of burnable materials tracked in the problem cutoff : float, optional @@ -323,25 +325,18 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): self._fast_yields[name] = fast @classmethod - def from_operator(cls, operator, cutoff=112.0, - thermal_energy=0.0253, fast_energy=500e3): + def from_operator(cls, operator, **kwargs): """Construct a helper from an operator - All keyword arguments are identical to their counterpart + All keyword arguments should be identical to their counterpart in the main ``__init__`` method Parameters ---------- operator : openmc.deplete.Operator Operator with a chain and burnable materials - cutoff : float, optional - Cutoff energy for tallying fast and thermal fissions - thermal_energy : float, optional - Energy to use when pulling thermal fission yields from - nuclides with multiple sets of yields - fast_energy : float, optional - Energy to use when pulling fast fission yields from - nuclides with multiple sets of yields + kwargs: + Additional keyword arguments to be used in construction Returns ------- @@ -349,8 +344,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): """ return cls(operator.chain.nuclides, len(operator.burnable_mats), - cutoff=cutoff, thermal_energy=thermal_energy, - fast_energy=fast_energy) + **kwargs) @staticmethod def _find_fallback_energy(name, energies, cutoff, check_under): @@ -481,8 +475,8 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): Parameters ---------- chain_nuclides : iterable of openmc.deplete.Nuclide - Nuclides tracked in the depletion chain. Not necessary - that all have yield data. + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. Attributes ---------- @@ -615,13 +609,18 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): return mat_yields @classmethod - def from_operator(cls, operator): + def from_operator(cls, operator, **kwargs): """Return a new helper with data from an operator + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + Parameters ---------- operator : openmc.deplete.TransportOperator Operator with a depletion chain + kwargs : + Additional keyword arguments to be used in construction Returns ------- From 6735657cb33e028e647104e19859f59d6f77cf10 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 09:45:09 -0500 Subject: [PATCH 32/45] Don't append to Tally filters on FY helpers As pointed out in https://github.com/openmc-dev/openmc/pull/1313/files#r315844775, appending to Tally.filters does not propagate to the C++ side. New filter lists are created using the previous lists, and then passed onto to the Tally.filters setter. --- openmc/deplete/helpers.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index ced3a2523a..cc8b63c0fa 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -378,9 +378,9 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): in parallel mode. """ super().generate_tallies(materials, mat_indexes) - energy_filter = EnergyFilter() - energy_filter.bins = (0.0, self._cutoff, self._upper_energy) - self._fission_rate_tally.filters.append(energy_filter) + energy_filter = EnergyFilter(bins=[0.0, self._cutoff, self._upper_energy]) + self._fission_rate_tally.filters = ( + self._fission_rate_tally.filters + [energy_filter]) def unpack(self): """Obtain fast and thermal fission fractions from tally""" @@ -510,18 +510,16 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): """ super().generate_tallies(materials, mat_indexes) fission_tally = self._fission_rate_tally + filters = fission_tally.filters + ene_filter = EnergyFilter(bins=[0, self._upper_energy]) + fission_tally.filters = filters + [ene_filter] + + func_filter = EnergyFunctionFilter() + func_filter.set_data((0, self._upper_energy), (0, self._upper_energy)) weighted_tally = Tally() - weighted_tally.filters = fission_tally.filters weighted_tally.scores = ['fission'] - - ene_bin = EnergyFilter() - ene_bin.bins = (0, self._upper_energy) - fission_tally.filters.append(ene_bin) - - ene_filter = EnergyFunctionFilter() - ene_filter.set_data((0, self._upper_energy), (0, self._upper_energy)) - weighted_tally.filters.append(ene_filter) + weighted_tally.filters = filters + [func_filter] self._weighted_tally = weighted_tally def update_tally_nuclides(self, nuclides): From b19a2db0eef3e5cab2b252e15b3f9bbd806ebec4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 10:05:59 -0500 Subject: [PATCH 33/45] Respond to reviewer comments on FY helpers - Change Operator attribute name for fission yield helper - Use copy.deepcopy for constant_yields, fast_yields, and thermal_yields - Change import location of openmc in tests - Use bisect.bisect_left to find replacement fission yields on FissionYieldCutoffHelper if thermal or fast energies not found --- openmc/deplete/abc.py | 5 +- openmc/deplete/chain.py | 11 ++-- openmc/deplete/helpers.py | 55 +++++++------------ openmc/deplete/operator.py | 20 +++---- .../unit_tests/test_deplete_fission_yields.py | 5 +- tests/unit_tests/test_deplete_nuclide.py | 1 - 6 files changed, 37 insertions(+), 60 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 7cccad3695..17ed33e854 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -396,10 +396,7 @@ class FissionYieldHelper(ABC): @property def constant_yields(self): - out = {} - for key, sub in self._constant_yields.items(): - out[key] = sub.copy() - return out + return deepcopy(self._constant_yields) @abstractmethod def weighted_yields(self, local_mat_index): diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 27b73134b3..ae1c33f194 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -663,11 +663,8 @@ class Chain(object): @fission_yields.setter def fission_yields(self, yields): - if yields is None: - self._fission_yields = None - return - if isinstance(yields, Mapping): - self._fission_yields = [yields] - return - check_type("fission_yields", yields, Iterable, Mapping) + if yields is not None: + if isinstance(yields, Mapping): + yields = [yields] + check_type("fission_yields", yields, Iterable, Mapping) self._fission_yields = yields diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index cc8b63c0fa..37c8075d48 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,9 +1,10 @@ """ Class for normalizing fission energy deposition """ +from copy import deepcopy from itertools import product from numbers import Real -import operator +import bisect from numpy import dot, zeros, newaxis @@ -200,8 +201,8 @@ class ConstantFissionYieldHelper(FissionYieldHelper): distances = [abs(energy - ene) for ene in nuc.yield_energies] min_index = min( range(len(nuc.yield_energies)), key=distances.__getitem__) - self._constant_yields[name] = ( - nuc.yield_data[nuc.yield_energies[min_index]]) + min_E = min(nuc.yield_energies, key=lambda e: abs(e - energy)) + self._constant_yields[name] = nuc.yield_data[min_E] @classmethod def from_operator(cls, operator, **kwargs): @@ -310,17 +311,21 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): yields = nuc.yield_data energies = nuc.yield_energies thermal = yields.get(thermal_energy) - if thermal is None: - # find first index >= cutoff - ix = self._find_fallback_energy( - name, energies, cutoff, True) - thermal = yields[energies[ix - 1]] fast = yields.get(fast_energy) - if fast is None: - # find first index <= cutoff - rev_ix = self._find_fallback_energy( - name, reversed(energies), cutoff, False) - fast = yields[energies[-rev_ix]] + if thermal is None or fast is None: + insert_ix = bisect.bisect_left(energies, cutoff) + # if zero, then all energies > cutoff + # if len(energies), then all energies <= cutoff + if insert_ix == 0 or insert_ix == len(energies): + tail = map("{:5.3e}".format, energies) + raise ValueError( + "Cannot find replacement fission yields for {} given " + "cutoff {:5.3e} eV. Yields provided at [{}] eV".format( + name, cutoff, ", ".join(tail))) + if thermal is None: + thermal = yields[energies[insert_ix - 1]] + if fast is None: + fast = yields[energies[insert_ix]] self._thermal_yields[name] = thermal self._fast_yields[name] = fast @@ -346,20 +351,6 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): return cls(operator.chain.nuclides, len(operator.burnable_mats), **kwargs) - @staticmethod - def _find_fallback_energy(name, energies, cutoff, check_under): - cutoff_func = operator.ge if check_under else operator.le - found = False - for ix, ene in enumerate(energies): - if cutoff_func(ene, cutoff): - found = True - break - if found and ix != 0: - return ix - domain = "thermal" if check_under else "fast" - raise ValueError("Could not find {} yields for {} " - "with cutoff {} eV".format(domain, name, cutoff)) - def generate_tallies(self, materials, mat_indexes): """Use C API to produce a fission rate tally in burnable materials @@ -435,17 +426,11 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): @property def thermal_yields(self): - out = {} - for key, sub in self._thermal_yields.items(): - out[key] = sub.copy() - return out + return deepcopy(self._thermal_yields) @property def fast_yields(self): - out = {} - for key, sub in self._fast_yields.items(): - out[key] = sub.copy() - return out + return deepcopy(self._fast_yields) class AveragedFissionYieldHelper(TalliedFissionYieldHelper): diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index b7a4d2c8f0..43523b0fe1 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -137,7 +137,7 @@ class Operator(TransportOperator): diff_burnable_mats : bool Whether to differentiate burnable materials with multiple instances """ - _fission_helpers_ = { + _fission_helpers = { "average": AveragedFissionYieldHelper, "constant": ConstantFissionYieldHelper, "cutoff": FissionYieldCutoffHelper, @@ -147,10 +147,10 @@ class Operator(TransportOperator): diff_burnable_mats=False, fission_q=None, dilute_initial=1.0e3, fission_yield_mode="constant", fission_yield_opts=None): - if fission_yield_mode not in self._fission_helpers_: + if fission_yield_mode not in self._fission_helpers: raise KeyError( "fission_yield_mode must be one of {}, not {}".format( - ", ".join(self._fission_helpers_), fission_yield_mode)) + ", ".join(self._fission_helpers), fission_yield_mode)) super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None @@ -207,10 +207,10 @@ class Operator(TransportOperator): self._energy_helper = ChainFissionHelper() # Select and create fission yield helper - fission_helper = self._fission_helpers_[fission_yield_mode] + fission_helper = self._fission_helpers[fission_yield_mode] fission_yield_opts = ( {} if fission_yield_opts is None else fission_yield_opts) - self._fsn_yield_helper = fission_helper.from_operator( + self._yield_helper = fission_helper.from_operator( self, **fission_yield_opts) def __call__(self, vec, power): @@ -242,7 +242,7 @@ class Operator(TransportOperator): nuclides = self._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._energy_helper.nuclides = nuclides - self._fsn_yield_helper.update_tally_nuclides(nuclides) + self._yield_helper.update_tally_nuclides(nuclides) # Run OpenMC openmc.capi.reset() @@ -448,8 +448,8 @@ class Operator(TransportOperator): self.chain.nuclides, self.reaction_rates.index_nuc, materials) # Tell fission yield helper what materials this process is # responsible for - self._fsn_yield_helper.generate_tallies( - materials, tuple(sorted(self._mat_index_map.values()))) + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) @@ -595,7 +595,7 @@ class Operator(TransportOperator): # Keep track of energy produced from all reactions in eV per source # particle self._energy_helper.reset() - self._fsn_yield_helper.unpack() + self._yield_helper.unpack() # Store fission yield dictionaries fission_yields = [] @@ -622,7 +622,7 @@ class Operator(TransportOperator): mat_index, nuc_ind, react_ind) # Compute fission yields for this material - fission_yields.append(self._fsn_yield_helper.weighted_yields(i)) + fission_yields.append(self._yield_helper.weighted_yields(i)) # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 547f34c227..2ef7c90b93 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -5,7 +5,6 @@ from unittest.mock import Mock import pytest import numpy - from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( FissionYieldCutoffHelper, ConstantFissionYieldHelper, @@ -75,10 +74,10 @@ def test_cutoff_construction(nuclide_bundle): "U235": nuclide_bundle.u235.yield_data[0.0253]} assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} # test failures in cutoff: super low, super high - with pytest.raises(ValueError, match="thermal yields"): + with pytest.raises(ValueError, match="replacement fission yields"): FissionYieldCutoffHelper( nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002) - with pytest.raises(ValueError, match="fast yields"): + with pytest.raises(ValueError, match="replacement fission yields"): FissionYieldCutoffHelper( nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6) diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 2cfb98e6ca..f1dfa4f770 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -4,7 +4,6 @@ import xml.etree.ElementTree as ET import numpy import pytest - from openmc.deplete import nuclide From f97796211306197194d21e5b79209f94d963a2a9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 10:44:28 -0500 Subject: [PATCH 34/45] Support passing nuclide name to Nuclide objects Allow the creation of a Nuclide with >>> n = Nuclide("U235") --- openmc/deplete/chain.py | 3 +-- openmc/deplete/nuclide.py | 14 +++++++++----- tests/unit_tests/test_deplete_chain.py | 15 +++++---------- tests/unit_tests/test_deplete_fission_yields.py | 9 +++------ tests/unit_tests/test_deplete_nuclide.py | 3 +-- 5 files changed, 19 insertions(+), 25 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index ae1c33f194..d590cd8d2c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -207,8 +207,7 @@ class Chain(object): for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)): data = decay_data[parent] - nuclide = Nuclide() - nuclide.name = parent + nuclide = Nuclide(parent) chain.nuclides.append(nuclide) chain.nuclide_dict[parent] = idx diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index d2b0acc169..d4f501d442 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -68,11 +68,16 @@ except AttributeError: class Nuclide(object): """Decay modes, reactions, and fission yields for a single nuclide. + Parameters + ---------- + name : str, optional + GND name of this nuclide, e.g. ``"He4"``, ``"Am242_m1"`` + Attributes ---------- - name : str + name : str or None Name of nuclide. - half_life : float + half_life : float or None Half life of nuclide in [s]. decay_energy : float Energy deposited from decay in [eV]. @@ -91,12 +96,11 @@ class Nuclide(object): treated as a nested dictionary ``{energy: {product: yield}}`` yield_energies : tuple of float or None Energies at which fission product yields exist - """ - def __init__(self): + def __init__(self, name=None): # Information about the nuclide - self.name = None + self.name = name self.half_life = None self.decay_energy = 0.0 diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 7a5a128450..cf6414b5cc 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -141,8 +141,7 @@ def test_export_to_xml(run_in_tmpdir): # Prevent different MPI ranks from conflicting filename = 'test{}.xml'.format(comm.rank) - A = nuclide.Nuclide() - A.name = "A" + A = nuclide.Nuclide("A") A.half_life = 2.36520e4 A.decay_modes = [ nuclide.DecayTuple("beta1", "B", 0.6), @@ -150,14 +149,12 @@ def test_export_to_xml(run_in_tmpdir): ] A.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] - B = nuclide.Nuclide() - B.name = "B" + B = nuclide.Nuclide("B") B.half_life = 3.29040e4 B.decay_modes = [nuclide.DecayTuple("beta", "A", 1.0)] B.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] - C = nuclide.Nuclide() - C.name = "C" + C = nuclide.Nuclide("C") C.reactions = [ nuclide.ReactionTuple("fission", None, 2.0e8, 1.0), nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7), @@ -292,8 +289,7 @@ def test_capture_branch_infer_ground(): chain = Chain.from_xml(chain_file) # Create nuclide to be added into the chain - xe136m = nuclide.Nuclide() - xe136m.name = "Xe136_m1" + xe136m = nuclide.Nuclide("Xe136_m1") chain.nuclides.append(xe136m) chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1 @@ -310,8 +306,7 @@ def test_capture_branch_no_rxn(): chain_file = Path(__file__).parents[1] / "chain_simple.xml" chain = Chain.from_xml(chain_file) - u5m = nuclide.Nuclide() - u5m.name = "U235_m1" + u5m = nuclide.Nuclide("U235_m1") chain.nuclides.append(u5m) chain.nuclide_dict[u5m.name] = len(chain.nuclides) - 1 diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 2ef7c90b93..4b71eba696 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -20,17 +20,14 @@ def nuclide_bundle(): 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} - u235 = Nuclide() - u235.name = "U235" + u235 = Nuclide("U235") u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} - u238 = Nuclide() - u238.name = "U238" + u238 = Nuclide("U238") u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) - xe135 = Nuclide() - xe135.name = "Xe135" + xe135 = Nuclide("Xe135") NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") return NuclideBundle(u235, u238, xe135) diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index f1dfa4f770..4310478c96 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -84,8 +84,7 @@ def test_from_xml(): def test_to_xml_element(): """Test writing nuclide data to an XML element.""" - C = nuclide.Nuclide() - C.name = "C" + C = nuclide.Nuclide("C") C.half_life = 0.123 C.decay_modes = [ nuclide.DecayTuple('beta-', 'B', 0.99), From 60f89794f684a30a5aed43868deef8e8b45fd847 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 10:55:43 -0500 Subject: [PATCH 35/45] Use dictionary to create FissionYieldDistribution Removes the FissionYieldDistribution.from_dict method as the most natural way to create such a distribution is with a dictionary. --- openmc/deplete/chain.py | 13 ++- openmc/deplete/nuclide.py | 79 ++++++------------- tests/unit_tests/test_deplete_chain.py | 2 +- .../unit_tests/test_deplete_fission_yields.py | 4 +- tests/unit_tests/test_deplete_nuclide.py | 6 +- 5 files changed, 38 insertions(+), 66 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index d590cd8d2c..ba959a3fcc 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -10,7 +10,6 @@ import math import re from collections import OrderedDict, defaultdict from collections.abc import Mapping, Iterable -from numbers import Real from warnings import warn from openmc.checkvalue import check_type, check_less_than @@ -127,8 +126,8 @@ class Chain(object): 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``. + 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 @@ -222,7 +221,8 @@ class Chain(object): if mode.daughter in decay_data: target = mode.daughter else: - print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter)) + 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 @@ -285,7 +285,7 @@ class Chain(object): yield_replace = 0.0 yields = defaultdict(float) for product, y in table.items(): - # Handle fission products that have no decay data available + # Handle fission products that have no decay data if product not in decay_data: daughter = replace_missing(product, decay_data) product = daughter @@ -297,8 +297,7 @@ class Chain(object): missing_fp.append((parent, E, yield_replace)) yield_data[E] = yields - nuclide.yield_data = ( - FissionYieldDistribution.from_dict(yield_data)) + nuclide.yield_data = FissionYieldDistribution(yield_data) # Display warnings if missing_daughter: diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index d4f501d442..a3c60019e3 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -3,16 +3,15 @@ Contains the per-nuclide components of a depletion chain. """ -from numbers import Real from collections import namedtuple -from collections.abc import Iterable, Mapping +from collections.abc import Mapping try: import lxml.etree as ET except ImportError: import xml.etree.ElementTree as ET -from numpy import asarray, empty +from numpy import empty from openmc.checkvalue import check_type @@ -132,7 +131,7 @@ class Nuclide(object): check_type("fission_yields", fission_yields, Mapping) if isinstance(fission_yields, FissionYieldDistribution): self._yield_data = fission_yields - self._yield_data = FissionYieldDistribution.from_dict(fission_yields) + self._yield_data = FissionYieldDistribution(fission_yields) @property def yield_energies(self): @@ -245,7 +244,7 @@ class FissionYieldDistribution(Mapping): Can be used as a dictionary mapping energies and products to fission yields:: - >>> fydist = FissionYieldDistribution.from_dict({ + >>> fydist = FissionYieldDistribution{ ... {0.0253: {"Xe135": 0.021}}) >>> fydist[0.0253]["Xe135"] 0.021 @@ -266,11 +265,10 @@ class FissionYieldDistribution(Mapping): Attributes ---------- energies : tuple - Energies for which fission yields exist. Converted for - indexing + Energies for which fission yields exist. Sorted by + increasing energy products : tuple - Fission products produced at all energies. Converted - for indexing + Fission products produced at all energies. Sorted by name. yield_matrix : numpy.ndarray Array ``(n_energy, n_products)`` where ``yield_matrix[g, j]`` is the fission yield of product @@ -278,22 +276,28 @@ class FissionYieldDistribution(Mapping): See Also -------- - * :meth:`from_xml_element`, :meth:`from_dict` - Construction methods + * :meth:`from_xml_element` - Construction methods * :class:`FissionYield` - Class used for storing yields at a given energy """ - def __init__(self, energies, products, yield_matrix): - check_type("energies", energies, Iterable, Real) + def __init__(self, fission_yields): + # mapping {energy: {product: value}} + energies = sorted(fission_yields) + + # Get a consistent set of products to produce a matrix of yields + shared_prod = set.union(*(set(x) for x in fission_yields.values())) + ordered_prod = sorted(shared_prod) + + yield_matrix = empty((len(energies), len(shared_prod))) + + 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) self.energies = tuple(energies) - check_type("products", products, Iterable, str) - self.products = tuple(products) - yield_matrix = asarray(yield_matrix, dtype=float) - if yield_matrix.shape != (len(self.energies), len(self.products)): - raise ValueError( - "Shape of yield matrix inconsistent. " - "Should be ({}, {}), is {}".format( - len(energies), len(products), - yield_matrix.shape)) + self.products = tuple(ordered_prod) self.yield_matrix = yield_matrix def __len__(self): @@ -329,38 +333,7 @@ class FissionYieldDistribution(Mapping): # Get a map of products to their corresponding yield all_yields[energy] = dict(zip(products, yields)) - return cls.from_dict(all_yields) - - @classmethod - def from_dict(cls, fission_yields): - """Construct a distribution from a dictionary of yields - - Parameters - ----------- - fission_yields : dict - Dictionary ``{energy: {product: yield}}`` - - Returns - ------- - FissionYieldDistribution - """ - # mapping {energy: {product: value}} - energies = sorted(fission_yields) - - # Get a consistent set of products to produce a matrix of yields - shared_prod = set.union(*(set(x) for x in fission_yields.values())) - ordered_prod = sorted(shared_prod) - - yield_matrix = empty((len(energies), len(shared_prod))) - - 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) - - return cls(energies, ordered_prod, yield_matrix) + return cls(all_yields) def to_xml_element(self, root): """Write fission yield data to an xml element diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index cf6414b5cc..b1ab1ae248 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -160,7 +160,7 @@ def test_export_to_xml(run_in_tmpdir): nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7), nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3) ] - C.yield_data = nuclide.FissionYieldDistribution.from_dict({ + C.yield_data = nuclide.FissionYieldDistribution({ 0.0253: {"A": 0.0292737, "B": 0.002566345}}) chain = Chain() diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 4b71eba696..81c3e2d420 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -21,11 +21,11 @@ def nuclide_bundle(): 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} u235 = Nuclide("U235") - u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) + u235.yield_data = FissionYieldDistribution(u5yield_dict) u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} u238 = Nuclide("U238") - u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) + u238.yield_data = FissionYieldDistribution(u8yield_dict) xe135 = Nuclide("Xe135") diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 4310478c96..3386e28984 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -71,7 +71,7 @@ def test_from_xml(): nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0), nuclide.ReactionTuple('fission', None, 193405400.0, 1.0), ] - expected_yield_data = nuclide.FissionYieldDistribution.from_dict({ + expected_yield_data = nuclide.FissionYieldDistribution({ 0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}}) assert u235.yield_data == expected_yield_data # test accessing the yield energies through the FissionYieldDistribution @@ -94,7 +94,7 @@ def test_to_xml_element(): nuclide.ReactionTuple('fission', None, 2.0e8, 1.0), nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0) ] - C.yield_data = nuclide.FissionYieldDistribution.from_dict( + C.yield_data = nuclide.FissionYieldDistribution( {0.0253: {"A": 0.0292737, "B": 0.002566345}}) element = C.to_xml_element() @@ -127,7 +127,7 @@ def test_fission_yield_distribution(): 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8, "Sm149": 2.69e-8}, 5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}, # drop Sm149 } - yield_dist = nuclide.FissionYieldDistribution.from_dict(yield_dict) + yield_dist = nuclide.FissionYieldDistribution(yield_dict) assert len(yield_dist) == len(yield_dict) assert yield_dist.energies == tuple(sorted(yield_dict.keys())) for exp_ene, exp_dist in yield_dict.items(): From 687d287b243aa0ac45e3eddd9a424519312d3202 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 11:34:54 -0500 Subject: [PATCH 36/45] Allow Nuclide.yield_data to be None Had to guard against some cases where it was assumed that yield_data would be a dictionary. For the most part, the fission yield helpers act on nuclides that have yield data already. The previous default state was to use an empty dictionary. This has been removed. A value of None for nuclide.yield_data indicates there is no yield data. --- openmc/deplete/abc.py | 2 ++ openmc/deplete/chain.py | 2 +- openmc/deplete/nuclide.py | 14 +++++++++----- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 17ed33e854..16d2d7abbd 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -387,6 +387,8 @@ class FissionYieldHelper(ABC): # Get all nuclides with fission yield data for nuc in chain_nuclides: + if nuc.yield_data is None: + continue if len(nuc.yield_data) == 1: self._constant_yields[nuc.name] = ( nuc.yield_data[nuc.yield_energies[0]]) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index ba959a3fcc..e38385ac52 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -400,7 +400,7 @@ class Chain(object): """ out = {} for nuc in self.nuclides: - if not nuc.yield_data: + if nuc.yield_data is None: continue yield_obj = nuc.yield_data[min(nuc.yield_energies)] out[nuc.name] = dict(yield_obj) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index a3c60019e3..068891a5f0 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -123,15 +123,19 @@ class Nuclide(object): @property def yield_data(self): if self._yield_data is None: - return {} + return None return self._yield_data @yield_data.setter def yield_data(self, fission_yields): - check_type("fission_yields", fission_yields, Mapping) - if isinstance(fission_yields, FissionYieldDistribution): - self._yield_data = fission_yields - self._yield_data = FissionYieldDistribution(fission_yields) + if fission_yields is None: + self._yield_data = None + else: + check_type("fission_yields", fission_yields, Mapping) + if isinstance(fission_yields, FissionYieldDistribution): + self._yield_data = fission_yields + else: + self._yield_data = FissionYieldDistribution(fission_yields) @property def yield_energies(self): From 8e2bb0bbbc99d6caf93b2ca1cdcc047d7fbe8346 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 17:20:19 -0500 Subject: [PATCH 37/45] Improve testing for FY helpers using C API Provide a module-scoped fixture that creates a temporary directory, adds simple settings, geometry, and material files, and initializes the openmc C API. This fixture also provides two empty openmc.capi.Material objects to help the helpers generate meaningful tallies. The main tests for AveragedFissionYieldHelper and FissionYieldCutoffHelper now go through the built tallies and examine the filters, nuclides, and scores. A helper function produces mocked-like tally data based on filters, nuclides, and scores found on a tally. Pu239 with 0.0253 eV, 500 keV, and 2 MeV yields is included in the nuclide_bundle fixture. This provides some additional heterogeneity in the results, as the 2 MeV data is not present on U235. --- openmc/deplete/helpers.py | 17 +- .../unit_tests/test_deplete_fission_yields.py | 297 +++++++++++++----- 2 files changed, 230 insertions(+), 84 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 37c8075d48..debc10403e 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -313,19 +313,24 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): thermal = yields.get(thermal_energy) fast = yields.get(fast_energy) if thermal is None or fast is None: - insert_ix = bisect.bisect_left(energies, cutoff) + cutoff_ix = bisect.bisect_left(energies, cutoff) # if zero, then all energies > cutoff # if len(energies), then all energies <= cutoff - if insert_ix == 0 or insert_ix == len(energies): + if cutoff_ix == 0 or cutoff_ix == len(energies): tail = map("{:5.3e}".format, energies) raise ValueError( "Cannot find replacement fission yields for {} given " "cutoff {:5.3e} eV. Yields provided at [{}] eV".format( name, cutoff, ", ".join(tail))) + # find closest energy to requested thermal, fast energies if thermal is None: - thermal = yields[energies[insert_ix - 1]] + min_E = min(energies[:cutoff_ix], + key=lambda e: abs(e - thermal_energy)) + thermal = yields[min_E] if fast is None: - fast = yields[energies[insert_ix]] + min_E = min(energies[cutoff_ix:], + key=lambda e: abs(e - fast_energy)) + fast = yields[min_E] self._thermal_yields[name] = thermal self._fast_yields[name] = fast @@ -369,7 +374,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): in parallel mode. """ super().generate_tallies(materials, mat_indexes) - energy_filter = EnergyFilter(bins=[0.0, self._cutoff, self._upper_energy]) + energy_filter = EnergyFilter([0.0, self._cutoff, self._upper_energy]) self._fission_rate_tally.filters = ( self._fission_rate_tally.filters + [energy_filter]) @@ -497,7 +502,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): fission_tally = self._fission_rate_tally filters = fission_tally.filters - ene_filter = EnergyFilter(bins=[0, self._upper_energy]) + ene_filter = EnergyFilter([0, self._upper_energy]) fission_tally.filters = filters + [ene_filter] func_filter = EnergyFunctionFilter() diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 81c3e2d420..05f62eaceb 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -1,17 +1,97 @@ """Test the FissionYieldHelpers""" +import os from collections import namedtuple from unittest.mock import Mock +import bisect import pytest import numpy +from openmc import capi from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( FissionYieldCutoffHelper, ConstantFissionYieldHelper, AveragedFissionYieldHelper) -MATERIALS = ["1", "2"] +@pytest.fixture(scope="module") +def materials(tmpdir_factory): + """Use C API to construct realistic materials for testing tallies""" + tmpdir = tmpdir_factory.mktemp("capi") + orig = tmpdir.chdir() + # Create proxy xml files to please openmc + with open("geometry.xml", "w") as stream: + stream.write(""" + + + + + + + + +""") + with open("settings.xml", "w") as stream: + stream.write(""" + + + eigenvalue + 100 + 10 + 0 + 1 + + + 0.0 0.0 0.0 + + + +""") + with open("materials.xml", "w") as stream: + stream.write(""" + + + + + + + + + +""") + try: + with capi.run_in_memory(): + yield [capi.Material(), capi.Material()] + finally: + print(os.path.abspath(os.curdir)) + os.remove(tmpdir / "settings.xml") + os.remove(tmpdir / "geometry.xml") + os.remove(tmpdir / "materials.xml") + os.remove(tmpdir / "summary.h5") + orig.chdir() + os.rmdir(tmpdir) + + +def proxy_tally_data(tally, fill=None): + """Construct an empty matrix built from a C tally + + The shape of tally.results will be + ``(n_bins, n_nuc * n_scores, 3)`` + """ + n_nucs = max(len(tally.nuclides), 1) + n_scores = max(len(tally.scores), 1) + n_bins = 1 + for tfilter in tally.filters: + if not hasattr(tfilter, "bins"): + continue + this_bins = len(tfilter.bins) + if isinstance(tfilter, capi.EnergyFilter): + this_bins -= 1 + n_bins *= max(this_bins, 1) + data = numpy.empty((n_bins, n_nucs * n_scores, 3)) + if fill is not None: + data.fill(fill) + return data @pytest.fixture(scope="module") @@ -29,47 +109,74 @@ def nuclide_bundle(): xe135 = Nuclide("Xe135") - NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") - return NuclideBundle(u235, u238, xe135) + pu239 = Nuclide("Pu239") + pu239.yield_data = FissionYieldDistribution({ + 0.0253: {"Xe135": 3.141e-3, "Sm149": 8.19e-10, "Gd155": 1.66e-9}, + 5.0e5: {"Xe135": 6.14e-3, "Sm149": 9.429e-10, "Gd155": 5.24e-9}, + 2e6: {"Xe135": 6.15e-3, "Sm149": 9.42e-10, "Gd155": 5.29e-9}}) + + NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135 pu239") + return NuclideBundle(u235, u238, xe135, pu239) @pytest.mark.parametrize( - "input_energy, u5_yield_energy", + "input_energy, yield_energy", ((0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) -def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): +def test_constant_helper(nuclide_bundle, input_energy, yield_energy): helper = ConstantFissionYieldHelper(nuclide_bundle, energy=input_energy) assert helper.energy == input_energy assert helper.constant_yields == { - "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], - "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal + "U235": nuclide_bundle.u235.yield_data[yield_energy], + "U238": nuclide_bundle.u238.yield_data[5.00e5], # only epithermal + "Pu239": nuclide_bundle.pu239.yield_data[yield_energy]} assert helper.constant_yields == helper.weighted_yields(1) def test_cutoff_construction(nuclide_bundle): + u235 = nuclide_bundle.u235 + u238 = nuclide_bundle.u238 + pu239 = nuclide_bundle.pu239 + # defaults helper = FissionYieldCutoffHelper(nuclide_bundle, 1) assert helper.constant_yields == { - "U238": nuclide_bundle.u238.yield_data[5.0e5]} + "U238": u238.yield_data[5.0e5]} assert helper.thermal_yields == { - "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} + "U235": u235.yield_data[0.0253], + "Pu239": pu239.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": u235.yield_data[5e5], + "Pu239": pu239.yield_data[5e5]} + # use 14 MeV yields helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6) assert helper.constant_yields == { - "U238": nuclide_bundle.u238.yield_data[5.0e5]} + "U238": u238.yield_data[5.0e5]} assert helper.thermal_yields == { - "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[14e6]} + "U235": u235.yield_data[0.0253], + "Pu239": pu239.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": u235.yield_data[14e6], + "Pu239": pu239.yield_data[2e6]} + # specify missing thermal yields -> use 0.0253 helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1) assert helper.thermal_yields == { - "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} + "U235": u235.yield_data[0.0253], + "Pu239": pu239.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": u235.yield_data[5e5], + "Pu239": pu239.yield_data[5e5]} + # request missing fast yields -> use epithermal helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4) assert helper.thermal_yields == { - "U235": nuclide_bundle.u235.yield_data[0.0253]} - assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]} + "U235": u235.yield_data[0.0253], + "Pu239": pu239.yield_data[0.0253]} + assert helper.fast_yields == { + "U235": u235.yield_data[5e5], + "Pu239": pu239.yield_data[5e5]} + # test failures in cutoff: super low, super high with pytest.raises(ValueError, match="replacement fission yields"): FissionYieldCutoffHelper( @@ -87,83 +194,117 @@ def test_cutoff_failure(key): FissionYieldCutoffHelper(None, None, **{key: -1}) -class ProxyMixin: - """Mixing that overloads the tally generation""" - def generate_tallies(self, materials, mat_indexes): - self._fission_rate_tally = Mock() - self._local_indexes = numpy.asarray(mat_indexes) - - -class CutoffProxy(ProxyMixin, FissionYieldCutoffHelper): - """Proxy that supplies a set of tallies""" - - # emulate some split between fast and thermal U235 fissions @pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) -def test_cutoff_helper(nuclide_bundle, therm_frac): - n_bmats = len(MATERIALS) - proxy = CutoffProxy(nuclide_bundle, n_bmats) - proxy.generate_tallies(MATERIALS, [0]) +def test_cutoff_helper(materials, nuclide_bundle, therm_frac): + helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials)) + helper.generate_tallies(materials, [0]) + non_zero_nucs = [n.name for n in nuclide_bundle] - tally_nucs = proxy.update_tally_nuclides(non_zero_nucs) - assert tally_nucs == ("U235",) + tally_nucs = helper.update_tally_nuclides(non_zero_nucs) + assert tally_nucs == ("Pu239", "U235",) + + # Check tallies + fission_tally = helper._fission_rate_tally + assert fission_tally is not None + filters = fission_tally.filters + assert len(filters) == 2 + assert isinstance(filters[0], capi.MaterialFilter) + assert len(filters[0].bins) == len(materials) + assert isinstance(filters[1], capi.EnergyFilter) + # lower, cutoff, and upper energy + assert len(filters[1].bins) == 3 + # Emulate building tallies # material x energy, tallied_nuclides, 3 - proxy_flux = 1e6 - tally_data = numpy.empty((n_bmats * 2, 1, 3)) - tally_data[0, 0, 1] = therm_frac * proxy_flux - tally_data[1, 0, 1] = (1 - therm_frac) * proxy_flux - proxy._fission_rate_tally.results = tally_data + tally_data = proxy_tally_data(fission_tally) + helper._fission_rate_tally = Mock() + helper_flux = 1e6 + tally_data[0, :, 1] = therm_frac * helper_flux + tally_data[1, :, 1] = (1 - therm_frac) * helper_flux + helper._fission_rate_tally.results = tally_data - proxy.unpack() + helper.unpack() # expected results of shape (n_mats, 2, n_tnucs) - expected_results = numpy.empty((1, 2, 1)) + expected_results = numpy.empty((1, 2, len(tally_nucs))) expected_results[:, 0] = therm_frac expected_results[:, 1] = 1 - therm_frac - assert proxy.results == pytest.approx(expected_results) + assert helper.results == pytest.approx(expected_results) - actual_yields = proxy.weighted_yields(0) + actual_yields = helper.weighted_yields(0) assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] - assert actual_yields["U235"] == ( - proxy.thermal_yields["U235"] * therm_frac - + proxy.fast_yields["U235"] * (1 - therm_frac)) + for nuc in tally_nucs: + assert actual_yields[nuc] == ( + helper.thermal_yields[nuc] * therm_frac + + helper.fast_yields[nuc] * (1 - therm_frac)) -class AverageProxy(ProxyMixin, AveragedFissionYieldHelper): - """Proxy for generating mock set of tallies""" - def generate_tallies(self, materials, mat_indexes): - super().generate_tallies(materials, mat_indexes) - self._weighted_tally = Mock() - - -@pytest.mark.parametrize("avg_energy", (0.01, 100, 15e6)) -def test_averaged_helper(nuclide_bundle, avg_energy): - proxy = AverageProxy(nuclide_bundle) - proxy.generate_tallies(MATERIALS, [0]) - tallied_nucs = proxy.update_tally_nuclides( +@pytest.mark.parametrize("avg_energy", (0.01, 6e5, 15e6)) +def test_averaged_helper(materials, nuclide_bundle, avg_energy): + helper = AveragedFissionYieldHelper(nuclide_bundle) + helper.generate_tallies(materials, [0]) + tallied_nucs = helper.update_tally_nuclides( [n.name for n in nuclide_bundle]) - assert tallied_nucs == ("U235", ) - # enforce some average energy - proxy_flux = 1e16 - fission_results = numpy.ones((len(MATERIALS), 1, 3)) * proxy_flux - weighted_results = fission_results * avg_energy - proxy._fission_rate_tally.results = fission_results - proxy._weighted_tally.results = weighted_results - proxy.unpack() - expected_results = numpy.ones((1, 1)) * avg_energy - assert proxy.results == pytest.approx(expected_results) + assert tallied_nucs == ("Pu239", "U235") - actual_yields = proxy.weighted_yields(0) + # check generated tallies + fission_tally = helper._fission_rate_tally + assert fission_tally is not None + fission_filters = fission_tally.filters + assert len(fission_filters) == 2 + assert isinstance(fission_filters[0], capi.MaterialFilter) + assert len(fission_filters[0].bins) == len(materials) + assert isinstance(fission_filters[1], capi.EnergyFilter) + assert len(fission_filters[1].bins) == 2 + assert fission_tally.scores == ["fission"] + assert fission_tally.nuclides == list(tallied_nucs) + + weighted_tally = helper._weighted_tally + assert weighted_tally is not None + weighted_filters = weighted_tally.filters + assert len(weighted_filters) == 2 + assert isinstance(weighted_filters[0], capi.MaterialFilter) + assert len(weighted_filters[0].bins) == len(materials) + assert isinstance(weighted_filters[1], capi.EnergyFunctionFilter) + assert len(weighted_filters[1].energy) == 2 + assert len(weighted_filters[1].y) == 2 + assert weighted_tally.scores == ["fission"] + assert weighted_tally.nuclides == list(tallied_nucs) + + helper_flux = 1e16 + fission_results = proxy_tally_data(fission_tally, helper_flux) + weighted_results = proxy_tally_data( + weighted_tally, helper_flux * avg_energy) + + helper._fission_rate_tally = Mock() + helper._weighted_tally = Mock() + helper._fission_rate_tally.results = fission_results + helper._weighted_tally.results = weighted_results + + helper.unpack() + expected_results = numpy.ones((1, len(tallied_nucs))) * avg_energy + assert helper.results == pytest.approx(expected_results) + + actual_yields = helper.weighted_yields(0) # constant U238 => no interpolation assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] # construct expected yields - if avg_energy < 0.0253: # take thermal U235 yields - exp_u235_yields = nuclide_bundle.u235.yield_data[0.0253] - elif avg_energy > 14e6: # take fastest U235 yields - exp_u235_yields = nuclide_bundle.u235.yield_data[14e6] - else: # reconstruct between thermal and epithermal - thermal = nuclide_bundle.u235.yield_data[0.0253] - epithermal = nuclide_bundle.u235.yield_data[5e5] - split = (avg_energy - 0.0253) / (5e5 - 0.0253) - exp_u235_yields = thermal * (1 - split) + epithermal * split + exp_u235_yields = interp_average_yields(nuclide_bundle.u235, avg_energy) assert actual_yields["U235"] == exp_u235_yields + exp_pu239_yields = interp_average_yields(nuclide_bundle.pu239, avg_energy) + assert actual_yields["Pu239"] == exp_pu239_yields + + +def interp_average_yields(nuc, avg_energy): + """Construct a set of yields by interpolation between neighbors""" + energies = nuc.yield_energies + yields = nuc.yield_data + if avg_energy < energies[0]: + return yields[energies[0]] + if avg_energy > energies[-1]: + return yields[energies[-1]] + thermal_ix = bisect.bisect_left(energies, avg_energy) + thermal_E, fast_E = energies[thermal_ix - 1:thermal_ix + 1] + assert thermal_E < avg_energy < fast_E + split = (avg_energy - thermal_E)/(fast_E - thermal_E) + return yields[thermal_E]*(1 - split) + yields[fast_E]*split From 4f46aa895c5abb84ef4827a75709e529d86d03f3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Aug 2019 17:46:51 -0500 Subject: [PATCH 38/45] Teach deplete.Nuclide.validate about new fission yields Reconfigure tests/unit_tests/test_deplete_chain.py and tests/unit_tests/test_deplete_nuclide.py to use the new dictionary-like representation of fission yields. --- openmc/deplete/nuclide.py | 4 ++-- tests/unit_tests/test_deplete_chain.py | 2 +- tests/unit_tests/test_deplete_nuclide.py | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index db503c8a59..a989bf7029 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -321,8 +321,8 @@ class Nuclide(object): valid = False if self.yield_data: - for energy, yield_list in self.yield_data.items(): - sum_yield = sum(y[1] for y in yield_list) + for energy, fission_yield in self.yield_data.items(): + sum_yield = fission_yield.yields.sum() stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance if stat: continue diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index ce45b84d27..bdaf6c5ab1 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -377,7 +377,7 @@ def test_validate(simple_chain): # Fix fission yields but keep to restore later old_yields = simple_chain["C"].yield_data - simple_chain["C"].yield_data = {0.0253: [("A", 1.4), ("B", 0.6)]} + simple_chain["C"].yield_data = {0.0253: {"A": 1.4, "B": 0.6}} assert simple_chain.validate(strict=True, tolerance=0.0) with pytest.warns(None) as record: diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index d81691932b..6f5356b129 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -177,8 +177,8 @@ def test_validate(): # fission yields nuc.yield_data = { - 0.0253: [("0", 1.5), ("1", 0.5)], - 1e6: [("0", 1.5), ("1", 0.5)], + 0.0253: {"0": 1.5, "1": 0.5}, + 1e6: {"0": 1.5, "1": 0.5}, } # nuclide is good and should have no warnings raise @@ -212,7 +212,7 @@ def test_validate(): # restore reactions, invalidate fission yields nuc.reactions.append(reaction) - nuc.yield_data[1e6].pop() + nuc.yield_data[1e6].yields *= 2 with pytest.raises(ValueError, match=r"fission yields.*1\.0*e"): nuc.validate(strict=True, quiet=False, tolerance=0.0) From 7b175961fa45bef717b30894f1659d0f2f436a57 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 11:24:39 -0500 Subject: [PATCH 39/45] Use constant FY if cutoff energy outside provided yields The FissionYieldCutoffHelper will always find a set of yields to use for all nuclides with yield data now, rather than raise an error. Previously, if the cutoff was outside the provided bounds, an error was raised because there wasn't a clear set of "fast" and "thermal" yields to use. However, this caused issues with nuclides that are missing a set of lower yields, like Th232 with yields at 5e5 and 6e6 per ENDF/B-VII.1 data. Now, if the cutoff energy is outside the bounds of provided yield data, the closet set of yields to the cutoff is taken to be constant. These nuclides will not be tallied during the transport routine. --- openmc/deplete/helpers.py | 23 +++--- .../unit_tests/test_deplete_fission_yields.py | 70 ++++++++++--------- 2 files changed, 50 insertions(+), 43 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index debc10403e..eb86b0ea4a 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -199,8 +199,6 @@ class ConstantFissionYieldHelper(FissionYieldHelper): continue # Specific energy not found, use closest energy distances = [abs(energy - ene) for ene in nuc.yield_energies] - min_index = min( - range(len(nuc.yield_energies)), key=distances.__getitem__) min_E = min(nuc.yield_energies, key=lambda e: abs(e - energy)) self._constant_yields[name] = nuc.yield_data[min_E] @@ -307,21 +305,24 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): self._cutoff = cutoff self._thermal_yields = {} self._fast_yields = {} + convert_to_constant = set() for name, nuc in self._chain_nuclides.items(): yields = nuc.yield_data energies = nuc.yield_energies thermal = yields.get(thermal_energy) fast = yields.get(fast_energy) if thermal is None or fast is None: + if cutoff <= energies[0]: + # use lowest energy yields as constant + self._constant_yields[name] = yields[energies[0]] + convert_to_constant.add(name) + continue + if cutoff >= energies[-1]: + # use highest energy yields as constant + self._constant_yields[name] = yields[energies[-1]] + convert_to_constant.add(name) + continue cutoff_ix = bisect.bisect_left(energies, cutoff) - # if zero, then all energies > cutoff - # if len(energies), then all energies <= cutoff - if cutoff_ix == 0 or cutoff_ix == len(energies): - tail = map("{:5.3e}".format, energies) - raise ValueError( - "Cannot find replacement fission yields for {} given " - "cutoff {:5.3e} eV. Yields provided at [{}] eV".format( - name, cutoff, ", ".join(tail))) # find closest energy to requested thermal, fast energies if thermal is None: min_E = min(energies[:cutoff_ix], @@ -333,6 +334,8 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): fast = yields[min_E] self._thermal_yields[name] = thermal self._fast_yields[name] = fast + for name in convert_to_constant: + self._chain_nuclides.pop(name) @classmethod def from_operator(cls, operator, **kwargs): diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 05f62eaceb..2842043d12 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -111,7 +111,6 @@ def nuclide_bundle(): pu239 = Nuclide("Pu239") pu239.yield_data = FissionYieldDistribution({ - 0.0253: {"Xe135": 3.141e-3, "Sm149": 8.19e-10, "Gd155": 1.66e-9}, 5.0e5: {"Xe135": 6.14e-3, "Sm149": 9.429e-10, "Gd155": 5.24e-9}, 2e6: {"Xe135": 6.15e-3, "Sm149": 9.42e-10, "Gd155": 5.29e-9}}) @@ -127,8 +126,8 @@ def test_constant_helper(nuclide_bundle, input_energy, yield_energy): assert helper.energy == input_energy assert helper.constant_yields == { "U235": nuclide_bundle.u235.yield_data[yield_energy], - "U238": nuclide_bundle.u238.yield_data[5.00e5], # only epithermal - "Pu239": nuclide_bundle.pu239.yield_data[yield_energy]} + "U238": nuclide_bundle.u238.yield_data[5.00e5], + "Pu239": nuclide_bundle.pu239.yield_data[5e5]} assert helper.constant_yields == helper.weighted_yields(1) @@ -140,50 +139,54 @@ def test_cutoff_construction(nuclide_bundle): # defaults helper = FissionYieldCutoffHelper(nuclide_bundle, 1) assert helper.constant_yields == { - "U238": u238.yield_data[5.0e5]} - assert helper.thermal_yields == { - "U235": u235.yield_data[0.0253], - "Pu239": pu239.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": u235.yield_data[5e5], + "U238": u238.yield_data[5.0e5], "Pu239": pu239.yield_data[5e5]} + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} # use 14 MeV yields helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6) assert helper.constant_yields == { - "U238": u238.yield_data[5.0e5]} - assert helper.thermal_yields == { - "U235": u235.yield_data[0.0253], - "Pu239": pu239.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": u235.yield_data[14e6], - "Pu239": pu239.yield_data[2e6]} + "U238": u238.yield_data[5.0e5], + "Pu239": pu239.yield_data[5e5]} + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[14e6]} # specify missing thermal yields -> use 0.0253 helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1) - assert helper.thermal_yields == { - "U235": u235.yield_data[0.0253], - "Pu239": pu239.yield_data[0.0253]} - assert helper.fast_yields == { - "U235": u235.yield_data[5e5], - "Pu239": pu239.yield_data[5e5]} + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} # request missing fast yields -> use epithermal helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4) + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} + + # higher cutoff energy -> obtain fast and "faster" yields + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, cutoff=1e6, + thermal_energy=5e5, fast_energy=14e6) + assert helper.constant_yields == {"U238": u238.yield_data[5e5]} assert helper.thermal_yields == { - "U235": u235.yield_data[0.0253], - "Pu239": pu239.yield_data[0.0253]} + "U235": u235.yield_data[5e5], "Pu239": pu239.yield_data[5e5]} assert helper.fast_yields == { - "U235": u235.yield_data[5e5], + "U235": u235.yield_data[14e6], "Pu239": pu239.yield_data[2e6]} + + # test super low and super high cutoff energies + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002) + assert helper.fast_yields == {} + assert helper.thermal_yields == {} + assert helper.constant_yields == { + "U235": u235.yield_data[0.0253], "U238": u238.yield_data[5e5], "Pu239": pu239.yield_data[5e5]} - # test failures in cutoff: super low, super high - with pytest.raises(ValueError, match="replacement fission yields"): - FissionYieldCutoffHelper( - nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002) - with pytest.raises(ValueError, match="replacement fission yields"): - FissionYieldCutoffHelper( - nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6) + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6) + assert helper.thermal_yields == {} + assert helper.fast_yields == {} + assert helper.constant_yields == { + "U235": u235.yield_data[14e6], "U238": u238.yield_data[5e5], + "Pu239": pu239.yield_data[2e6]} @pytest.mark.parametrize("key", ("cutoff", "thermal_energy", "fast_energy")) @@ -197,7 +200,8 @@ def test_cutoff_failure(key): # emulate some split between fast and thermal U235 fissions @pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) def test_cutoff_helper(materials, nuclide_bundle, therm_frac): - helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials)) + helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials), + cutoff=1e6, fast_energy=14e6) helper.generate_tallies(materials, [0]) non_zero_nucs = [n.name for n in nuclide_bundle] From 91c08ef3232b5789532f788ee679f37027ad5664 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 15:20:56 -0500 Subject: [PATCH 40/45] Improve/expand special methods for FissionYield Remove copy method. __contains__ and __getitem__ now take advantage of the ordering of products. This can improve time searching for fission products and retrieving yields. radd and rmul methods defer to their left counterparts, e.g. x * fy => fy * x. Return NotImplemented types if addition is not done with another set of fission yields and multiplication is not done with a scalar. Improve/expand special methods for FissionYield Remove copy method. __contains__ and __getitem__ now take advantage of the ordering of products. This can improve time searching for fission products and retrieving yields. radd and rmul methods defer to their left counterparts, e.g. x * fy => fy * x. Return NotImplemented types if addition is not done with another set of fission yields and multiplication is not done with a scalar. --- openmc/deplete/nuclide.py | 39 ++++++++++++++++-------- tests/unit_tests/test_deplete_nuclide.py | 39 +++++++++++++++++------- 2 files changed, 55 insertions(+), 23 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index a989bf7029..2ec1d3a66d 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -3,9 +3,11 @@ Contains the per-nuclide components of a depletion chain. """ +import bisect from collections.abc import Mapping from collections import namedtuple, defaultdict from warnings import warn +from numbers import Real try: import lxml.etree as ET except ImportError: @@ -339,6 +341,7 @@ class Nuclide(object): return valid + class FissionYieldDistribution(Mapping): """Energy-dependent fission product yields for a single nuclide @@ -504,10 +507,15 @@ class FissionYield(Mapping): self.products = products self.yields = yields + def __contains__(self, product): + ix = bisect.bisect_left(self.products, product) + return ix != len(self.products) and self.products[ix] == product + def __getitem__(self, product): - if product not in self.products: + ix = bisect.bisect_left(self.products, product) + if ix == len(self.products) or self.products[ix] != product: raise KeyError(product) - return self.yields[self.products.index(product)] + return self.yields[ix] def __len__(self): return len(self.products) @@ -516,35 +524,42 @@ class FissionYield(Mapping): return iter(self.products) def items(self): + """Return pairs of product, yield""" return zip(self.products, self.yields) def __add__(self, other): - new = self.copy() + if not isinstance(other, FissionYield): + return NotImplemented + new = FissionYield(self.products, self.yields.copy()) new += other return new def __iadd__(self, other): """Increment value from other fission yield""" + if not isinstance(other, FissionYield): + return NotImplemented self.yields += other.yields return self + def __radd__(self, other): + return self + other + def __imul__(self, scalar): + if not isinstance(scalar, Real): + return NotImplemented self.yields *= scalar return self def __mul__(self, scalar): - new = self.copy() + if not isinstance(scalar, Real): + return NotImplemented + new = FissionYield(self.products, self.yields.copy()) new *= scalar return new def __rmul__(self, scalar): - new = self.copy() - new *= scalar - return new + return self * scalar def __repr__(self): - return "<{}: {}>".format(self.__class__.__name__, repr(dict(self))) - - def copy(self): - """Return an identical yield object, with unique yields""" - return FissionYield(self.products, self.yields.copy()) + return "<{} containing {} products and yields>".format( + self.__class__.__name__, len(self)) diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 6f5356b129..7d68b3a3d8 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -134,25 +134,42 @@ def test_fission_yield_distribution(): act_dist = yield_dict[exp_ene] for exp_prod, exp_yield in exp_dist.items(): assert act_dist[exp_prod] == exp_yield - exp_yield_matrix = numpy.array([ + exp_yield = numpy.array([ [4.08e-12, 1.71e-12, 7.85e-4], [1.32e-12, 0.0, 1.12e-3], [5.83e-8, 2.69e-8, 4.54e-3]]) - assert numpy.array_equal(yield_dist.yield_matrix, exp_yield_matrix) + assert numpy.array_equal(yield_dist.yield_matrix, exp_yield) # Test the operations / special methods for fission yield - orig_yield_obj = yield_dist[0.0253] + orig_yields = yield_dist[0.0253] + assert len(orig_yields) == len(yield_dict[0.0253]) + for key, value in yield_dict[0.0253].items(): + assert key in orig_yields + assert orig_yields[key] == value # __getitem__ return yields as a view into yield matrix - assert orig_yield_obj.yields.base is yield_dist.yield_matrix - copied_yield = orig_yield_obj.copy() - # copied yields own their own memory -> not a view - assert copied_yield.yields.base is None + assert orig_yields.yields.base is yield_dist.yield_matrix # Fission yield feature uses scaled and incremented - mod_yields = orig_yield_obj * 2 - assert numpy.array_equal(orig_yield_obj.yields * 2, mod_yields.yields) - mod_yields += orig_yield_obj - assert numpy.array_equal(orig_yield_obj.yields * 3, mod_yields.yields) + mod_yields = orig_yields * 2 + assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields) + mod_yields += orig_yields + assert numpy.array_equal(orig_yields.yields * 3, mod_yields.yields) + + # Failure modes for adding, multiplying yields + similar = numpy.empty_like(orig_yields.yields) + with pytest.raises(TypeError): + orig_yields + similar + with pytest.raises(TypeError): + similar + orig_yields + with pytest.raises(TypeError): + orig_yields += similar + with pytest.raises(TypeError): + orig_yields * similar + with pytest.raises(TypeError): + similar * orig_yields + with pytest.raises(TypeError): + orig_yields *= similar + def test_validate(): From 61232248dffef6ec14a1f23387949af6cfe452f5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 15:23:45 -0500 Subject: [PATCH 41/45] Update documentation for FissionYieldDistribution --- openmc/deplete/nuclide.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 2ec1d3a66d..8896180a88 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -355,16 +355,12 @@ class FissionYieldDistribution(Mapping): Parameters ---------- - energies : iterable of float - Energies for which fission yield data exist. Must be ordered - by increasing energy - products : iterable of str - Fission products produced by this parent at all energies. - Must be ordered alphabetically - yield_matrix : numpy.ndarray or iterable of iterable of float - Array of shape ``(n_energy, n_products)`` where - ``yield_matrix[g][j]`` is the yield of - ``ordered_products[j]`` due to a fission in energy region ``g``. + fission_yields : dict + Dictionary of energies and fission product yields for that energy. + Expected to be of the form ``{float: {str: float}}``. The first + float is the energy, typically in eV, that represents this + distribution. The underlying dictionary maps fission products + to their respective yields. Attributes ---------- @@ -416,6 +412,11 @@ class FissionYieldDistribution(Mapping): def __iter__(self): return iter(self.energies) + def __repr__(self): + return "<{} with {} products at {} energies>".format( + self.__class__.__name__, self.yield_matrix.shape[1], + len(self.energies)) + @classmethod def from_xml_element(cls, element): """Construct a distribution from a depletion chain xml file From 1e5cfe09613cf8260840d56900408702e662eba8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 15:32:45 -0500 Subject: [PATCH 42/45] Return list of str from FY helpers update_tally_nuclides --- openmc/deplete/abc.py | 8 ++++---- tests/unit_tests/test_deplete_fission_yields.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 16d2d7abbd..fb638e5c21 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -458,13 +458,13 @@ class FissionYieldHelper(ABC): Returns ------- - nuclides : tuple of str + nuclides : list of str Union of nuclides that the :class:`openmc.deplete.Operator` says have non-zero densities at this stage and those that have yield data. Sorted by nuclide name """ - return tuple(sorted(self._chain_set & set(nuclides))) + return sorted(self._chain_set & set(nuclides)) @classmethod def from_operator(cls, operator, **kwargs): @@ -550,7 +550,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): Returns ------- - nuclides : tuple of str + nuclides : list of str Union of input nuclides and those that have multiple sets of yield data. Sorted by nuclide name @@ -562,7 +562,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): assert self._fission_rate_tally is not None, ( "Run generate_tallies first") overlap = set(self._chain_nuclides).intersection(set(nuclides)) - nuclides = tuple(sorted(overlap)) + nuclides = sorted(overlap) self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] self._fission_rate_tally.nuclides = nuclides return nuclides diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 2842043d12..d3faec35bf 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -206,7 +206,7 @@ def test_cutoff_helper(materials, nuclide_bundle, therm_frac): non_zero_nucs = [n.name for n in nuclide_bundle] tally_nucs = helper.update_tally_nuclides(non_zero_nucs) - assert tally_nucs == ("Pu239", "U235",) + assert tally_nucs == ["Pu239", "U235",] # Check tallies fission_tally = helper._fission_rate_tally @@ -249,7 +249,7 @@ def test_averaged_helper(materials, nuclide_bundle, avg_energy): helper.generate_tallies(materials, [0]) tallied_nucs = helper.update_tally_nuclides( [n.name for n in nuclide_bundle]) - assert tallied_nucs == ("Pu239", "U235") + assert tallied_nucs == ["Pu239", "U235"] # check generated tallies fission_tally = helper._fission_rate_tally From 991a36df8dce07c15a0a94f2e813ea9f309ea9be Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 22 Aug 2019 10:13:28 -0500 Subject: [PATCH 43/45] Pass strings to os functions in test_deplete_fission_yields os.remove and os.chdir support path-like objects only for python 3.6+ --- tests/unit_tests/test_deplete_fission_yields.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index d3faec35bf..3a63d718b9 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -63,13 +63,12 @@ def materials(tmpdir_factory): with capi.run_in_memory(): yield [capi.Material(), capi.Material()] finally: - print(os.path.abspath(os.curdir)) - os.remove(tmpdir / "settings.xml") - os.remove(tmpdir / "geometry.xml") - os.remove(tmpdir / "materials.xml") - os.remove(tmpdir / "summary.h5") + # Convert to strings as os.remove in py 3.5 doesn't support Paths + for file_path in ("settings.xml", "geometry.xml", "materials.xml", + "summary.h5"): + os.remove(str(tmpdir / file_path)) orig.chdir() - os.rmdir(tmpdir) + os.rmdir(str(tmpdir)) def proxy_tally_data(tally, fill=None): From 399b77e38fea9f4488049c7c918038e04cfc3971 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 30 Aug 2019 09:36:28 -0500 Subject: [PATCH 44/45] Clean up internals for AveragedFissionYieldHelper Respond to some other reviewer comments for #1313 --- openmc/deplete/chain.py | 1 + openmc/deplete/helpers.py | 14 ++++++-------- openmc/deplete/operator.py | 2 +- tests/unit_tests/test_deplete_fission_yields.py | 2 +- 4 files changed, 9 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 9f6b9ecf7f..c2b5b75b31 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -731,6 +731,7 @@ class Chain(object): yields = [yields] check_type("fission_yields", yields, Iterable, Mapping) self._fission_yields = yields + def validate(self, strict=True, quiet=False, tolerance=1e-4): """Search for possible inconsistencies diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index eb86b0ea4a..56644a3bba 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -425,11 +425,9 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): return yields rates = self.results[local_mat_index] # iterate over thermal then fast yields, prefer __mul__ to __rmul__ - for therm_frac, nuc in zip(rates[0], self._tally_nucs): - yields[nuc.name] = self._thermal_yields[nuc.name] * therm_frac - - for fast_frac, nuc in zip(rates[1], self._tally_nucs): - yields[nuc.name] += self._fast_yields[nuc.name] * fast_frac + for therm_frac, fast_frac, nuc in zip(rates[0], rates[1], self._tally_nucs): + yields[nuc.name] = (self._thermal_yields[nuc.name] * therm_frac + + self._fast_yields[nuc.name] * fast_frac) return yields @property @@ -444,9 +442,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): class AveragedFissionYieldHelper(TalliedFissionYieldHelper): r"""Class that computes fission yields based on average fission energy - Computes average energy at which fission events occured - reactions for all nuclides with multiple sets of fission yields - by + Computes average energy at which fission events occured with .. math:: @@ -481,6 +477,8 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): have shape ``(n_mats, n_tnucs)``, where ``n_mats`` is the number of materials where fission reactions were tallied and ``n_tnucs`` is the number of nuclides with multiple sets of fission yields. + Data in the array are the average energy of fission events for + tallied nuclides across burnable materials. """ def __init__(self, chain_nuclides): diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 43523b0fe1..01bf92dcc9 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -86,7 +86,7 @@ class Operator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. - fission_yield_mode : ("constant", "cutoff", "average") + fission_yield_mode : {"constant", "cutoff", "average"} Key indicating what fission product yield scheme to use. The key determines what fission energy helper is used: diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 3a63d718b9..62d2f65866 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -205,7 +205,7 @@ def test_cutoff_helper(materials, nuclide_bundle, therm_frac): non_zero_nucs = [n.name for n in nuclide_bundle] tally_nucs = helper.update_tally_nuclides(non_zero_nucs) - assert tally_nucs == ["Pu239", "U235",] + assert tally_nucs == ["Pu239", "U235"] # Check tallies fission_tally = helper._fission_rate_tally From c486b4187ffbf6ae3c36e2ab9a7e3489c0cbdf70 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 30 Aug 2019 09:55:58 -0500 Subject: [PATCH 45/45] Use python API for setup in test_deplete_fission_yields --- .../unit_tests/test_deplete_fission_yields.py | 68 +++++++------------ 1 file changed, 24 insertions(+), 44 deletions(-) diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 62d2f65866..bedd702a90 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -6,7 +6,8 @@ from unittest.mock import Mock import bisect import pytest -import numpy +import numpy as np +import openmc from openmc import capi from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( @@ -19,46 +20,25 @@ def materials(tmpdir_factory): """Use C API to construct realistic materials for testing tallies""" tmpdir = tmpdir_factory.mktemp("capi") orig = tmpdir.chdir() - # Create proxy xml files to please openmc - with open("geometry.xml", "w") as stream: - stream.write(""" - - - - - - - - -""") - with open("settings.xml", "w") as stream: - stream.write(""" - - - eigenvalue - 100 - 10 - 0 - 1 - - - 0.0 0.0 0.0 - - - -""") - with open("materials.xml", "w") as stream: - stream.write(""" - - - - - - - - - -""") + # Create proxy problem to please openmc + mfuel = openmc.Material(name="test_fuel") + mfuel.volume = 1.0 + for nuclide in ["U235", "U238", "Xe135", "Pu239"]: + mfuel.add_nuclide(nuclide, 1.0) + openmc.Materials([mfuel]).export_to_xml() + # Geometry + box = openmc.rectangular_prism(1.0, 1.0, boundary_type="reflective") + cell = openmc.Cell(fill=mfuel, region=box) + root = openmc.Universe(cells=[cell]) + openmc.Geometry(root).export_to_xml() + # settings + settings = openmc.Settings() + settings.particles = 100 + settings.inactive = 0 + settings.batches = 10 + settings.verbosity = 1 + settings.export_to_xml() + try: with capi.run_in_memory(): yield [capi.Material(), capi.Material()] @@ -87,7 +67,7 @@ def proxy_tally_data(tally, fill=None): if isinstance(tfilter, capi.EnergyFilter): this_bins -= 1 n_bins *= max(this_bins, 1) - data = numpy.empty((n_bins, n_nucs * n_scores, 3)) + data = np.empty((n_bins, n_nucs * n_scores, 3)) if fill is not None: data.fill(fill) return data @@ -229,7 +209,7 @@ def test_cutoff_helper(materials, nuclide_bundle, therm_frac): helper.unpack() # expected results of shape (n_mats, 2, n_tnucs) - expected_results = numpy.empty((1, 2, len(tally_nucs))) + expected_results = np.empty((1, 2, len(tally_nucs))) expected_results[:, 0] = therm_frac expected_results[:, 1] = 1 - therm_frac assert helper.results == pytest.approx(expected_results) @@ -285,7 +265,7 @@ def test_averaged_helper(materials, nuclide_bundle, avg_energy): helper._weighted_tally.results = weighted_results helper.unpack() - expected_results = numpy.ones((1, len(tallied_nucs))) * avg_energy + expected_results = np.ones((1, len(tallied_nucs))) * avg_energy assert helper.results == pytest.approx(expected_results) actual_yields = helper.weighted_yields(0)