Use energy deposition score for depletion

Introduce a new subclass of EnergyHelper, EnergyScoreHelper,
that computes the system energy using the energy-deposition score.
This energy is fed back to the Operator to normalize reaction rates.

The energy from the tally is only stored on the helper on the
MPI process 0 as to avoid scaling the system energy by the number
of processes. During the Operator unpacking, the energy reported
by each process is reduced, as the previous implementations took
fission reaction rates from the "local materials", e.g. the
materials each process is responsible for depleting.
The tally results are shared across all processes and only
contains a single quantity, the tallied energy deposition across
all materials.

This mode is controlled by passing the "energy_mode" argument passed
to the Operator. The two options are "fission-q" [default and previous
behavior] and "energy-deposition" [new features]. If energy_mode indicates using
the energy deposition score, the user-supplied fission-q dictionary is not used.

(cherry picked from commit d566f3080f)

docs/source/pythonapi/deplete.rst

@@ -95,6 +95,7 @@ total system energy.
    helpers.ChainFissionHelper
    helpers.ConstantFissionYieldHelper
    helpers.DirectReactionRateHelper
+   helpers.EnergyScoreHelper
    helpers.FissionYieldCutoffHelper
 
 The following classes are abstract classes that can be used to extend the

openmc/deplete/helpers.py

@@ -9,6 +9,7 @@ from collections import defaultdict
 
 from numpy import dot, zeros, newaxis
 
+from . import comm
 from openmc.checkvalue import check_type, check_greater_than
 from openmc.lib import (
     Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter)
@@ -157,6 +158,47 @@ class ChainFissionHelper(EnergyHelper):
         self._energy += dot(fission_rates, self._fission_q_vector)
 
 
+class EnergyScoreHelper(EnergyHelper):
+    """Class responsible for obtaining system energy via a tally score
+
+    Attributes
+    ----------
+    nuclides : list of str
+        List of nuclides with reaction rates. Not needed, but provided
+        for a consistent API across other :class:`EnergyHelper`
+    energy : float
+        System energy [eV] computed from the tally. Will be zero for
+        all MPI processes that are not the "master" process to avoid
+        artificially increasing the tallied energy.
+
+    """
+
+    def __init__(self):
+        super().__init__()
+        self._tally = None
+
+    def prepare(self, *args, **kwargs):
+        """Create a tally for system energy production
+
+        Input arguments are not used, as the only information needed
+        is :attr:`score`
+
+        """
+        self._tally = Tally()
+        self._tally.scores = ["energy-deposition"]
+
+    def reset(self):
+        """Obtain system energy from tally
+
+        Only the master process, ``comm.rank == 0`` will
+        have a non-zero :attr:`energy` taken from the tally.
+        This avoids accidentally scaling the system power by
+        the number of MPI processes
+        """
+        super().reset()
+        if not comm.rank:
+            self._energy = self._tally.results[0, 0, 1]
+
 # ------------------------------------
 # Helper for collapsing fission yields
 # ------------------------------------

openmc/deplete/operator.py

@@ -13,6 +13,7 @@ from itertools import chain
 import os
 import time
 import xml.etree.ElementTree as ET
+from warnings import warn
 
 import h5py
 import numpy as np
@@ -27,7 +28,7 @@ from .reaction_rates import ReactionRates
 from .results_list import ResultsList
 from .helpers import (
     DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
-    FissionYieldCutoffHelper, AveragedFissionYieldHelper)
+    FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper)
 
 
 def _distribute(items):
@@ -78,9 +79,16 @@ class Operator(TransportOperator):
     diff_burnable_mats : bool, optional
         Whether to differentiate burnable materials with multiple instances.
         Default: False.
+    energy_mode : {"energy-deposition", "fission-q"}
+        Indicator for computing system energy. ``"energy-deposition"`` will
+        compute with a single energy deposition tally, taking fission energy
+        release data and heating into consideration. ``"fission-q"`` will
+        use the fission Q values from the depletion chain, not taking
+        heating into consideration.
     fission_q : dict, optional
         Dictionary of nuclides and their fission Q values [eV]. If not given,
-        values will be pulled from the ``chain_file``.
+        values will be pulled from the ``chain_file``. Only applicable
+        if ``"energy_mode" == "fission-q"``
     dilute_initial : float, optional
         Initial atom density [atoms/cm^3] to add for nuclides that are zero
         in initial condition to ensure they exist in the decay chain.
@@ -144,13 +152,22 @@ class Operator(TransportOperator):
     }
 
     def __init__(self, geometry, settings, chain_file=None, prev_results=None,
-                 diff_burnable_mats=False, fission_q=None,
-                 dilute_initial=1.0e3, fission_yield_mode="constant",
-                 fission_yield_opts=None):
+                 diff_burnable_mats=False, energy_mode="fission-q",
+                 fission_q=None, dilute_initial=1.0e3,
+                 fission_yield_mode="constant", fission_yield_opts=None):
         if fission_yield_mode not in self._fission_helpers:
             raise KeyError(
                 "fission_yield_mode must be one of {}, not {}".format(
                     ", ".join(self._fission_helpers), fission_yield_mode))
+        if energy_mode == "energy-deposition":
+            if fission_q is not None:
+                warn("Fission Q dictionary not used if energy deposition "
+                     "is used")
+                fission_q = None
+        elif energy_mode != "fission-q":
+            raise ValueError(
+                "energy_mode {} not supported. Must be energy-deposition "
+                "or fission-q".format(energy_mode))
         super().__init__(chain_file, fission_q, dilute_initial, prev_results)
         self.round_number = False
         self.prev_res = None
@@ -204,7 +221,9 @@ class Operator(TransportOperator):
         # Get classes to assist working with tallies
         self._rate_helper = DirectReactionRateHelper(
             self.reaction_rates.n_nuc, self.reaction_rates.n_react)
-        self._energy_helper = ChainFissionHelper()
+        self._energy_helper = (
+            ChainFissionHelper() if energy_mode == "fission-q"
+            else EnergyScoreHelper())
 
         # Select and create fission yield helper
         fission_helper = self._fission_helpers[fission_yield_mode]