diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index 37c8075d4..debc10403 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -313,19 +313,24 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
thermal = yields.get(thermal_energy)
fast = yields.get(fast_energy)
if thermal is None or fast is None:
- insert_ix = bisect.bisect_left(energies, cutoff)
+ cutoff_ix = bisect.bisect_left(energies, cutoff)
# if zero, then all energies > cutoff
# if len(energies), then all energies <= cutoff
- if insert_ix == 0 or insert_ix == len(energies):
+ if cutoff_ix == 0 or cutoff_ix == len(energies):
tail = map("{:5.3e}".format, energies)
raise ValueError(
"Cannot find replacement fission yields for {} given "
"cutoff {:5.3e} eV. Yields provided at [{}] eV".format(
name, cutoff, ", ".join(tail)))
+ # find closest energy to requested thermal, fast energies
if thermal is None:
- thermal = yields[energies[insert_ix - 1]]
+ min_E = min(energies[:cutoff_ix],
+ key=lambda e: abs(e - thermal_energy))
+ thermal = yields[min_E]
if fast is None:
- fast = yields[energies[insert_ix]]
+ min_E = min(energies[cutoff_ix:],
+ key=lambda e: abs(e - fast_energy))
+ fast = yields[min_E]
self._thermal_yields[name] = thermal
self._fast_yields[name] = fast
@@ -369,7 +374,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
in parallel mode.
"""
super().generate_tallies(materials, mat_indexes)
- energy_filter = EnergyFilter(bins=[0.0, self._cutoff, self._upper_energy])
+ energy_filter = EnergyFilter([0.0, self._cutoff, self._upper_energy])
self._fission_rate_tally.filters = (
self._fission_rate_tally.filters + [energy_filter])
@@ -497,7 +502,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
fission_tally = self._fission_rate_tally
filters = fission_tally.filters
- ene_filter = EnergyFilter(bins=[0, self._upper_energy])
+ ene_filter = EnergyFilter([0, self._upper_energy])
fission_tally.filters = filters + [ene_filter]
func_filter = EnergyFunctionFilter()
diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py
index 81c3e2d42..05f62eace 100644
--- a/tests/unit_tests/test_deplete_fission_yields.py
+++ b/tests/unit_tests/test_deplete_fission_yields.py
@@ -1,17 +1,97 @@
"""Test the FissionYieldHelpers"""
+import os
from collections import namedtuple
from unittest.mock import Mock
+import bisect
import pytest
import numpy
+from openmc import capi
from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution
from openmc.deplete.helpers import (
FissionYieldCutoffHelper, ConstantFissionYieldHelper,
AveragedFissionYieldHelper)
-MATERIALS = ["1", "2"]
+@pytest.fixture(scope="module")
+def materials(tmpdir_factory):
+ """Use C API to construct realistic materials for testing tallies"""
+ tmpdir = tmpdir_factory.mktemp("capi")
+ orig = tmpdir.chdir()
+ # Create proxy xml files to please openmc
+ with open("geometry.xml", "w") as stream:
+ stream.write("""
+
+
+ |
+
+
+
+
+
+""")
+ with open("settings.xml", "w") as stream:
+ stream.write("""
+
+
+ eigenvalue
+ 100
+ 10
+ 0
+ 1
+
+
+ 0.0 0.0 0.0
+
+
+
+""")
+ with open("materials.xml", "w") as stream:
+ stream.write("""
+
+
+
+
+
+
+
+
+
+""")
+ try:
+ with capi.run_in_memory():
+ yield [capi.Material(), capi.Material()]
+ finally:
+ print(os.path.abspath(os.curdir))
+ os.remove(tmpdir / "settings.xml")
+ os.remove(tmpdir / "geometry.xml")
+ os.remove(tmpdir / "materials.xml")
+ os.remove(tmpdir / "summary.h5")
+ orig.chdir()
+ os.rmdir(tmpdir)
+
+
+def proxy_tally_data(tally, fill=None):
+ """Construct an empty matrix built from a C tally
+
+ The shape of tally.results will be
+ ``(n_bins, n_nuc * n_scores, 3)``
+ """
+ n_nucs = max(len(tally.nuclides), 1)
+ n_scores = max(len(tally.scores), 1)
+ n_bins = 1
+ for tfilter in tally.filters:
+ if not hasattr(tfilter, "bins"):
+ continue
+ this_bins = len(tfilter.bins)
+ if isinstance(tfilter, capi.EnergyFilter):
+ this_bins -= 1
+ n_bins *= max(this_bins, 1)
+ data = numpy.empty((n_bins, n_nucs * n_scores, 3))
+ if fill is not None:
+ data.fill(fill)
+ return data
@pytest.fixture(scope="module")
@@ -29,47 +109,74 @@ def nuclide_bundle():
xe135 = Nuclide("Xe135")
- NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135")
- return NuclideBundle(u235, u238, xe135)
+ pu239 = Nuclide("Pu239")
+ pu239.yield_data = FissionYieldDistribution({
+ 0.0253: {"Xe135": 3.141e-3, "Sm149": 8.19e-10, "Gd155": 1.66e-9},
+ 5.0e5: {"Xe135": 6.14e-3, "Sm149": 9.429e-10, "Gd155": 5.24e-9},
+ 2e6: {"Xe135": 6.15e-3, "Sm149": 9.42e-10, "Gd155": 5.29e-9}})
+
+ NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135 pu239")
+ return NuclideBundle(u235, u238, xe135, pu239)
@pytest.mark.parametrize(
- "input_energy, u5_yield_energy",
+ "input_energy, yield_energy",
((0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5)))
-def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy):
+def test_constant_helper(nuclide_bundle, input_energy, yield_energy):
helper = ConstantFissionYieldHelper(nuclide_bundle, energy=input_energy)
assert helper.energy == input_energy
assert helper.constant_yields == {
- "U235": nuclide_bundle.u235.yield_data[u5_yield_energy],
- "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal
+ "U235": nuclide_bundle.u235.yield_data[yield_energy],
+ "U238": nuclide_bundle.u238.yield_data[5.00e5], # only epithermal
+ "Pu239": nuclide_bundle.pu239.yield_data[yield_energy]}
assert helper.constant_yields == helper.weighted_yields(1)
def test_cutoff_construction(nuclide_bundle):
+ u235 = nuclide_bundle.u235
+ u238 = nuclide_bundle.u238
+ pu239 = nuclide_bundle.pu239
+
# defaults
helper = FissionYieldCutoffHelper(nuclide_bundle, 1)
assert helper.constant_yields == {
- "U238": nuclide_bundle.u238.yield_data[5.0e5]}
+ "U238": u238.yield_data[5.0e5]}
assert helper.thermal_yields == {
- "U235": nuclide_bundle.u235.yield_data[0.0253]}
- assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]}
+ "U235": u235.yield_data[0.0253],
+ "Pu239": pu239.yield_data[0.0253]}
+ assert helper.fast_yields == {
+ "U235": u235.yield_data[5e5],
+ "Pu239": pu239.yield_data[5e5]}
+
# use 14 MeV yields
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6)
assert helper.constant_yields == {
- "U238": nuclide_bundle.u238.yield_data[5.0e5]}
+ "U238": u238.yield_data[5.0e5]}
assert helper.thermal_yields == {
- "U235": nuclide_bundle.u235.yield_data[0.0253]}
- assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[14e6]}
+ "U235": u235.yield_data[0.0253],
+ "Pu239": pu239.yield_data[0.0253]}
+ assert helper.fast_yields == {
+ "U235": u235.yield_data[14e6],
+ "Pu239": pu239.yield_data[2e6]}
+
# specify missing thermal yields -> use 0.0253
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1)
assert helper.thermal_yields == {
- "U235": nuclide_bundle.u235.yield_data[0.0253]}
- assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]}
+ "U235": u235.yield_data[0.0253],
+ "Pu239": pu239.yield_data[0.0253]}
+ assert helper.fast_yields == {
+ "U235": u235.yield_data[5e5],
+ "Pu239": pu239.yield_data[5e5]}
+
# request missing fast yields -> use epithermal
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4)
assert helper.thermal_yields == {
- "U235": nuclide_bundle.u235.yield_data[0.0253]}
- assert helper.fast_yields == {"U235": nuclide_bundle.u235.yield_data[5e5]}
+ "U235": u235.yield_data[0.0253],
+ "Pu239": pu239.yield_data[0.0253]}
+ assert helper.fast_yields == {
+ "U235": u235.yield_data[5e5],
+ "Pu239": pu239.yield_data[5e5]}
+
# test failures in cutoff: super low, super high
with pytest.raises(ValueError, match="replacement fission yields"):
FissionYieldCutoffHelper(
@@ -87,83 +194,117 @@ def test_cutoff_failure(key):
FissionYieldCutoffHelper(None, None, **{key: -1})
-class ProxyMixin:
- """Mixing that overloads the tally generation"""
- def generate_tallies(self, materials, mat_indexes):
- self._fission_rate_tally = Mock()
- self._local_indexes = numpy.asarray(mat_indexes)
-
-
-class CutoffProxy(ProxyMixin, FissionYieldCutoffHelper):
- """Proxy that supplies a set of tallies"""
-
-
# emulate some split between fast and thermal U235 fissions
@pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8))
-def test_cutoff_helper(nuclide_bundle, therm_frac):
- n_bmats = len(MATERIALS)
- proxy = CutoffProxy(nuclide_bundle, n_bmats)
- proxy.generate_tallies(MATERIALS, [0])
+def test_cutoff_helper(materials, nuclide_bundle, therm_frac):
+ helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials))
+ helper.generate_tallies(materials, [0])
+
non_zero_nucs = [n.name for n in nuclide_bundle]
- tally_nucs = proxy.update_tally_nuclides(non_zero_nucs)
- assert tally_nucs == ("U235",)
+ tally_nucs = helper.update_tally_nuclides(non_zero_nucs)
+ assert tally_nucs == ("Pu239", "U235",)
+
+ # Check tallies
+ fission_tally = helper._fission_rate_tally
+ assert fission_tally is not None
+ filters = fission_tally.filters
+ assert len(filters) == 2
+ assert isinstance(filters[0], capi.MaterialFilter)
+ assert len(filters[0].bins) == len(materials)
+ assert isinstance(filters[1], capi.EnergyFilter)
+ # lower, cutoff, and upper energy
+ assert len(filters[1].bins) == 3
+
# Emulate building tallies
# material x energy, tallied_nuclides, 3
- proxy_flux = 1e6
- tally_data = numpy.empty((n_bmats * 2, 1, 3))
- tally_data[0, 0, 1] = therm_frac * proxy_flux
- tally_data[1, 0, 1] = (1 - therm_frac) * proxy_flux
- proxy._fission_rate_tally.results = tally_data
+ tally_data = proxy_tally_data(fission_tally)
+ helper._fission_rate_tally = Mock()
+ helper_flux = 1e6
+ tally_data[0, :, 1] = therm_frac * helper_flux
+ tally_data[1, :, 1] = (1 - therm_frac) * helper_flux
+ helper._fission_rate_tally.results = tally_data
- proxy.unpack()
+ helper.unpack()
# expected results of shape (n_mats, 2, n_tnucs)
- expected_results = numpy.empty((1, 2, 1))
+ expected_results = numpy.empty((1, 2, len(tally_nucs)))
expected_results[:, 0] = therm_frac
expected_results[:, 1] = 1 - therm_frac
- assert proxy.results == pytest.approx(expected_results)
+ assert helper.results == pytest.approx(expected_results)
- actual_yields = proxy.weighted_yields(0)
+ actual_yields = helper.weighted_yields(0)
assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5]
- assert actual_yields["U235"] == (
- proxy.thermal_yields["U235"] * therm_frac
- + proxy.fast_yields["U235"] * (1 - therm_frac))
+ for nuc in tally_nucs:
+ assert actual_yields[nuc] == (
+ helper.thermal_yields[nuc] * therm_frac
+ + helper.fast_yields[nuc] * (1 - therm_frac))
-class AverageProxy(ProxyMixin, AveragedFissionYieldHelper):
- """Proxy for generating mock set of tallies"""
- def generate_tallies(self, materials, mat_indexes):
- super().generate_tallies(materials, mat_indexes)
- self._weighted_tally = Mock()
-
-
-@pytest.mark.parametrize("avg_energy", (0.01, 100, 15e6))
-def test_averaged_helper(nuclide_bundle, avg_energy):
- proxy = AverageProxy(nuclide_bundle)
- proxy.generate_tallies(MATERIALS, [0])
- tallied_nucs = proxy.update_tally_nuclides(
+@pytest.mark.parametrize("avg_energy", (0.01, 6e5, 15e6))
+def test_averaged_helper(materials, nuclide_bundle, avg_energy):
+ helper = AveragedFissionYieldHelper(nuclide_bundle)
+ helper.generate_tallies(materials, [0])
+ tallied_nucs = helper.update_tally_nuclides(
[n.name for n in nuclide_bundle])
- assert tallied_nucs == ("U235", )
- # enforce some average energy
- proxy_flux = 1e16
- fission_results = numpy.ones((len(MATERIALS), 1, 3)) * proxy_flux
- weighted_results = fission_results * avg_energy
- proxy._fission_rate_tally.results = fission_results
- proxy._weighted_tally.results = weighted_results
- proxy.unpack()
- expected_results = numpy.ones((1, 1)) * avg_energy
- assert proxy.results == pytest.approx(expected_results)
+ assert tallied_nucs == ("Pu239", "U235")
- actual_yields = proxy.weighted_yields(0)
+ # check generated tallies
+ fission_tally = helper._fission_rate_tally
+ assert fission_tally is not None
+ fission_filters = fission_tally.filters
+ assert len(fission_filters) == 2
+ assert isinstance(fission_filters[0], capi.MaterialFilter)
+ assert len(fission_filters[0].bins) == len(materials)
+ assert isinstance(fission_filters[1], capi.EnergyFilter)
+ assert len(fission_filters[1].bins) == 2
+ assert fission_tally.scores == ["fission"]
+ assert fission_tally.nuclides == list(tallied_nucs)
+
+ weighted_tally = helper._weighted_tally
+ assert weighted_tally is not None
+ weighted_filters = weighted_tally.filters
+ assert len(weighted_filters) == 2
+ assert isinstance(weighted_filters[0], capi.MaterialFilter)
+ assert len(weighted_filters[0].bins) == len(materials)
+ assert isinstance(weighted_filters[1], capi.EnergyFunctionFilter)
+ assert len(weighted_filters[1].energy) == 2
+ assert len(weighted_filters[1].y) == 2
+ assert weighted_tally.scores == ["fission"]
+ assert weighted_tally.nuclides == list(tallied_nucs)
+
+ helper_flux = 1e16
+ fission_results = proxy_tally_data(fission_tally, helper_flux)
+ weighted_results = proxy_tally_data(
+ weighted_tally, helper_flux * avg_energy)
+
+ helper._fission_rate_tally = Mock()
+ helper._weighted_tally = Mock()
+ helper._fission_rate_tally.results = fission_results
+ helper._weighted_tally.results = weighted_results
+
+ helper.unpack()
+ expected_results = numpy.ones((1, len(tallied_nucs))) * avg_energy
+ assert helper.results == pytest.approx(expected_results)
+
+ actual_yields = helper.weighted_yields(0)
# constant U238 => no interpolation
assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5]
# construct expected yields
- if avg_energy < 0.0253: # take thermal U235 yields
- exp_u235_yields = nuclide_bundle.u235.yield_data[0.0253]
- elif avg_energy > 14e6: # take fastest U235 yields
- exp_u235_yields = nuclide_bundle.u235.yield_data[14e6]
- else: # reconstruct between thermal and epithermal
- thermal = nuclide_bundle.u235.yield_data[0.0253]
- epithermal = nuclide_bundle.u235.yield_data[5e5]
- split = (avg_energy - 0.0253) / (5e5 - 0.0253)
- exp_u235_yields = thermal * (1 - split) + epithermal * split
+ exp_u235_yields = interp_average_yields(nuclide_bundle.u235, avg_energy)
assert actual_yields["U235"] == exp_u235_yields
+ exp_pu239_yields = interp_average_yields(nuclide_bundle.pu239, avg_energy)
+ assert actual_yields["Pu239"] == exp_pu239_yields
+
+
+def interp_average_yields(nuc, avg_energy):
+ """Construct a set of yields by interpolation between neighbors"""
+ energies = nuc.yield_energies
+ yields = nuc.yield_data
+ if avg_energy < energies[0]:
+ return yields[energies[0]]
+ if avg_energy > energies[-1]:
+ return yields[energies[-1]]
+ thermal_ix = bisect.bisect_left(energies, avg_energy)
+ thermal_E, fast_E = energies[thermal_ix - 1:thermal_ix + 1]
+ assert thermal_E < avg_energy < fast_E
+ split = (avg_energy - thermal_E)/(fast_E - thermal_E)
+ return yields[thermal_E]*(1 - split) + yields[fast_E]*split