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

@@ -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"""

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)

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