mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
291 lines
11 KiB
Python
291 lines
11 KiB
Python
"""Test the FissionYieldHelpers"""
|
|
|
|
import os
|
|
from collections import namedtuple
|
|
from unittest.mock import Mock
|
|
import bisect
|
|
|
|
import pytest
|
|
import numpy as np
|
|
import openmc
|
|
from openmc import lib
|
|
from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution
|
|
from openmc.deplete.helpers import (
|
|
FissionYieldCutoffHelper, ConstantFissionYieldHelper,
|
|
AveragedFissionYieldHelper)
|
|
|
|
|
|
@pytest.fixture(scope="module")
|
|
def materials(tmpdir_factory):
|
|
"""Use C API to construct realistic materials for testing tallies"""
|
|
tmpdir = tmpdir_factory.mktemp("lib")
|
|
orig = tmpdir.chdir()
|
|
# Create proxy problem to please openmc
|
|
mfuel = openmc.Material(name="test_fuel")
|
|
mfuel.volume = 1.0
|
|
for nuclide in ["U235", "U238", "Xe135", "Pu239"]:
|
|
mfuel.add_nuclide(nuclide, 1.0)
|
|
openmc.Materials([mfuel]).export_to_xml()
|
|
# Geometry
|
|
box = openmc.model.RectangularPrism(1.0, 1.0, boundary_type="reflective")
|
|
cell = openmc.Cell(fill=mfuel, region=-box)
|
|
root = openmc.Universe(cells=[cell])
|
|
openmc.Geometry(root).export_to_xml()
|
|
# settings
|
|
settings = openmc.Settings()
|
|
settings.particles = 100
|
|
settings.inactive = 0
|
|
settings.batches = 10
|
|
settings.verbosity = 1
|
|
settings.export_to_xml()
|
|
|
|
try:
|
|
with lib.run_in_memory():
|
|
yield [lib.Material(), lib.Material()]
|
|
finally:
|
|
for file_path in ("settings.xml", "geometry.xml", "materials.xml",
|
|
"summary.h5"):
|
|
os.remove(tmpdir / file_path)
|
|
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.mean will be ``(n_bins, n_nuc * n_scores)``
|
|
"""
|
|
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, lib.EnergyFilter):
|
|
this_bins -= 1
|
|
n_bins *= max(this_bins, 1)
|
|
data = np.empty((n_bins, n_nucs * n_scores))
|
|
if fill is not None:
|
|
data.fill(fill)
|
|
return data
|
|
|
|
|
|
@pytest.fixture(scope="module")
|
|
def nuclide_bundle():
|
|
u5yield_dict = {
|
|
0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12},
|
|
5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12},
|
|
1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}}
|
|
u235 = Nuclide("U235")
|
|
u235.yield_data = FissionYieldDistribution(u5yield_dict)
|
|
|
|
u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}}
|
|
u238 = Nuclide("U238")
|
|
u238.yield_data = FissionYieldDistribution(u8yield_dict)
|
|
|
|
xe135 = Nuclide("Xe135")
|
|
|
|
pu239 = Nuclide("Pu239")
|
|
pu239.yield_data = FissionYieldDistribution({
|
|
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, yield_energy",
|
|
((0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5)))
|
|
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[yield_energy],
|
|
"U238": nuclide_bundle.u238.yield_data[5.00e5],
|
|
"Pu239": nuclide_bundle.pu239.yield_data[5e5]}
|
|
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": u238.yield_data[5.0e5],
|
|
"Pu239": pu239.yield_data[5e5]}
|
|
assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]}
|
|
assert helper.fast_yields == {"U235": u235.yield_data[5e5]}
|
|
|
|
# use 14 MeV yields
|
|
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6)
|
|
assert helper.constant_yields == {
|
|
"U238": u238.yield_data[5.0e5],
|
|
"Pu239": pu239.yield_data[5e5]}
|
|
assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]}
|
|
assert helper.fast_yields == {"U235": u235.yield_data[14e6]}
|
|
|
|
# specify missing thermal yields -> use 0.0253
|
|
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1)
|
|
assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]}
|
|
assert helper.fast_yields == {"U235": u235.yield_data[5e5]}
|
|
|
|
# request missing fast yields -> use epithermal
|
|
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4)
|
|
assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]}
|
|
assert helper.fast_yields == {"U235": u235.yield_data[5e5]}
|
|
|
|
# higher cutoff energy -> obtain fast and "faster" yields
|
|
helper = FissionYieldCutoffHelper(nuclide_bundle, 1, cutoff=1e6,
|
|
thermal_energy=5e5, fast_energy=14e6)
|
|
assert helper.constant_yields == {"U238": u238.yield_data[5e5]}
|
|
assert helper.thermal_yields == {
|
|
"U235": u235.yield_data[5e5], "Pu239": pu239.yield_data[5e5]}
|
|
assert helper.fast_yields == {
|
|
"U235": u235.yield_data[14e6], "Pu239": pu239.yield_data[2e6]}
|
|
|
|
# test super low and super high cutoff energies
|
|
helper = FissionYieldCutoffHelper(
|
|
nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002)
|
|
assert helper.fast_yields == {}
|
|
assert helper.thermal_yields == {}
|
|
assert helper.constant_yields == {
|
|
"U235": u235.yield_data[0.0253], "U238": u238.yield_data[5e5],
|
|
"Pu239": pu239.yield_data[5e5]}
|
|
|
|
helper = FissionYieldCutoffHelper(
|
|
nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6)
|
|
assert helper.thermal_yields == {}
|
|
assert helper.fast_yields == {}
|
|
assert helper.constant_yields == {
|
|
"U235": u235.yield_data[14e6], "U238": u238.yield_data[5e5],
|
|
"Pu239": pu239.yield_data[2e6]}
|
|
|
|
|
|
@pytest.mark.parametrize("key", ("cutoff", "thermal_energy", "fast_energy"))
|
|
def test_cutoff_failure(key):
|
|
with pytest.raises(TypeError, match=key):
|
|
FissionYieldCutoffHelper(None, None, **{key: None})
|
|
with pytest.raises(ValueError, match=key):
|
|
FissionYieldCutoffHelper(None, None, **{key: -1})
|
|
|
|
|
|
# emulate some split between fast and thermal U235 fissions
|
|
@pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8))
|
|
def test_cutoff_helper(materials, nuclide_bundle, therm_frac):
|
|
helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials),
|
|
cutoff=1e6, fast_energy=14e6)
|
|
helper.generate_tallies(materials, [0])
|
|
|
|
non_zero_nucs = [n.name for n in nuclide_bundle]
|
|
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], lib.MaterialFilter)
|
|
assert len(filters[0].bins) == len(materials)
|
|
assert isinstance(filters[1], lib.EnergyFilter)
|
|
# lower, cutoff, and upper energy
|
|
assert len(filters[1].bins) == 3
|
|
|
|
# Emulate building tallies
|
|
# material x energy, tallied_nuclides, 3
|
|
tally_data = proxy_tally_data(fission_tally)
|
|
helper._fission_rate_tally = Mock()
|
|
helper_flux = 1e6
|
|
tally_data[0] = therm_frac * helper_flux
|
|
tally_data[1] = (1 - therm_frac) * helper_flux
|
|
helper._fission_rate_tally.mean = tally_data
|
|
|
|
helper.unpack()
|
|
# expected results of shape (n_mats, 2, n_tnucs)
|
|
expected_results = np.empty((1, 2, len(tally_nucs)))
|
|
expected_results[:, 0] = therm_frac
|
|
expected_results[:, 1] = 1 - therm_frac
|
|
assert helper.results == pytest.approx(expected_results)
|
|
|
|
actual_yields = helper.weighted_yields(0)
|
|
assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5]
|
|
for nuc in tally_nucs:
|
|
assert actual_yields[nuc] == (
|
|
helper.thermal_yields[nuc] * therm_frac
|
|
+ helper.fast_yields[nuc] * (1 - therm_frac))
|
|
|
|
|
|
@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 == ["Pu239", "U235"]
|
|
|
|
# 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], lib.MaterialFilter)
|
|
assert len(fission_filters[0].bins) == len(materials)
|
|
assert isinstance(fission_filters[1], lib.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], lib.MaterialFilter)
|
|
assert len(weighted_filters[0].bins) == len(materials)
|
|
assert isinstance(weighted_filters[1], lib.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.mean = fission_results
|
|
helper._weighted_tally.mean = weighted_results
|
|
|
|
helper.unpack()
|
|
expected_results = np.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
|
|
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
|