diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 1be3e27ff9..4a332af590 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -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 diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py new file mode 100644 index 0000000000..4b79d5b3d3 --- /dev/null +++ b/openmc/deplete/helpers.py @@ -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) diff --git a/openmc/deplete/tally_helpers.py b/openmc/deplete/tally_helpers.py deleted file mode 100644 index d61bfb90b6..0000000000 --- a/openmc/deplete/tally_helpers.py +++ /dev/null @@ -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)