Refactor tally_helpers to Rates and Energy Helpers

The TallyHelperBase class has been split into two
abstract classes: ReactionRateHelper and FissionEnergyHelper.
The former is responsible for creating reaction rates
for all materials, given nuclides with non-zero densities from
the Operator.
The latter is responsible for producing recoverable fission
energy for each nuclide in each material tracked.

Two concrete classes are included to make the Operator functional.
These can be imported from openmc/deplete/helpers.py.
First, the DirectRxnRateHelper preserves the old mode for
directly tallying one-group reaction rates using the OpenMC C API.
The second, ChainFissHelper, generates the fission energy
vector from the chain data from before.
This commit is contained in:
Andrew Johnson 2019-07-05 13:09:30 -05:00
parent b2f1eef449
commit e70c5d539f
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
3 changed files with 232 additions and 112 deletions

View file

@ -7,11 +7,14 @@ to run a full depletion simulation.
from collections import namedtuple
import os
from pathlib import Path
from abc import ABCMeta, abstractmethod
from abc import ABC, abstractmethod
from xml.etree import ElementTree as ET
from warnings import warn
from numpy import zeros, nonzero
from openmc.data import DataLibrary
from openmc.checkvalue import check_type
from .chain import Chain
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
@ -34,7 +37,7 @@ except AttributeError:
pass
class TransportOperator(metaclass=ABCMeta):
class TransportOperator(ABC):
"""Abstract class defining a transport operator
Each depletion integrator is written to work with a generic transport
@ -160,3 +163,113 @@ class TransportOperator(metaclass=ABCMeta):
def finalize(self):
pass
class ReactionRateHelper(ABC):
"""Abstract class for generating reaction rates for operators
Responsible for generating reaction rate tallies for burnable
materials, given nuclides and scores from the operator.
Reaction rates are passed back to the operator for be used in
an :class:`openmc.deplete.OperatorResult` instance
"""
def __init__(self):
self._nuclides = None
self._rate_tally = None
self._results_cache = None
@abstractmethod
def generate_tallies(self, materials, scores):
"""Use the capi to build tallies needed for reaction rates"""
pass
@property
def nuclides(self):
"""List of nuclides with requested reaction rates"""
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
self._rate_tally.nuclides = nuclides
def _reset_results_cache(self, nnucs, nreact):
"""Cache for results for a given material"""
if self._results_cache is None or self._results_cache.shape != (nnucs, nreact):
self._results_cache = zeros((nnucs, nreact))
else:
self._results_cache.fill(0.0)
return self._results_cache
@abstractmethod
def get_material_rates(self, mat_id, nuc_index, react_index):
"""Return 2D array of [nuclide, reaction] reaction rates
``nuc_index`` and ``react_index`` are orderings of nuclides
and reactions such that the ordering is consistent between
reaction tallies and energy deposition tallies"""
pass
def divide_by_adens(self, number):
"""Normalize reaction rates by number of nuclides
Acts on the current material examined by
:meth:`get_material_rates`
Parameters
----------
number : iterable of float
Number density [#/b/cm] of each nuclide tracked in the calculation.
Ordered identically to :attr:`nuclides`
Returns
-------
results : :class:`numpy.ndarray`
2D array ``[n_nuclides, n_rxns]`` of reaction rates normalized by
the number of nuclides
"""
mask = nonzero(number)
results = self._results_cache
for col in range(results.shape[1]):
results[mask, col] /= number[mask]
return results
class FissionEnergyHelper(ABC):
"""Abstract class for normalizing fission reactions to a given level
"""
def __init__(self):
self._nuclides = None
self._fission_E = None
@abstractmethod
def prepare(self, chain_nucs, rate_index, materials):
"""Perform work needed to obtain fission energy per material
``chain_nucs`` is all nuclides tracked in the depletion chain,
while ``rate_index`` should be a mapping from nuclide name
to index in the reaction rate vector used in
:meth:`get_fission_energy`.
``materials`` should be a list of all materials tracked
on the operator to which this object is attached"""
pass
@abstractmethod
def get_fission_energy(self, fission_rates, mat_index):
"""Return fission energy in this material given fission rates"""
pass
@property
def nuclides(self):
"""List of nuclides with requested reaction rates"""
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides

117
openmc/deplete/helpers.py Normal file
View file

@ -0,0 +1,117 @@
"""
Class for normalizing fission energy deposition
"""
from itertools import product
from numpy import dot, zeros
from openmc.capi import Tally, MaterialFilter
from .abc import ReactionRateHelper, FissionEnergyHelper
# -------------------------------------
# Helpers for generating reaction rates
# -------------------------------------
class DirectRxnRateHelper(ReactionRateHelper):
"""Class that generates tallies for one-group rates"""
def generate_tallies(self, materials, scores):
"""Produce one-group reaction rate tally
Uses the :mod:`openmc.capi` to generate a tally
of relevant reactions across all burnable materials.
Parameters
----------
materials : iterable of :class:`openmc.Material`
Burnable materials in the problem. Used to
construct a :class:`openmc.MaterialFilter`
scores : iterable of str
Reaction identifiers, e.g. ``"(n, fission)"``,
``"(n, gamma)"``, needed for the reaction rate tally.
"""
self._rate_tally = Tally()
self._rate_tally.scores = scores
self._rate_tally.filters = [MaterialFilter(materials)]
def get_material_rates(self, mat_id, nuc_index, react_index):
"""Return an array of reaction rates for a material
Parameters
----------
mat_id : int
Unique id for the requested material
nuc_index : iterable of int
Index for each nuclide in :attr:`nuclides` in the
desired reaction rate matrix
react_index : iterable of int
Index for each reaction scored in the tally
Returns
-------
rates : :class:`numpy.ndarray`
2D matrix ``(len(nuc_index), len(react_index))`` with the
reaction rates in this material
"""
results = self._reset_results_cache(len(nuc_index), len(react_index))
full_tally_res = self._rate_tally.results[mat_id, :, 1]
for i_tally, (i_nuc, i_react) in enumerate(
product(nuc_index, react_index)):
results[i_nuc, i_react] = full_tally_res[i_tally]
return results
# ------------------------------------
# Helpers for obtaining fission energy
# ------------------------------------
class ChainFissHelper(FissionEnergyHelper):
"""Fission Q-values are pulled from chain"""
def prepare(self, chain_nucs, rate_index, _materials):
"""Populate the fission Q value vector from a chain.
Paramters
---------
chain_nucs : iterable of :class:`openmc.deplete.Nuclide`
Nuclides used in this depletion chain. Do not need
to be ordered
rate_index : dict of str to int
Dictionary mapping names of nuclides, e.g. ``"U235"``,
to a corresponding index in the desired fission Q
vector.
_materials : list of str
Unused. Materials to be tracked for this helper.
"""
if (self._fission_E is not None
and self._fission_E.shape == (len(rate_index),)):
return
fiss_E = zeros(len(rate_index))
for nuclide in chain_nucs:
if nuclide.name in rate_index:
for rx in nuclide.reactions:
if rx.type == "fission":
fiss_E[rate_index[nuclide.name]] = rx.Q
break
self._fission_E = fiss_E
def get_fission_energy(self, fiss_rates, _mat_index):
"""Return a vector of the isotopic fission energy for this material
parameters
----------
fission_rates : numpy.ndarray
fission reaction rate for each isotope in the specified
material. should be ordered corresponding to initial
``rate_index`` used in :meth:`set_fission_q`
_mat_index : int
index for the material requested. Unused, as all
isotopes in all materials have the same Q value.
"""
return dot(fiss_rates, self._fission_E)

View file

@ -1,110 +0,0 @@
"""
Class for normalizing fission energy deposition
"""
from abc import ABC, abstractmethod
from numpy import dot, zeros
from openmc.checkvalue import check_type
from openmc.capi import Tally, MaterialFilter
class TallyHelperBase(ABC):
"""Base class for working with tallies for depletion"""
def __init__(self):
self._fiss_q = None
self._rx_tally = None
self._nuclides = []
@abstractmethod
def set_fission_q(self, chain_nucs, rate_index):
"""Populate the energy released per fission Q value array"""
pass
@property
def reaction_tally(self):
if self._rx_tally is None:
raise AttributeError(
"Reaction tally for {} not set.".format(
self.__class__.__name__
)
)
return self._rx_tally
def generate_tallies(self, materials, scores):
self._rx_tally = Tally()
self._rx_tally.scores = scores
self._rx_tally.filters = [MaterialFilter(materials)]
@property
def nuclides(self):
"""List of nuclides with requested reaction rates"""
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
self._rx_tally.nuclides = nuclides
@abstractmethod
def get_fission_energy(self, fission_rates, mat_index):
"""return a vector of the isotopic fission energy for this material
parameters
----------
fission_rates: numpy.ndarray
fission reaction rate for each isotope in the specified
material. should be ordered corresponding to initial
``rate_index`` used in :meth:`set_fission_q`
mat_index: int
index for the material requested.
"""
class ChainFissTallyHelper(TallyHelperBase):
"""Fission Q-values are pulled from chain"""
def set_fission_q(self, chain_nucs, rate_index):
"""Populate the fission Q value vector from a chain.
Paramters
---------
chain_nucs: iterable of :class:`openmc.deplete.Nuclide`
Nuclides used in this depletion chain. Do not need
to be ordered
rate_index: dict of str to int
Dictionary mapping names of nuclides, e.g. ``"U235"``,
to a corresponding index in the desired fission Q
vector.
"""
if (self._fiss_q is not None
and self._fiss_q.shape == (len(rate_index), )):
return
fq = zeros(len(rate_index))
for nuclide in chain_nucs:
if nuclide.name in rate_index:
for rx in nuclide.reactions:
if rx.type == "fission":
fq[rate_index[nuclide.name]] = rx.Q
break
self._fiss_q = fq
def get_fission_energy(self, fiss_rates, _mat_index):
"""Return a vector of the isotopic fission energy for this material
parameters
----------
fission_rates: numpy.ndarray
fission reaction rate for each isotope in the specified
material. should be ordered corresponding to initial
``rate_index`` used in :meth:`set_fission_q`
mat_index: int
index for the material requested. Unused, as all
isotopes in all materials have the same Q value.
"""
return dot(fiss_rates, self._fiss_q)