From e70c5d539f33f48049d8af1ef1e303e34cf41a10 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 5 Jul 2019 13:09:30 -0500 Subject: [PATCH] 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. --- openmc/deplete/abc.py | 117 +++++++++++++++++++++++++++++++- openmc/deplete/helpers.py | 117 ++++++++++++++++++++++++++++++++ openmc/deplete/tally_helpers.py | 110 ------------------------------ 3 files changed, 232 insertions(+), 112 deletions(-) create mode 100644 openmc/deplete/helpers.py delete mode 100644 openmc/deplete/tally_helpers.py 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)