diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 4706bfabb..b3e2dae29 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -13,7 +13,7 @@ from collections.abc import Mapping from numbers import Real from warnings import warn -from openmc.checkvalue import check_type, check_less_than, check_greater_than +from openmc.checkvalue import check_type, check_greater_than from openmc.data import gnd_name, zam # Try to use lxml if it is available. It preserves the order of attributes and @@ -104,6 +104,19 @@ def replace_missing(product, decay_data): return product +_SECONDARY_PARTICLES = { + "(n,p)": ["H1"], "(n,d)": ["H2"], "(n,t)": ["H3"], "(n,3He)": ["He3"], + "(n,a)": ["He4"], "(n,2nd)": ["H2"], "(n,na)": ["He4"], "(n,3na)": ["He4"], + "(n,n3a)": ["He4"] * 3, "(n,2na)": ["He4"], "(n,np)": ["H1"], + "(n,n2a)": ["He4"] * 2, "(n,2n2a)": ["He4"] * 2, "(n,nd)": ["H2"], + "(n,nt)": ["H3"], "(n,nHe-3)": ["He3"], "(n,nd2a)": ["H2", "He4"], + "(n,nt2a)": ["H3", "He4", "He4"], "(n,2np)": ["H1"], "(n,3np)": ["H1"], + "(n,n2p)": ["H1"] * 2, "(n,2a)": ["He4"] * 2, "(n,3a)": ["He4"] * 3, + "(n,2p)": ["H1"] * 2, "(n,pa)": ["H1", "He4"], + "(n,t2a)": ["H3", "He4", "He4"], "(n,d2a)": ["H2", "He4", "He4"], + "(n,pd)": ["H1", "H2"], "(n,pt)": ["H1", "H3"], "(n,da)": ["H2", "He4"]} + + class Chain(object): """Full representation of a depletion chain. @@ -123,7 +136,6 @@ 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. - """ def __init__(self): @@ -455,70 +467,99 @@ class Chain(object): dict.update(matrix_dok, matrix) return matrix_dok.tocsr() - def get_capture_branches(self): - """Return a dictionary with capture branching ratios + def get_branch_ratios(self, reaction="(n,gamma)"): + """Return a dictionary with reaction branching ratios + + Parameters + ---------- + reaction : str, optional + Reaction name like ``"(n,gamma)"`` [default], or + ``"(n,alpha)"``. Returns ------- - capt : - nested dict of parent nuclide keys with capture targets and - branching ratios:: + branches : dict + nested dict of parent nuclide keys with reaction targets and + branching ratios. Consider the capture, ``"(n,gamma)"``, + reaction for Am241:: {"Am241": {"Am242": 0.91, "Am242_m1": 0.09}} See Also -------- - :meth:`set_capture_branches` - + :meth:`set_branch_ratios` """ capt = {} for nuclide in self.nuclides: nuc_capt = {} for rx in nuclide.reactions: - if rx.type == "(n,gamma)" and rx.branching_ratio != 1.0: + if rx.type == reaction and rx.branching_ratio != 1.0: nuc_capt[rx.target] = rx.branching_ratio if len(nuc_capt) > 0: capt[nuclide.name] = nuc_capt return capt - def set_capture_branches(self, branch_ratios, strict=True): - """Set the capture branching ratios - - To provide a buffer around floating point precisions, - the sum of all branching ratios from a single parent - cannot be greater than 1.00001. + def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", + strict=True, tolerance=1e-5): + """Set the branching ratios for a given reactions Parameters ---------- branch_ratios : dict of {str: {str: float}} Capture branching ratios to be inserted. First layer keys are names of parent nuclides, e.g. - ``"Am241"``. The capture branching ratios for these + ``"Am241"``. The branching ratios for these parents will be modified. Corresponding values are dictionaries of ``{target: branching_ratio}`` - strict : bool - If this evalutes to ``True``, then all parents and - products must exist in the :class:`Chain`. A - :class:`KeyError` will be raised at the first - nuclide that does not exist. Otherwise, print - a warning message for missing parents and/or - products. + reaction : str, optional + Reaction name like ``"(n,gamma)"`` [default], or + ``"(n, alpha)"``. + strict : bool, optional + Error control. If this evalutes to ``True``, then errors will + be raised if inconsistencies are found. Otherwise, warnings + will be raised for most issues. + tolerance : float, optional + Tolerance on the sum of all branching ratios for a + single parent. Will be checked with:: + + 1 - tol < sum_br < 1 + tol + + Raises + ------ + IndexError + If no isotopes were found on the chain that have the requested + reaction + KeyError + If ``strict`` evaluates to ``False`` and a parent isotope in + ``branch_ratios`` does not exist on the chain + AttributeError + If ``strict`` evaluates to ``False`` and a parent isotope in + ``branch_ratios`` does not have the requested reaction + ValueError + If ``strict`` evalutes to ``False`` and the sum of one parents + branch ratios is outside 1 +/- ``tolerance`` See Also -------- - :meth:`get_capture_branches` + :meth:`get_branch_ratios` """ # Store some useful information through the validation stage sums = {} - capt_ix_map = {} + rxn_ix_map = {} grounds = {} + tolerance = abs(tolerance) + missing_parents = set() missing_products = {} - no_capture = set() + missing_reaction = set() + bad_sums = {} + + # Secondary products, like alpha particles, should not be modified + secondary = _SECONDARY_PARTICLES.get(reaction, []) # Check for validity before manipulation @@ -544,11 +585,11 @@ class Chain(object): if prod_flag: continue - # Make sure this nuclide has capture reactions + # Make sure this nuclide has the reaction indexes = [] for ix, rx in enumerate(self[parent].reactions): - if rx.type == "(n,gamma)": + if rx.type == reaction and rx.target not in secondary: indexes.append(ix) if "_m" not in rx.target: grounds[parent] = rx.target @@ -556,24 +597,39 @@ class Chain(object): if len(indexes) == 0: if strict: raise AttributeError( - "Nuclide {} does not have capture reactions in " - "this {}".format(parent, self.__class__.__name__)) - no_capture.add(parent) + "Nuclide {} does not have {} reactions".format( + parent, reaction)) + missing_reaction.add(parent) continue - capt_ix_map[parent] = indexes - this_sum = sum(sub.values()) - check_less_than(parent + " ratios", this_sum, 1.00001) - sums[parent] = this_sum + # sum of branching ratios can be lower than 1 if no ground + # target is given, but never greater + if (this_sum >= 1 + tolerance or (grounds[parent] in sub + and this_sum <= 1 - tolerance)): + if strict: + msg = ("Sum of {} branching ratios for {} " + "({:7.3f}) outside tolerance of 1 +/- " + "{:5.3e}".format( + reaction, parent, this_sum, tolerance)) + raise ValueError(msg) + bad_sums[parent] = this_sum + else: + rxn_ix_map[parent] = indexes + sums[parent] = this_sum + + if len(rxn_ix_map) == 0: + raise IndexError( + "No {} reactions found in this {}".format( + reaction, self.__class__.__name__)) if len(missing_parents) > 0: warn("The following nuclides were not found in {}: {}".format( self.__class__.__name__, ", ".join(sorted(missing_parents)))) - if len(no_capture) > 0: - warn("The following nuclides did not have capture reactions: " - "{}".format(", ".join(sorted(no_capture)))) + if len(missing_reaction) > 0: + warn("The following nuclides did not have {} reactions: " + "{}".format(reaction, ", ".join(sorted(missing_reaction)))) if len(missing_products) > 0: tail = ("{} -> {}".format(k, v) @@ -582,28 +638,35 @@ class Chain(object): "parents were unmodified: \n{}".format( self.__class__.__name__, ", ".join(tail))) + if len(bad_sums) > 0: + tail = ("{}: {:5.3f}".format(k, s) + for k, s in sorted(bad_sums.items())) + warn("The following parent nuclides were given {} branch ratios " + "with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format( + reaction, tolerance, "\n".join(tail))) + # Insert new ReactionTuples with updated branch ratios - for parent_name, capt_index in capt_ix_map.items(): + for parent_name, rxn_index in rxn_ix_map.items(): parent = self[parent_name] new_ratios = branch_ratios[parent_name] - capt_index = capt_ix_map[parent_name] + rxn_index = rxn_ix_map[parent_name] # Assume Q value is independent of target state - capt_Q = parent.reactions[capt_index[0]].Q + rxn_Q = parent.reactions[rxn_index[0]].Q - # Remove existing capture reactions + # Remove existing reactions - for ix in reversed(capt_index): + for ix in reversed(rxn_index): parent.reactions.pop(ix) all_meta = True for tgt, br in new_ratios.items(): - all_meta = all_meta and ("_m" in tgt) + all_meta = all_meta and ("_m" in tgt) parent.reactions.append(ReactionTuple( - "(n,gamma)", tgt, capt_Q, br)) + reaction, tgt, rxn_Q, br)) if all_meta and sums[parent_name] != 1.0: ground_br = 1.0 - sums[parent_name] @@ -613,7 +676,7 @@ class Chain(object): ground_tgt = gnd_name(pz, pa + 1, 0) new_ratios[ground_tgt] = ground_br parent.reactions.append(ReactionTuple( - "(n,gamma)", ground_tgt, capt_Q, ground_br)) + reaction, ground_tgt, rxn_Q, ground_br)) def validate(self, strict=True, quiet=False, tolerance=1e-4): """Search for possible inconsistencies diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 294e56637..1fd974e05 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -248,29 +248,29 @@ def test_set_fiss_q(): def test_get_set_chain_br(simple_chain): """Test minor modifications to capture branch ratios""" expected = {"C": {"A": 0.7, "B": 0.3}} - assert simple_chain.get_capture_branches() == expected + assert simple_chain.get_branch_ratios() == expected # safely modify new_chain = Chain.from_xml("chain_test.xml") new_br = {"C": {"A": 0.5, "B": 0.5}, "A": {"C": 0.99, "B": 0.01}} - new_chain.set_capture_branches(new_br) - assert new_chain.get_capture_branches() == new_br + new_chain.set_branch_ratios(new_br) + assert new_chain.get_branch_ratios() == new_br # write, re-read new_chain.export_to_xml("chain_mod.xml") - assert Chain.from_xml("chain_mod.xml").get_capture_branches() == new_br + assert Chain.from_xml("chain_mod.xml").get_branch_ratios() == new_br # Test non-strict [warn, not error] setting bad_br = {"B": {"X": 0.6, "A": 0.4}, "X": {"A": 0.5, "C": 0.5}} bad_br.update(new_br) - new_chain.set_capture_branches(bad_br, strict=False) - assert new_chain.get_capture_branches() == new_br + new_chain.set_branch_ratios(bad_br, strict=False) + assert new_chain.get_branch_ratios() == new_br # Ensure capture reactions are removed rem_br = {"A": {"C": 1.0}} - new_chain.set_capture_branches(rem_br) + new_chain.set_branch_ratios(rem_br) # A is not in returned dict because there is no branch - assert "A" not in new_chain.get_capture_branches() + assert "A" not in new_chain.get_branch_ratios() def test_capture_branch_infer_ground(): @@ -289,9 +289,9 @@ def test_capture_branch_infer_ground(): chain.nuclides.append(xe136m) chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1 - chain.set_capture_branches(infer_br) + chain.set_branch_ratios(infer_br, "(n,gamma)") - assert chain.get_capture_branches() == set_br + assert chain.get_branch_ratios("(n,gamma)") == set_br def test_capture_branch_no_rxn(): @@ -307,9 +307,8 @@ def test_capture_branch_no_rxn(): chain.nuclides.append(u5m) chain.nuclide_dict[u5m.name] = len(chain.nuclides) - 1 - phrase = "U234 does not have capture reactions" - with pytest.raises(AttributeError, match=phrase): - chain.set_capture_branches(u4br) + with pytest.raises(AttributeError, match="U234"): + chain.set_branch_ratios(u4br) def test_capture_branch_failures(simple_chain): @@ -318,17 +317,65 @@ def test_capture_branch_failures(simple_chain): # Parent isotope not present br = {"X": {"A": 0.6, "B": 0.7}} with pytest.raises(KeyError, match="X"): - simple_chain.set_capture_branches(br) + simple_chain.set_branch_ratios(br) # Product isotope not present br = {"C": {"X": 0.4, "A": 0.2, "B": 0.4}} with pytest.raises(KeyError, match="X"): - simple_chain.set_capture_branches(br) + simple_chain.set_branch_ratios(br) # Sum of ratios > 1.0 br = {"C": {"A": 1.0, "B": 1.0}} - with pytest.raises(ValueError, match="C ratios"): - simple_chain.set_capture_branches(br) + with pytest.raises(ValueError, match=r"Sum of \(n,gamma\).*for C"): + simple_chain.set_branch_ratios(br, "(n,gamma)") + + +def test_set_alpha_branches(): + """Test setting of alpha reaction branching ratios""" + # Build a mock chain + chain = Chain() + + parent = nuclide.Nuclide() + parent.name = "A" + + he4 = nuclide.Nuclide() + he4.name = "He4" + + ground_tgt = nuclide.Nuclide() + ground_tgt.name = "B" + + meta_tgt = nuclide.Nuclide() + meta_tgt.name = "B_m1" + + for ix, nuc in enumerate((parent, ground_tgt, meta_tgt, he4)): + chain.nuclides.append(nuc) + chain.nuclide_dict[nuc.name] = ix + + # add reactions to parent + parent.reactions.append(nuclide.ReactionTuple( + "(n,a)", ground_tgt.name, 1.0, 0.6)) + parent.reactions.append(nuclide.ReactionTuple( + "(n,a)", meta_tgt.name, 1.0, 0.4)) + parent.reactions.append(nuclide.ReactionTuple( + "(n,a)", he4.name, 1.0, 1.0)) + + expected_ref = {"A": {"B": 0.6, "B_m1": 0.4}} + + assert chain.get_branch_ratios("(n,a)") == expected_ref + + # alter and check again + + altered = {"A": {"B": 0.5, "B_m1": 0.5}} + + chain.set_branch_ratios(altered, "(n,a)") + assert chain.get_branch_ratios("(n,a)") == altered + + # make sure that alpha particle still produced + for r in parent.reactions: + if r.target == he4.name: + break + else: + raise ValueError("Helium has been removed and should not have been") def test_validate(simple_chain):