mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #1313 from drewejohnson/chain-fission-yields
Use energy-dependent fission yields when constructing depletion matrix
This commit is contained in:
Paul Romano
GitHub
commit
d4782d5719
12 changed files with 1621 additions and 166 deletions
|
|
@ -6,7 +6,7 @@
|
|||
|
||||
.. module:: openmc.deplete
|
||||
|
||||
Several functions are provided that implement different time-integration
|
||||
Several classes are provided that implement different time-integration
|
||||
algorithms for depletion calculations, which are described in detail in Colin
|
||||
Josey's thesis, `Development and analysis of high order neutron
|
||||
transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
|
||||
|
|
@ -25,7 +25,7 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
|
|||
SICELIIntegrator
|
||||
SILEQIIntegrator
|
||||
|
||||
Each of these functions expects a "transport operator" to be passed. An operator
|
||||
Each of these classes expects a "transport operator" to be passed. An operator
|
||||
specific to OpenMC is available using the following class:
|
||||
|
||||
.. autosummary::
|
||||
|
|
@ -65,6 +65,8 @@ for a depletion chain:
|
|||
DecayTuple
|
||||
Nuclide
|
||||
ReactionTuple
|
||||
FissionYieldDistribution
|
||||
FissionYield
|
||||
|
||||
The following classes are used during a depletion simulation and store auxiliary
|
||||
data, such as number densities and reaction rates for each material.
|
||||
|
|
@ -75,13 +77,25 @@ data, such as number densities and reaction rates for each material.
|
|||
:template: myclass.rst
|
||||
|
||||
AtomNumber
|
||||
ChainFissionHelper
|
||||
DirectReactionRateHelper
|
||||
OperatorResult
|
||||
ReactionRates
|
||||
Results
|
||||
ResultsList
|
||||
|
||||
The following classes are used to help the :class:`openmc.deplete.Operator`
|
||||
compute quantities like effective fission yields, reaction rates, and
|
||||
total system energy.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
helpers.AveragedFissionYieldHelper
|
||||
helpers.ChainFissionHelper
|
||||
helpers.ConstantFissionYieldHelper
|
||||
helpers.DirectReactionRateHelper
|
||||
helpers.FissionYieldCutoffHelper
|
||||
|
||||
The following classes are abstract classes that can be used to extend the
|
||||
:mod:`openmc.deplete` capabilities:
|
||||
|
|
@ -91,8 +105,10 @@ The following classes are abstract classes that can be used to extend the
|
|||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
ReactionRateHelper
|
||||
EnergyHelper
|
||||
FissionYieldHelper
|
||||
ReactionRateHelper
|
||||
TalliedFissionYieldHelper
|
||||
TransportOperator
|
||||
|
||||
Custom integrators can be developed by subclassing from the following abstract
|
||||
|
|
|
|||
|
|
@ -1,77 +1,78 @@
|
|||
"""Functions to form the special matrix for depletion"""
|
||||
|
||||
|
||||
def celi_f1(chain, rates):
|
||||
return (5 / 12 * chain.form_matrix(rates[0])
|
||||
+ 1 / 12 * chain.form_matrix(rates[1]))
|
||||
def celi_f1(chain, rates, fission_yields=None):
|
||||
return (5 / 12 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ 1 / 12 * chain.form_matrix(rates[1], fission_yields))
|
||||
|
||||
|
||||
def celi_f2(chain, rates):
|
||||
return (1 / 12 * chain.form_matrix(rates[0])
|
||||
+ 5 / 12 * chain.form_matrix(rates[1]))
|
||||
def celi_f2(chain, rates, fission_yields=None):
|
||||
return (1 / 12 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ 5 / 12 * chain.form_matrix(rates[1], fission_yields))
|
||||
|
||||
|
||||
def cf4_f1(chain, rates):
|
||||
return 1 / 2 * chain.form_matrix(rates)
|
||||
def cf4_f1(chain, rates, fission_yields=None):
|
||||
return 1 / 2 * chain.form_matrix(rates, fission_yields)
|
||||
|
||||
|
||||
def cf4_f2(chain, rates):
|
||||
return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1])
|
||||
def cf4_f2(chain, rates, fission_yields=None):
|
||||
return (-1 / 2 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ chain.form_matrix(rates[1], fission_yields))
|
||||
|
||||
|
||||
def cf4_f3(chain, rates):
|
||||
return (1 / 4 * chain.form_matrix(rates[0])
|
||||
+ 1 / 6 * chain.form_matrix(rates[1])
|
||||
+ 1 / 6 * chain.form_matrix(rates[2])
|
||||
- 1 / 12 * chain.form_matrix(rates[3]))
|
||||
def cf4_f3(chain, rates, fission_yields=None):
|
||||
return (1 / 4 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ 1 / 6 * chain.form_matrix(rates[1], fission_yields)
|
||||
+ 1 / 6 * chain.form_matrix(rates[2], fission_yields)
|
||||
- 1 / 12 * chain.form_matrix(rates[3], fission_yields))
|
||||
|
||||
|
||||
def cf4_f4(chain, rates):
|
||||
return (-1 / 12 * chain.form_matrix(rates[0])
|
||||
+ 1 / 6 * chain.form_matrix(rates[1])
|
||||
+ 1 / 6 * chain.form_matrix(rates[2])
|
||||
+ 1 / 4 * chain.form_matrix(rates[3]))
|
||||
def cf4_f4(chain, rates, fission_yields=None):
|
||||
return (-1 / 12 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ 1 / 6 * chain.form_matrix(rates[1], fission_yields)
|
||||
+ 1 / 6 * chain.form_matrix(rates[2], fission_yields)
|
||||
+ 1 / 4 * chain.form_matrix(rates[3], fission_yields))
|
||||
|
||||
|
||||
def rk4_f1(chain, rates):
|
||||
return 1 / 2 * chain.form_matrix(rates)
|
||||
def rk4_f1(chain, rates, fission_yields=None):
|
||||
return 1 / 2 * chain.form_matrix(rates, fission_yields)
|
||||
|
||||
|
||||
def rk4_f4(chain, rates):
|
||||
return (1 / 6 * chain.form_matrix(rates[0])
|
||||
+ 1 / 3 * chain.form_matrix(rates[1])
|
||||
+ 1 / 3 * chain.form_matrix(rates[2])
|
||||
+ 1 / 6 * chain.form_matrix(rates[3]))
|
||||
def rk4_f4(chain, rates, fission_yields=None):
|
||||
return (1 / 6 * chain.form_matrix(rates[0], fission_yields)
|
||||
+ 1 / 3 * chain.form_matrix(rates[1], fission_yields)
|
||||
+ 1 / 3 * chain.form_matrix(rates[2], fission_yields)
|
||||
+ 1 / 6 * chain.form_matrix(rates[3], fission_yields))
|
||||
|
||||
|
||||
def leqi_f1(chain, inputs):
|
||||
f1 = chain.form_matrix(inputs[0])
|
||||
f2 = chain.form_matrix(inputs[1])
|
||||
def leqi_f1(chain, inputs, fission_yields):
|
||||
f1 = chain.form_matrix(inputs[0], fission_yields)
|
||||
f2 = chain.form_matrix(inputs[1], fission_yields)
|
||||
dt_l, dt = inputs[2], inputs[3]
|
||||
return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2
|
||||
|
||||
|
||||
def leqi_f2(chain, inputs):
|
||||
f1 = chain.form_matrix(inputs[0])
|
||||
f2 = chain.form_matrix(inputs[1])
|
||||
def leqi_f2(chain, inputs, fission_yields=None):
|
||||
f1 = chain.form_matrix(inputs[0], fission_yields)
|
||||
f2 = chain.form_matrix(inputs[1], fission_yields)
|
||||
dt_l, dt = inputs[2], inputs[3]
|
||||
return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2
|
||||
|
||||
|
||||
def leqi_f3(chain, inputs):
|
||||
f1 = chain.form_matrix(inputs[0])
|
||||
f2 = chain.form_matrix(inputs[1])
|
||||
f3 = chain.form_matrix(inputs[2])
|
||||
def leqi_f3(chain, inputs, fission_yields=None):
|
||||
f1 = chain.form_matrix(inputs[0], fission_yields)
|
||||
f2 = chain.form_matrix(inputs[1], fission_yields)
|
||||
f3 = chain.form_matrix(inputs[2], fission_yields)
|
||||
dt_l, dt = inputs[3], inputs[4]
|
||||
return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
|
||||
+ (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2)
|
||||
/ (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3)
|
||||
|
||||
|
||||
def leqi_f4(chain, inputs):
|
||||
f1 = chain.form_matrix(inputs[0])
|
||||
f2 = chain.form_matrix(inputs[1])
|
||||
f3 = chain.form_matrix(inputs[2])
|
||||
def leqi_f4(chain, inputs, fission_yields=None):
|
||||
f1 = chain.form_matrix(inputs[0], fission_yields)
|
||||
f2 = chain.form_matrix(inputs[1], fission_yields)
|
||||
f3 = chain.form_matrix(inputs[2], fission_yields)
|
||||
dt_l, dt = inputs[3], inputs[4]
|
||||
return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
|
||||
+ (dt ** 2 + 2 * dt * dt_l + dt_l ** 2)
|
||||
|
|
|
|||
|
|
@ -13,10 +13,11 @@ from copy import deepcopy
|
|||
from warnings import warn
|
||||
from numbers import Real, Integral
|
||||
|
||||
from numpy import nonzero, empty
|
||||
from numpy import nonzero, empty, asarray
|
||||
from uncertainties import ufloat
|
||||
|
||||
from openmc.data import DataLibrary, JOULE_PER_EV
|
||||
from openmc.capi import MaterialFilter, Tally
|
||||
from openmc.checkvalue import check_type, check_greater_than
|
||||
from .results import Results
|
||||
from .chain import Chain
|
||||
|
|
@ -364,6 +365,217 @@ class EnergyHelper(ABC):
|
|||
self._nuclides = nuclides
|
||||
|
||||
|
||||
class FissionYieldHelper(ABC):
|
||||
"""Abstract class for processing energy dependent fission yields
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.Nuclide
|
||||
Nuclides tracked in the depletion chain. All nuclides are
|
||||
not required to have fission yield data.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
constant_yields : dict of str to :class:`openmc.deplete.FissionYield`
|
||||
Fission yields for all nuclides that only have one set of
|
||||
fission yield data. Can be accessed as ``{parent: {product: yield}}``
|
||||
"""
|
||||
|
||||
def __init__(self, chain_nuclides):
|
||||
self._chain_nuclides = {}
|
||||
self._constant_yields = {}
|
||||
|
||||
# Get all nuclides with fission yield data
|
||||
for nuc in chain_nuclides:
|
||||
if nuc.yield_data is None:
|
||||
continue
|
||||
if len(nuc.yield_data) == 1:
|
||||
self._constant_yields[nuc.name] = (
|
||||
nuc.yield_data[nuc.yield_energies[0]])
|
||||
elif len(nuc.yield_data) > 1:
|
||||
self._chain_nuclides[nuc.name] = nuc
|
||||
self._chain_set = set(self._chain_nuclides) | set(self._constant_yields)
|
||||
|
||||
@property
|
||||
def constant_yields(self):
|
||||
return deepcopy(self._constant_yields)
|
||||
|
||||
@abstractmethod
|
||||
def weighted_yields(self, local_mat_index):
|
||||
"""Return fission yields for a specific material
|
||||
|
||||
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`
|
||||
|
||||
Returns
|
||||
-------
|
||||
library : dict
|
||||
Dictionary of ``{parent: {product: fyield}}``
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def unpack():
|
||||
"""Unpack tally data prior to compute fission yields.
|
||||
|
||||
Called after a :meth:`openmc.deplete.Operator.__call__`
|
||||
routine during the normalization of reaction rates.
|
||||
|
||||
Not necessary for all subclasses to implement, unless tallies
|
||||
are used.
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def generate_tallies(materials, mat_indexes):
|
||||
"""Construct tallies necessary for computing fission yields
|
||||
|
||||
Called during the operator set up phase prior to depleting.
|
||||
Not necessary for subclasses to implement
|
||||
|
||||
Parameters
|
||||
----------
|
||||
materials : iterable of C-API materials
|
||||
Materials to be used in :class:`openmc.capi.MaterialFilter`
|
||||
mat_indexes : iterable of int
|
||||
Indices of tallied materials that will have their fission
|
||||
yields computed by this helper. Necessary as the
|
||||
:class:`openmc.deplete.Operator` that uses this helper
|
||||
may only burn a subset of all materials when running
|
||||
in parallel mode.
|
||||
"""
|
||||
|
||||
def update_tally_nuclides(self, nuclides):
|
||||
"""Return nuclides with non-zero densities and yield data
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : iterable of str
|
||||
Nuclides with non-zero densities from the
|
||||
:class:`openmc.deplete.Operator`
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : list of str
|
||||
Union of nuclides that the :class:`openmc.deplete.Operator`
|
||||
says have non-zero densities at this stage and those that
|
||||
have yield data. Sorted by nuclide name
|
||||
|
||||
"""
|
||||
return sorted(self._chain_set & set(nuclides))
|
||||
|
||||
@classmethod
|
||||
def from_operator(cls, operator, **kwargs):
|
||||
"""Create a new instance by pulling data from the operator
|
||||
|
||||
All keyword arguments should be identical to their counterpart
|
||||
in the main ``__init__`` method
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator with a depletion chain
|
||||
kwargs: optional
|
||||
Additional keyword arguments to be used in constuction
|
||||
"""
|
||||
return cls(operator.chain.nuclides, **kwargs)
|
||||
|
||||
|
||||
class TalliedFissionYieldHelper(FissionYieldHelper):
|
||||
"""Abstract class for computing fission yields with tallies
|
||||
|
||||
Generates a basic fission rate tally in all burnable materials with
|
||||
:meth:`generate_tallies`, and set nuclides to be tallied with
|
||||
:meth:`update_tally_nuclides`. Subclasses will need to implement
|
||||
:meth:`unpack` and :meth:`weighted_yields`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.Nuclide
|
||||
Nuclides tracked in the depletion chain. Not necessary
|
||||
that all have yield data.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
constant_yields : dict of str to :class:`openmc.deplete.FissionYield`
|
||||
Fission yields for all nuclides that only have one set of
|
||||
fission yield data. Can be accessed as ``{parent: {product: yield}}``
|
||||
results : None or numpy.ndarray
|
||||
Tally results shaped in a manner useful to this helper.
|
||||
"""
|
||||
|
||||
_upper_energy = 20.0e6 # upper energy for tallies
|
||||
|
||||
def __init__(self, chain_nuclides):
|
||||
super().__init__(chain_nuclides)
|
||||
self._local_indexes = None
|
||||
self._fission_rate_tally = None
|
||||
self._tally_nucs = []
|
||||
self.results = None
|
||||
|
||||
def generate_tallies(self, materials, mat_indexes):
|
||||
"""Construct the fission rate tally
|
||||
|
||||
Parameters
|
||||
----------
|
||||
materials : iterable of :class:`openmc.capi.Material`
|
||||
Materials to be used in :class:`openmc.capi.MaterialFilter`
|
||||
mat_indexes : iterable of int
|
||||
Indices of tallied materials that will have their fission
|
||||
yields computed by this helper. Necessary as the
|
||||
:class:`openmc.deplete.Operator` that uses this helper
|
||||
may only burn a subset of all materials when running
|
||||
in parallel mode.
|
||||
"""
|
||||
self._local_indexes = asarray(mat_indexes)
|
||||
|
||||
# Tally group-wise fission reaction rates
|
||||
self._fission_rate_tally = Tally()
|
||||
self._fission_rate_tally.scores = ['fission']
|
||||
|
||||
self._fission_rate_tally.filters = [MaterialFilter(materials)]
|
||||
|
||||
def update_tally_nuclides(self, nuclides):
|
||||
"""Tally nuclides with non-zero density and multiple yields
|
||||
|
||||
Must be run after :meth:`generate_tallies`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : iterable of str
|
||||
Potential nuclides to be tallied, such as those with
|
||||
non-zero density at this stage.
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : list of str
|
||||
Union of input nuclides and those that have multiple sets
|
||||
of yield data. Sorted by nuclide name
|
||||
|
||||
Raises
|
||||
------
|
||||
AttributeError
|
||||
If tallies not generated
|
||||
"""
|
||||
assert self._fission_rate_tally is not None, (
|
||||
"Run generate_tallies first")
|
||||
overlap = set(self._chain_nuclides).intersection(set(nuclides))
|
||||
nuclides = sorted(overlap)
|
||||
self._tally_nucs = [self._chain_nuclides[n] for n in nuclides]
|
||||
self._fission_rate_tally.nuclides = nuclides
|
||||
return nuclides
|
||||
|
||||
@abstractmethod
|
||||
def unpack(self):
|
||||
"""Unpack tallies after a transport run.
|
||||
|
||||
Abstract because each subclass will need to arrange its
|
||||
tally data.
|
||||
"""
|
||||
|
||||
|
||||
class Integrator(ABC):
|
||||
"""Abstract class for solving the time-integration for depletion
|
||||
|
||||
|
|
|
|||
|
|
@ -9,12 +9,13 @@ from itertools import chain
|
|||
import math
|
||||
import re
|
||||
from collections import OrderedDict, defaultdict
|
||||
from collections.abc import Mapping
|
||||
from collections.abc import Mapping, Iterable
|
||||
from numbers import Real
|
||||
from warnings import warn
|
||||
|
||||
from openmc.checkvalue import check_type, check_greater_than
|
||||
from openmc.data import gnd_name, zam
|
||||
from .nuclide import FissionYieldDistribution
|
||||
|
||||
# Try to use lxml if it is available. It preserves the order of attributes and
|
||||
# provides a pretty-printer by default. If not available,
|
||||
|
|
@ -136,12 +137,22 @@ class Chain(object):
|
|||
Reactions that are tracked in the depletion chain
|
||||
nuclide_dict : OrderedDict of str to int
|
||||
Maps a nuclide name to an index in nuclides.
|
||||
fission_yields : None or iterable of dict
|
||||
List of effective fission yields for materials. Each dictionary
|
||||
should be of the form ``{parent: {product: yield}}`` with
|
||||
types ``{str: {str: float}}``, where ``yield`` is the fission product
|
||||
yield for isotope ``parent`` producing isotope ``product``.
|
||||
A single entry indicates yields are constant across all materials.
|
||||
Otherwise, an entry can be added for each material to be burned.
|
||||
Ordering should be identical to how the operator orders reaction
|
||||
rates for burnable materials.
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.nuclides = []
|
||||
self.reactions = []
|
||||
self.nuclide_dict = OrderedDict()
|
||||
self._fission_yields = None
|
||||
|
||||
def __contains__(self, nuclide):
|
||||
return nuclide in self.nuclide_dict
|
||||
|
|
@ -209,8 +220,7 @@ class Chain(object):
|
|||
for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)):
|
||||
data = decay_data[parent]
|
||||
|
||||
nuclide = Nuclide()
|
||||
nuclide.name = parent
|
||||
nuclide = Nuclide(parent)
|
||||
|
||||
chain.nuclides.append(nuclide)
|
||||
chain.nuclide_dict[parent] = idx
|
||||
|
|
@ -225,7 +235,8 @@ class Chain(object):
|
|||
if mode.daughter in decay_data:
|
||||
target = mode.daughter
|
||||
else:
|
||||
print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter))
|
||||
print('missing {} {} {}'.format(
|
||||
parent, ','.join(mode.modes), mode.daughter))
|
||||
target = replace_missing(mode.daughter, decay_data)
|
||||
|
||||
# Write branching ratio, taking care to ensure sum is unity
|
||||
|
|
@ -279,15 +290,16 @@ class Chain(object):
|
|||
fpy = fpy_data[parent]
|
||||
|
||||
if fpy.energies is not None:
|
||||
nuclide.yield_energies = fpy.energies
|
||||
yield_energies = fpy.energies
|
||||
else:
|
||||
nuclide.yield_energies = [0.0]
|
||||
yield_energies = [0.0]
|
||||
|
||||
for E, table in zip(nuclide.yield_energies, fpy.independent):
|
||||
yield_data = {}
|
||||
for E, table in zip(yield_energies, fpy.independent):
|
||||
yield_replace = 0.0
|
||||
yields = defaultdict(float)
|
||||
for product, y in table.items():
|
||||
# Handle fission products that have no decay data available
|
||||
# Handle fission products that have no decay data
|
||||
if product not in decay_data:
|
||||
daughter = replace_missing(product, decay_data)
|
||||
product = daughter
|
||||
|
|
@ -297,10 +309,9 @@ class Chain(object):
|
|||
|
||||
if yield_replace > 0.0:
|
||||
missing_fp.append((parent, E, yield_replace))
|
||||
yield_data[E] = yields
|
||||
|
||||
nuclide.yield_data[E] = []
|
||||
for k in sorted(yields, key=openmc.data.zam):
|
||||
nuclide.yield_data[E].append((k, yields[k]))
|
||||
nuclide.yield_data = FissionYieldDistribution(yield_data)
|
||||
|
||||
# Display warnings
|
||||
if missing_daughter:
|
||||
|
|
@ -387,23 +398,56 @@ class Chain(object):
|
|||
clean_indentation(root_elem)
|
||||
tree.write(str(filename), encoding='utf-8')
|
||||
|
||||
def form_matrix(self, rates):
|
||||
def get_thermal_fission_yields(self):
|
||||
"""Return fission yields at lowest incident neutron energy
|
||||
|
||||
Used as the default set of fission yields for :meth:`form_matrix`
|
||||
if ``fission_yields`` are not provided
|
||||
|
||||
Returns
|
||||
-------
|
||||
fission_yields : dict
|
||||
Dictionary of ``{parent: {product: f_yield}}``
|
||||
where ``parent`` and ``product`` are both string
|
||||
names of nuclides with yield data and ``f_yield``
|
||||
is a float for the fission yield.
|
||||
"""
|
||||
out = {}
|
||||
for nuc in self.nuclides:
|
||||
if nuc.yield_data is None:
|
||||
continue
|
||||
yield_obj = nuc.yield_data[min(nuc.yield_energies)]
|
||||
out[nuc.name] = dict(yield_obj)
|
||||
return out
|
||||
|
||||
def form_matrix(self, rates, fission_yields=None):
|
||||
"""Forms depletion matrix.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rates : numpy.ndarray
|
||||
2D array indexed by (nuclide, reaction)
|
||||
fission_yields : dict, optional
|
||||
Option to use a custom set of fission yields. Expected
|
||||
to be of the form ``{parent : {product : f_yield}}``
|
||||
with string nuclide names for ``parent`` and ``product``,
|
||||
and ``f_yield`` as the respective fission yield
|
||||
|
||||
Returns
|
||||
-------
|
||||
scipy.sparse.csr_matrix
|
||||
Sparse matrix representing depletion.
|
||||
|
||||
See Also
|
||||
--------
|
||||
:meth:`get_thermal_fission_yields`
|
||||
"""
|
||||
matrix = defaultdict(float)
|
||||
reactions = set()
|
||||
|
||||
if fission_yields is None:
|
||||
fission_yields = self.get_thermal_fission_yields()
|
||||
|
||||
for i, nuc in enumerate(self.nuclides):
|
||||
|
||||
if nuc.n_decay_modes != 0:
|
||||
|
|
@ -448,11 +492,7 @@ 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..
|
||||
energy, data = sorted(nuc.yield_data.items())[0]
|
||||
for product, y in data:
|
||||
for product, y in fission_yields[nuc.name].items():
|
||||
yield_val = y * path_rate
|
||||
if yield_val != 0.0:
|
||||
k = self.nuclide_dict[product]
|
||||
|
|
@ -678,6 +718,20 @@ class Chain(object):
|
|||
parent.reactions.append(ReactionTuple(
|
||||
reaction, ground_tgt, rxn_Q, ground_br))
|
||||
|
||||
@property
|
||||
def fission_yields(self):
|
||||
if self._fission_yields is None:
|
||||
self._fission_yields = [self.get_thermal_fission_yields()]
|
||||
return self._fission_yields
|
||||
|
||||
@fission_yields.setter
|
||||
def fission_yields(self, yields):
|
||||
if yields is not None:
|
||||
if isinstance(yields, Mapping):
|
||||
yields = [yields]
|
||||
check_type("fission_yields", yields, Iterable, Mapping)
|
||||
self._fission_yields = yields
|
||||
|
||||
def validate(self, strict=True, quiet=False, tolerance=1e-4):
|
||||
"""Search for possible inconsistencies
|
||||
|
||||
|
|
|
|||
|
|
@ -30,7 +30,9 @@ def deplete(chain, x, rates, dt, matrix_func=None):
|
|||
dt : float
|
||||
Time in [s] to deplete for
|
||||
maxtrix_func : Callable, optional
|
||||
Function of two variables: ``chain`` and ``rates``.
|
||||
Function to form the depletion matrix after calling
|
||||
``matrix_func(chain, rates, fission_yields)``, where
|
||||
``fission_yields = {parent: {product: yield_frac}}``
|
||||
Expected to return the depletion matrix required by
|
||||
:func:`CRAM48`.
|
||||
|
||||
|
|
@ -40,9 +42,19 @@ def deplete(chain, x, rates, dt, matrix_func=None):
|
|||
Updated atom number vectors for each material
|
||||
"""
|
||||
|
||||
fission_yields = chain.fission_yields
|
||||
if len(fission_yields) == 1:
|
||||
fission_yields = repeat(fission_yields[0])
|
||||
elif len(fission_yields) != len(x):
|
||||
raise ValueError(
|
||||
"Number of material fission yield distributions {} is not equal "
|
||||
"to the number of compositions {}".format(len(fission_yields),
|
||||
len(x)))
|
||||
|
||||
# Use multiprocessing pool to distribute work
|
||||
with Pool() as pool:
|
||||
iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func))
|
||||
iters = zip(repeat(chain), x, rates, repeat(dt),
|
||||
fission_yields, repeat(matrix_func))
|
||||
x_result = list(pool.starmap(_cram_wrapper, iters))
|
||||
|
||||
return x_result
|
||||
|
|
@ -67,7 +79,7 @@ def timed_deplete(*args, **kwargs):
|
|||
return time.time() - start, results
|
||||
|
||||
|
||||
def _cram_wrapper(chain, n0, rates, dt, matrix_func=None):
|
||||
def _cram_wrapper(chain, n0, rates, dt, fission_yields, matrix_func=None):
|
||||
"""Wraps depletion matrix creation / CRAM solve for multiprocess execution
|
||||
|
||||
Parameters
|
||||
|
|
@ -82,6 +94,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None):
|
|||
Time to integrate to.
|
||||
maxtrix_func : function, optional
|
||||
Function to form the depletion matrix
|
||||
fission_yields : dict
|
||||
Single-energy fission yields of the form
|
||||
``{parent: {product: fission_yield}}``
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -90,9 +105,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None):
|
|||
"""
|
||||
|
||||
if matrix_func is None:
|
||||
A = chain.form_matrix(rates)
|
||||
A = chain.form_matrix(rates, fission_yields)
|
||||
else:
|
||||
A = matrix_func(chain, rates)
|
||||
A = matrix_func(chain, rates, fission_yields)
|
||||
return CRAM48(A, n0, dt)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,12 +1,24 @@
|
|||
"""
|
||||
Class for normalizing fission energy deposition
|
||||
"""
|
||||
from copy import deepcopy
|
||||
from itertools import product
|
||||
from numbers import Real
|
||||
import bisect
|
||||
|
||||
from numpy import dot, zeros
|
||||
from numpy import dot, zeros, newaxis
|
||||
|
||||
from openmc.capi import Tally, MaterialFilter
|
||||
from .abc import ReactionRateHelper, EnergyHelper
|
||||
from openmc.checkvalue import check_type, check_greater_than
|
||||
from openmc.capi import (
|
||||
Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter)
|
||||
from .abc import (
|
||||
ReactionRateHelper, EnergyHelper, FissionYieldHelper,
|
||||
TalliedFissionYieldHelper)
|
||||
|
||||
__all__ = (
|
||||
"DirectReactionRateHelper", "ChainFissionHelper",
|
||||
"ConstantFissionYieldHelper", "FissionYieldCutoffHelper",
|
||||
"AveragedFissionYieldHelper")
|
||||
|
||||
# -------------------------------------
|
||||
# Helpers for generating reaction rates
|
||||
|
|
@ -142,3 +154,465 @@ 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 ConstantFissionYieldHelper(FissionYieldHelper):
|
||||
"""Class that uses a single set of fission yields on each isotope
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.Nuclide
|
||||
Nuclides tracked in the depletion chain. All nuclides are
|
||||
not required to have fission yield data.
|
||||
energy : float, optional
|
||||
Key in :attr:`openmc.deplete.Nuclide.yield_data` corresponding
|
||||
to the desired set of fission yield data. Typically one of
|
||||
``{0.0253, 500000, 14000000}`` corresponding to 0.0253 eV,
|
||||
500 keV, and 14 MeV yield libraries. If the specific key is not
|
||||
found, will fall back to closest energy present.
|
||||
Default: 0.0253 eV for thermal yields
|
||||
|
||||
Attributes
|
||||
----------
|
||||
constant_yields : dict of str to :class:`openmc.deplete.FissionYield`
|
||||
Fission yields for all nuclides that only have one set of
|
||||
fission yield data. Can be accessed as ``{parent: {product: yield}}``
|
||||
energy : float
|
||||
Energy of fission yield libraries.
|
||||
"""
|
||||
|
||||
def __init__(self, chain_nuclides, energy=0.0253):
|
||||
check_type("energy", energy, Real)
|
||||
check_greater_than("energy", energy, 0.0, equality=True)
|
||||
self._energy = energy
|
||||
super().__init__(chain_nuclides)
|
||||
# Iterate over all nuclides with > 1 set of yields
|
||||
for name, nuc in self._chain_nuclides.items():
|
||||
yield_data = nuc.yield_data.get(energy)
|
||||
if yield_data is not None:
|
||||
self._constant_yields[name] = yield_data
|
||||
continue
|
||||
# Specific energy not found, use closest energy
|
||||
distances = [abs(energy - ene) for ene in nuc.yield_energies]
|
||||
min_E = min(nuc.yield_energies, key=lambda e: abs(e - energy))
|
||||
self._constant_yields[name] = nuc.yield_data[min_E]
|
||||
|
||||
@classmethod
|
||||
def from_operator(cls, operator, **kwargs):
|
||||
"""Return a new ConstantFissionYieldHelper using operator data
|
||||
|
||||
All keyword arguments should be identical to their counterpart
|
||||
in the main ``__init__`` method
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
operator with a depletion chain
|
||||
kwargs:
|
||||
Additional keyword arguments to be used in construction
|
||||
|
||||
Returns
|
||||
-------
|
||||
ConstantFissionYieldHelper
|
||||
"""
|
||||
return cls(operator.chain.nuclides, **kwargs)
|
||||
|
||||
@property
|
||||
def energy(self):
|
||||
return self._energy
|
||||
|
||||
def weighted_yields(self, _local_mat_index=None):
|
||||
"""Return fission yields for all nuclides requested
|
||||
|
||||
Parameters
|
||||
----------
|
||||
_local_mat_index : int, optional
|
||||
Current material index. Not used since all yields are
|
||||
constant
|
||||
|
||||
Returns
|
||||
-------
|
||||
library : dict
|
||||
Dictionary of ``{parent: {product: fyield}}``
|
||||
"""
|
||||
return self.constant_yields
|
||||
|
||||
|
||||
class FissionYieldCutoffHelper(TalliedFissionYieldHelper):
|
||||
"""Helper that computes fission yields based on a cutoff energy
|
||||
|
||||
Tally fission rates above and below the cutoff energy.
|
||||
Assume that all fissions below cutoff energy have use thermal fission
|
||||
product yield distributions, while all fissions above use a faster
|
||||
set of yield distributions.
|
||||
|
||||
Uses a limit of 20 MeV for tallying fission.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.Nuclide
|
||||
Nuclides tracked in the depletion chain. All nuclides are
|
||||
not required to have fission yield data.
|
||||
n_bmats : int, optional
|
||||
Number of burnable materials tracked in the problem
|
||||
cutoff : float, optional
|
||||
Cutoff energy in [eV] below which all fissions will be
|
||||
use thermal yields. All other fissions will use a
|
||||
faster set of yields. Default: 112 [eV]
|
||||
thermal_energy : float, optional
|
||||
Energy of yield data corresponding to thermal yields.
|
||||
Default: 0.0253 [eV]
|
||||
fast_energy : float, optional
|
||||
Energy of yield data corresponding to fast yields.
|
||||
Default: 500 [kev]
|
||||
|
||||
Attributes
|
||||
----------
|
||||
n_bmats : int
|
||||
Number of burnable materials tracked in the problem.
|
||||
Must be set prior to generating tallies
|
||||
thermal_yields : dict
|
||||
Dictionary of the form ``{parent: {product: yield}}``
|
||||
with thermal yields
|
||||
fast_yields : dict
|
||||
Dictionary of the form ``{parent: {product: yield}}``
|
||||
with fast yields
|
||||
results : numpy.ndarray
|
||||
Array of fission rate fractions with shape
|
||||
``(n_mats, 2, n_nucs)``. ``results[:, 0]``
|
||||
corresponds to the fraction of all fissions
|
||||
that occured below ``cutoff``. The number
|
||||
of materials in the first axis corresponds
|
||||
to the number of materials burned by the
|
||||
:class:`openmc.deplete.Operator`
|
||||
"""
|
||||
|
||||
def __init__(self, chain_nuclides, n_bmats, cutoff=112.0,
|
||||
thermal_energy=0.0253, fast_energy=500.0e3):
|
||||
check_type("cutoff", cutoff, Real)
|
||||
check_type("thermal_energy", thermal_energy, Real)
|
||||
check_type("fast_energy", fast_energy, Real)
|
||||
check_greater_than("thermal_energy", thermal_energy, 0.0, equality=True)
|
||||
check_greater_than("cutoff", cutoff, thermal_energy, equality=False)
|
||||
check_greater_than("fast_energy", fast_energy, cutoff, equality=False)
|
||||
self.n_bmats = n_bmats
|
||||
super().__init__(chain_nuclides)
|
||||
self._cutoff = cutoff
|
||||
self._thermal_yields = {}
|
||||
self._fast_yields = {}
|
||||
convert_to_constant = set()
|
||||
for name, nuc in self._chain_nuclides.items():
|
||||
yields = nuc.yield_data
|
||||
energies = nuc.yield_energies
|
||||
thermal = yields.get(thermal_energy)
|
||||
fast = yields.get(fast_energy)
|
||||
if thermal is None or fast is None:
|
||||
if cutoff <= energies[0]:
|
||||
# use lowest energy yields as constant
|
||||
self._constant_yields[name] = yields[energies[0]]
|
||||
convert_to_constant.add(name)
|
||||
continue
|
||||
if cutoff >= energies[-1]:
|
||||
# use highest energy yields as constant
|
||||
self._constant_yields[name] = yields[energies[-1]]
|
||||
convert_to_constant.add(name)
|
||||
continue
|
||||
cutoff_ix = bisect.bisect_left(energies, cutoff)
|
||||
# find closest energy to requested thermal, fast energies
|
||||
if thermal is None:
|
||||
min_E = min(energies[:cutoff_ix],
|
||||
key=lambda e: abs(e - thermal_energy))
|
||||
thermal = yields[min_E]
|
||||
if fast is None:
|
||||
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
|
||||
for name in convert_to_constant:
|
||||
self._chain_nuclides.pop(name)
|
||||
|
||||
@classmethod
|
||||
def from_operator(cls, operator, **kwargs):
|
||||
"""Construct a helper from an operator
|
||||
|
||||
All keyword arguments should be identical to their counterpart
|
||||
in the main ``__init__`` method
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.Operator
|
||||
Operator with a chain and burnable materials
|
||||
kwargs:
|
||||
Additional keyword arguments to be used in construction
|
||||
|
||||
Returns
|
||||
-------
|
||||
FissionYieldCutoffHelper
|
||||
|
||||
"""
|
||||
return cls(operator.chain.nuclides, len(operator.burnable_mats),
|
||||
**kwargs)
|
||||
|
||||
def generate_tallies(self, materials, mat_indexes):
|
||||
"""Use C API to produce a fission rate tally in burnable materials
|
||||
|
||||
Include a :class:`openmc.capi.EnergyFilter` to tally fission rates
|
||||
above and below cutoff energy.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
materials : iterable of :class:`openmc.capi.Material`
|
||||
Materials to be used in :class:`openmc.capi.MaterialFilter`
|
||||
mat_indexes : iterable of int
|
||||
Indices of tallied materials that will have their fission
|
||||
yields computed by this helper. Necessary as the
|
||||
:class:`openmc.deplete.Operator` that uses this helper
|
||||
may only burn a subset of all materials when running
|
||||
in parallel mode.
|
||||
"""
|
||||
super().generate_tallies(materials, mat_indexes)
|
||||
energy_filter = EnergyFilter([0.0, self._cutoff, self._upper_energy])
|
||||
self._fission_rate_tally.filters = (
|
||||
self._fission_rate_tally.filters + [energy_filter])
|
||||
|
||||
def unpack(self):
|
||||
"""Obtain fast and thermal fission fractions from tally"""
|
||||
if not self._tally_nucs:
|
||||
self.results = None
|
||||
return
|
||||
fission_rates = self._fission_rate_tally.results[..., 1].reshape(
|
||||
self.n_bmats, 2, len(self._tally_nucs))
|
||||
self.results = fission_rates[self._local_indexes]
|
||||
total_fission = self.results.sum(axis=1)
|
||||
nz_mat, nz_nuc = total_fission.nonzero()
|
||||
self.results[nz_mat, :, nz_nuc] /= total_fission[nz_mat, newaxis, nz_nuc]
|
||||
|
||||
def weighted_yields(self, local_mat_index):
|
||||
"""Return fission yields for a specific material
|
||||
|
||||
For nuclides with both yield data above and below
|
||||
the cutoff energy, the effective yield for nuclide ``A``
|
||||
will be a weighted sum of fast and thermal yields. The
|
||||
weights will be the fraction of ``A`` fission events
|
||||
in the above and below the cutoff energy.
|
||||
|
||||
If ``A`` has fission product distribution ``F``
|
||||
for fast fissions and ``T`` for thermal fissions, and
|
||||
70% of ``A`` fissions are considered thermal, then
|
||||
the effective fission product yield distributions
|
||||
for ``A`` is ``0.7 * T + 0.3 * F``
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mat_index : int
|
||||
Index for specific burnable material. Effective
|
||||
yields will be produced using
|
||||
``self.results[local_mat_index]``
|
||||
|
||||
Returns
|
||||
-------
|
||||
library : dict
|
||||
Dictionary of ``{parent: {product: fyield}}``
|
||||
"""
|
||||
yields = self.constant_yields
|
||||
if not self._tally_nucs:
|
||||
return yields
|
||||
rates = self.results[local_mat_index]
|
||||
# iterate over thermal then fast yields, prefer __mul__ to __rmul__
|
||||
for therm_frac, fast_frac, nuc in zip(rates[0], rates[1], self._tally_nucs):
|
||||
yields[nuc.name] = (self._thermal_yields[nuc.name] * therm_frac
|
||||
+ self._fast_yields[nuc.name] * fast_frac)
|
||||
return yields
|
||||
|
||||
@property
|
||||
def thermal_yields(self):
|
||||
return deepcopy(self._thermal_yields)
|
||||
|
||||
@property
|
||||
def fast_yields(self):
|
||||
return deepcopy(self._fast_yields)
|
||||
|
||||
|
||||
class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
|
||||
r"""Class that computes fission yields based on average fission energy
|
||||
|
||||
Computes average energy at which fission events occured with
|
||||
|
||||
.. math::
|
||||
|
||||
\bar{E} = \frac{
|
||||
\int_0^\infty E\sigma_f(E)\phi(E)dE
|
||||
}{
|
||||
\int_0^\infty\sigma_f(E)\phi(E)dE
|
||||
}
|
||||
|
||||
If the average energy for a nuclide is below the lowest energy
|
||||
with yield data, that set of fission yields is taken.
|
||||
Conversely, if the average energy is above the highest energy
|
||||
with yield data, that set of fission yields is used.
|
||||
For the case where the average energy is between two sets
|
||||
of yields, the effective fission yield computed by
|
||||
linearly interpolating between yields provided at the
|
||||
nearest energies above and below the average.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_nuclides : iterable of openmc.deplete.Nuclide
|
||||
Nuclides tracked in the depletion chain. All nuclides are
|
||||
not required to have fission yield data.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
constant_yields : dict of str to :class:`openmc.deplete.FissionYield`
|
||||
Fission yields for all nuclides that only have one set of
|
||||
fission yield data. Can be accessed as ``{parent: {product: yield}}``
|
||||
results : None or numpy.ndarray
|
||||
If tallies have been generated and unpacked, then the array will
|
||||
have shape ``(n_mats, n_tnucs)``, where ``n_mats`` is the number
|
||||
of materials where fission reactions were tallied and ``n_tnucs``
|
||||
is the number of nuclides with multiple sets of fission yields.
|
||||
Data in the array are the average energy of fission events for
|
||||
tallied nuclides across burnable materials.
|
||||
"""
|
||||
|
||||
def __init__(self, chain_nuclides):
|
||||
super().__init__(chain_nuclides)
|
||||
self._weighted_tally = None
|
||||
|
||||
def generate_tallies(self, materials, mat_indexes):
|
||||
"""Construct tallies to determine average energy of fissions
|
||||
|
||||
Parameters
|
||||
----------
|
||||
materials : iterable of :class:`openmc.capi.Material`
|
||||
Materials to be used in :class:`openmc.capi.MaterialFilter`
|
||||
mat_indexes : iterable of int
|
||||
Indices of tallied materials that will have their fission
|
||||
yields computed by this helper. Necessary as the
|
||||
:class:`openmc.deplete.Operator` that uses this helper
|
||||
may only burn a subset of all materials when running
|
||||
in parallel mode.
|
||||
"""
|
||||
super().generate_tallies(materials, mat_indexes)
|
||||
fission_tally = self._fission_rate_tally
|
||||
filters = fission_tally.filters
|
||||
|
||||
ene_filter = EnergyFilter([0, self._upper_energy])
|
||||
fission_tally.filters = filters + [ene_filter]
|
||||
|
||||
func_filter = EnergyFunctionFilter()
|
||||
func_filter.set_data((0, self._upper_energy), (0, self._upper_energy))
|
||||
weighted_tally = Tally()
|
||||
weighted_tally.scores = ['fission']
|
||||
weighted_tally.filters = filters + [func_filter]
|
||||
self._weighted_tally = weighted_tally
|
||||
|
||||
def update_tally_nuclides(self, nuclides):
|
||||
"""Tally nuclides with non-zero density and multiple yields
|
||||
|
||||
Must be run after :meth:`generate_tallies`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : iterable of str
|
||||
Potential nuclides to be tallied, such as those with
|
||||
non-zero density at this stage.
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : tuple of str
|
||||
Union of input nuclides and those that have multiple sets
|
||||
of yield data. Sorted by nuclide name
|
||||
|
||||
Raises
|
||||
------
|
||||
AttributeError
|
||||
If tallies not generated
|
||||
"""
|
||||
tally_nucs = super().update_tally_nuclides(nuclides)
|
||||
self._weighted_tally.nuclides = tally_nucs
|
||||
return tally_nucs
|
||||
|
||||
def unpack(self):
|
||||
"""Unpack tallies and populate :attr:`results` with average energies"""
|
||||
if not self._tally_nucs:
|
||||
self.results = None
|
||||
return
|
||||
fission_results = (
|
||||
self._fission_rate_tally.results[self._local_indexes, :, 1])
|
||||
self.results = (
|
||||
self._weighted_tally.results[self._local_indexes, :, 1]).copy()
|
||||
nz_mat, nz_nuc = fission_results.nonzero()
|
||||
self.results[nz_mat, nz_nuc] /= fission_results[nz_mat, nz_nuc]
|
||||
|
||||
def weighted_yields(self, local_mat_index):
|
||||
"""Return fission yields for a specific material
|
||||
|
||||
Use the computed average energy of fission
|
||||
events to determine fission yields. If average
|
||||
energy is between two sets of yields, linearly
|
||||
interpolate bewteen the two.
|
||||
Otherwise take the closet set of yields.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mat_index : int
|
||||
Index for specific burnable material. Effective
|
||||
yields will be produced using
|
||||
``self.results[local_mat_index]``
|
||||
|
||||
Returns
|
||||
-------
|
||||
library : dict
|
||||
Dictionary of ``{parent: {product: fyield}}``
|
||||
"""
|
||||
if not self._tally_nucs:
|
||||
return self.constant_yields
|
||||
mat_yields = {}
|
||||
average_energies = self.results[local_mat_index]
|
||||
for avg_e, nuc in zip(average_energies, self._tally_nucs):
|
||||
nuc_energies = nuc.yield_energies
|
||||
if avg_e <= nuc_energies[0]:
|
||||
mat_yields[nuc.name] = nuc.yield_data[nuc_energies[0]]
|
||||
continue
|
||||
if avg_e >= nuc_energies[-1]:
|
||||
mat_yields[nuc.name] = nuc.yield_data[nuc_energies[-1]]
|
||||
continue
|
||||
# in-between two energies
|
||||
# linear search since there are usually ~3 energies
|
||||
for ix, ene in enumerate(nuc_energies[:-1]):
|
||||
if nuc_energies[ix + 1] > avg_e:
|
||||
break
|
||||
lower, upper = nuc_energies[ix:ix + 2]
|
||||
fast_frac = (avg_e - lower) / (upper - lower)
|
||||
mat_yields[nuc.name] = (
|
||||
nuc.yield_data[lower] * (1 - fast_frac)
|
||||
+ nuc.yield_data[upper] * fast_frac)
|
||||
mat_yields.update(self.constant_yields)
|
||||
return mat_yields
|
||||
|
||||
@classmethod
|
||||
def from_operator(cls, operator, **kwargs):
|
||||
"""Return a new helper with data from an operator
|
||||
|
||||
All keyword arguments should be identical to their counterpart
|
||||
in the main ``__init__`` method
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator with a depletion chain
|
||||
kwargs :
|
||||
Additional keyword arguments to be used in construction
|
||||
|
||||
Returns
|
||||
-------
|
||||
AveragedFissionYieldHelper
|
||||
"""
|
||||
return cls(operator.chain.nuclides)
|
||||
|
|
|
|||
|
|
@ -3,13 +3,20 @@
|
|||
Contains the per-nuclide components of a depletion chain.
|
||||
"""
|
||||
|
||||
import bisect
|
||||
from collections.abc import Mapping
|
||||
from collections import namedtuple, defaultdict
|
||||
from warnings import warn
|
||||
from numbers import Real
|
||||
try:
|
||||
import lxml.etree as ET
|
||||
except ImportError:
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
from numpy import empty
|
||||
|
||||
from openmc.checkvalue import check_type
|
||||
|
||||
|
||||
DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio')
|
||||
DecayTuple.__doc__ = """\
|
||||
|
|
@ -62,11 +69,16 @@ except AttributeError:
|
|||
class Nuclide(object):
|
||||
"""Decay modes, reactions, and fission yields for a single nuclide.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
name : str, optional
|
||||
GND name of this nuclide, e.g. ``"He4"``, ``"Am242_m1"``
|
||||
|
||||
Attributes
|
||||
----------
|
||||
name : str
|
||||
name : str or None
|
||||
Name of nuclide.
|
||||
half_life : float
|
||||
half_life : float or None
|
||||
Half life of nuclide in [s].
|
||||
decay_energy : float
|
||||
Energy deposited from decay in [eV].
|
||||
|
|
@ -80,17 +92,16 @@ class Nuclide(object):
|
|||
reactions : list of openmc.deplete.ReactionTuple
|
||||
Reaction information. Each element of the list is a named tuple with
|
||||
attribute 'type', 'target', 'Q', and 'branching_ratio'.
|
||||
yield_data : dict of float to list
|
||||
Maps tabulated energy to list of (product, yield) for all
|
||||
neutron-induced fission products.
|
||||
yield_energies : list of float
|
||||
Energies at which fission product yiels exist
|
||||
|
||||
yield_data : FissionYieldDistribution or None
|
||||
Fission product yields at tabulated energies for this nuclide. Can be
|
||||
treated as a nested dictionary ``{energy: {product: yield}}``
|
||||
yield_energies : tuple of float or None
|
||||
Energies at which fission product yields exist
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
def __init__(self, name=None):
|
||||
# Information about the nuclide
|
||||
self.name = None
|
||||
self.name = name
|
||||
self.half_life = None
|
||||
self.decay_energy = 0.0
|
||||
|
||||
|
|
@ -101,8 +112,7 @@ class Nuclide(object):
|
|||
self.reactions = []
|
||||
|
||||
# Neutron fission yields, if present
|
||||
self.yield_data = {}
|
||||
self.yield_energies = []
|
||||
self._yield_data = None
|
||||
|
||||
@property
|
||||
def n_decay_modes(self):
|
||||
|
|
@ -112,6 +122,29 @@ class Nuclide(object):
|
|||
def n_reaction_paths(self):
|
||||
return len(self.reactions)
|
||||
|
||||
@property
|
||||
def yield_data(self):
|
||||
if self._yield_data is None:
|
||||
return None
|
||||
return self._yield_data
|
||||
|
||||
@yield_data.setter
|
||||
def yield_data(self, fission_yields):
|
||||
if fission_yields is None:
|
||||
self._yield_data = None
|
||||
else:
|
||||
check_type("fission_yields", fission_yields, Mapping)
|
||||
if isinstance(fission_yields, FissionYieldDistribution):
|
||||
self._yield_data = fission_yields
|
||||
else:
|
||||
self._yield_data = FissionYieldDistribution(fission_yields)
|
||||
|
||||
@property
|
||||
def yield_energies(self):
|
||||
if self._yield_data is None:
|
||||
return None
|
||||
return self.yield_data.energies
|
||||
|
||||
@classmethod
|
||||
def from_xml(cls, element, fission_q=None):
|
||||
"""Read nuclide from an XML element.
|
||||
|
|
@ -166,13 +199,7 @@ class Nuclide(object):
|
|||
|
||||
fpy_elem = element.find('neutron_fission_yields')
|
||||
if fpy_elem is not None:
|
||||
for yields_elem in fpy_elem.iter('fission_yields'):
|
||||
E = float(yields_elem.get('energy'))
|
||||
products = yields_elem.find('products').text.split()
|
||||
yields = [float(y) for y in
|
||||
yields_elem.find('data').text.split()]
|
||||
nuc.yield_data[E] = list(zip(products, yields))
|
||||
nuc.yield_energies = list(sorted(nuc.yield_data.keys()))
|
||||
nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem)
|
||||
|
||||
return nuc
|
||||
|
||||
|
|
@ -212,15 +239,7 @@ class Nuclide(object):
|
|||
fpy_elem = ET.SubElement(elem, 'neutron_fission_yields')
|
||||
energy_elem = ET.SubElement(fpy_elem, 'energies')
|
||||
energy_elem.text = ' '.join(str(E) for E in self.yield_energies)
|
||||
|
||||
for E in self.yield_energies:
|
||||
yields_elem = ET.SubElement(fpy_elem, 'fission_yields')
|
||||
yields_elem.set('energy', str(E))
|
||||
|
||||
products_elem = ET.SubElement(yields_elem, 'products')
|
||||
products_elem.text = ' '.join(x[0] for x in self.yield_data[E])
|
||||
data_elem = ET.SubElement(yields_elem, 'data')
|
||||
data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E])
|
||||
self.yield_data.to_xml_element(fpy_elem)
|
||||
|
||||
return elem
|
||||
|
||||
|
|
@ -304,8 +323,8 @@ class Nuclide(object):
|
|||
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)
|
||||
for energy, fission_yield in self.yield_data.items():
|
||||
sum_yield = fission_yield.yields.sum()
|
||||
stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance
|
||||
if stat:
|
||||
continue
|
||||
|
|
@ -321,3 +340,227 @@ class Nuclide(object):
|
|||
valid = False
|
||||
|
||||
return valid
|
||||
|
||||
|
||||
class FissionYieldDistribution(Mapping):
|
||||
"""Energy-dependent fission product yields for a single nuclide
|
||||
|
||||
Can be used as a dictionary mapping energies and products to fission
|
||||
yields::
|
||||
|
||||
>>> fydist = FissionYieldDistribution{
|
||||
... {0.0253: {"Xe135": 0.021}})
|
||||
>>> fydist[0.0253]["Xe135"]
|
||||
0.021
|
||||
|
||||
Parameters
|
||||
----------
|
||||
fission_yields : dict
|
||||
Dictionary of energies and fission product yields for that energy.
|
||||
Expected to be of the form ``{float: {str: float}}``. The first
|
||||
float is the energy, typically in eV, that represents this
|
||||
distribution. The underlying dictionary maps fission products
|
||||
to their respective yields.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
energies : tuple
|
||||
Energies for which fission yields exist. Sorted by
|
||||
increasing energy
|
||||
products : tuple
|
||||
Fission products produced at all energies. Sorted by name.
|
||||
yield_matrix : numpy.ndarray
|
||||
Array ``(n_energy, n_products)`` where
|
||||
``yield_matrix[g, j]`` is the fission yield of product
|
||||
``j`` for energy group ``g``.
|
||||
|
||||
See Also
|
||||
--------
|
||||
* :meth:`from_xml_element` - Construction methods
|
||||
* :class:`FissionYield` - Class used for storing yields at a given energy
|
||||
"""
|
||||
|
||||
def __init__(self, fission_yields):
|
||||
# mapping {energy: {product: value}}
|
||||
energies = sorted(fission_yields)
|
||||
|
||||
# Get a consistent set of products to produce a matrix of yields
|
||||
shared_prod = set.union(*(set(x) for x in fission_yields.values()))
|
||||
ordered_prod = sorted(shared_prod)
|
||||
|
||||
yield_matrix = empty((len(energies), len(shared_prod)))
|
||||
|
||||
for g_index, energy in enumerate(energies):
|
||||
prod_map = fission_yields[energy]
|
||||
for prod_ix, product in enumerate(ordered_prod):
|
||||
yield_val = prod_map.get(product)
|
||||
yield_matrix[g_index, prod_ix] = (
|
||||
0.0 if yield_val is None else yield_val)
|
||||
self.energies = tuple(energies)
|
||||
self.products = tuple(ordered_prod)
|
||||
self.yield_matrix = yield_matrix
|
||||
|
||||
def __len__(self):
|
||||
return len(self.energies)
|
||||
|
||||
def __getitem__(self, energy):
|
||||
if energy not in self.energies:
|
||||
raise KeyError(energy)
|
||||
return FissionYield(
|
||||
self.products, self.yield_matrix[self.energies.index(energy)])
|
||||
|
||||
def __iter__(self):
|
||||
return iter(self.energies)
|
||||
|
||||
def __repr__(self):
|
||||
return "<{} with {} products at {} energies>".format(
|
||||
self.__class__.__name__, self.yield_matrix.shape[1],
|
||||
len(self.energies))
|
||||
|
||||
@classmethod
|
||||
def from_xml_element(cls, element):
|
||||
"""Construct a distribution from a depletion chain xml file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
element : xml.etree.ElementTree.Element
|
||||
XML element to pull fission yield data from
|
||||
|
||||
Returns
|
||||
-------
|
||||
FissionYieldDistribution
|
||||
"""
|
||||
all_yields = {}
|
||||
for elem_index, yield_elem in enumerate(element.iter("fission_yields")):
|
||||
energy = float(yield_elem.get("energy"))
|
||||
products = yield_elem.find("products").text.split()
|
||||
yields = map(float, yield_elem.find("data").text.split())
|
||||
# Get a map of products to their corresponding yield
|
||||
all_yields[energy] = dict(zip(products, yields))
|
||||
|
||||
return cls(all_yields)
|
||||
|
||||
def to_xml_element(self, root):
|
||||
"""Write fission yield data to an xml element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
root : xml.etree.ElementTree.Element
|
||||
Element to write distribution data to
|
||||
"""
|
||||
for energy, yield_obj in self.items():
|
||||
yield_element = ET.SubElement(root, "fission_yields")
|
||||
yield_element.set("energy", str(energy))
|
||||
product_elem = ET.SubElement(yield_element, "products")
|
||||
product_elem.text = " ".join(map(str, yield_obj.products))
|
||||
data_elem = ET.SubElement(yield_element, "data")
|
||||
data_elem.text = " ".join(map(str, yield_obj.yields))
|
||||
|
||||
|
||||
class FissionYield(Mapping):
|
||||
"""Mapping for fission yields of a parent at a specific energy
|
||||
|
||||
Separated to support nested dictionary-like behavior for
|
||||
:class:`FissionYieldDistribution`, and allowing math operations
|
||||
on a single vector of yields. Can in turn be used like a
|
||||
dictionary to fetch fission yields.
|
||||
Supports multiplication of a scalar to scale the fission
|
||||
yields and addition of another set of yields.
|
||||
|
||||
Does not support resizing / inserting new products that do
|
||||
not exist.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
products : tuple of str
|
||||
Products for this specific distribution
|
||||
yields : numpy.ndarray
|
||||
Fission product yields for each product in ``products``
|
||||
|
||||
Attributes
|
||||
----------
|
||||
products : tuple of str
|
||||
Products for this specific distribution
|
||||
yields : numpy.ndarray
|
||||
Fission product yields for each product in ``products``
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy
|
||||
>>> fy_vector = FissionYield(
|
||||
... ("Xe135", "I129", "Sm149"),
|
||||
... numpy.array((0.002, 0.001, 0.0003)))
|
||||
>>> fy_vector["Xe135"]
|
||||
0.002
|
||||
>>> new = fy_vector.copy()
|
||||
>>> fy_vector *= 2
|
||||
>>> fy_vector["Xe135"]
|
||||
0.004
|
||||
>>> new["Xe135"]
|
||||
0.002
|
||||
>>> (new + fy_vector)["Sm149"]
|
||||
0.0009
|
||||
>>> dict(new)
|
||||
{"Xe135": 0.002, "I129": 0.001, "Sm149": 0.0003}
|
||||
"""
|
||||
|
||||
def __init__(self, products, yields):
|
||||
self.products = products
|
||||
self.yields = yields
|
||||
|
||||
def __contains__(self, product):
|
||||
ix = bisect.bisect_left(self.products, product)
|
||||
return ix != len(self.products) and self.products[ix] == product
|
||||
|
||||
def __getitem__(self, product):
|
||||
ix = bisect.bisect_left(self.products, product)
|
||||
if ix == len(self.products) or self.products[ix] != product:
|
||||
raise KeyError(product)
|
||||
return self.yields[ix]
|
||||
|
||||
def __len__(self):
|
||||
return len(self.products)
|
||||
|
||||
def __iter__(self):
|
||||
return iter(self.products)
|
||||
|
||||
def items(self):
|
||||
"""Return pairs of product, yield"""
|
||||
return zip(self.products, self.yields)
|
||||
|
||||
def __add__(self, other):
|
||||
if not isinstance(other, FissionYield):
|
||||
return NotImplemented
|
||||
new = FissionYield(self.products, self.yields.copy())
|
||||
new += other
|
||||
return new
|
||||
|
||||
def __iadd__(self, other):
|
||||
"""Increment value from other fission yield"""
|
||||
if not isinstance(other, FissionYield):
|
||||
return NotImplemented
|
||||
self.yields += other.yields
|
||||
return self
|
||||
|
||||
def __radd__(self, other):
|
||||
return self + other
|
||||
|
||||
def __imul__(self, scalar):
|
||||
if not isinstance(scalar, Real):
|
||||
return NotImplemented
|
||||
self.yields *= scalar
|
||||
return self
|
||||
|
||||
def __mul__(self, scalar):
|
||||
if not isinstance(scalar, Real):
|
||||
return NotImplemented
|
||||
new = FissionYield(self.products, self.yields.copy())
|
||||
new *= scalar
|
||||
return new
|
||||
|
||||
def __rmul__(self, scalar):
|
||||
return self * scalar
|
||||
|
||||
def __repr__(self):
|
||||
return "<{} containing {} products and yields>".format(
|
||||
self.__class__.__name__, len(self))
|
||||
|
|
|
|||
|
|
@ -25,7 +25,9 @@ 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, ConstantFissionYieldHelper,
|
||||
FissionYieldCutoffHelper, AveragedFissionYieldHelper)
|
||||
|
||||
|
||||
def _distribute(items):
|
||||
|
|
@ -84,6 +86,21 @@ class Operator(TransportOperator):
|
|||
in initial condition to ensure they exist in the decay chain.
|
||||
Only done for nuclides with reaction rates.
|
||||
Defaults to 1.0e3.
|
||||
fission_yield_mode : {"constant", "cutoff", "average"}
|
||||
Key indicating what fission product yield scheme to use. The
|
||||
key determines what fission energy helper is used:
|
||||
|
||||
* "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper`
|
||||
* "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper`
|
||||
* "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper`
|
||||
|
||||
The documentation on these classes describe their methodology
|
||||
and differences. Default: ``"constant"``
|
||||
fission_yield_opts : dict of str to option, optional
|
||||
Optional arguments to pass to the helper determined by
|
||||
``fission_yield_mode``. Will be passed directly on to the
|
||||
helper. Passing a value of None will use the defaults for
|
||||
the associated helper.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -120,9 +137,20 @@ class Operator(TransportOperator):
|
|||
diff_burnable_mats : bool
|
||||
Whether to differentiate burnable materials with multiple instances
|
||||
"""
|
||||
_fission_helpers = {
|
||||
"average": AveragedFissionYieldHelper,
|
||||
"constant": ConstantFissionYieldHelper,
|
||||
"cutoff": FissionYieldCutoffHelper,
|
||||
}
|
||||
|
||||
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
|
||||
diff_burnable_mats=False, fission_q=None,
|
||||
dilute_initial=1.0e3):
|
||||
dilute_initial=1.0e3, fission_yield_mode="constant",
|
||||
fission_yield_opts=None):
|
||||
if fission_yield_mode not in self._fission_helpers:
|
||||
raise KeyError(
|
||||
"fission_yield_mode must be one of {}, not {}".format(
|
||||
", ".join(self._fission_helpers), fission_yield_mode))
|
||||
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
|
||||
self.round_number = False
|
||||
self.prev_res = None
|
||||
|
|
@ -178,6 +206,13 @@ class Operator(TransportOperator):
|
|||
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
|
||||
self._energy_helper = ChainFissionHelper()
|
||||
|
||||
# Select and create fission yield helper
|
||||
fission_helper = self._fission_helpers[fission_yield_mode]
|
||||
fission_yield_opts = (
|
||||
{} if fission_yield_opts is None else fission_yield_opts)
|
||||
self._yield_helper = fission_helper.from_operator(
|
||||
self, **fission_yield_opts)
|
||||
|
||||
def __call__(self, vec, power):
|
||||
"""Runs a simulation.
|
||||
|
||||
|
|
@ -204,8 +239,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._yield_helper.update_tally_nuclides(nuclides)
|
||||
|
||||
# Run OpenMC
|
||||
openmc.capi.reset()
|
||||
|
|
@ -409,6 +446,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._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 +595,10 @@ class Operator(TransportOperator):
|
|||
# Keep track of energy produced from all reactions in eV per source
|
||||
# particle
|
||||
self._energy_helper.reset()
|
||||
self._yield_helper.unpack()
|
||||
|
||||
# Store fission yield dictionaries
|
||||
fission_yields = []
|
||||
|
||||
# Create arrays to store fission Q values, reaction rates, and nuclide
|
||||
# numbers, zeroed out in material iteration
|
||||
|
|
@ -576,6 +621,9 @@ class Operator(TransportOperator):
|
|||
tally_rates = self._rate_helper.get_material_rates(
|
||||
mat_index, nuc_ind, react_ind)
|
||||
|
||||
# Compute fission yields for this material
|
||||
fission_yields.append(self._yield_helper.weighted_yields(i))
|
||||
|
||||
# Accumulate energy from fission
|
||||
self._energy_helper.update(tally_rates[:, fission_ind], mat_index)
|
||||
|
||||
|
|
@ -589,6 +637,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 = fission_yields
|
||||
|
||||
return OperatorResult(k_combined, rates)
|
||||
|
||||
def _get_nuclides_with_data(self):
|
||||
|
|
|
|||
|
|
@ -73,6 +73,49 @@ SCHEMES = {
|
|||
}
|
||||
|
||||
|
||||
class TestChain(object):
|
||||
"""Empty chain to assist with unit testing depletion routines
|
||||
|
||||
Only really provides the form_matrix function, but acts like
|
||||
a real Chain
|
||||
"""
|
||||
|
||||
fission_yields = [None]
|
||||
|
||||
@staticmethod
|
||||
def get_thermal_fission_yields():
|
||||
return None
|
||||
|
||||
def form_matrix(self, rates, _fission_yields=None):
|
||||
"""Forms the f(y) matrix in y' = f(y)y.
|
||||
|
||||
Nominally a depletion matrix, this is abstracted on the off chance
|
||||
that the function f has nothing to do with depletion at all.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rates : numpy.ndarray
|
||||
Slice of reaction rates for a single material
|
||||
_fission_yields : optional
|
||||
Not used
|
||||
|
||||
Returns
|
||||
-------
|
||||
scipy.sparse.csr_matrix
|
||||
Sparse matrix representing f(y).
|
||||
"""
|
||||
|
||||
y_1 = rates[0, 0]
|
||||
y_2 = rates[1, 0]
|
||||
|
||||
a11 = np.sin(y_2)
|
||||
a12 = np.cos(y_1)
|
||||
a21 = -np.cos(y_2)
|
||||
a22 = np.sin(y_1)
|
||||
|
||||
return sp.csr_matrix(np.array([[a11, a12], [a21, a22]]))
|
||||
|
||||
|
||||
class DummyOperator(TransportOperator):
|
||||
"""This is a dummy operator class with no statistical uncertainty.
|
||||
|
||||
|
|
@ -88,6 +131,7 @@ class DummyOperator(TransportOperator):
|
|||
"""
|
||||
def __init__(self, previous_results=None):
|
||||
self.prev_res = previous_results
|
||||
self.chain = TestChain()
|
||||
self.output_dir = "."
|
||||
|
||||
def __call__(self, vec, power, print_out=False):
|
||||
|
|
@ -120,37 +164,6 @@ class DummyOperator(TransportOperator):
|
|||
# Create a fake rates object
|
||||
return OperatorResult(ufloat(0.0, 0.0), reaction_rates)
|
||||
|
||||
@property
|
||||
def chain(self):
|
||||
return self
|
||||
|
||||
def form_matrix(self, rates):
|
||||
"""Forms the f(y) matrix in y' = f(y)y.
|
||||
|
||||
Nominally a depletion matrix, this is abstracted on the off chance
|
||||
that the function f has nothing to do with depletion at all.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rates : numpy.ndarray
|
||||
Slice of reaction rates for a single material
|
||||
|
||||
Returns
|
||||
-------
|
||||
scipy.sparse.csr_matrix
|
||||
Sparse matrix representing f(y).
|
||||
"""
|
||||
|
||||
y_1 = rates[0, 0]
|
||||
y_2 = rates[1, 0]
|
||||
|
||||
a11 = np.sin(y_2)
|
||||
a12 = np.cos(y_1)
|
||||
a21 = -np.cos(y_2)
|
||||
a22 = np.sin(y_1)
|
||||
|
||||
return sp.csr_matrix(np.array([[a11, a12], [a21, a22]]))
|
||||
|
||||
@property
|
||||
def volume(self):
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -3,10 +3,11 @@
|
|||
from collections.abc import Mapping
|
||||
import os
|
||||
from pathlib import Path
|
||||
from itertools import product
|
||||
|
||||
import numpy as np
|
||||
from openmc.data import zam, ATOMIC_SYMBOL
|
||||
from openmc.deplete import comm, Chain, reaction_rates, nuclide
|
||||
from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram
|
||||
import pytest
|
||||
|
||||
from tests import cdtemp
|
||||
|
|
@ -128,9 +129,10 @@ def test_from_xml(simple_chain):
|
|||
assert [r.branching_ratio for r in nuc.reactions] == [1.0, 0.7, 0.3]
|
||||
|
||||
# Yield tests
|
||||
assert nuc.yield_energies == [0.0253]
|
||||
assert nuc.yield_energies == (0.0253,)
|
||||
assert list(nuc.yield_data) == [0.0253]
|
||||
assert nuc.yield_data[0.0253] == [("A", 0.0292737), ("B", 0.002566345)]
|
||||
assert nuc.yield_data[0.0253].products == ("A", "B")
|
||||
assert (nuc.yield_data[0.0253].yields == [0.0292737, 0.002566345]).all()
|
||||
|
||||
|
||||
def test_export_to_xml(run_in_tmpdir):
|
||||
|
|
@ -139,8 +141,7 @@ def test_export_to_xml(run_in_tmpdir):
|
|||
# Prevent different MPI ranks from conflicting
|
||||
filename = 'test{}.xml'.format(comm.rank)
|
||||
|
||||
A = nuclide.Nuclide()
|
||||
A.name = "A"
|
||||
A = nuclide.Nuclide("A")
|
||||
A.half_life = 2.36520e4
|
||||
A.decay_modes = [
|
||||
nuclide.DecayTuple("beta1", "B", 0.6),
|
||||
|
|
@ -148,21 +149,19 @@ def test_export_to_xml(run_in_tmpdir):
|
|||
]
|
||||
A.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)]
|
||||
|
||||
B = nuclide.Nuclide()
|
||||
B.name = "B"
|
||||
B = nuclide.Nuclide("B")
|
||||
B.half_life = 3.29040e4
|
||||
B.decay_modes = [nuclide.DecayTuple("beta", "A", 1.0)]
|
||||
B.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)]
|
||||
|
||||
C = nuclide.Nuclide()
|
||||
C.name = "C"
|
||||
C = nuclide.Nuclide("C")
|
||||
C.reactions = [
|
||||
nuclide.ReactionTuple("fission", None, 2.0e8, 1.0),
|
||||
nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7),
|
||||
nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3)
|
||||
]
|
||||
C.yield_energies = [0.0253]
|
||||
C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
|
||||
C.yield_data = nuclide.FissionYieldDistribution({
|
||||
0.0253: {"A": 0.0292737, "B": 0.002566345}})
|
||||
|
||||
chain = Chain()
|
||||
chain.nuclides = [A, B, C]
|
||||
|
|
@ -220,6 +219,13 @@ def test_form_matrix(simple_chain):
|
|||
assert mat[1, 2] == mat12
|
||||
assert mat[2, 2] == mat22
|
||||
|
||||
# Pass equivalent fission yields directly
|
||||
# Ensure identical matrix is formed
|
||||
f_yields = {"C": {"A": 0.0292737, "B": 0.002566345}}
|
||||
new_mat = chain.form_matrix(react[0], f_yields)
|
||||
for r, c in product(range(3), range(3)):
|
||||
assert new_mat[r, c] == mat[r, c]
|
||||
|
||||
|
||||
def test_getitem():
|
||||
""" Test nuc_by_ind converter function. """
|
||||
|
|
@ -283,8 +289,7 @@ def test_capture_branch_infer_ground():
|
|||
chain = Chain.from_xml(chain_file)
|
||||
|
||||
# Create nuclide to be added into the chain
|
||||
xe136m = nuclide.Nuclide()
|
||||
xe136m.name = "Xe136_m1"
|
||||
xe136m = nuclide.Nuclide("Xe136_m1")
|
||||
|
||||
chain.nuclides.append(xe136m)
|
||||
chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1
|
||||
|
|
@ -301,8 +306,7 @@ def test_capture_branch_no_rxn():
|
|||
chain_file = Path(__file__).parents[1] / "chain_simple.xml"
|
||||
chain = Chain.from_xml(chain_file)
|
||||
|
||||
u5m = nuclide.Nuclide()
|
||||
u5m.name = "U235_m1"
|
||||
u5m = nuclide.Nuclide("U235_m1")
|
||||
|
||||
chain.nuclides.append(u5m)
|
||||
chain.nuclide_dict[u5m.name] = len(chain.nuclides) - 1
|
||||
|
|
@ -378,6 +382,32 @@ def test_set_alpha_branches():
|
|||
raise ValueError("Helium has been removed and should not have been")
|
||||
|
||||
|
||||
def test_simple_fission_yields(simple_chain):
|
||||
"""Check the default fission yields that can be used to form the matrix
|
||||
"""
|
||||
fission_yields = simple_chain.get_thermal_fission_yields()
|
||||
assert fission_yields == {"C": {"A": 0.0292737, "B": 0.002566345}}
|
||||
|
||||
|
||||
def test_fission_yield_attribute(simple_chain):
|
||||
"""Test the fission_yields property"""
|
||||
thermal_yields = simple_chain.get_thermal_fission_yields()
|
||||
# generate default with property
|
||||
assert simple_chain.fission_yields[0] == thermal_yields
|
||||
empty_chain = Chain()
|
||||
empty_chain.fission_yields = thermal_yields
|
||||
assert empty_chain.fission_yields[0] == thermal_yields
|
||||
empty_chain.fission_yields = [thermal_yields] * 2
|
||||
assert empty_chain.fission_yields[0] == thermal_yields
|
||||
assert empty_chain.fission_yields[1] == thermal_yields
|
||||
|
||||
# test failure with deplete function
|
||||
# number fission yields != number of materials
|
||||
dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1)
|
||||
with pytest.raises(
|
||||
ValueError, match="fission yield.*not equal.*compositions"):
|
||||
cram.deplete(empty_chain, dummy_conc, None, 0.5)
|
||||
|
||||
def test_validate(simple_chain):
|
||||
"""Test the validate method"""
|
||||
|
||||
|
|
@ -394,7 +424,7 @@ def test_validate(simple_chain):
|
|||
|
||||
# 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)]}
|
||||
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:
|
||||
|
|
|
|||
293
tests/unit_tests/test_deplete_fission_yields.py
Normal file
293
tests/unit_tests/test_deplete_fission_yields.py
Normal file
|
|
@ -0,0 +1,293 @@
|
|||
"""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 capi
|
||||
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("capi")
|
||||
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.rectangular_prism(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 capi.run_in_memory():
|
||||
yield [capi.Material(), capi.Material()]
|
||||
finally:
|
||||
# Convert to strings as os.remove in py 3.5 doesn't support Paths
|
||||
for file_path in ("settings.xml", "geometry.xml", "materials.xml",
|
||||
"summary.h5"):
|
||||
os.remove(str(tmpdir / file_path))
|
||||
orig.chdir()
|
||||
os.rmdir(str(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 = np.empty((n_bins, n_nucs * n_scores, 3))
|
||||
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], 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
|
||||
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
|
||||
|
||||
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], 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 = 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
|
||||
|
|
@ -2,8 +2,8 @@
|
|||
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
import numpy
|
||||
import pytest
|
||||
|
||||
from openmc.deplete import nuclide
|
||||
|
||||
|
||||
|
|
@ -71,18 +71,20 @@ def test_from_xml():
|
|||
nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0),
|
||||
nuclide.ReactionTuple('fission', None, 193405400.0, 1.0),
|
||||
]
|
||||
assert u235.yield_energies == [0.0253]
|
||||
assert u235.yield_data == {
|
||||
0.0253: [('Te134', 0.062155), ('Zr100', 0.0497641),
|
||||
('Xe138', 0.0481413)]
|
||||
}
|
||||
expected_yield_data = nuclide.FissionYieldDistribution({
|
||||
0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}})
|
||||
assert u235.yield_data == expected_yield_data
|
||||
# test accessing the yield energies through the FissionYieldDistribution
|
||||
assert u235.yield_energies == (0.0253,)
|
||||
assert u235.yield_energies is u235.yield_data.energies
|
||||
with pytest.raises(AttributeError): # not settable
|
||||
u235.yield_energies = [0.0253, 5e5]
|
||||
|
||||
|
||||
def test_to_xml_element():
|
||||
"""Test writing nuclide data to an XML element."""
|
||||
|
||||
C = nuclide.Nuclide()
|
||||
C.name = "C"
|
||||
C = nuclide.Nuclide("C")
|
||||
C.half_life = 0.123
|
||||
C.decay_modes = [
|
||||
nuclide.DecayTuple('beta-', 'B', 0.99),
|
||||
|
|
@ -92,8 +94,8 @@ def test_to_xml_element():
|
|||
nuclide.ReactionTuple('fission', None, 2.0e8, 1.0),
|
||||
nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0)
|
||||
]
|
||||
C.yield_energies = [0.0253]
|
||||
C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
|
||||
C.yield_data = nuclide.FissionYieldDistribution(
|
||||
{0.0253: {"A": 0.0292737, "B": 0.002566345}})
|
||||
element = C.to_xml_element()
|
||||
|
||||
assert element.get("half_life") == "0.123"
|
||||
|
|
@ -118,6 +120,57 @@ def test_to_xml_element():
|
|||
assert element.find('neutron_fission_yields') is not None
|
||||
|
||||
|
||||
def test_fission_yield_distribution():
|
||||
"""Test an energy-dependent yield distribution"""
|
||||
yield_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, "Sm149": 2.69e-8},
|
||||
5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}, # drop Sm149
|
||||
}
|
||||
yield_dist = nuclide.FissionYieldDistribution(yield_dict)
|
||||
assert len(yield_dist) == len(yield_dict)
|
||||
assert yield_dist.energies == tuple(sorted(yield_dict.keys()))
|
||||
for exp_ene, exp_dist in yield_dict.items():
|
||||
act_dist = yield_dict[exp_ene]
|
||||
for exp_prod, exp_yield in exp_dist.items():
|
||||
assert act_dist[exp_prod] == exp_yield
|
||||
exp_yield = numpy.array([
|
||||
[4.08e-12, 1.71e-12, 7.85e-4],
|
||||
[1.32e-12, 0.0, 1.12e-3],
|
||||
[5.83e-8, 2.69e-8, 4.54e-3]])
|
||||
assert numpy.array_equal(yield_dist.yield_matrix, exp_yield)
|
||||
|
||||
# Test the operations / special methods for fission yield
|
||||
orig_yields = yield_dist[0.0253]
|
||||
assert len(orig_yields) == len(yield_dict[0.0253])
|
||||
for key, value in yield_dict[0.0253].items():
|
||||
assert key in orig_yields
|
||||
assert orig_yields[key] == value
|
||||
# __getitem__ return yields as a view into yield matrix
|
||||
assert orig_yields.yields.base is yield_dist.yield_matrix
|
||||
|
||||
# Fission yield feature uses scaled and incremented
|
||||
mod_yields = orig_yields * 2
|
||||
assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields)
|
||||
mod_yields += orig_yields
|
||||
assert numpy.array_equal(orig_yields.yields * 3, mod_yields.yields)
|
||||
|
||||
# Failure modes for adding, multiplying yields
|
||||
similar = numpy.empty_like(orig_yields.yields)
|
||||
with pytest.raises(TypeError):
|
||||
orig_yields + similar
|
||||
with pytest.raises(TypeError):
|
||||
similar + orig_yields
|
||||
with pytest.raises(TypeError):
|
||||
orig_yields += similar
|
||||
with pytest.raises(TypeError):
|
||||
orig_yields * similar
|
||||
with pytest.raises(TypeError):
|
||||
similar * orig_yields
|
||||
with pytest.raises(TypeError):
|
||||
orig_yields *= similar
|
||||
|
||||
|
||||
def test_validate():
|
||||
|
||||
nuc = nuclide.Nuclide()
|
||||
|
|
@ -141,8 +194,8 @@ def test_validate():
|
|||
# fission yields
|
||||
|
||||
nuc.yield_data = {
|
||||
0.0253: [("0", 1.5), ("1", 0.5)],
|
||||
1e6: [("0", 1.5), ("1", 0.5)],
|
||||
0.0253: {"0": 1.5, "1": 0.5},
|
||||
1e6: {"0": 1.5, "1": 0.5},
|
||||
}
|
||||
|
||||
# nuclide is good and should have no warnings raise
|
||||
|
|
@ -176,7 +229,7 @@ def test_validate():
|
|||
|
||||
# restore reactions, invalidate fission yields
|
||||
nuc.reactions.append(reaction)
|
||||
nuc.yield_data[1e6].pop()
|
||||
nuc.yield_data[1e6].yields *= 2
|
||||
|
||||
with pytest.raises(ValueError, match=r"fission yields.*1\.0*e"):
|
||||
nuc.validate(strict=True, quiet=False, tolerance=0.0)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue