mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Add FissionYieldHelper for energy-dependent fission yields
Operator now has a new helper for tallying fission rates
and computing fission yields using energy dependency.
The FissionYieldHelper creates fission rate tallies across
burnable materials with an energy filter. The energy filter
corresponds to each energy group for which there are fission
yields.
During the operator's unpacking of tallies, this helper
computes a relative fission contribution by normalizing
fission rates, such that the sum of all group fission rates
in each material and each nuclide is unity.
A single fission yield for each parent nuclide is computed
by weighting the energy-dependent fission yields by
these relative contributions.
A nested dictionary {parent: {product: yield}} is added to
a library, and then attached to the depletion chain
at the end of the unpacking. These libraries are used
during the depletion process instead of the only-thermal
yields from before.
A test was added that creates a proxy helper that doesn't
use the C-API, but unpacks the "fission tally" and
computes the fission yield library in an identical manner
as the actual class. Since the proxy helper is a subclass
of the one used in real transport, the functionality
is consistent.
This commit is contained in:
Andrew Johnson
parent
8cd71be864
commit
599d473888
4 changed files with 309 additions and 8 deletions
|
|
@ -452,9 +452,6 @@ class Chain(object):
|
|||
k = self.nuclide_dict[target]
|
||||
matrix[k, i] += path_rate * br
|
||||
else:
|
||||
# Assume that we should always use thermal fission
|
||||
# yields. At some point it would be nice to account
|
||||
# for the energy-dependence..
|
||||
for product, y in fission_yields[nuc.name].items():
|
||||
yield_val = y * path_rate
|
||||
if yield_val != 0.0:
|
||||
|
|
|
|||
|
|
@ -1,11 +1,13 @@
|
|||
"""
|
||||
Class for normalizing fission energy deposition
|
||||
"""
|
||||
from collections import defaultdict
|
||||
from itertools import product
|
||||
|
||||
from numpy import dot, zeros
|
||||
from numpy import dot, zeros, newaxis, divide, asarray
|
||||
|
||||
from openmc.capi import Tally, MaterialFilter
|
||||
from openmc.capi import (
|
||||
Tally, MaterialFilter, EnergyFilter)
|
||||
from .abc import ReactionRateHelper, EnergyHelper
|
||||
|
||||
# -------------------------------------
|
||||
|
|
@ -142,3 +144,190 @@ class ChainFissionHelper(EnergyHelper):
|
|||
isotopes in all materials have the same Q value.
|
||||
"""
|
||||
self._energy += dot(fission_rates, self._fission_q_vector)
|
||||
|
||||
|
||||
# ------------------------------------
|
||||
# Helper for collapsing fission yields
|
||||
# ------------------------------------
|
||||
|
||||
|
||||
class FissionYieldHelper(object):
|
||||
"""Class for using energy-dependent fission yields in depletion chain
|
||||
|
||||
Creates a tally across all burnable materials to score the fission
|
||||
rate in nuclides with yield data. An energy filter is used to
|
||||
compute this rates in a group structure corresponding to the
|
||||
fission yield data. This tally data is used to compute the
|
||||
relative number of fission events in each energy region,
|
||||
which serve as the weights for each energy-dependent fission
|
||||
yield distribution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.nuclide
|
||||
Nuclides tracked in the depletion chain. Not necessary
|
||||
that all have yield data.
|
||||
n_bmats : int
|
||||
Number of burnable materials tracked in the problem
|
||||
|
||||
Attributes
|
||||
----------
|
||||
energy_bounds : tuple of float
|
||||
Sorted energy bounds from the tally filter
|
||||
results : numpy.ndarray
|
||||
Array of tally results for this process with shape
|
||||
``(n_local_mat, n_energy, n_nucs)``
|
||||
libraries : list of dict
|
||||
List of fission yield dictionaries of the form
|
||||
``{parent: {product: yield}}``. Populated in
|
||||
:meth:`compute_yields` and reset during
|
||||
:meth:`unpack`.
|
||||
"""
|
||||
|
||||
def __init__(self, chain_nuclides, n_bmats):
|
||||
self.libraries = []
|
||||
self._chain_nuclides = {}
|
||||
self._chain_set = set()
|
||||
self._tally_map = {}
|
||||
# TODO Support user-requested minimum energy?
|
||||
yield_energies = {0.0}
|
||||
|
||||
# Get nuclides with fission yield data, names
|
||||
# and all energy points
|
||||
# Names are provided from operator tally nuclides
|
||||
for nuc in chain_nuclides:
|
||||
if len(nuc.yield_data) == 0:
|
||||
continue
|
||||
self._chain_nuclides[nuc.name] = nuc
|
||||
self._chain_set.add(nuc.name)
|
||||
yield_energies.update(nuc.yield_energies)
|
||||
|
||||
# Create energy grid
|
||||
self._energy_bounds = tuple(sorted(yield_energies))
|
||||
self.n_bmats = n_bmats
|
||||
|
||||
self._reaction_tally = None
|
||||
self.results = None
|
||||
self.local_indexes = None
|
||||
|
||||
@property
|
||||
def energy_bounds(self):
|
||||
return self._energy_bounds
|
||||
|
||||
def generate_tallies(self, materials, mat_indexes):
|
||||
"""Construct the fission rate tally
|
||||
|
||||
Parameters
|
||||
----------
|
||||
materials : iterable of C-API materials
|
||||
Materials to be used in :class:`openmc.capi.MaterialFilter`
|
||||
mat_indexes : iterable of int
|
||||
Indexes for materials in ``materials`` tracked on this
|
||||
process
|
||||
"""
|
||||
# Tally group-wise fission reaction rates
|
||||
self._reaction_tally = Tally()
|
||||
self._reaction_tally.scores = ['fission']
|
||||
|
||||
# Tally energy-weighted group-wise fission reaction rate
|
||||
# Used to evaluated linear interpolation between fission yield points
|
||||
self._weighted_reaction_tally = Tally()
|
||||
self._weighted_reaction_tally.scores = ['fission']
|
||||
|
||||
filters = [
|
||||
MaterialFilter(materials), EnergyFilter(self._energy_bounds)]
|
||||
|
||||
self._reaction_tally.filters = filters
|
||||
self.local_indexes = asarray(mat_indexes)
|
||||
|
||||
def set_fissionable_nuclides(self, nuclides):
|
||||
"""List of string of nuclides with data to be tallied
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : iterable of str
|
||||
nuclides with non-zero densities that are candidates
|
||||
for the fission tally. Not necessary that all are nuclides
|
||||
with fission yields, but at least one must be
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : tuple of str
|
||||
Nuclides ordered as they appear in the tally and in
|
||||
the nuclide column of :attr:`results`
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If no nuclides in ``nuclides`` are tracked on this
|
||||
object
|
||||
"""
|
||||
# Set of all nuclides with positive density
|
||||
# and fission yield data
|
||||
nuc_set = self._chain_set & set(nuclides)
|
||||
if len(nuc_set) == 0:
|
||||
raise ValueError(
|
||||
"No overlap between chain nuclides with fission yields and "
|
||||
"requested tally nuclides")
|
||||
nuclides = tuple(sorted(nuc_set))
|
||||
self._tally_index = [self._chain_nuclides[n] for n in nuclides]
|
||||
self._reaction_tally.nuclides = nuclides
|
||||
return nuclides
|
||||
|
||||
def unpack(self):
|
||||
"""Unpack fission rate tallies to produce :attr:`results`
|
||||
|
||||
Resets :attr:`libraries` under the assumption this is called
|
||||
during the :class:`openmc.deplete.Operator` unpackign process
|
||||
"""
|
||||
# clear old libraries
|
||||
self.libraries = []
|
||||
|
||||
# get view into tally results
|
||||
# new shape: [material, energy, parent nuclide]
|
||||
result_view = self._reaction_tally.results[..., 1].reshape(
|
||||
self.n_bmats, len(self._energy_bounds) - 1,
|
||||
len(self._reaction_tally.nuclides))
|
||||
|
||||
# Get results specific to this process
|
||||
self.results = result_view[self.local_indexes, ...]
|
||||
|
||||
# scale fission yields proportional to total fission rate
|
||||
fsn_rate = self.results.sum(axis=1)
|
||||
# TODO Guard against divide by zero
|
||||
self.results /= fsn_rate[:, newaxis, :]
|
||||
|
||||
def compute_yields(self, local_mat_index):
|
||||
"""Compute single fission yields using :attr:`results`
|
||||
|
||||
Produces a new library in :attr:`self.libraries`
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mat_index : int
|
||||
Index for material tracked on this process that
|
||||
exists in :attr:`local_mat_index` and fits within
|
||||
the first axis in :attr:`results`
|
||||
"""
|
||||
tally_results = self.results[local_mat_index]
|
||||
|
||||
# Dictionary {parent_nuclide : [product, yield_vector]}
|
||||
initial_library = {}
|
||||
for i_energy, energy in enumerate(self._energy_bounds[1:]):
|
||||
for i_nuc, fiss_frac in enumerate(tally_results[i_energy]):
|
||||
parent = self._tally_index[i_nuc]
|
||||
yield_data = parent.yield_data.get(energy)
|
||||
if yield_data is None:
|
||||
continue
|
||||
if parent not in initial_library:
|
||||
initial_library[parent] = yield_data * fiss_frac
|
||||
continue
|
||||
initial_library[parent] += yield_data * fiss_frac
|
||||
|
||||
# convert to dictionary that can be passed to Chain.form_matrix
|
||||
# {parent: {product: yield}}
|
||||
library = {}
|
||||
for k, yield_obj in initial_library.items():
|
||||
library[k.name] = dict(zip(yield_obj.products, yield_obj.yields))
|
||||
|
||||
self.libraries.append(library)
|
||||
|
|
|
|||
|
|
@ -25,7 +25,8 @@ from .abc import TransportOperator, OperatorResult
|
|||
from .atom_number import AtomNumber
|
||||
from .reaction_rates import ReactionRates
|
||||
from .results_list import ResultsList
|
||||
from .helpers import DirectReactionRateHelper, ChainFissionHelper
|
||||
from .helpers import (
|
||||
DirectReactionRateHelper, ChainFissionHelper, FissionYieldHelper)
|
||||
|
||||
|
||||
def _distribute(items):
|
||||
|
|
@ -177,6 +178,8 @@ class Operator(TransportOperator):
|
|||
self._rate_helper = DirectReactionRateHelper(
|
||||
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
|
||||
self._energy_helper = ChainFissionHelper()
|
||||
self._fsn_yield_helper = FissionYieldHelper(
|
||||
self.chain.nuclides, len(self.burnable_mats))
|
||||
|
||||
def __call__(self, vec, power):
|
||||
"""Runs a simulation.
|
||||
|
|
@ -204,8 +207,10 @@ class Operator(TransportOperator):
|
|||
|
||||
# Update material compositions and tally nuclides
|
||||
self._update_materials()
|
||||
self._rate_helper.nuclides = self._get_tally_nuclides()
|
||||
self._energy_helper.nuclides = self._rate_helper.nuclides
|
||||
nuclides = self._get_tally_nuclides()
|
||||
self._rate_helper.nuclides = nuclides
|
||||
self._energy_helper.nuclides = nuclides
|
||||
self._fsn_yield_helper.set_fissionable_nuclides(nuclides)
|
||||
|
||||
# Run OpenMC
|
||||
openmc.capi.reset()
|
||||
|
|
@ -409,6 +414,10 @@ class Operator(TransportOperator):
|
|||
self._rate_helper.generate_tallies(materials, self.chain.reactions)
|
||||
self._energy_helper.prepare(
|
||||
self.chain.nuclides, self.reaction_rates.index_nuc, materials)
|
||||
# Tell fission yield helper what materials this process is
|
||||
# responsible for
|
||||
self._fsn_yield_helper.generate_tallies(
|
||||
materials, tuple(sorted(self._mat_index_map.values())))
|
||||
|
||||
# Return number density vector
|
||||
return list(self.number.get_mat_slice(np.s_[:]))
|
||||
|
|
@ -554,6 +563,7 @@ class Operator(TransportOperator):
|
|||
# Keep track of energy produced from all reactions in eV per source
|
||||
# particle
|
||||
self._energy_helper.reset()
|
||||
self._fsn_yield_helper.unpack()
|
||||
|
||||
# Create arrays to store fission Q values, reaction rates, and nuclide
|
||||
# numbers, zeroed out in material iteration
|
||||
|
|
@ -576,6 +586,9 @@ class Operator(TransportOperator):
|
|||
tally_rates = self._rate_helper.get_material_rates(
|
||||
mat_index, nuc_ind, react_ind)
|
||||
|
||||
# Compute fission yields for this material
|
||||
self._fsn_yield_helper.compute_yields(i)
|
||||
|
||||
# Accumulate energy from fission
|
||||
self._energy_helper.update(tally_rates[:, fission_ind], mat_index)
|
||||
|
||||
|
|
@ -589,6 +602,9 @@ class Operator(TransportOperator):
|
|||
# Scale reaction rates to obtain units of reactions/sec
|
||||
rates *= power / energy
|
||||
|
||||
# Store new fission yields on the chain
|
||||
self.chain.fission_yields = self._fsn_yield_helper.libraries
|
||||
|
||||
return OperatorResult(k_combined, rates)
|
||||
|
||||
def _get_nuclides_with_data(self):
|
||||
|
|
|
|||
99
tests/unit_tests/test_deplete_helpers.py
Normal file
99
tests/unit_tests/test_deplete_helpers.py
Normal file
|
|
@ -0,0 +1,99 @@
|
|||
"""Test the Operator helpers"""
|
||||
|
||||
from unittest.mock import Mock
|
||||
|
||||
import pytest
|
||||
import numpy
|
||||
from numpy.testing import assert_array_equal
|
||||
|
||||
from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution
|
||||
from openmc.deplete.helpers import FissionYieldHelper
|
||||
|
||||
|
||||
class FissionYieldHelperProxy(FissionYieldHelper):
|
||||
|
||||
def __init__(self, chain_nuclides, n_bmats):
|
||||
super().__init__(chain_nuclides, n_bmats)
|
||||
self._reaction_tally = Mock()
|
||||
self.local_indexes = numpy.array([0])
|
||||
|
||||
def generate_tallies(self, *args, **kwargs):
|
||||
# Avoid calls to the C-API
|
||||
pass
|
||||
|
||||
|
||||
def test_fission_yield_helper():
|
||||
"""Test the collection of fission yield data using approximated tallies
|
||||
"""
|
||||
u5yield_dict = {
|
||||
0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12},
|
||||
1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}}
|
||||
u235 = Nuclide()
|
||||
u235.name = "U235"
|
||||
u235.yield_energies = (0.0253, 1.40e7)
|
||||
u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict)
|
||||
|
||||
u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}}
|
||||
u238 = Nuclide()
|
||||
u238.name = "U238"
|
||||
u238.yield_energies = (5.00e5, )
|
||||
u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict)
|
||||
|
||||
xe135 = Nuclide()
|
||||
xe135.name = "Xe135"
|
||||
|
||||
n_bmats = 2
|
||||
helper = FissionYieldHelperProxy([u235, u238, xe135], n_bmats)
|
||||
|
||||
assert helper.energy_bounds == (0, 0.0253, 5.00e5, 1.40e7)
|
||||
|
||||
# test that tally must be created with nuclides with yields
|
||||
with pytest.raises(ValueError, match="No overlap"):
|
||||
helper.set_fissionable_nuclides(["Xe135", ])
|
||||
|
||||
act_nucs = helper.set_fissionable_nuclides(["U235", "U238", "Xe135"])
|
||||
assert act_nucs == ("U235", "U238")
|
||||
|
||||
# Emulate getting tally data from transport run
|
||||
# Test as if this Helper is responsible for one of two materials
|
||||
# Tally results ordered [n_mat * n_ene, n_fiss_nuc, 3]
|
||||
|
||||
tally_res = numpy.zeros((3 * n_bmats, 2, 3))
|
||||
u5_fiss_rates = numpy.array([1.0, 1.5, 2.0])
|
||||
u8_fiss_rates = numpy.array([0.0, 1.0, 1.0])
|
||||
tally_res[:3, 0, 1] = u5_fiss_rates
|
||||
tally_res[:3, 1, 1] = u8_fiss_rates
|
||||
|
||||
helper._reaction_tally.results = tally_res
|
||||
helper.unpack() # compute yield fractions
|
||||
|
||||
# Compare fraction fission rate from helper.reset
|
||||
exp_results = numpy.empty((1, 3, 2))
|
||||
# Fraction of fission events in each energy range
|
||||
u5_vec = u5_fiss_rates / u5_fiss_rates.sum()
|
||||
u8_vec = u8_fiss_rates / u8_fiss_rates.sum()
|
||||
exp_results[..., 0] = u5_vec
|
||||
exp_results[..., 1] = u8_vec
|
||||
assert_array_equal(helper.results, exp_results)
|
||||
|
||||
# Compute and compare the library of fission yields
|
||||
exp_lib = {
|
||||
"U235": {
|
||||
"Xe135": (u5yield_dict[0.0253]["Xe135"] * u5_vec[0]
|
||||
+ u5yield_dict[1.40e7]["Xe135"] * u5_vec[2]),
|
||||
"Gd155": (u5yield_dict[0.0253]["Gd155"] * u5_vec[0]
|
||||
+ u5yield_dict[1.40e7]["Gd155"] * u5_vec[2]),
|
||||
"Sm149": u5yield_dict[0.0253]["Sm149"] * u5_vec[0],
|
||||
},
|
||||
"U238": {
|
||||
"Xe135": u8yield_dict[5.00e5]["Xe135"] * u8_vec[1],
|
||||
"Gd155": u8yield_dict[5.00e5]["Gd155"] * u8_vec[1],
|
||||
}
|
||||
}
|
||||
|
||||
assert len(helper.libraries) == 0
|
||||
helper.compute_yields(0)
|
||||
assert len(helper.libraries) == 1
|
||||
act_library = helper.libraries[0]
|
||||
for parent, sub_yields in exp_lib.items():
|
||||
assert act_library[parent] == pytest.approx(sub_yields)
|
||||
Loading…
Add table
Add a link
Reference in a new issue