crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1317 from drewejohnson/chain-set-branches-general

Browse source 
Chain set branches general
This commit is contained in:
Paul Romano 2019-08-25 08:58:20 -05:00 • committed by GitHub
commit 83028b8dcd
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 174 additions and 64 deletions
Download patch file Download diff file Expand all files Collapse all files

157
openmc/deplete/chain.py
View file

@ -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

81
tests/unit_tests/test_deplete_chain.py
View file

@ -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):