From c44154b5f1ae0e89d3c48aba4bb5b22e4c8c67bf Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 1 Jul 2019 08:33:55 -0500
Subject: [PATCH 01/14] Tally helper for helping operator with tallies, q
 values

openmc.deplete.tally_helers.ChainFissTallyHelper is responsible
for populating the fission_q vector, building the reaction rate
tallies, and updating the energy produced by fission per material.

This abstraction will be helpful as more methods for pulling
fission q values are introduced, namely as issue #1238, indirect
energy release, gets resolved
---
 openmc/deplete/operator.py      | 49 ++++++++++----------------------
 openmc/deplete/tally_helpers.py | 50 +++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 34 deletions(-)
 create mode 100644 openmc/deplete/tally_helpers.py

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index ffbffbc73..3976fe7c9 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -24,6 +24,7 @@ from . import comm
 from .abc import TransportOperator, OperatorResult
 from .atom_number import AtomNumber
 from .reaction_rates import ReactionRates
+from .tally_helpers import ChainFissTallyHelper
 
 
 def _distribute(items):
@@ -152,6 +153,10 @@ class Operator(TransportOperator):
         self.reaction_rates = ReactionRates(
             self.local_mats, self._burnable_nucs, self.chain.reactions)
 
+        # Get class to assist working with tallies
+        self._tally_helper = ChainFissTallyHelper(len(self.local_mats))
+
+
     def __call__(self, vec, power, print_out=True):
         """Runs a simulation.
 
@@ -180,7 +185,7 @@ class Operator(TransportOperator):
 
         # Update material compositions and tally nuclides
         self._update_materials()
-        self._tally.nuclides = self._get_tally_nuclides()
+        self._tally_helper.reaction_tally.nuclides = self._get_tally_nuclides()
 
         # Run OpenMC
         openmc.capi.reset()
@@ -373,7 +378,9 @@ class Operator(TransportOperator):
         openmc.capi.init(intracomm=comm)
 
         # Generate tallies in memory
-        self._generate_tallies()
+        materials = [openmc.capi.materials[int(i)]
+                     for i in self.burnable_mats]
+        self._tally_helper.generate_tallies(materials, self.chain.reactions)
 
         # Return number density vector
         return list(self.number.get_mat_slice(np.s_[:]))
@@ -482,27 +489,6 @@ class Operator(TransportOperator):
         nuc_list = comm.bcast(nuc_list)
         return [nuc for nuc in nuc_list if nuc in self.chain]
 
-    def _generate_tallies(self):
-        """Generates depletion tallies.
-
-        Using information from the depletion chain as well as the nuclides
-        currently in the problem, this function automatically generates a
-        tally.xml for the simulation.
-
-        """
-        # Create tallies for depleting regions
-        materials = [openmc.capi.materials[int(i)]
-                     for i in self.burnable_mats]
-        mat_filter = openmc.capi.MaterialFilter(materials)
-
-        # Set up a tally that has a material filter covering each depletable
-        # material and scores corresponding to all reactions that cause
-        # transmutation. The nuclides for the tally are set later when eval() is
-        # called.
-        self._tally = openmc.capi.Tally()
-        self._tally.scores = self.chain.reactions
-        self._tally.filters = [mat_filter]
-
     def _unpack_tallies_and_normalize(self, power):
         """Unpack tallies from OpenMC and return an operator result
 
@@ -529,7 +515,7 @@ class Operator(TransportOperator):
 
         # Extract tally bins
         materials = self.burnable_mats
-        nuclides = self._tally.nuclides
+        nuclides = self._tally_helper.reaction_tally.nuclides
 
         # Form fast map
         nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
@@ -544,19 +530,14 @@ class Operator(TransportOperator):
 
         # Create arrays to store fission Q values, reaction rates, and nuclide
         # numbers
-        fission_Q = np.zeros(rates.n_nuc)
         rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
         number = np.zeros(rates.n_nuc)
 
         fission_ind = rates.index_rx["fission"]
 
-        for nuclide in self.chain.nuclides:
-            if nuclide.name in rates.index_nuc:
-                for rx in nuclide.reactions:
-                    if rx.type == 'fission':
-                        ind = rates.index_nuc[nuclide.name]
-                        fission_Q[ind] = rx.Q
-                        break
+        self._tally_helper.set_fission_q(self.chain.nuclides, rates.index_nuc)
+
+        tally_results = self._tally_helper.reaction_tally.results
 
         # Extract results
         for i, mat in enumerate(self.local_mats):
@@ -564,7 +545,7 @@ class Operator(TransportOperator):
             slab = materials.index(mat)
 
             # Get material results hyperslab
-            results = self._tally.results[slab, :, 1]
+            results = tally_results[slab, :, 1]
 
             # Zero out reaction rates and nuclide numbers
             rates_expanded[:] = 0.0
@@ -579,7 +560,7 @@ class Operator(TransportOperator):
                     j += 1
 
             # Accumulate energy from fission
-            energy += np.dot(rates_expanded[:, fission_ind], fission_Q)
+            energy += self._tally_helper.get_fiss_energy(rates_expanded[:, fission_ind], i)
 
             # Divide by total number and store
             for i_nuc_results in nuc_ind:
diff --git a/openmc/deplete/tally_helpers.py b/openmc/deplete/tally_helpers.py
new file mode 100644
index 000000000..47c78ba8e
--- /dev/null
+++ b/openmc/deplete/tally_helpers.py
@@ -0,0 +1,50 @@
+"""
+Class for normalizing fission energy deposition
+"""
+
+from numpy import dot, zeros
+
+from openmc.capi import Tally, MaterialFilter
+
+
+class ChainFissTallyHelper(object):
+    """Fission Q-values are pulled from chain"""
+
+    def __init__(self, n_materials):
+        self._fiss_q = None
+        self.n_materials = n_materials
+        self._rx_tally = None
+
+    def set_fission_q(self, chain_nucs, rate_index):
+        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
+
+    @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)]
+
+    def get_fiss_energy(self, fiss_rates, _mat_index):
+        return dot(fiss_rates, self._fiss_q)

From 888086ac297d46e244b9b61d3e01c7de33bf6803 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 1 Jul 2019 09:22:01 -0500
Subject: [PATCH 02/14] Add base class for TallyHelpers

Interface is designed such that one could make a tally
helper work with unique materials, mainly for pulling
fission Q values on a per material basis
---
 openmc/deplete/operator.py      |  9 +--
 openmc/deplete/tally_helpers.py | 98 ++++++++++++++++++++++++++-------
 2 files changed, 84 insertions(+), 23 deletions(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 3976fe7c9..9d3f6035e 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -154,7 +154,7 @@ class Operator(TransportOperator):
             self.local_mats, self._burnable_nucs, self.chain.reactions)
 
         # Get class to assist working with tallies
-        self._tally_helper = ChainFissTallyHelper(len(self.local_mats))
+        self._tally_helper = ChainFissTallyHelper()
 
 
     def __call__(self, vec, power, print_out=True):
@@ -185,7 +185,7 @@ class Operator(TransportOperator):
 
         # Update material compositions and tally nuclides
         self._update_materials()
-        self._tally_helper.reaction_tally.nuclides = self._get_tally_nuclides()
+        self._tally_helper.nuclides = self._get_tally_nuclides()
 
         # Run OpenMC
         openmc.capi.reset()
@@ -515,7 +515,7 @@ class Operator(TransportOperator):
 
         # Extract tally bins
         materials = self.burnable_mats
-        nuclides = self._tally_helper.reaction_tally.nuclides
+        nuclides = self._tally_helper.nuclides
 
         # Form fast map
         nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
@@ -560,7 +560,8 @@ class Operator(TransportOperator):
                     j += 1
 
             # Accumulate energy from fission
-            energy += self._tally_helper.get_fiss_energy(rates_expanded[:, fission_ind], i)
+            energy += self._tally_helper.get_fission_energy(
+                rates_expanded[:, fission_ind], i)
 
             # Divide by total number and store
             for i_nuc_results in nuc_ind:
diff --git a/openmc/deplete/tally_helpers.py b/openmc/deplete/tally_helpers.py
index 47c78ba8e..d61bfb90b 100644
--- a/openmc/deplete/tally_helpers.py
+++ b/openmc/deplete/tally_helpers.py
@@ -1,35 +1,26 @@
 """
 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 ChainFissTallyHelper(object):
-    """Fission Q-values are pulled from chain"""
+class TallyHelperBase(ABC):
+    """Base class for working with tallies for depletion"""
 
-    def __init__(self, n_materials):
+    def __init__(self):
         self._fiss_q = None
-        self.n_materials = n_materials
         self._rx_tally = None
+        self._nuclides = []
 
+    @abstractmethod
     def set_fission_q(self, chain_nucs, rate_index):
-        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
+        """Populate the energy released per fission Q value array"""
+        pass
 
     @property
     def reaction_tally(self):
@@ -46,5 +37,74 @@ class ChainFissTallyHelper(object):
         self._rx_tally.scores = scores
         self._rx_tally.filters = [MaterialFilter(materials)]
 
-    def get_fiss_energy(self, fiss_rates, _mat_index):
+    @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)

From e70c5d539f33f48049d8af1ef1e303e34cf41a10 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Fri, 5 Jul 2019 13:09:30 -0500
Subject: [PATCH 03/14] 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 1be3e27ff..4a332af59 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 000000000..4b79d5b3d
--- /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 d61bfb90b..000000000
--- 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)

From 4832708ad79cc7c9dc7e5a40e893dc6ae54ed071 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Fri, 5 Jul 2019 13:28:32 -0500
Subject: [PATCH 04/14] Apply new reaction rate, fission energy helpers to
 Operator

rate_helper responsible for passing reaction rates
onto the operator, while energy_helper is responsible
for computing the actual fission energy produced in each
material
---
 openmc/deplete/abc.py      |  3 ++-
 openmc/deplete/operator.py | 47 +++++++++++++++-----------------------
 2 files changed, 20 insertions(+), 30 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 4a332af59..13be50812 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -198,7 +198,8 @@ class ReactionRateHelper(ABC):
 
     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):
+        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)
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index f31978908..fab2a9f48 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -24,7 +24,7 @@ from . import comm
 from .abc import TransportOperator, OperatorResult
 from .atom_number import AtomNumber
 from .reaction_rates import ReactionRates
-from .tally_helpers import ChainFissTallyHelper
+from .helpers import DirectRxnRateHelper, ChainFissHelper
 
 
 def _distribute(items):
@@ -154,7 +154,8 @@ class Operator(TransportOperator):
             self.local_mats, self._burnable_nucs, self.chain.reactions)
 
         # Get class to assist working with tallies
-        self._tally_helper = ChainFissTallyHelper()
+        self._rate_helper = DirectRxnRateHelper()
+        self._energy_helper = ChainFissHelper()
 
 
     def __call__(self, vec, power, print_out=True):
@@ -185,7 +186,8 @@ class Operator(TransportOperator):
 
         # Update material compositions and tally nuclides
         self._update_materials()
-        self._tally_helper.nuclides = self._get_tally_nuclides()
+        self._rate_helper.nuclides = self._get_tally_nuclides()
+        self._energy_helper.nuclides = self._rate_helper.nuclides
 
         # Run OpenMC
         openmc.capi.reset()
@@ -380,7 +382,9 @@ class Operator(TransportOperator):
         # Generate tallies in memory
         materials = [openmc.capi.materials[int(i)]
                      for i in self.burnable_mats]
-        self._tally_helper.generate_tallies(materials, self.chain.reactions)
+        self._rate_helper.generate_tallies(materials, self.chain.reactions)
+        self._energy_helper.prepare(
+            self.chain.nuclides, self.reaction_rates.index_nuc, materials)
 
         # Return number density vector
         return list(self.number.get_mat_slice(np.s_[:]))
@@ -516,7 +520,7 @@ class Operator(TransportOperator):
 
         # Extract tally bins
         materials = self.burnable_mats
-        nuclides = self._tally_helper.nuclides
+        nuclides = self._rate_helper.nuclides
 
         # Form fast map
         nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
@@ -530,47 +534,32 @@ class Operator(TransportOperator):
         energy = 0.0
 
         # Create arrays to store fission Q values, reaction rates, and nuclide
-        # numbers
-        rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
-        number = np.zeros(rates.n_nuc)
+        # numbers, zeroed out in material iteration
+        number = np.empty(rates.n_nuc)
 
         fission_ind = rates.index_rx["fission"]
 
-        self._tally_helper.set_fission_q(self.chain.nuclides, rates.index_nuc)
-
-        tally_results = self._tally_helper.reaction_tally.results
-
         # Extract results
         for i, mat in enumerate(self.local_mats):
             # Get tally index
             slab = materials.index(mat)
 
-            # Get material results hyperslab
-            results = tally_results[slab, :, 1]
-
             # Zero out reaction rates and nuclide numbers
-            rates_expanded[:] = 0.0
             number[:] = 0.0
 
-            # Expand into our memory layout
-            j = 0
+            # Get new number densities
             for nuc, i_nuc_results in zip(nuclides, nuc_ind):
                 number[i_nuc_results] = self.number[mat, nuc]
-                for react in react_ind:
-                    rates_expanded[i_nuc_results, react] = results[j]
-                    j += 1
+
+            tally_rates = self._rate_helper.get_material_rates(
+                i, nuc_ind, react_ind)
 
             # Accumulate energy from fission
-            energy += self._tally_helper.get_fission_energy(
-                rates_expanded[:, fission_ind], i)
+            energy += self._energy_helper.get_fission_energy(
+                tally_rates[:, fission_ind], i)
 
             # Divide by total number and store
-            for i_nuc_results in nuc_ind:
-                if number[i_nuc_results] != 0.0:
-                    for react in react_ind:
-                        rates_expanded[i_nuc_results, react] /= number[i_nuc_results]
-
-            rates[i, :, :] = rates_expanded
+            rates[i, :, :] = self._rate_helper.divide_by_adens(number)
 
         # Reduce energy produced from all processes
         energy = comm.allreduce(energy)

From 9cbbc6e6b9ab405c473c3c6b756dafab76485222 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Fri, 5 Jul 2019 13:37:12 -0500
Subject: [PATCH 05/14] Document abstract and concrete operator helpers

---
 docs/source/pythonapi/deplete.rst | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index 532e06655..091b25da7 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -75,10 +75,24 @@ data, such as number densities and reaction rates for each material.
    :template: myclass.rst
 
    AtomNumber
+   ChainFissHelper
+   DirectRxnRateHelper
    OperatorResult
    ReactionRates
    Results
    ResultsList
+
+
+The following classes are abstract classes that can be used to extend the
+:mod:`openmc.deplete` capabilities:
+
+.. autosummary::
+   :toctree:generated
+   :nosignatures:
+   :template: myclass.rst
+
+   ReactionRateHelper
+   FissionEnergyHelper
    TransportOperator
 
 Each of the integrator functions also relies on a number of "helper" functions

From 0aebbc13791008e418adcac5f65405cd3be70117 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 15 Jul 2019 11:13:45 -0500
Subject: [PATCH 06/14] Document nuclides attributes for Operator helpers

Remove some pass statements for abstract methods as well
---
 openmc/deplete/abc.py | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 13be50812..565a04089 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -173,6 +173,12 @@ class ReactionRateHelper(ABC):
 
     Reaction rates are passed back to the operator for be used in
     an :class:`openmc.deplete.OperatorResult` instance
+
+    Attributes
+    ----------
+    nuclides : list of str
+        All nuclides with desired reaction rates. Ordered to be
+        consistent with :class:`openmc.deplete.Operator`
     """
 
     def __init__(self):
@@ -183,7 +189,6 @@ class ReactionRateHelper(ABC):
     @abstractmethod
     def generate_tallies(self, materials, scores):
         """Use the capi to build tallies needed for reaction rates"""
-        pass
 
     @property
     def nuclides(self):
@@ -222,6 +227,8 @@ class ReactionRateHelper(ABC):
 
         Parameters
         ----------
+        energy : float
+            Energy produced in this region [W]
         number : iterable of float
             Number density [#/b/cm] of each nuclide tracked in the calculation.
             Ordered identically to :attr:`nuclides`
@@ -242,6 +249,12 @@ class ReactionRateHelper(ABC):
 
 class FissionEnergyHelper(ABC):
     """Abstract class for normalizing fission reactions to a given level
+
+    Attributes
+    ----------
+    nuclides : list of str
+        All nuclides with desired reaction rates. Ordered to be
+        consistent with :class:`openmc.deplete.Operator`
     """
 
     def __init__(self):
@@ -258,12 +271,10 @@ class FissionEnergyHelper(ABC):
         :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):

From f7628035130ff48e306729efb92048693acc27a5 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 15 Jul 2019 11:14:16 -0500
Subject: [PATCH 07/14] DirectRxnRateHelper -> DirectReactionRateHelper

---
 openmc/deplete/helpers.py  | 2 +-
 openmc/deplete/operator.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index 4b79d5b3d..f039eeb70 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -13,7 +13,7 @@ from .abc import ReactionRateHelper, FissionEnergyHelper
 # -------------------------------------
 
 
-class DirectRxnRateHelper(ReactionRateHelper):
+class DirectReactionRateHelper(ReactionRateHelper):
     """Class that generates tallies for one-group rates"""
 
     def generate_tallies(self, materials, scores):
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index fab2a9f48..586a512ca 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -24,7 +24,7 @@ from . import comm
 from .abc import TransportOperator, OperatorResult
 from .atom_number import AtomNumber
 from .reaction_rates import ReactionRates
-from .helpers import DirectRxnRateHelper, ChainFissHelper
+from .helpers import DirectReactionRateHelper, ChainFissHelper
 
 
 def _distribute(items):
@@ -154,7 +154,7 @@ class Operator(TransportOperator):
             self.local_mats, self._burnable_nucs, self.chain.reactions)
 
         # Get class to assist working with tallies
-        self._rate_helper = DirectRxnRateHelper()
+        self._rate_helper = DirectReactionRateHelper()
         self._energy_helper = ChainFissHelper()
 
 

From d23fa7cf657c1bf902c078ad039cbddba33dab69 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 15 Jul 2019 11:39:05 -0500
Subject: [PATCH 08/14] Improve documentation for Operator helpers

---
 openmc/deplete/abc.py     | 45 ++++++++++++++++++++++++++++-----------
 openmc/deplete/helpers.py | 10 ++++-----
 2 files changed, 37 insertions(+), 18 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 565a04089..62db6915e 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -214,10 +214,13 @@ class ReactionRateHelper(ABC):
     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
+        Parameters
+        ----------
+        nuc_index : list of str
+            Ordering of desired nuclides
+        react_index : list of str
+            Ordering of reactions
+        """
 
     def divide_by_adens(self, number):
         """Normalize reaction rates by number of nuclides
@@ -230,12 +233,12 @@ class ReactionRateHelper(ABC):
         energy : float
             Energy produced in this region [W]
         number : iterable of float
-            Number density [#/b/cm] of each nuclide tracked in the calculation.
+            Number density [atoms/b/cm] of each nuclide tracked in the calculation.
             Ordered identically to :attr:`nuclides`
 
         Returns
         -------
-        results : :class:`numpy.ndarray`
+        results : `numpy.ndarray`
             2D array ``[n_nuclides, n_rxns]`` of reaction rates normalized by
             the number of nuclides
         """
@@ -265,16 +268,32 @@ class FissionEnergyHelper(ABC):
     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"""
+        Parameters
+        ----------
+        chain_nucs : list of str
+            All nuclides to be tracked in this problem
+        rate_index : dict of str to int
+            Mapping from nuclide name to index in the
+            reaction rate vector used in :meth:`get_energy`.
+        materials : list of str
+            All materials tracked on the operator helped by this
+            object
+        """
 
     @abstractmethod
     def get_fission_energy(self, fission_rates, mat_index):
-        """Return fission energy in this material given fission rates"""
+        """Return fission energy in this material given fission rates
+
+        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:`prepare`
+        mat_index : int
+            Index for the material requested.
+        """
+
 
     @property
     def nuclides(self):
diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index f039eeb70..5e8b4146e 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -102,16 +102,16 @@ class ChainFissHelper(FissionEnergyHelper):
         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
+        """Return fission energy for this material
 
-        parameters
+        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`
+            material. Should be ordered corresponding to initial
+            ``rate_index`` used in :meth:`prepare`
         _mat_index : int
-            index for the material requested. Unused, as all
+            index for the material requested. Unused, as identical
             isotopes in all materials have the same Q value.
         """
         return dot(fiss_rates, self._fission_E)

From 3b8d9758bf77b5a89d94e3935ea138d60c8ebffa Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 15 Jul 2019 11:53:01 -0500
Subject: [PATCH 09/14] Improve array filling/accessing for Operator.__call__

When possible, replaced ``x[:, :, :] = y`` with  ``x.fill(y)``.
Simple performance tests showed the approaches to be identical
with respect to process time. The latter maintains more
flexibility; whatever the shape of x becomes, no changes
need to be made, so long as x has a ``fill`` method.

Replaced ``rates[i, :, :] = ...`` with ``rates[i]`` per
reviewer comments
---
 openmc/deplete/operator.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 586a512ca..73fad5093 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -513,7 +513,7 @@ class Operator(TransportOperator):
 
         """
         rates = self.reaction_rates
-        rates[:, :, :] = 0.0
+        rates.fill(0.0)
 
         # Get k and uncertainty
         k_combined = openmc.capi.keff()
@@ -527,7 +527,6 @@ class Operator(TransportOperator):
         react_ind = [rates.index_rx[react] for react in self.chain.reactions]
 
         # Compute fission power
-        # TODO : improve this calculation
 
         # Keep track of energy produced from all reactions in eV per source
         # particle
@@ -545,7 +544,7 @@ class Operator(TransportOperator):
             slab = materials.index(mat)
 
             # Zero out reaction rates and nuclide numbers
-            number[:] = 0.0
+            number.fill(0.0)
 
             # Get new number densities
             for nuc, i_nuc_results in zip(nuclides, nuc_ind):
@@ -559,7 +558,7 @@ class Operator(TransportOperator):
                 tally_rates[:, fission_ind], i)
 
             # Divide by total number and store
-            rates[i, :, :] = self._rate_helper.divide_by_adens(number)
+            rates[i] = self._rate_helper.divide_by_adens(number)
 
         # Reduce energy produced from all processes
         energy = comm.allreduce(energy)

From 15c2bc2de3d51715cb073016c5d8c58daebfe82a Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 18 Jul 2019 12:39:05 -0500
Subject: [PATCH 10/14] Use correct material index in operator unpacking

Fix a bug introduced where the material iteration index,

        for i, mat in enumerate(self.local_mats):

was used for the global __material__ index. The value of i
was used to extract reaction rate tallies for that material.
This causes an issue with multiple operators on MPI processes,
where an Operator`s material i does not equal the i-th burnable
material.
---
 openmc/deplete/operator.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 73fad5093..a3c7c2ba8 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -551,7 +551,7 @@ class Operator(TransportOperator):
                 number[i_nuc_results] = self.number[mat, nuc]
 
             tally_rates = self._rate_helper.get_material_rates(
-                i, nuc_ind, react_ind)
+                slab, nuc_ind, react_ind)
 
             # Accumulate energy from fission
             energy += self._energy_helper.get_fission_energy(

From d0d4a05c7599bc7a150e7fb76118b34309d2df96 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 18 Jul 2019 13:13:13 -0500
Subject: [PATCH 11/14] Add a dictionary to map local to global burnable
 materials

Removes the need to call self.burnable_mats.index
for every local material at every depletion step
---
 openmc/deplete/operator.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index a3c7c2ba8..03c9813f5 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -138,6 +138,7 @@ class Operator(TransportOperator):
         openmc.reset_auto_ids()
         self.burnable_mats, volume, nuclides = self._get_burnable_mats()
         self.local_mats = _distribute(self.burnable_mats)
+        self._mat_index_map = {}
 
         # Determine which nuclides have incident neutron data
         self.nuclides_with_data = self._get_nuclides_with_data()
@@ -386,6 +387,10 @@ class Operator(TransportOperator):
         self._energy_helper.prepare(
             self.chain.nuclides, self.reaction_rates.index_nuc, materials)
 
+        # Generate map from local materials => material index
+        self._mat_index_map = {
+            lm: self.burnable_mats.index(lm) for lm in self.local_mats}
+
         # Return number density vector
         return list(self.number.get_mat_slice(np.s_[:]))
 
@@ -519,7 +524,6 @@ class Operator(TransportOperator):
         k_combined = openmc.capi.keff()
 
         # Extract tally bins
-        materials = self.burnable_mats
         nuclides = self._rate_helper.nuclides
 
         # Form fast map
@@ -541,7 +545,7 @@ class Operator(TransportOperator):
         # Extract results
         for i, mat in enumerate(self.local_mats):
             # Get tally index
-            slab = materials.index(mat)
+            mat_index = self._mat_index_map[mat]
 
             # Zero out reaction rates and nuclide numbers
             number.fill(0.0)
@@ -551,7 +555,7 @@ class Operator(TransportOperator):
                 number[i_nuc_results] = self.number[mat, nuc]
 
             tally_rates = self._rate_helper.get_material_rates(
-                slab, nuc_ind, react_ind)
+                mat_index, nuc_ind, react_ind)
 
             # Accumulate energy from fission
             energy += self._energy_helper.get_fission_energy(
@@ -561,7 +565,10 @@ class Operator(TransportOperator):
             rates[i] = self._rate_helper.divide_by_adens(number)
 
         # Reduce energy produced from all processes
+        print("Energy - {}: {:9.7e}".format(comm.rank, energy))
         energy = comm.allreduce(energy)
+        if comm.rank == 0:
+            print("Energy: {:9.7e}".format(energy))
 
         # Determine power in eV/s
         power /= JOULE_PER_EV

From 808a41b4e395c0c86ff0cfb957b82fc1d5dfa9d9 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 18 Jul 2019 14:35:39 -0500
Subject: [PATCH 12/14] Improve energy calculation with EnergyHelper.update

EnergyHelper (was FissionEnergyHelper) has a new
``update`` method that replaces get_energy. This method
is no longer abstract and defaults to not performing any actions.
This works in conjunction with a new private ``_energy`` atribute
and ``reset`` method to streamline the energy normalization procedure.
For the current concrete ChainFissionHelper class, the reset method
sets _energy to zero, while the update method updates the energy
with the dot product fission_rates X fission_q_vector.
The total energy produced is made accessible, but not publically
writable with an energy property. After cycling through all
local materials on an operator, the total system energy is computed
by

        energy = comm.allreduce(self._energy_helper.energy)

For a concrete class where the system energy is computed via tallies
and not with Q-values, the update methods could continue to do
nothing. The total system energy could be stored either in the
reset method or returned directly from the energy attribute/property.
---
 openmc/deplete/abc.py      | 43 ++++++++++++++++++++++++++++----------
 openmc/deplete/helpers.py  | 24 ++++++++++++---------
 openmc/deplete/operator.py | 16 ++++++--------
 3 files changed, 52 insertions(+), 31 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 62db6915e..b171d9356 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -250,23 +250,44 @@ class ReactionRateHelper(ABC):
         return results
 
 
-class FissionEnergyHelper(ABC):
-    """Abstract class for normalizing fission reactions to a given level
+class EnergyHelper(ABC):
+    """Abstract class for obtaining energy produced
+
+    The ultimate goal of this helper is to provide instances of
+    :class:`openmc.deplete.Operator` with the total energy produced
+    in a transport simulation. This information, provided with the
+    power requested by the user and reaction rates from a
+    :class:`ReactionRateHelper` will scale reaction rates to the
+    correct values.
 
     Attributes
     ----------
     nuclides : list of str
         All nuclides with desired reaction rates. Ordered to be
         consistent with :class:`openmc.deplete.Operator`
+    energy : float
+        Total energy [eV/s] produced in a transport simulation.
+        Updated in the material iteration with :meth:`update`.
     """
 
     def __init__(self):
         self._nuclides = None
-        self._fission_E = None
+        self._energy = 0.0
+
+    @property
+    def energy(self):
+        return self._energy
+
+    def reset(self):
+        """Reset energy produced prior to unpacking tallies"""
+        self._energy = 0.0
 
     @abstractmethod
     def prepare(self, chain_nucs, rate_index, materials):
-        """Perform work needed to obtain fission energy per material
+        """Perform work needed to obtain energy produced
+
+        This method is called prior to the transport simulations
+        in :meth:`openmc.deplete.Operator.initial_condition`.
 
         Parameters
         ----------
@@ -274,15 +295,15 @@ class FissionEnergyHelper(ABC):
             All nuclides to be tracked in this problem
         rate_index : dict of str to int
             Mapping from nuclide name to index in the
-            reaction rate vector used in :meth:`get_energy`.
+            `fission_rates` for :meth:`update`.
         materials : list of str
             All materials tracked on the operator helped by this
-            object
+            object. Should correspond to
+            :attr:`openmc.deplete.Operator.burnable_materials`
         """
 
-    @abstractmethod
-    def get_fission_energy(self, fission_rates, mat_index):
-        """Return fission energy in this material given fission rates
+    def update(self, fission_rates, mat_index):
+        """Update the energy produced
 
         Parameters
         ----------
@@ -291,10 +312,10 @@ class FissionEnergyHelper(ABC):
             material. Should be ordered corresponding to initial
             ``rate_index`` used in :meth:`prepare`
         mat_index : int
-            Index for the material requested.
+            Index for the specific material in the list of all burnable
+            materials.
         """
 
-
     @property
     def nuclides(self):
         """List of nuclides with requested reaction rates"""
diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index 5e8b4146e..da85312a8 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -6,7 +6,7 @@ from itertools import product
 from numpy import dot, zeros
 
 from openmc.capi import Tally, MaterialFilter
-from .abc import ReactionRateHelper, FissionEnergyHelper
+from .abc import ReactionRateHelper, EnergyHelper
 
 # -------------------------------------
 # Helpers for generating reaction rates
@@ -68,9 +68,13 @@ class DirectReactionRateHelper(ReactionRateHelper):
 # ------------------------------------
 
 
-class ChainFissHelper(FissionEnergyHelper):
+class ChainFissionHelper(EnergyHelper):
     """Fission Q-values are pulled from chain"""
 
+    def __init__(self):
+        super().__init__()
+        self._fission_q_vector = None
+
     def prepare(self, chain_nucs, rate_index, _materials):
         """Populate the fission Q value vector from a chain.
 
@@ -86,23 +90,23 @@ class ChainFissHelper(FissionEnergyHelper):
         _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),)):
+        if (self._fission_q_vector is not None
+                and self._fission_q_vector.shape == (len(rate_index),)):
             return
 
-        fiss_E = zeros(len(rate_index))
+        fission_qs = 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
+                        fission_qs[rate_index[nuclide.name]] = rx.Q
                         break
 
-        self._fission_E = fiss_E
+        self._fission_q_vector = fission_qs
 
-    def get_fission_energy(self, fiss_rates, _mat_index):
-        """Return fission energy for this material
+    def update(self, fission_rates, _mat_index):
+        """Update energy produced with fission rates in a material
 
         Parameters
         ----------
@@ -114,4 +118,4 @@ class ChainFissHelper(FissionEnergyHelper):
             index for the material requested. Unused, as identical
             isotopes in all materials have the same Q value.
         """
-        return dot(fiss_rates, self._fission_E)
+        self._energy += dot(fission_rates, self._fission_q_vector)
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 03c9813f5..815602247 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -24,7 +24,7 @@ from . import comm
 from .abc import TransportOperator, OperatorResult
 from .atom_number import AtomNumber
 from .reaction_rates import ReactionRates
-from .helpers import DirectReactionRateHelper, ChainFissHelper
+from .helpers import DirectReactionRateHelper, ChainFissionHelper
 
 
 def _distribute(items):
@@ -154,9 +154,9 @@ class Operator(TransportOperator):
         self.reaction_rates = ReactionRates(
             self.local_mats, self._burnable_nucs, self.chain.reactions)
 
-        # Get class to assist working with tallies
+        # Get classes to assist working with tallies
         self._rate_helper = DirectReactionRateHelper()
-        self._energy_helper = ChainFissHelper()
+        self._energy_helper = ChainFissionHelper()
 
 
     def __call__(self, vec, power, print_out=True):
@@ -534,7 +534,7 @@ class Operator(TransportOperator):
 
         # Keep track of energy produced from all reactions in eV per source
         # particle
-        energy = 0.0
+        self._energy_helper.reset()
 
         # Create arrays to store fission Q values, reaction rates, and nuclide
         # numbers, zeroed out in material iteration
@@ -558,17 +558,13 @@ class Operator(TransportOperator):
                 mat_index, nuc_ind, react_ind)
 
             # Accumulate energy from fission
-            energy += self._energy_helper.get_fission_energy(
-                tally_rates[:, fission_ind], i)
+            self._energy_helper.update(tally_rates[:, fission_ind], mat_index)
 
             # Divide by total number and store
             rates[i] = self._rate_helper.divide_by_adens(number)
 
         # Reduce energy produced from all processes
-        print("Energy - {}: {:9.7e}".format(comm.rank, energy))
-        energy = comm.allreduce(energy)
-        if comm.rank == 0:
-            print("Energy: {:9.7e}".format(energy))
+        energy = comm.allreduce(self._energy_helper.energy)
 
         # Determine power in eV/s
         power /= JOULE_PER_EV

From 9f41dce1d253b876bf5417dc33771422a5bf549f Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Tue, 23 Jul 2019 08:29:37 -0500
Subject: [PATCH 13/14] Return energy in J/s/source neutron for EnergyHelper

Addressing comments in review for #1278

- Documentation cleanup
- Better naming convention regarding ReactionRateHelper
  results cache
- The power in Operator tally unpacking and normalizing is
  no longer converted to eV/s, since the EnergyHelper.energy
  property is now returned in J/s/source neutron
---
 docs/source/pythonapi/deplete.rst |  6 ++---
 openmc/deplete/abc.py             | 34 +++++++++++++-------------
 openmc/deplete/helpers.py         | 40 ++++++++++++++++++++++---------
 openmc/deplete/operator.py        | 16 +++++--------
 4 files changed, 55 insertions(+), 41 deletions(-)

diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index 091b25da7..70d48f931 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -75,8 +75,8 @@ data, such as number densities and reaction rates for each material.
    :template: myclass.rst
 
    AtomNumber
-   ChainFissHelper
-   DirectRxnRateHelper
+   ChainFissionHelper
+   DirectReactionRateHelper
    OperatorResult
    ReactionRates
    Results
@@ -92,7 +92,7 @@ The following classes are abstract classes that can be used to extend the
    :template: myclass.rst
 
    ReactionRateHelper
-   FissionEnergyHelper
+   EnergyHelper
    TransportOperator
 
 Each of the integrator functions also relies on a number of "helper" functions
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index b171d9356..ecff16868 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -11,9 +11,9 @@ from abc import ABC, abstractmethod
 from xml.etree import ElementTree as ET
 from warnings import warn
 
-from numpy import zeros, nonzero
+from numpy import nonzero, empty
 
-from openmc.data import DataLibrary
+from openmc.data import DataLibrary, JOULE_PER_EV
 from openmc.checkvalue import check_type
 from .chain import Chain
 
@@ -188,7 +188,7 @@ class ReactionRateHelper(ABC):
 
     @abstractmethod
     def generate_tallies(self, materials, scores):
-        """Use the capi to build tallies needed for reaction rates"""
+        """Use the C API to build tallies needed for reaction rates"""
 
     @property
     def nuclides(self):
@@ -201,13 +201,15 @@ class ReactionRateHelper(ABC):
         self._nuclides = nuclides
         self._rate_tally.nuclides = nuclides
 
-    def _reset_results_cache(self, nnucs, nreact):
-        """Cache for results for a given material"""
+    def _get_results_cache(self, nnucs, nreact):
+        """Cache for results for a given material
+
+        Creates an empty array of shape ``(nnucs, nreact)``
+        if the shape does not match the current cache.
+        """
         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)
+            self._results_cache = empty((nnucs, nreact))
         return self._results_cache
 
     @abstractmethod
@@ -216,6 +218,8 @@ class ReactionRateHelper(ABC):
 
         Parameters
         ----------
+        mat_id : int
+            Unique ID for the requested material
         nuc_index : list of str
             Ordering of desired nuclides
         react_index : list of str
@@ -230,17 +234,13 @@ class ReactionRateHelper(ABC):
 
         Parameters
         ----------
-        energy : float
-            Energy produced in this region [W]
-        number : iterable of float
-            Number density [atoms/b/cm] of each nuclide tracked in the calculation.
             Ordered identically to :attr:`nuclides`
 
         Returns
         -------
-        results : `numpy.ndarray`
-            2D array ``[n_nuclides, n_rxns]`` of reaction rates normalized by
-            the number of nuclides
+        results : numpy.ndarray
+            Array of reactions rates of shape ``(n_nuclides, n_rxns)``
+            normalized by the number of nuclides
         """
 
         mask = nonzero(number)
@@ -266,7 +266,7 @@ class EnergyHelper(ABC):
         All nuclides with desired reaction rates. Ordered to be
         consistent with :class:`openmc.deplete.Operator`
     energy : float
-        Total energy [eV/s] produced in a transport simulation.
+        Total energy [J/s/source neutron] produced in a transport simulation.
         Updated in the material iteration with :meth:`update`.
     """
 
@@ -276,7 +276,7 @@ class EnergyHelper(ABC):
 
     @property
     def energy(self):
-        return self._energy
+        return self._energy * JOULE_PER_EV
 
     def reset(self):
         """Reset energy produced prior to unpacking tallies"""
diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index da85312a8..edf0feb06 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -14,7 +14,14 @@ from .abc import ReactionRateHelper, EnergyHelper
 
 
 class DirectReactionRateHelper(ReactionRateHelper):
-    """Class that generates tallies for one-group rates"""
+    """Class that generates tallies for one-group rates
+
+    Attributes
+    ----------
+    nuclides : list of str
+        All nuclides with desired reaction rates. Ordered to be
+        consistent with :class:`openmc.deplete.Operator`
+    """
 
     def generate_tallies(self, materials, scores):
         """Produce one-group reaction rate tally
@@ -41,7 +48,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
         Parameters
         ----------
         mat_id : int
-            Unique id for the requested material
+            Unique ID for the requested material
         nuc_index : iterable of int
             Index for each nuclide in :attr:`nuclides` in the
             desired reaction rate matrix
@@ -50,11 +57,12 @@ class DirectReactionRateHelper(ReactionRateHelper):
 
         Returns
         -------
-        rates : :class:`numpy.ndarray`
-            2D matrix ``(len(nuc_index), len(react_index))`` with the
+        rates : numpy.ndarray
+            Array with shape ``(n_nuclides, n_rxns)`` with the
             reaction rates in this material
         """
-        results = self._reset_results_cache(len(nuc_index), len(react_index))
+        results = self._get_results_cache(len(nuc_index), len(react_index))
+        results.fill(0.0)
         full_tally_res = self._rate_tally.results[mat_id, :, 1]
         for i_tally, (i_nuc, i_react) in enumerate(
                 product(nuc_index, react_index)):
@@ -63,13 +71,23 @@ class DirectReactionRateHelper(ReactionRateHelper):
         return results
 
 
-# ------------------------------------
-# Helpers for obtaining fission energy
-# ------------------------------------
+# ----------------------------
+# Helpers for obtaining energy
+# ----------------------------
 
 
 class ChainFissionHelper(EnergyHelper):
-    """Fission Q-values are pulled from chain"""
+    """Computes energy using fission Q values from depletion chain
+
+    Attributes
+    ----------
+    nuclides : list of str
+        All nuclides with desired reaction rates. Ordered to be
+        consistent with :class:`openmc.deplete.Operator`
+    energy : float
+        Total energy [J/s/source neutron] produced in a transport simulation.
+        Updated in the material iteration with :meth:`update`.
+    """
 
     def __init__(self):
         super().__init__()
@@ -78,8 +96,8 @@ class ChainFissionHelper(EnergyHelper):
     def prepare(self, chain_nucs, rate_index, _materials):
         """Populate the fission Q value vector from a chain.
 
-        Paramters
-        ---------
+        Parameters
+        ----------
         chain_nucs : iterable of :class:`openmc.deplete.Nuclide`
             Nuclides used in this depletion chain. Do not need
             to be ordered
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 815602247..bb3c96e4f 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -19,7 +19,6 @@ import numpy as np
 
 import openmc
 import openmc.capi
-from openmc.data import JOULE_PER_EV
 from . import comm
 from .abc import TransportOperator, OperatorResult
 from .atom_number import AtomNumber
@@ -138,7 +137,11 @@ class Operator(TransportOperator):
         openmc.reset_auto_ids()
         self.burnable_mats, volume, nuclides = self._get_burnable_mats()
         self.local_mats = _distribute(self.burnable_mats)
-        self._mat_index_map = {}
+
+        # Generate map from local materials => material index
+        self._mat_index_map = {
+            lm: self.burnable_mats.index(lm) for lm in self.local_mats}
+
 
         # Determine which nuclides have incident neutron data
         self.nuclides_with_data = self._get_nuclides_with_data()
@@ -158,7 +161,6 @@ class Operator(TransportOperator):
         self._rate_helper = DirectReactionRateHelper()
         self._energy_helper = ChainFissionHelper()
 
-
     def __call__(self, vec, power, print_out=True):
         """Runs a simulation.
 
@@ -387,10 +389,6 @@ class Operator(TransportOperator):
         self._energy_helper.prepare(
             self.chain.nuclides, self.reaction_rates.index_nuc, materials)
 
-        # Generate map from local materials => material index
-        self._mat_index_map = {
-            lm: self.burnable_mats.index(lm) for lm in self.local_mats}
-
         # Return number density vector
         return list(self.number.get_mat_slice(np.s_[:]))
 
@@ -564,11 +562,9 @@ class Operator(TransportOperator):
             rates[i] = self._rate_helper.divide_by_adens(number)
 
         # Reduce energy produced from all processes
+        # J / s / source neutron
         energy = comm.allreduce(self._energy_helper.energy)
 
-        # Determine power in eV/s
-        power /= JOULE_PER_EV
-
         # Scale reaction rates to obtain units of reactions/sec
         rates *= power / energy
 

From 6ae45931a67df84592e016da8ea85a2b1c111b80 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Tue, 23 Jul 2019 11:28:58 -0500
Subject: [PATCH 14/14] Re-document number parameter in
 ReactionRateHelper.divide_by_adens

---
 openmc/deplete/abc.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index ecff16868..c5c086854 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -234,6 +234,8 @@ class ReactionRateHelper(ABC):
 
         Parameters
         ----------
+        number : iterable of float
+            Number density [atoms/b-cm] of each nuclide tracked in the calculation.
             Ordered identically to :attr:`nuclides`
 
         Returns