diff --git a/docs/source/usersguide/scripts.rst b/docs/source/usersguide/scripts.rst
index 76c7c2f22..388130947 100644
--- a/docs/source/usersguide/scripts.rst
+++ b/docs/source/usersguide/scripts.rst
@@ -106,6 +106,37 @@ Compton profile data using an existing data library from `Geant4
`_. Note that OpenMC includes this data file by default
so it should not be necessary in practice to generate it yourself.
+
+.. _scripts_depletion_chain:
+
+-------------------------------
+``openmc-make-depletion-chain``
+-------------------------------
+
+This script generates a depletion chain file called ``chain_endfb71.xml``
+using ENDF/B-VII.1 nuclear data. If the :envvar:`OPENMC_ENDF_DATA` variable
+is not set, and ``"neutron"``, ``"decay"``, ``"nfy"`` directories
+do not exist, then ENDF/B-VII.1 data will be downloaded.
+
+.. _scripts_depletion_chain_casl:
+
+------------------------------------
+``openmc-make-depletion-chain-casl``
+------------------------------------
+
+This script generates a depletion chain called ``chain_casl.xml``
+using ENDF/B-VII.1 nuclear data for a simplified chain.
+The nuclides were chosen by CASL-ORIGEN, which can be found in
+Appendix A of Kang Seog Kim, `"Specification for the VERA Depletion
+Benchmark Suite" `_,
+CASL-U-2015-1014-000, Rev. 0, ORNL/TM-2016/53, 2016.
+``Te129`` has been added into this chain due to its link to
+``I129`` production.
+
+If the :envvar:`OPENMC_ENDF_DATA` variable is not set,
+and ``"neutron"``, ``"decay"``, ``"nfy"`` directories
+to not exist, then ENDF/B-VII.1 data will be downloaded.
+
.. _scripts_stopping:
-------------------------------
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 3cb3dbf01..04fdf658e 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -10,9 +10,10 @@ import math
import re
from collections import OrderedDict, defaultdict
from collections.abc import Mapping
+from numbers import Real
from warnings import warn
-from openmc.checkvalue import check_type
+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
@@ -677,3 +678,54 @@ class Chain(object):
new_ratios[ground_tgt] = ground_br
parent.reactions.append(ReactionTuple(
reaction, ground_tgt, rxn_Q, ground_br))
+
+ def validate(self, strict=True, quiet=False, tolerance=1e-4):
+ """Search for possible inconsistencies
+
+ The following checks are performed for all nuclides present:
+
+ 1) For all non-fission reactions, does the sum of branching
+ ratios equal about one?
+ 2) For fission reactions, does the sum of fission yield
+ fractions equal about two?
+
+ Parameters
+ ----------
+ strict : bool, optional
+ Raise exceptions at the first inconsistency if true.
+ Otherwise mark a warning
+ quiet : bool, optional
+ Flag to suppress warnings and return immediately at
+ the first inconsistency. Used only if
+ ``strict`` does not evaluate to ``True``.
+ tolerance : float, optional
+ Absolute tolerance for comparisons. Used to compare computed
+ value ``x`` to intended value ``y`` as::
+
+ valid = (y - tolerance <= x <= y + tolerance)
+
+ Returns
+ -------
+ valid : bool
+ True if no inconsistencies were found
+
+ Raises
+ ------
+ ValueError
+ If ``strict`` evaluates to ``True`` and an inconistency was
+ found
+
+ See Also
+ --------
+ openmc.deplete.Nuclide.validate
+ """
+ check_type("tolerance", tolerance, Real)
+ check_greater_than("tolerance", tolerance, 0.0, True)
+ valid = True
+ # Sort through nuclides by name
+ for name in sorted(self.nuclide_dict):
+ stat = self[name].validate(strict, quiet, tolerance)
+ if quiet and not stat:
+ return stat
+ valid = valid and stat
+ return valid
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 88d5ac6ea..538b69790 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -3,7 +3,8 @@
Contains the per-nuclide components of a depletion chain.
"""
-from collections import namedtuple
+from collections import namedtuple, defaultdict
+from warnings import warn
try:
import lxml.etree as ET
except ImportError:
@@ -222,3 +223,101 @@ class Nuclide(object):
data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E])
return elem
+
+ def validate(self, strict=True, quiet=False, tolerance=1e-4):
+ """Search for possible inconsistencies
+
+ The following checks are performed:
+
+ 1) for all non-fission reactions and decay modes,
+ does the sum of branching ratios equal about one?
+ 2) for fission reactions, does the sum of fission yield
+ fractions equal about two?
+
+ Parameters
+ ----------
+ strict : bool, optional
+ Raise exceptions at the first inconsistency if true.
+ Otherwise mark a warning
+ quiet : bool, optional
+ Flag to suppress warnings and return immediately at
+ the first inconsistency. Used only if
+ ``strict`` does not evaluate to ``True``.
+ tolerance : float, optional
+ Absolute tolerance for comparisons. Used to compare computed
+ value ``x`` to intended value ``y`` as::
+
+ valid = (y - tolerance <= x <= y + tolerance)
+
+ Returns
+ -------
+ valid : bool
+ True if no inconsistencies were found
+
+ Raises
+ ------
+ ValueError
+ If ``strict`` evaluates to ``True`` and an inconistency was
+ found
+
+ See Also
+ --------
+ openmc.deplete.Chain.validate
+ """
+
+ msg_func = ("Nuclide {name} has {prop} that sum to {actual} "
+ "instead of {expected} +/- {tol:7.4e}").format
+ valid = True
+
+ # check decay modes
+ if self.decay_modes:
+ sum_br = sum(m.branching_ratio for m in self.decay_modes)
+ stat = 1.0 - tolerance <= sum_br <= 1.0 + tolerance
+ if not stat:
+ msg = msg_func(
+ name=self.name, actual=sum_br, expected=1.0, tol=tolerance,
+ prop="decay mode branch ratios")
+ if strict:
+ raise ValueError(msg)
+ elif quiet:
+ return False
+ warn(msg)
+ valid = False
+
+ if self.reactions:
+ type_map = defaultdict(set)
+ for reaction in self.reactions:
+ type_map[reaction.type].add(reaction)
+ for rxn_type, reactions in type_map.items():
+ sum_rxn = sum(rx.branching_ratio for rx in reactions)
+ stat = 1.0 - tolerance <= sum_rxn <= 1.0 + tolerance
+ if stat:
+ continue
+ msg = msg_func(
+ name=self.name, actual=sum_br, expected=1.0, tol=tolerance,
+ prop="{} reaction branch ratios".format(rxn_type))
+ if strict:
+ raise ValueError(msg)
+ elif quiet:
+ return False
+ warn(msg)
+ 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)
+ stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance
+ if stat:
+ continue
+ msg = msg_func(
+ name=self.name, actual=sum_yield,
+ expected=2.0, tol=tolerance,
+ prop="fission yields (E = {:7.4e} eV)".format(energy))
+ if strict:
+ raise ValueError(msg)
+ elif quiet:
+ return False
+ warn(msg)
+ valid = False
+
+ return valid
diff --git a/scripts/casl_chain.py b/scripts/casl_chain.py
index 86389e3ce..94fc9dd96 100755
--- a/scripts/casl_chain.py
+++ b/scripts/casl_chain.py
@@ -6,6 +6,8 @@
# Note 32 of the 255 nuclides appeare twice as they are both activation
# nuclides (category 1) and fission product nuclides (category 3).
+# Te129 has been added due to it's link to I129 production.
+
CASL_CHAIN = {
# Nuclide: (Stable, CAT, IFPY, Special yield treatment)
# Stable: True if nuclide has no decay reactions
@@ -187,6 +189,7 @@ CASL_CHAIN = {
'Sb127': (False, 3, 2, None),
'Te127': (False, 3, -1, None),
'Te127_m1': (False, 3, -1, None),
+ 'Te129': (False, 3, 1, None),
'Te129_m1': (False, 3, 2, None),
'Te132': (False, 3, 2, None),
'I127': (True, 3, 1, None),
diff --git a/scripts/openmc-make-depletion-chain b/scripts/openmc-make-depletion-chain
index 7c08f984e..51092a722 100755
--- a/scripts/openmc-make-depletion-chain
+++ b/scripts/openmc-make-depletion-chain
@@ -1,7 +1,7 @@
#!/usr/bin/env python3
-import glob
import os
+from pathlib import Path
from zipfile import ZipFile
from openmc._utils import download
@@ -15,15 +15,28 @@ URLS = [
]
def main():
- for url in URLS:
- basename = download(url)
- with ZipFile(basename, 'r') as zf:
- print('Extracting {}...'.format(basename))
- zf.extractall()
+ endf_dir = os.environ.get("OPENMC_ENDF_DATA")
+ if endf_dir is not None:
+ endf_dir = Path(endf_dir)
+ elif all(os.path.isdir(lib) for lib in ("neutrons", "decay", "nfy")):
+ endf_dir = Path(".")
+ else:
+ for url in URLS:
+ basename = download(url)
+ with ZipFile(basename, 'r') as zf:
+ print('Extracting {}...'.format(basename))
+ zf.extractall()
+ endf_dir = Path(".")
- decay_files = glob.glob(os.path.join('decay', '*.endf'))
- nfy_files = glob.glob(os.path.join('nfy', '*.endf'))
- neutron_files = glob.glob(os.path.join('neutrons', '*.endf'))
+ decay_files = tuple((endf_dir / "decay").glob("*endf"))
+ neutron_files = tuple((endf_dir / "neutrons").glob("*endf"))
+ nfy_files = tuple((endf_dir / "nfy").glob("*endf"))
+
+ # check files exist
+ for flist, ftype in [(decay_files, "decay"), (neutron_files, "neutron"),
+ (nfy_files, "neutron fission product yield")]:
+ if not flist:
+ raise IOError("No {} endf files found in {}".format(ftype, endf_dir))
chain = openmc.deplete.Chain.from_endf(decay_files, nfy_files, neutron_files)
chain.export_to_xml('chain_endfb71.xml')
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index 3cf726ed1..1fd974e05 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -330,7 +330,6 @@ def test_capture_branch_failures(simple_chain):
simple_chain.set_branch_ratios(br, "(n,gamma)")
-
def test_set_alpha_branches():
"""Test setting of alpha reaction branching ratios"""
# Build a mock chain
@@ -377,3 +376,55 @@ def test_set_alpha_branches():
break
else:
raise ValueError("Helium has been removed and should not have been")
+
+
+def test_validate(simple_chain):
+ """Test the validate method"""
+
+ # current chain is invalid
+ # fission yields do not sum to 2.0
+ with pytest.raises(ValueError, match="Nuclide C.*fission yields"):
+ simple_chain.validate(strict=True, tolerance=0.0)
+
+ with pytest.warns(UserWarning) as record:
+ assert not simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
+ assert not simple_chain.validate(strict=False, quiet=True, tolerance=0.0)
+ assert len(record) == 1
+ assert "Nuclide C" in record[0].message.args[0]
+
+ # 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)]}
+
+ assert simple_chain.validate(strict=True, tolerance=0.0)
+ with pytest.warns(None) as record:
+ assert simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
+ assert len(record) == 0
+
+ # Mess up "earlier" nuclide's reactions
+ decay_mode = simple_chain["A"].decay_modes.pop()
+
+ with pytest.raises(ValueError, match="Nuclide A.*decay mode"):
+ simple_chain.validate(strict=True, tolerance=0.0)
+
+ # restore old fission yields
+ simple_chain["C"].yield_data = old_yields
+
+ with pytest.warns(UserWarning) as record:
+ assert not simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
+ assert len(record) == 2
+ assert "Nuclide A" in record[0].message.args[0]
+ assert "Nuclide C" in record[1].message.args[0]
+
+ # restore decay modes
+ simple_chain["A"].decay_modes.append(decay_mode)
+
+
+def test_validate_inputs():
+ c = Chain()
+
+ with pytest.raises(TypeError, match="tolerance"):
+ c.validate(tolerance=None)
+
+ with pytest.raises(ValueError, match="tolerance"):
+ c.validate(tolerance=-1)
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index f2a101d2a..38848e9c1 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 pytest
+
from openmc.deplete import nuclide
@@ -114,3 +116,92 @@ 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_validate():
+
+ nuc = nuclide.Nuclide()
+ nuc.name = "Test"
+
+ # decay modes: type, target, branching_ratio
+
+ nuc.decay_modes = [
+ nuclide.DecayTuple("type 0", "0", 0.5),
+ nuclide.DecayTuple("type 1", "1", 0.5),
+ ]
+
+ # reactions: type, target, Q, branching_ratio
+ nuc.reactions = [
+ nuclide.ReactionTuple("0", "0", 1000, 0.3),
+ nuclide.ReactionTuple("0", "1", 1000, 0.3),
+ nuclide.ReactionTuple("1", "2", 1000, 1.0),
+ nuclide.ReactionTuple("0", "3", 1000, 0.4),
+ ]
+
+ # fission yields
+
+ nuc.yield_data = {
+ 0.0253: [("0", 1.5), ("1", 0.5)],
+ 1e6: [("0", 1.5), ("1", 0.5)],
+ }
+
+ # nuclide is good and should have no warnings raise
+ with pytest.warns(None) as record:
+ assert nuc.validate(strict=True, quiet=False, tolerance=0.0)
+ assert len(record) == 0
+
+ # invalidate decay modes
+ decay = nuc.decay_modes.pop()
+ with pytest.raises(ValueError, match="decay mode"):
+ nuc.validate(strict=True, quiet=False, tolerance=0.0)
+
+ with pytest.warns(UserWarning) as record:
+ assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
+ assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
+ assert len(record) == 1
+ assert "decay mode" in record[0].message.args[0]
+
+ # restore decay modes, invalidate reactions
+ nuc.decay_modes.append(decay)
+ reaction = nuc.reactions.pop()
+
+ with pytest.raises(ValueError, match="0 reaction"):
+ nuc.validate(strict=True, quiet=False, tolerance=0.0)
+
+ with pytest.warns(UserWarning) as record:
+ assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
+ assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
+ assert len(record) == 1
+ assert "0 reaction" in record[0].message.args[0]
+
+ # restore reactions, invalidate fission yields
+ nuc.reactions.append(reaction)
+ nuc.yield_data[1e6].pop()
+
+ with pytest.raises(ValueError, match=r"fission yields.*1\.0*e"):
+ nuc.validate(strict=True, quiet=False, tolerance=0.0)
+
+ with pytest.warns(UserWarning) as record:
+ assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
+ assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
+ assert len(record) == 1
+ assert "1.0" in record[0].message.args[0]
+
+ # invalidate everything, check that error is raised at decay modes
+
+ decay = nuc.decay_modes.pop()
+ reaction = nuc.reactions.pop()
+
+ with pytest.raises(ValueError, match="decay mode"):
+ nuc.validate(strict=True, quiet=False, tolerance=0.0)
+
+ # check for warnings
+ # should be one warning for decay modes, reactions, fission yields
+
+ with pytest.warns(UserWarning) as record:
+ assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
+ assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
+ assert len(record) == 3
+ assert "decay mode" in record[0].message.args[0]
+ assert "0 reaction" in record[1].message.args[0]
+ assert "1.0" in record[2].message.args[0]