From d1878a205e8020bc72d06c0f5f64503fdd808929 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 20 Mar 2020 18:07:17 -0400 Subject: [PATCH 01/24] Replace "Nothing" depletion targets with None This change is more declarative and mitigates any changes that there are errors if a target is Nothing, None, or their lower cased varieties. This is helpful as the depletion chain gets reduced down and we want to keep total removal rates for isotopes that are present, but their products don't exist in the new reduced chain. Chain.form_matrix has been updated to handle targets that are None for non-fission reactions and for decay reactions. Documentation for ReactionTuple and DecayTuple have been updated to indicate that the target attribute could be None. --- openmc/deplete/chain.py | 23 +++++++++++------------ openmc/deplete/nuclide.py | 22 ++++++++++++++++------ 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 814de12fc4..8c31914acb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -500,7 +500,7 @@ class Chain(object): # Gain for _, target, branching_ratio in nuc.decay_modes: # Allow for total annihilation for debug purposes - if target != 'Nothing': + if target is not None: branch_val = branching_ratio * decay_constant if branch_val != 0.0: @@ -525,17 +525,16 @@ class Chain(object): matrix[i, i] -= path_rate # Gain term; allow for total annihilation for debug purposes - if target != 'Nothing': - if r_type != 'fission': - if path_rate != 0.0: - k = self.nuclide_dict[target] - matrix[k, i] += path_rate * br - else: - for product, y in fission_yields[nuc.name].items(): - yield_val = y * path_rate - if yield_val != 0.0: - k = self.nuclide_dict[product] - matrix[k, i] += yield_val + if r_type != 'fission': + if target is not None and path_rate != 0.0: + k = self.nuclide_dict[target] + matrix[k, i] += path_rate * br + else: + for product, y in fission_yields[nuc.name].items(): + yield_val = y * path_rate + if yield_val != 0.0: + k = self.nuclide_dict[product] + matrix[k, i] += yield_val # Clear set of reactions reactions.clear() diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 0aae5adb72..7235745aa2 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -30,8 +30,10 @@ Parameters ---------- type : str Type of the decay mode, e.g., 'beta-' -target : str - Nuclide resulting from decay +target : str or None + Nuclide resulting from decay. A value of ``None`` implies the + target does not exist in the currently configured depletion + chain branching_ratio : float Branching ratio of the decay mode @@ -53,8 +55,11 @@ Parameters ---------- type : str Type of the reaction, e.g., 'fission' -target : str - nuclide resulting from reaction +target : str or None + Nuclide resulting from reaction. A value of ``None`` + implies either no single target, e.g. from fission, + or that the target nuclide is not considered + in the current depletion chain Q : float Q value of the reaction in [eV] branching_ratio : float @@ -179,6 +184,8 @@ class Nuclide(object): for decay_elem in element.iter('decay'): d_type = decay_elem.get('type') target = decay_elem.get('target') + if target is not None and target.lower() == "nothing": + target = None branching_ratio = float(decay_elem.get('branching_ratio')) nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio)) @@ -192,6 +199,8 @@ class Nuclide(object): # just set null values if r_type != 'fission': target = reaction_elem.get('target') + if target is not None and target.lower() == "nothing": + target = None else: target = None if fission_q is not None: @@ -226,7 +235,7 @@ class Nuclide(object): for mode, daughter, br in self.decay_modes: mode_elem = ET.SubElement(elem, 'decay') mode_elem.set('type', mode) - mode_elem.set('target', daughter) + mode_elem.set('target', daughter or "Nothing") mode_elem.set('branching_ratio', str(br)) elem.set('reactions', str(len(self.reactions))) @@ -234,7 +243,7 @@ class Nuclide(object): rx_elem = ET.SubElement(elem, 'reaction') rx_elem.set('type', rx) rx_elem.set('Q', str(Q)) - if rx != 'fission': + if rx != 'fission' or daughter is not None: rx_elem.set('target', daughter) if br != 1.0: rx_elem.set('branching_ratio', str(br)) @@ -461,6 +470,7 @@ class FissionYieldDistribution(Mapping): data_elem.text = " ".join(map(str, yield_obj.yields)) + class FissionYield(Mapping): """Mapping for fission yields of a parent at a specific energy From 58b644cd3b7dfc2872895b73f7795a9e9d96561b Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 18 Mar 2020 20:08:03 -0400 Subject: [PATCH 02/24] Minor fission yield touch ups --- openmc/deplete/nuclide.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 7235745aa2..960ba6eff5 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -406,9 +406,7 @@ class FissionYieldDistribution(Mapping): for g_index, energy in enumerate(energies): prod_map = fission_yields[energy] for prod_ix, product in enumerate(ordered_prod): - yield_val = prod_map.get(product) - yield_matrix[g_index, prod_ix] = ( - 0.0 if yield_val is None else yield_val) + yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0) self.energies = tuple(energies) self.products = tuple(ordered_prod) self.yield_matrix = yield_matrix @@ -444,7 +442,7 @@ class FissionYieldDistribution(Mapping): FissionYieldDistribution """ all_yields = {} - for elem_index, yield_elem in enumerate(element.iter("fission_yields")): + for yield_elem in element.iter("fission_yields"): energy = float(yield_elem.get("energy")) products = yield_elem.find("products").text.split() yields = map(float, yield_elem.find("data").text.split()) From 2c2d9d7f312e8b7dfcea10de2ca9aee8e8361f3c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 18 Mar 2020 20:09:19 -0400 Subject: [PATCH 03/24] Provide FissionYieldDistribution.restrict_products --- openmc/deplete/nuclide.py | 32 +++++++++++++++++++++++- tests/unit_tests/test_deplete_nuclide.py | 15 +++++++++++ 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 960ba6eff5..7ba72b6d56 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -13,7 +13,7 @@ try: except ImportError: import xml.etree.ElementTree as ET -from numpy import empty +from numpy import empty, searchsorted from openmc.checkvalue import check_type @@ -467,6 +467,36 @@ class FissionYieldDistribution(Mapping): data_elem = ET.SubElement(yield_element, "data") data_elem.text = " ".join(map(str, yield_obj.yields)) + def restrict_products(self, possible_products): + """Return a new distribution with select products + + Parameters + ---------- + possible_products : iterable of str + Candidate pool of fission products. Existing products + not contained here will not exist in the new instance + + Returns + ------- + FissionYieldDistribution or None + A value of None indicates no values in + ``possible_products`` exist in :attr:`products` + + """ + + overlap = set(self.products).intersection(possible_products) + if not overlap: + return None + + products = sorted(overlap) + indices = searchsorted(self.products, products) + + # coerce back to dictionary to pass back to __init__ + new_yields = {} + for ene, yields in zip(self.energies, self.yield_matrix.copy()): + new_yields[ene] = dict(zip(products, yields[indices])) + + return type(self)(new_yields) class FissionYield(Mapping): diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index b8a258df00..92b004845a 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -176,6 +176,21 @@ def test_fission_yield_distribution(): with pytest.raises(TypeError): orig_yields *= similar + # Test restriction of fission products + strict_restrict = yield_dist.restrict_products(["Xe135", "Sm149"]) + with_extras = yield_dist.restrict_products( + ["Xe135", "Sm149", "H1", "U235"]) + + assert strict_restrict.products == ("Sm149", "Xe135", ) + assert strict_restrict.energies == yield_dist.energies + assert with_extras.products == ("Sm149", "Xe135", ) + assert with_extras.energies == yield_dist.energies + for ene, new_yields in strict_restrict.items(): + for product in strict_restrict.products: + assert new_yields[product] == yield_dist[ene][product] + assert with_extras[ene][product] == yield_dist[ene][product] + + assert yield_dist.restrict_products(["U235"]) is None def test_validate(): From 789926d107ce447886034a4e6f01bed0587a8091 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 18 Mar 2020 20:10:49 -0400 Subject: [PATCH 04/24] Provide Chain.reduce for following paths --- openmc/deplete/chain.py | 142 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 141 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 8c31914acb..d820fe520b 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -10,7 +10,7 @@ import math import re from collections import OrderedDict, defaultdict from collections.abc import Mapping, Iterable -from numbers import Real +from numbers import Real, Integral from warnings import warn from openmc.checkvalue import check_type, check_greater_than @@ -820,3 +820,143 @@ class Chain(object): return stat valid = valid and stat return valid + + def reduce(self, initial_isotopes, level=None): + """Reduce the size of the chain by follwing transmutation paths + + Parameters + ---------- + initial_isotopes : iterable of str + Start the search based on the contents of these isotopes + level : int, optional + Depth of transmuation path to follow. Must be greater than + or equal to zero. A value of zero returns a chain with + ``initial_isotopes``. The default value of None implies + that all isotopes that appear in the transmutation paths + of the initial isotopes and their progeny should be + explored + + Returns + ------- + Chain + + """ + check_type("initial_isotopes", initial_isotopes, Iterable, str) + if level is None: + level = math.inf + else: + check_type("level", level, Integral) + check_greater_than("level", level, 0, equality=True) + + all_isotopes = self._follow(set(initial_isotopes), level) + + # Avoid re-sorting for fission yields + name_sort = sorted(all_isotopes) + + nuclides = [] + nuclide_dict = {} + reactions = set() + + for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)): + previous = self[iso] + new_nuclide = Nuclide(previous.name) + new_nuclide.half_life = previous.half_life + new_nuclide.decay_energy = new_nuclide.decay_energy + + new_decay = [] + for mode in previous.decay_modes: + if mode.target in all_isotopes: + new_decay.append(mode) + else: + new_decay.append(DecayTuple( + mode.type, None, mode.branching_ratio)) + new_nuclide.decay_modes = new_decay + + new_reactions = [] + for rxn in previous.reactions: + if rxn.target in all_isotopes: + new_reactions.append(rxn) + reactions.add(rxn.type) + elif rxn.type == "fission": + new_yields = new_nuclide.yield_data = ( + previous.yield_data.restrict_products(name_sort)) + if new_yields is not None: + new_reactions.append(rxn) + reactions.add("fission") + # Maintain total destruction rates but set no target + else: + new_reactions.append(ReactionTuple( + rxn.type, None, rxn.Q, rxn.branching_ratio)) + reactions.add(rxn.type) + + new_nuclide.reactions = new_reactions + + nuclides.append(new_nuclide) + nuclide_dict[iso] = idx + + new_chain = type(self)() + new_chain.nuclides = nuclides + new_chain.nuclide_dict = nuclide_dict + + # Doesn't appear that the ordering matters for the reactions, + # just the contents + new_chain.reactions = sorted(reactions) + + return new_chain + + def _follow(self, isotopes, level): + """Return all isotopes present up to depth level""" + found = set(isotopes) + remaining = set(self.nuclide_dict) + if not found.issubset(remaining): + raise IndexError( + "The following isotopes were not found in the chain: " + "{}".format(", ".join(found - remaining))) + + if level == 0: + return found + + remaining.difference_update(found) + + depth = 0 + next_iso = set() + + while depth < level and remaining: + # Exhaust all isotopes at this level + while isotopes: + iso = isotopes.pop() + found.add(iso) + nuclide = self[iso] + + # Follow all transmutation paths for this nuclide + for rxn in nuclide.reactions + nuclide.decay_modes: + if rxn.type == "fission" or rxn.target is None: + continue + # Skip if we've already come across this isotope + elif (rxn.target in next_iso + or rxn.target in found or rxn.target in isotopes): + continue + next_iso.add(rxn.target) + + if nuclide.yield_data is not None: + for product in nuclide.yield_data.products: + if (product in next_iso + or product in found or product in isotopes): + continue + next_iso.add(product) + + if not next_iso: + # No additional isotopes to process, nor to update the + # current set of discovered isotopes + return found + + # Prepare for next dig + depth += 1 + isotopes.update(next_iso) + remaining.difference_update(next_iso) + next_iso.clear() + + # Process isotope that would have started next depth + found.update(isotopes) + + return found From b01a3e011be15f096f265fec7c19620861d08985 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 18 Mar 2020 20:11:47 -0400 Subject: [PATCH 05/24] Add unit tests for Chain.reduce --- tests/unit_tests/test_deplete_chain.py | 96 +++++++++++++++++++++++++- 1 file changed, 94 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 08c60d3ed1..b7a1e6c0f2 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -6,7 +6,6 @@ from pathlib import Path from itertools import product import numpy as np -from openmc.data import zam, ATOMIC_SYMBOL from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram import pytest @@ -405,7 +404,8 @@ def test_fission_yield_attribute(simple_chain): dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1) with pytest.raises( ValueError, match="fission yield.*not equal.*compositions"): - cram.deplete(empty_chain, dummy_conc, None, 0.5) + cram.deplete(empty_chain, dummy_conc, None, 0.5) + def test_validate(simple_chain): """Test the validate method""" @@ -457,3 +457,95 @@ def test_validate_inputs(): with pytest.raises(ValueError, match="tolerance"): c.validate(tolerance=-1) + + +@pytest.fixture +def gnd_simple_chain(): + chainfile = Path(__file__).parents[1] / "chain_simple.xml" + return Chain.from_xml(chainfile) + + +def test_reduce(gnd_simple_chain): + ref_U5 = gnd_simple_chain["U235"] + ref_iodine = gnd_simple_chain["I135"] + ref_U5_yields = ref_U5.yield_data + + no_depth = gnd_simple_chain.reduce(["U235", "I135"], 0) + # We should get a chain just containing U235 and I135 + assert len(no_depth) == 2 + assert set(no_depth.reactions) == set(gnd_simple_chain.reactions) + + u5_round0 = no_depth["U235"] + assert u5_round0.n_decay_modes == ref_U5.n_decay_modes + for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes): + assert newmode.target is None + assert newmode.type == refmode.type, newmode + assert newmode.branching_ratio == refmode.branching_ratio, newmode + + assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths + for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions): + assert newrxn.target is None, newrxn + assert newrxn.type == refrxn.type, newrxn + assert newrxn.Q == refrxn.Q, newrxn + assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn + + assert u5_round0.yield_data is not None + assert u5_round0.yield_data.products == ("I135", ) + assert u5_round0.yield_data.yield_matrix == ( + ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")] + ) + + bareI5 = no_depth["I135"] + assert bareI5.n_decay_modes == ref_iodine.n_decay_modes + for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes): + assert newmode.target is None + assert newmode.type == refmode.type, newmode + assert newmode.branching_ratio == refmode.branching_ratio, newmode + + assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths + for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions): + assert newrxn.target is None, newrxn + assert newrxn.type == refrxn.type, newrxn + assert newrxn.Q == refrxn.Q, newrxn + assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn + + follow_u5 = gnd_simple_chain.reduce(["U235"], 1) + u5_round1 = follow_u5["U235"] + assert u5_round1.decay_modes == ref_U5.decay_modes + assert u5_round1.reactions == ref_U5.reactions + assert u5_round1.yield_data is not None + assert ( + u5_round1.yield_data.yield_matrix == ref_U5_yields.yield_matrix + ).all() + + # Per the chain_simple.xml + # I135 -> Xe135 -> Cs135 + # I135 -> Xe136 + # No limit on depth + iodine_chain = gnd_simple_chain.reduce(["I135"]) + truncated_iodine = gnd_simple_chain.reduce(["I135"], 1) + assert len(iodine_chain) == 4 + assert len(truncated_iodine) == 3 + assert set(iodine_chain.nuclide_dict) == { + "I135", "Xe135", "Xe136", "Cs135"} + assert set(truncated_iodine.nuclide_dict) == {"I135", "Xe135", "Xe136"} + assert iodine_chain.reactions == ["(n,gamma)"] + assert iodine_chain["I135"].decay_modes == ref_iodine.decay_modes + assert iodine_chain["I135"].reactions == ref_iodine.reactions + for mode in truncated_iodine["Xe135"].decay_modes: + assert mode.target is None + + # Test that no FissionYieldDistribution is made if there are no + # fission products + u5_noyields = gnd_simple_chain.reduce(["U235"], 0)["U235"] + assert u5_noyields.yield_data is None + + # Check early termination if the eventual full chain + # is specified by using the iodine isotopes + new_iodine = gnd_simple_chain.reduce(set(iodine_chain.nuclide_dict)) + assert set(iodine_chain.nuclide_dict) == set(new_iodine.nuclide_dict) + + # Failure if some requested isotopes not in chain + + with pytest.raises(IndexError, match=".*not found.*Xx999"): + gnd_simple_chain.reduce(["U235", "Xx999"]) From d67d0cd0f9deb79d5627fa79bc985854d89f5d03 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 20 Mar 2020 18:11:54 -0400 Subject: [PATCH 06/24] Allow Operator to reduce depletion chain from geometry A new parameter, reduce_chain, can be used to instruct the operator to create a smaller depletion chain. The geometry is used to find all nuclides in burnable materials. The argument can be a boolean, or a non-negative integer. A value of True indicates no depth for following the paths, while an integer can be used to directly specify the depth. --- openmc/deplete/operator.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 8eb033783f..2572150012 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -113,6 +113,12 @@ class Operator(TransportOperator): ``fission_yield_mode``. Will be passed directly on to the helper. Passing a value of None will use the defaults for the associated helper. + reduce_chain : bool or int, optional + If ``True`` or an integer, create a reduced depletion chain by + following transmuation paths for isotopes in burnable materials. + A value of ``True`` implies to follow all paths to completion, + while a non-negative integer indicates the depth. See + :meth:`openmc.deplete.Chain.reduce` Attributes ---------- @@ -158,7 +164,8 @@ class Operator(TransportOperator): def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, energy_mode="fission-q", fission_q=None, dilute_initial=1.0e3, - fission_yield_mode="constant", fission_yield_opts=None): + fission_yield_mode="constant", fission_yield_opts=None, + reduce_chain=False): if fission_yield_mode not in self._fission_helpers: raise KeyError( "fission_yield_mode must be one of {}, not {}".format( @@ -179,6 +186,17 @@ class Operator(TransportOperator): self.geometry = geometry self.diff_burnable_mats = diff_burnable_mats + # Reduce the chain before we create more materials + if reduce_chain is not False: + all_isotopes = set() + for material in geometry.get_all_materials().values(): + if not material.depletable: + continue + for name, _dens_percent, _dens_type in material.nuclides: + all_isotopes.add(name) + level = None if reduce_chain is True else reduce_chain + self.chain = self.chain.reduce(all_isotopes, level) + # Differentiate burnable materials with multiple instances if self.diff_burnable_mats: self._differentiate_burnable_mats() From eab604263772de069db9b44548a8d54e8fd1d332 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 21 Mar 2020 10:07:40 -0400 Subject: [PATCH 07/24] Update TransportOperator.get_results_info docstring The volume entry returned by the canonical Operator is a dictionary of string material ids to volumes. This form is expected by the Result objects in transfer_volumes and when writing attributes to the hdf5 file. --- openmc/deplete/abc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 31479f9952..ddc48bc4de 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -190,7 +190,7 @@ class TransportOperator(ABC): Returns ------- - volume : list of float + volume : dict of str to float Volumes corresponding to materials in burn_list nuc_list : list of str A list of all nuclide names. Used for sorting the simulation. From ceaecef45a776d4bb95c8d96b3ed2ae3037ec114 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 21 Mar 2020 10:11:56 -0400 Subject: [PATCH 08/24] Update Result.volume docstring: dict of str -> float --- openmc/deplete/results.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index b44739bc5f..0b6d92d7c9 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -36,7 +36,7 @@ class Results(object): Number of nuclides. rates : list of ReactionRates The reaction rates for each substep. - volume : OrderedDict of int to float + volume : OrderedDict of str to float Dictionary mapping mat id to volume. mat_to_ind : OrderedDict of str to int A dictionary mapping mat ID as string to index. From a86c3481d682743cb6d78735a56b6cffd46385c0 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 21 Mar 2020 10:15:32 -0400 Subject: [PATCH 09/24] ResultList.get_atoms supports atom densities, days Two new optional arguments to ResultsList.get_atoms are introduced that control the units on the return values. Time can be returned in second or in days if time_units is "s" or "d", respectively. Concentrations can be returned in atoms/cm^3 or atoms/b/cm if nuc_units is one of those values. The default is to return time in seconds and number of atoms, retaining back compatibility. --- openmc/deplete/results_list.py | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 1b5b81ac4c..9e7cd9ffb9 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -2,7 +2,7 @@ import h5py import numpy as np from .results import Results, _VERSION_RESULTS -from openmc.checkvalue import check_filetype_version +from openmc.checkvalue import check_filetype_version, check_value __all__ = ["ResultsList"] @@ -40,7 +40,7 @@ class ResultsList(list): new.append(Results.from_hdf5(fh, i)) return new - def get_atoms(self, mat, nuc): + def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"): """Get number of nuclides over time from a single material .. note:: @@ -57,15 +57,24 @@ class ResultsList(list): Material name to evaluate nuc : str Nuclide name to evaluate + nuc_units : {"atoms", "atoms/b/cm", "atoms/cm^3"}, optional + Units for the returned concentration. Default is ``"atoms"`` + time_units : {"s", "d"}, optional + Units for the returned time array. Default is ``"s"`` to + return the value in seconds Returns ------- time : numpy.ndarray - Array of times in [s] + Array of times in units of ``time_units`` concentration : numpy.ndarray - Total number of atoms for specified nuclide + Concentration of specified nuclide in units of ``nuc_units`` """ + check_value("time_units", time_units, {"s", "d"}) + check_value("nuc_units", nuc_units, + {"atoms", "atoms/b/cm", "atoms/cm^3"}) + time = np.empty_like(self, dtype=float) concentration = np.empty_like(self, dtype=float) @@ -74,6 +83,17 @@ class ResultsList(list): time[i] = result.time[0] concentration[i] = result[0, mat, nuc] + # Unit conversions + if time_units == "d": + time /= (60 * 60 * 24) + + if nuc_units != "atoms": + # Divide by volume to get density + concentration /= self[0].volume[mat] + if nuc_units == "atoms/b/cm": + # 1 barn = 1e-24 cm^2 + concentration *= 1e-24 + return time, concentration def get_reaction_rate(self, mat, nuc, rx): From 4677714866be47e2557cb602faaae1df1f568a8a Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 21 Mar 2020 10:18:06 -0400 Subject: [PATCH 10/24] Add tests for new ResultsList.get_atoms units --- tests/unit_tests/test_deplete_resultslist.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index c3ceda5d56..f95bc31539 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -19,12 +19,24 @@ def test_get_atoms(res): """Tests evaluating single nuclide concentration.""" t, n = res.get_atoms("1", "Xe135") - t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - n_ref = [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14] + t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0]) + n_ref = np.array( + [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14]) np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(n, n_ref) + # Check alternate units + volume = res[0].volume["1"] + + t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atoms/cm^3", time_units="d") + + assert t_days == pytest.approx(t_ref / (60 * 60 * 24)) + assert n_cm3 == pytest.approx(n_ref / volume) + + _t, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atoms/b/cm") + assert n_bcm == pytest.approx(n_cm3 * 1e-24) + def test_get_reaction_rate(res): """Tests evaluating reaction rate.""" From a2691a32d5a61e379ecc061511e21bc9e83dd35f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 23 Mar 2020 13:43:36 -0400 Subject: [PATCH 11/24] Use atoms/b-cm, atoms/cm3 in ResultsList.get_atoms Co-Authored-By: Paul Romano --- openmc/deplete/results_list.py | 6 +++--- tests/unit_tests/test_deplete_resultslist.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 9e7cd9ffb9..eee214d8ff 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -57,7 +57,7 @@ class ResultsList(list): Material name to evaluate nuc : str Nuclide name to evaluate - nuc_units : {"atoms", "atoms/b/cm", "atoms/cm^3"}, optional + nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional Units for the returned concentration. Default is ``"atoms"`` time_units : {"s", "d"}, optional Units for the returned time array. Default is ``"s"`` to @@ -73,7 +73,7 @@ class ResultsList(list): """ check_value("time_units", time_units, {"s", "d"}) check_value("nuc_units", nuc_units, - {"atoms", "atoms/b/cm", "atoms/cm^3"}) + {"atoms", "atom/b-cm", "atom/cm3"}) time = np.empty_like(self, dtype=float) concentration = np.empty_like(self, dtype=float) @@ -90,7 +90,7 @@ class ResultsList(list): if nuc_units != "atoms": # Divide by volume to get density concentration /= self[0].volume[mat] - if nuc_units == "atoms/b/cm": + if nuc_units == "atom/b-cm": # 1 barn = 1e-24 cm^2 concentration *= 1e-24 diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index f95bc31539..d3609135b3 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -29,12 +29,12 @@ def test_get_atoms(res): # Check alternate units volume = res[0].volume["1"] - t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atoms/cm^3", time_units="d") + t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atom/cm3", time_units="d") assert t_days == pytest.approx(t_ref / (60 * 60 * 24)) assert n_cm3 == pytest.approx(n_ref / volume) - _t, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atoms/b/cm") + _t, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atom/b-cm") assert n_bcm == pytest.approx(n_cm3 * 1e-24) From 2394468a49ef23ce660fcc865ed8a3445fed84e3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 23 Mar 2020 13:56:35 -0400 Subject: [PATCH 12/24] Support time in h, min from ResultsList.get_atoms --- openmc/deplete/results_list.py | 10 +++++++--- tests/unit_tests/test_deplete_resultslist.py | 6 +++++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index eee214d8ff..db1f0ff9db 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -59,9 +59,9 @@ class ResultsList(list): Nuclide name to evaluate nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional Units for the returned concentration. Default is ``"atoms"`` - time_units : {"s", "d"}, optional + time_units : {"s", "min", "h", "d"}, optional Units for the returned time array. Default is ``"s"`` to - return the value in seconds + return the value in seconds. Returns ------- @@ -71,7 +71,7 @@ class ResultsList(list): Concentration of specified nuclide in units of ``nuc_units`` """ - check_value("time_units", time_units, {"s", "d"}) + check_value("time_units", time_units, {"s", "d", "min", "h"}) check_value("nuc_units", nuc_units, {"atoms", "atom/b-cm", "atom/cm3"}) @@ -86,6 +86,10 @@ class ResultsList(list): # Unit conversions if time_units == "d": time /= (60 * 60 * 24) + elif time_units == "h": + time /= (60 * 60) + elif time_units == "min": + time /= 60 if nuc_units != "atoms": # Divide by volume to get density diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index d3609135b3..9b40a1ac91 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -34,8 +34,12 @@ def test_get_atoms(res): assert t_days == pytest.approx(t_ref / (60 * 60 * 24)) assert n_cm3 == pytest.approx(n_ref / volume) - _t, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atom/b-cm") + t_min, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atom/b-cm", time_units="min") assert n_bcm == pytest.approx(n_cm3 * 1e-24) + assert t_min == pytest.approx(t_ref / 60) + + t_hour, _n = res.get_atoms("1", "Xe135", time_units="h") + assert t_hour == pytest.approx(t_ref / (60 * 60)) def test_get_reaction_rate(res): From 8ba369d6cc44cbdc1fe2eaf12d043eafcd7512ce Mon Sep 17 00:00:00 2001 From: =shimwell Date: Mon, 23 Mar 2020 20:15:59 +0000 Subject: [PATCH 13/24] added elements from formula --- openmc/material.py | 96 ++++++++++++++++++++++++++++++- tests/unit_tests/test_material.py | 55 ++++++++++++++++++ 2 files changed, 150 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 3952fe4e99..8d5168d512 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1,9 +1,10 @@ -from collections import OrderedDict, defaultdict +from collections import OrderedDict, defaultdict, Counter from collections.abc import Iterable from copy import deepcopy from numbers import Real, Integral from pathlib import Path import warnings +import re from xml.etree import ElementTree as ET import numpy as np @@ -588,6 +589,99 @@ class Material(IDManagerMixin): enrichment_type): self.add_nuclide(*nuclide) + def add_elements_from_formula(self, formula, percent_type='ao', enrichment=None, + enrichment_target=None, enrichment_type=None): + """Add a elements from a chemical formula to the material. + + Parameters + ---------- + formula : str + Formula to add, e.g., 'C2O', 'C6H12O6', or (NH4)2SO4. + Note this is case sensitive, elements must start with an uppercase + character. Any numbers will be converted to integers (rounded down). + percent : float + Atom or weight percent + percent_type : {'ao', 'wo'}, optional + 'ao' for atom percent and 'wo' for weight percent. Defaults to atom + percent. + enrichment : float, optional + Enrichment of an enrichment_taget nuclide in percent (ao or wo). + If enrichment_taget is not supplied then it is enrichment for U235 + in weight percent. For example, input 4.95 for 4.95 weight percent + enriched U. + Default is None (natural composition). + enrichment_target: str, optional + Single nuclide name to enrich from a natural composition (e.g., 'O16') + enrichment_type: {'ao', 'wo'}, optional + 'ao' for enrichment as atom percent and 'wo' for weight percent. + Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment + + Notes + ----- + General enrichment procedure is allowed only for elements composed of + two isotopes. If `enrichment_target` is given without `enrichment` + natural composition is added to the material. + + """ + cv.check_type('formula', formula, str) + + # Tokenizes the formula and check validity of tokens + tokens = re.findall(r"([A-Z][a-z]*)(\d*)|(\()|(\))(\d*)", formula) + for row in tokens: + for token in row: + if token.isalpha(): + if token not in list(openmc.data.ATOMIC_NUMBER.keys())[1:]: + msg = 'Formula entry {} not an element symbol.' \ + .format(token) + raise ValueError(msg) + elif token not in ['(', ')', ''] and not token.isdigit(): + msg = 'Formula must be made from a sequence of ' \ + 'element symbols, integers, and backets. ' \ + '{} is not an allowable entry.'.format(token) + raise ValueError(msg) + + # Checks that the number of opening and closing brackets are equal + if formula.count('(') != formula.count(')'): + msg = 'Number of opening and closing brackets is not equal ' \ + 'in the input formula {}.'.format(formula) + raise ValueError(msg) + + # Checks that every part of the original formula has been tokenized + for row in tokens: + for token in row: + formula = formula.replace(token, '', 1) + if len(formula) != 0: + msg = 'Part of formula was not successfully parsed as an ' \ + 'element symbol, bracket or integer. {} was not parsed.' \ + .format(formula) + raise ValueError(msg) + + # Works through the tokens building a stack + mat_stack = [Counter()] + for symbol, multi1, opening_bracket, closing_bracket, multi2 in tokens: + if symbol: + mat_stack[-1][symbol] += int(multi1 or 1) + if opening_bracket: + mat_stack.append(Counter()) + if closing_bracket: + stack_top = mat_stack.pop() + for i in stack_top: + mat_stack[-1][i] += int(multi2 or 1) * stack_top[i] + + # Normalizing percentages + percents = mat_stack[0].values() + norm_percents = [float(i) / sum(percents) for i in percents] + elements = mat_stack[0].keys() + + # Adds each element and percent to the material + for element, percent in zip(elements, norm_percents): + if enrichment_target is not None and element == re.sub(r'\d+$', '', enrichment_target): + self.add_element(element, percent, percent_type, enrichment, + enrichment_target, enrichment_type) + else: + self.add_element(element, percent, percent_type) + + def add_s_alpha_beta(self, name, fraction=1.0): r"""Add an :math:`S(\alpha,\beta)` table to the material diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 6f51135ab2..7153fded5a 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -58,6 +58,61 @@ def test_elements_by_name(): assert a._nuclides == b._nuclides assert b._nuclides == c._nuclides +def test_adding_elements_by_formula(): + """Test adding elements from a formula""" + # testing the correct nuclides are added to the Material + m = openmc.Material() + m.add_elements_from_formula('Li4SiO4') + ref_dens = {'Li6': 0.033728, 'Li7': 0.410715, + 'Si28': 0.102477, 'Si29': 0.0052035, 'Si30': 0.0034301, + 'O16': 0.443386, 'O17': 0.000168} + nuc_dens = m.get_nuclide_atom_densities() + for nuclide in ref_dens: + assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + + # testing the correct nuclides are added to the Material when enriched + m = openmc.Material() + m.add_elements_from_formula('Li4SiO4', + enrichment=60., + enrichment_target='Li6') + ref_dens = {'Li6': 0.2666, 'Li7': 0.1777, + 'Si28': 0.102477, 'Si29': 0.0052035, 'Si30': 0.0034301, + 'O16': 0.443386, 'O17': 0.000168} + nuc_dens = m.get_nuclide_atom_densities() + for nuclide in ref_dens: + assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + + # testing the use of brackets + m = openmc.Material() + m.add_elements_from_formula('Mg2(NO3)2') + ref_dens = {'Mg24': 0.157902, 'Mg25': 0.02004, 'Mg26': 0.022058, + 'N14': 0.199267, 'N15': 0.000732, + 'O16': 0.599772, 'O17': 0.000227} + nuc_dens = m.get_nuclide_atom_densities() + for nuclide in ref_dens: + assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + + # testing lowercase elements results in a value error + m = openmc.Material() + with pytest.raises(ValueError): + m.add_elements_from_formula('li4SiO4') + + # testing lowercase elements results in a value error + m = openmc.Material() + with pytest.raises(ValueError): + m.add_elements_from_formula('Li4Sio4') + + # testing incorrect character in formula results in a value error + m = openmc.Material() + with pytest.raises(ValueError): + m.add_elements_from_formula('Li4$SiO4') + + # testing unequal opening and closing brackets + m = openmc.Material() + with pytest.raises(ValueError): + m.add_elements_from_formula('Fe(H2O)4(OH)2)') + + def test_density(): m = openmc.Material() for unit in ['g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3']: From 4e22c6a936c17f60d2c726818acf8836b76615d7 Mon Sep 17 00:00:00 2001 From: =shimwell Date: Mon, 23 Mar 2020 20:47:43 +0000 Subject: [PATCH 14/24] removed unused argument description --- openmc/material.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 8d5168d512..bdcf6adf17 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -599,8 +599,6 @@ class Material(IDManagerMixin): Formula to add, e.g., 'C2O', 'C6H12O6', or (NH4)2SO4. Note this is case sensitive, elements must start with an uppercase character. Any numbers will be converted to integers (rounded down). - percent : float - Atom or weight percent percent_type : {'ao', 'wo'}, optional 'ao' for atom percent and 'wo' for weight percent. Defaults to atom percent. From 605b3eade124c5a87cbd1d5e6e7abbdc75be0597 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 24 Mar 2020 11:34:20 +0000 Subject: [PATCH 15/24] code review imporements Co-Authored-By: Andrew Johnson --- openmc/material.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index bdcf6adf17..092389cbce 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -603,11 +603,10 @@ class Material(IDManagerMixin): 'ao' for atom percent and 'wo' for weight percent. Defaults to atom percent. enrichment : float, optional - Enrichment of an enrichment_taget nuclide in percent (ao or wo). - If enrichment_taget is not supplied then it is enrichment for U235 + Enrichment of an enrichment_target nuclide in percent (ao or wo). + If enrichment_target is not supplied then it is enrichment for U235 in weight percent. For example, input 4.95 for 4.95 weight percent - enriched U. - Default is None (natural composition). + enriched U. Default is None (natural composition). enrichment_target: str, optional Single nuclide name to enrich from a natural composition (e.g., 'O16') enrichment_type: {'ao', 'wo'}, optional @@ -628,7 +627,7 @@ class Material(IDManagerMixin): for row in tokens: for token in row: if token.isalpha(): - if token not in list(openmc.data.ATOMIC_NUMBER.keys())[1:]: + if token == "n" or token not in openmc.data.ATOMIC_NUMBER: msg = 'Formula entry {} not an element symbol.' \ .format(token) raise ValueError(msg) @@ -663,8 +662,8 @@ class Material(IDManagerMixin): mat_stack.append(Counter()) if closing_bracket: stack_top = mat_stack.pop() - for i in stack_top: - mat_stack[-1][i] += int(multi2 or 1) * stack_top[i] + for symbol, value in stack_top.items(): + mat_stack[-1][symbol] += int(multi2 or 1) * value # Normalizing percentages percents = mat_stack[0].values() @@ -678,8 +677,6 @@ class Material(IDManagerMixin): enrichment_target, enrichment_type) else: self.add_element(element, percent, percent_type) - - def add_s_alpha_beta(self, name, fraction=1.0): r"""Add an :math:`S(\alpha,\beta)` table to the material From 5fe6576a9cba1679c491344efa250133eafec5a5 Mon Sep 17 00:00:00 2001 From: =shimwell Date: Tue, 24 Mar 2020 11:57:17 +0000 Subject: [PATCH 16/24] added test to check decimal values fail --- tests/unit_tests/test_material.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 7153fded5a..43aeb62789 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -92,6 +92,11 @@ def test_adding_elements_by_formula(): for nuclide in ref_dens: assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + # testing non integer multiplier results in a value error + m = openmc.Material() + with pytest.raises(ValueError): + m.add_elements_from_formula('Li4.2SiO4') + # testing lowercase elements results in a value error m = openmc.Material() with pytest.raises(ValueError): From 430b42609896f4872f50b9ea1274114a67eec407 Mon Sep 17 00:00:00 2001 From: =shimwell Date: Tue, 24 Mar 2020 12:06:56 +0000 Subject: [PATCH 17/24] added check for decimal points in formula --- openmc/material.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 092389cbce..7cc7e3db60 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -598,7 +598,7 @@ class Material(IDManagerMixin): formula : str Formula to add, e.g., 'C2O', 'C6H12O6', or (NH4)2SO4. Note this is case sensitive, elements must start with an uppercase - character. Any numbers will be converted to integers (rounded down). + character. Multiplier numbers must be integers. percent_type : {'ao', 'wo'}, optional 'ao' for atom percent and 'wo' for weight percent. Defaults to atom percent. @@ -622,6 +622,11 @@ class Material(IDManagerMixin): """ cv.check_type('formula', formula, str) + if '.' in formula: + msg = 'Non-integer multiplier values are not accepted. The ' \ + 'input formula {} contains a "." character.'.format(formula) + raise ValueError(msg) + # Tokenizes the formula and check validity of tokens tokens = re.findall(r"([A-Z][a-z]*)(\d*)|(\()|(\))(\d*)", formula) for row in tokens: From ffa0cfc03e953fdadc2e4d9fd61319a3b18c1555 Mon Sep 17 00:00:00 2001 From: =shimwell Date: Tue, 24 Mar 2020 13:33:42 +0000 Subject: [PATCH 18/24] added tests for number of elements --- tests/unit_tests/test_material.py | 34 +++++++++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 43aeb62789..8e9ebb362d 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -3,7 +3,7 @@ import openmc.model import openmc.stats import openmc.examples import pytest - +from collections import defaultdict def test_attributes(uo2): assert uo2.name == 'UO2' @@ -60,9 +60,23 @@ def test_elements_by_name(): def test_adding_elements_by_formula(): """Test adding elements from a formula""" - # testing the correct nuclides are added to the Material + # testing the correct nuclides and elements are added to a material m = openmc.Material() m.add_elements_from_formula('Li4SiO4') + # checking the ratio of elements is 4:1:4 for Li:Si:O + elem = defaultdict(float) + for nuclide, adens in m.get_nuclide_atom_densities().values(): + if nuclide.startswith("Li"): + elem["Li"] += adens + if nuclide.startswith("Si"): + elem["Si"] += adens + if nuclide.startswith("O"): + elem["O"] += adens + total_number_of_atoms = 9 + assert elem["Li"] == pytest.approx(4./total_number_of_atoms) + assert elem["Si"] == pytest.approx(1./total_number_of_atoms) + assert elem["O"] == pytest.approx(4/total_number_of_atoms) + # testing the correct nuclides are added to the Material ref_dens = {'Li6': 0.033728, 'Li7': 0.410715, 'Si28': 0.102477, 'Si29': 0.0052035, 'Si30': 0.0034301, 'O16': 0.443386, 'O17': 0.000168} @@ -85,6 +99,22 @@ def test_adding_elements_by_formula(): # testing the use of brackets m = openmc.Material() m.add_elements_from_formula('Mg2(NO3)2') + + # checking the ratio of elements is 2:2:6 for Mg:N:O + elem = defaultdict(float) + for nuclide, adens in m.get_nuclide_atom_densities().values(): + if nuclide.startswith("Mg"): + elem["Mg"] += adens + if nuclide.startswith("N"): + elem["N"] += adens + if nuclide.startswith("O"): + elem["O"] += adens + total_number_of_atoms = 10 + assert elem["Mg"] == pytest.approx(2./total_number_of_atoms) + assert elem["N"] == pytest.approx(2./total_number_of_atoms) + assert elem["O"] == pytest.approx(6/total_number_of_atoms) + + # testing the correct nuclides are added when brackets are used ref_dens = {'Mg24': 0.157902, 'Mg25': 0.02004, 'Mg26': 0.022058, 'N14': 0.199267, 'N15': 0.000732, 'O16': 0.599772, 'O17': 0.000227} From 3669709061f6dc51591bdb46ee6dbf6170b500d7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 25 Mar 2020 06:29:21 -0400 Subject: [PATCH 19/24] Apply suggestions from code review * Use set operations (|=) rather than in-place methods (update) * Remove second arguments to some asserts in Chain.reduce testing Co-Authored-By: Paul Romano --- openmc/deplete/chain.py | 8 ++++---- tests/unit_tests/test_deplete_chain.py | 26 ++++++++++++------------ tests/unit_tests/test_deplete_nuclide.py | 4 ++-- 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index d820fe520b..fce5051b0b 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -906,7 +906,7 @@ class Chain(object): def _follow(self, isotopes, level): """Return all isotopes present up to depth level""" - found = set(isotopes) + found = isotopes.copy() remaining = set(self.nuclide_dict) if not found.issubset(remaining): raise IndexError( @@ -916,7 +916,7 @@ class Chain(object): if level == 0: return found - remaining.difference_update(found) + remaining -= found depth = 0 next_iso = set() @@ -952,8 +952,8 @@ class Chain(object): # Prepare for next dig depth += 1 - isotopes.update(next_iso) - remaining.difference_update(next_iso) + isotopes |= next_iso + remaining -= next_iso next_iso.clear() # Process isotope that would have started next depth diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index b7a1e6c0f2..ba3be4e37a 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -479,18 +479,18 @@ def test_reduce(gnd_simple_chain): assert u5_round0.n_decay_modes == ref_U5.n_decay_modes for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes): assert newmode.target is None - assert newmode.type == refmode.type, newmode - assert newmode.branching_ratio == refmode.branching_ratio, newmode + assert newmode.type == refmode.type + assert newmode.branching_ratio == refmode.branching_ratio assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions): - assert newrxn.target is None, newrxn - assert newrxn.type == refrxn.type, newrxn - assert newrxn.Q == refrxn.Q, newrxn - assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn + assert newrxn.target is None + assert newrxn.type == refrxn.type + assert newrxn.Q == refrxn.Q + assert newrxn.branching_ratio == refrxn.branching_ratio assert u5_round0.yield_data is not None - assert u5_round0.yield_data.products == ("I135", ) + assert u5_round0.yield_data.products == ("I135",) assert u5_round0.yield_data.yield_matrix == ( ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")] ) @@ -499,15 +499,15 @@ def test_reduce(gnd_simple_chain): assert bareI5.n_decay_modes == ref_iodine.n_decay_modes for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes): assert newmode.target is None - assert newmode.type == refmode.type, newmode - assert newmode.branching_ratio == refmode.branching_ratio, newmode + assert newmode.type == refmode.type + assert newmode.branching_ratio == refmode.branching_ratio assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions): - assert newrxn.target is None, newrxn - assert newrxn.type == refrxn.type, newrxn - assert newrxn.Q == refrxn.Q, newrxn - assert newrxn.branching_ratio == refrxn.branching_ratio, newrxn + assert newrxn.target is None + assert newrxn.type == refrxn.type + assert newrxn.Q == refrxn.Q + assert newrxn.branching_ratio == refrxn.branching_ratio follow_u5 = gnd_simple_chain.reduce(["U235"], 1) u5_round1 = follow_u5["U235"] diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 92b004845a..126206cc70 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -181,9 +181,9 @@ def test_fission_yield_distribution(): with_extras = yield_dist.restrict_products( ["Xe135", "Sm149", "H1", "U235"]) - assert strict_restrict.products == ("Sm149", "Xe135", ) + assert strict_restrict.products == ("Sm149", "Xe135") assert strict_restrict.energies == yield_dist.energies - assert with_extras.products == ("Sm149", "Xe135", ) + assert with_extras.products == ("Sm149", "Xe135") assert with_extras.energies == yield_dist.energies for ene, new_yields in strict_restrict.items(): for product in strict_restrict.products: From c2562322ab7a788bdc09ed1846a2934582756060 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 25 Mar 2020 06:48:43 -0400 Subject: [PATCH 20/24] Expand Chain.reduce docstring: methodology and return The methodology is explained through a simple example mostly borrowed from the chain_simple.xml test chain. A few different cases are presented, showing the effect of initial isotopes and search depth. A description of the returned chain is also included. --- openmc/deplete/chain.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index fce5051b0b..a23ba2bc1c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -822,7 +822,30 @@ class Chain(object): return valid def reduce(self, initial_isotopes, level=None): - """Reduce the size of the chain by follwing transmutation paths + """Reduce the size of the chain by following transmutation paths + + As an example, consider a simple chain with the following + isotopes and transmutation paths:: + + U235 (n,gamma) U236 + (n,fission) (Xe135, I135, Cs135) + I135 (beta decay) Xe135 (beta decay) Cs135 + Xe135 (n,gamma) Xe136 + + Calling ``chain.reduce(["I135"])`` will produce a depletion + chain that contains only isotopes that would originate from + I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not + be included, but multiple isotopes can be used to start + the search. + + The ``level`` value controls the depth of the search. + ``chain.reduce(["U235"], level=1)`` would return a chain + with all isotopes except Xe136, since it is two transmutations + removed from U235 in this case. + + While targets will not be included in the new chain, the + total destruction rate and decay rate of included isotopes + will be preserved. Parameters ---------- @@ -839,6 +862,8 @@ class Chain(object): Returns ------- Chain + Depletion chain containing isotopes that would appear + after following up to ``level`` reactions and decay paths """ check_type("initial_isotopes", initial_isotopes, Iterable, str) From d3368db5ee030bdedf5e6b039ddda2c126aaeafe Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 25 Mar 2020 07:10:32 -0400 Subject: [PATCH 21/24] Use reduce_chain and reduce_chain_level with Operator Control if the chain is to be reduced, and how far down the paths to follow. Defaults to not reducing the depletion chain. If the chain is reduced, the default depth is None, implying no limit --- openmc/deplete/operator.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 2572150012..5fa6745c8e 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -113,12 +113,13 @@ class Operator(TransportOperator): ``fission_yield_mode``. Will be passed directly on to the helper. Passing a value of None will use the defaults for the associated helper. - reduce_chain : bool or int, optional - If ``True`` or an integer, create a reduced depletion chain by - following transmuation paths for isotopes in burnable materials. - A value of ``True`` implies to follow all paths to completion, - while a non-negative integer indicates the depth. See - :meth:`openmc.deplete.Chain.reduce` + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. Default is False. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. Attributes ---------- @@ -165,7 +166,7 @@ class Operator(TransportOperator): diff_burnable_mats=False, energy_mode="fission-q", fission_q=None, dilute_initial=1.0e3, fission_yield_mode="constant", fission_yield_opts=None, - reduce_chain=False): + reduce_chain=False, reduce_chain_level=None): if fission_yield_mode not in self._fission_helpers: raise KeyError( "fission_yield_mode must be one of {}, not {}".format( @@ -187,15 +188,14 @@ class Operator(TransportOperator): self.diff_burnable_mats = diff_burnable_mats # Reduce the chain before we create more materials - if reduce_chain is not False: + if reduce_chain: all_isotopes = set() for material in geometry.get_all_materials().values(): if not material.depletable: continue for name, _dens_percent, _dens_type in material.nuclides: all_isotopes.add(name) - level = None if reduce_chain is True else reduce_chain - self.chain = self.chain.reduce(all_isotopes, level) + self.chain = self.chain.reduce(all_isotopes, reduce_chain_level) # Differentiate burnable materials with multiple instances if self.diff_burnable_mats: From f882cf666bb02a2650342bda9842fcf804f91cbc Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 25 Mar 2020 07:12:40 -0400 Subject: [PATCH 22/24] Remove unused imports from operator.py time, itertools.chain, and h5py not used --- openmc/deplete/operator.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 5fa6745c8e..ce7e709d9c 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem. import sys import copy from collections import OrderedDict -from itertools import chain import os -import time import xml.etree.ElementTree as ET from warnings import warn -import h5py import numpy as np from uncertainties import ufloat From 4234fe45e71a7a3971e06c0733c253e96222bb8f Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 25 Mar 2020 14:58:02 +0000 Subject: [PATCH 23/24] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/material.py | 4 ++-- tests/unit_tests/test_material.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 7cc7e3db60..e34123e934 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -607,9 +607,9 @@ class Material(IDManagerMixin): If enrichment_target is not supplied then it is enrichment for U235 in weight percent. For example, input 4.95 for 4.95 weight percent enriched U. Default is None (natural composition). - enrichment_target: str, optional + enrichment_target : str, optional Single nuclide name to enrich from a natural composition (e.g., 'O16') - enrichment_type: {'ao', 'wo'}, optional + enrichment_type : {'ao', 'wo'}, optional 'ao' for enrichment as atom percent and 'wo' for weight percent. Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 8e9ebb362d..e366401827 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -58,7 +58,8 @@ def test_elements_by_name(): assert a._nuclides == b._nuclides assert b._nuclides == c._nuclides -def test_adding_elements_by_formula(): + +def test_add_elements_by_formula(): """Test adding elements from a formula""" # testing the correct nuclides and elements are added to a material m = openmc.Material() @@ -354,4 +355,3 @@ def test_mix_materials(): assert m3.density == pytest.approx(dens3) assert m4.density == pytest.approx(dens4) assert m5.density == pytest.approx(dens5) - From c013f6373075c64205511a182e1632bd87dea1a3 Mon Sep 17 00:00:00 2001 From: =shimwell Date: Wed, 25 Mar 2020 15:03:39 +0000 Subject: [PATCH 24/24] pep8 format applied to imports --- tests/unit_tests/test_material.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index e366401827..c47faee8cd 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -1,9 +1,12 @@ +from collections import defaultdict + +import pytest + import openmc +import openmc.examples import openmc.model import openmc.stats -import openmc.examples -import pytest -from collections import defaultdict + def test_attributes(uo2): assert uo2.name == 'UO2'