mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 20:45:35 -04:00
Improve testing for FY helpers using C API
Provide a module-scoped fixture that creates a temporary directory, adds simple settings, geometry, and material files, and initializes the openmc C API. This fixture also provides two empty openmc.capi.Material objects to help the helpers generate meaningful tallies. The main tests for AveragedFissionYieldHelper and FissionYieldCutoffHelper now go through the built tallies and examine the filters, nuclides, and scores. A helper function produces mocked-like tally data based on filters, nuclides, and scores found on a tally. Pu239 with 0.0253 eV, 500 keV, and 2 MeV yields is included in the nuclide_bundle fixture. This provides some additional heterogeneity in the results, as the 2 MeV data is not present on U235.
This commit is contained in:
Andrew Johnson
parent
687d287b24
commit
8e2bb0bbbc
2 changed files with 230 additions and 84 deletions
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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("""
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<geometry>
|
||||
<cell id="1" material="1" name="fuel" region="1 -2 3 -4" universe="1" />
|
||||
<surface boundary="reflective" coeffs="-0.63" id="1" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="0.63" id="2" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.63" id="3" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="0.63" id="4" type="y-plane" />
|
||||
</geometry>
|
||||
""")
|
||||
with open("settings.xml", "w") as stream:
|
||||
stream.write("""
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<settings>
|
||||
<run_mode>eigenvalue</run_mode>
|
||||
<particles>100</particles>
|
||||
<batches>10</batches>
|
||||
<inactive>0</inactive>
|
||||
<verbosity>1</verbosity>
|
||||
<source strength="1.0">
|
||||
<space type="point">
|
||||
<parameters>0.0 0.0 0.0</parameters>
|
||||
</space>
|
||||
</source>
|
||||
</settings>
|
||||
""")
|
||||
with open("materials.xml", "w") as stream:
|
||||
stream.write("""
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
<material depletable="true" id="1" name="U" volume="71.67537585">
|
||||
<density units="g/cc" value="10.4" />
|
||||
<nuclide ao="1.0" name="U235" />
|
||||
<nuclide ao="1.0" name="U238" />
|
||||
<nuclide ao="1.0" name="Xe135" />
|
||||
<nuclide ao="1.0" name="Pu239" />
|
||||
</material>
|
||||
</materials>""")
|
||||
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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue