Merge pull request #1350 from drewejohnson/energy-dep-local-heating

Use MT301,901 for energy deposition for depletion

docs/source/methods/energy_deposition.rst

@@ -0,0 +1,149 @@
+.. _methods_heating:
+
+=============================
+Heating and Energy Deposition
+=============================
+
+As particles traverse a problem, some portion of their energy is deposited at 
+collision sites. This energy is deposited when charged particles, including
+electrons and recoil nuclei, undergo electromagnetic interactions with
+surrounding electons and ions. The information describing how much energy
+is deposited for a specific reaction is referred to as 
+"heating numbers" and can be computed using a program like NJOY with the
+``heatr`` module.
+
+These heating rate is the product of reaction-specific coefficients and
+a reaction cross section
+
+.. math::
+
+    H(E) = \phi(E)\sum_i\rho_i\sum_rk_{i, r}(E),
+
+and has units energy per time, typically eV / s.
+Here, :math:`k_{i, r}` are the KERMA (Kinetic Energy Release in Materials)
+[Mack97]_ coefficients for reaction :math:`r` of isotope :math:`i`.
+The KERMA coefficients have units energy :math:`\times` cross-section, e.g.
+eV-barn, and can be used much like a reaction cross section for the purpose
+of tallying energy deposition.
+
+KERMA coefficients can be computed using the energy-balance method with
+a nuclear data processing code like NJOY, which performs the following
+iteration over all reactions :math:`r` for all isotopes :math:`i`
+requested
+
+.. math::
+
+    k_{i, r}(E) = \left(E + Q_{i, r} - \bar{E}_{i, r, n}
+    - \bar{E}_{i, r, \gamma}\right)\sigma_{i, r}(E),
+
+removing the energy of neutral particles (neutrons and photons) that are
+transported away from the reaction site :math:`\bar{E}`, and the reaction
+:math:`Q` value.
+
+-------
+Fission
+-------
+
+During a fission event, there are potentially many secondary particles, and all
+must be considered. The total energy released in a fission event is typically
+broken up into the following categories:
+
+- :math:`E_{fr}` - kinetic energy of fission fragments
+- :math:`E_{n,p}` - energy of prompt fission neutrons
+- :math:`E_{n,d}` - energy of delayed fission neutrons
+- :math:`E_{\gamma,p}` - energy of prompt fission photons
+- :math:`E_{\gamma,d}` - energy of delayed fission photons
+- :math:`E_{\beta}` - energy of released :math:`\beta` particles
+- :math:`E_{\nu}` - energy of neutrinos
+
+These components are defined in MF=1,MT=458 data in a standard ENDF/B-6 formatted
+file. All these quantities may depend upon incident neutron energy,
+but this dependence is not shown to make the following demonstrations cleaner.
+As neutrinos scarcely interact with matter, the recoverable energy from
+fission is defined as
+
+.. math::
+
+    E_r\equiv E_{fr} + E_{n,p} + E_{n, d} + E_{\gamma, p}
+    + E_{\gamma, d} + E_{\beta}
+
+Furthermore, the energy of the secondary neutrons and photons is given as
+:math:`E_{n, p}` and :math:`E_{\gamma, p}`, respectively.
+
+NJOY computes the fission KERMA coefficient using this energy-balance method to be
+
+.. math::
+
+    k_{i, f}(E) = \left[E + Q(E) - \bar{E}(E)\right]\sigma_{i, f}(E)
+    = \left[E_{fr} + E_{\gamma, p}\right]\sigma_{i, j}(E)
+
+.. note::
+
+    The energy from delayed neutrons and photons and beta particles is intentionally
+    left out from the NJOY calculations.
+
+---------------------
+OpenMC Implementation
+---------------------
+
+For fissile isotopes, OpenMC makes modifications to the heating reaction to
+include all relevant components of fission energy release. These modifications
+are made to the total heating reaction, MT=301. Breaking the total heating
+KERMA into a fission and non-fission section, one can write
+
+.. math::
+
+    k_i(E) = k_{i, nf}(E) + \left[E_{fr}(E) + E_{\gamma, p}\right]\sigma_{i, f}(E)
+
+OpenMC seeks to modify the total heating data to include energy from
+:math:`\beta` particles and, conditionally, delayed photons. This conditional
+inclusion depends on the simulation mode: neutron transport, or coupled
+neutron-photon transport. The heating due to fission is removed using MT=318
+data, and then re-built using the desired components of fission energy release
+from MF=1,MT=458 data.
+
+Neutron Transport
+-----------------
+
+For this case, OpenMC instructs ``heatr`` to produce heating coefficients
+assuming that energy from photons, :math:`E_{\gamma, p}` and
+:math:`E_{\gamma, d}`, is deposited at the fission site.
+Let :math:`N901` represent the total heating number returned from this ``heatr`` 
+run with :math:`N918` reflecting fission heating computed from NJOY.
+:math:`M901` represent the following modification
+
+.. math::
+
+    M901_{i}(E)\equiv N901_{i}(E) - N918_{i}(E)
+      + \left[E_{i, fr} + E_{i, \beta} + E_{i, \gamma, p}
+      + E_{i, \gamma, d}\right]\sigma_{i, f}(E).
+
+This modified heating data is stored as the MT=901 reaction and will be scored
+if ``901`` is included in :attr:`openmc.Tally.scores`.
+
+Coupled neutron-photon transport
+--------------------------------
+
+Here, OpenMC instructs ``heatr`` to assume that energy from photons is not
+deposited locally. However, the definitions provided in the NJOY manual
+indicate that, regardless of this mode, the prompt photon energy is still
+included in :math:`k_{i, f}`, and therefore must be manually removed.
+Let :math:`N301` represent the total heating number returned from this
+``heatr`` run and :math:`M301` be
+
+.. math::
+
+    M301_{i}(E)\equiv N301_{i}(E) - N318_{i}(E)
+      + \left[E_{i, fr}(E) + E_{i, \beta}(E)\right]\sigma_{i, f}(E).
+
+This modified heating data is stored as the MT=301 reaction and will be scored
+if ``heating`` is included in :attr:`openmc.Tally.scores`.
+
+----------
+References
+----------
+
+.. [Mack97] Abdou, M.A., Maynard, C.W., and Wright, R.Q. MACK: computer
+   program to calculate neutron energy release parameters (fluence-to-kerma
+   factors) and multigroup neutron reaction cross sections from nuclear data
+   in ENDF Format. Oak Ridge National Laboratory report ORNL-TM-3994.

docs/source/methods/index.rst

@@ -18,3 +18,4 @@ Theory and Methodology
     eigenvalue
     parallelization
     cmfd
+    energy_deposition

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

docs/source/usersguide/tallies.rst

@@ -261,12 +261,13 @@ The following tables show all valid scores:
     |                      |produced by NJOY's HEATR module while for photons, |
     |                      |this is tallied from either direct photon energy   |
     |                      |deposition (analog estimator) or pre-generated     |
-    |                      |photon heating number.                             |
+    |                      |photon heating number. See :ref:`methods_heating`  |
     +----------------------+---------------------------------------------------+
     |heating-local         |Total nuclear heating in units of eV per source    |
     |                      |particle assuming energy from secondary photons is |
     |                      |deposited locally. Note that this score should only|
-    |                      |be used for incident neutrons.                     |
+    |                      |be used for incident neutrons. See                 |
+    |                      |:ref:`methods_heating`.                            |
     +----------------------+---------------------------------------------------+
     |kappa-fission         |The recoverable energy production rate due to      |
     |                      |fission. The recoverable energy is defined as the  |

openmc/deplete/dummy_comm.py

@@ -1,3 +1,5 @@
+import sys
+
 class DummyCommunicator(object):
     rank = 0
     size = 1
@@ -25,3 +27,6 @@ class DummyCommunicator(object):
 
     def scatter(self, sendobj, root=0):
         return sendobj[0]
+
+    def Abort(self, exit_code_or_msg):
+        sys.exit(exit_code_or_msg)

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,56 @@ 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
+
+    Parameters
+    ----------
+    score : string
+        Valid score to use when obtaining system energy from OpenMC.
+        Defaults to "heating-local"
+
+    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.
+    score : str
+        Score used to obtain system energy
+
+    """
+
+    def __init__(self, score="heating-local"):
+        super().__init__()
+        self.score = score
+        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 = [self.score]
+
+    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 comm.rank == 0:
+            self._energy = self._tally.results[0, 0, 1]
+
 # ------------------------------------
 # Helper for collapsing fission yields
 # ------------------------------------

openmc/deplete/operator.py

@@ -7,12 +7,14 @@ densities is all done in-memory instead of through the filesystem.
 
 """
 
+import sys
 import copy
 from collections import OrderedDict
 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 +29,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 +80,15 @@ 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
     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,11 @@ 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()
+        if energy_mode == "fission-q":
+            self._energy_helper = ChainFissionHelper()
+        else:
+            score = "heating" if settings.photon_transport else "heating-local"
+            self._energy_helper = EnergyScoreHelper(score)
 
         # Select and create fission yield helper
         fission_helper = self._fission_helpers[fission_yield_mode]
@@ -216,6 +237,10 @@ class Operator(TransportOperator):
     def __call__(self, vec, power):
         """Runs a simulation.
 
+        Simulation will abort under the following circumstances:
+
+            1) No energy is computed using OpenMC tallies.
+
         Parameters
         ----------
         vec : list of numpy.ndarray
@@ -235,8 +260,6 @@ class Operator(TransportOperator):
         # Update status
         self.number.set_density(vec)
 
-        time_start = time.time()
-
         # Update material compositions and tally nuclides
         self._update_materials()
         nuclides = self._get_tally_nuclides()
@@ -248,8 +271,6 @@ class Operator(TransportOperator):
         openmc.lib.reset()
         openmc.lib.run()
 
-        time_openmc = time.time()
-
         # Extract results
         op_result = self._unpack_tallies_and_normalize(power)
 
@@ -634,6 +655,15 @@ class Operator(TransportOperator):
         # J / s / source neutron
         energy = comm.allreduce(self._energy_helper.energy)
 
+        # Guard against divide by zero
+        if energy == 0:
+            if comm.rank == 0:
+                sys.stderr.flush()
+                print(" No energy reported from OpenMC tallies. Do your HDF5 "
+                      "files have heating data?\n", file=sys.stderr, flush=True)
+            comm.barrier()
+            comm.Abort(1)
+
         # Scale reaction rates to obtain units of reactions/sec
         rates *= power / energy
 

openmc/deplete/results_list.py

@@ -137,7 +137,7 @@ class ResultsList(list):
     def get_depletion_time(self):
         """Return an array of the average time to deplete a material
 
-        ..note::
+        .. note::
 
             Will have one fewer row than number of other methods,
             like :meth:`get_eigenvalues`, because no depletion
@@ -145,7 +145,6 @@ class ResultsList(list):
 
         Returns
         -------
-
         times : :class:`numpy.ndarray`
             Vector of average time to deplete a single material
             across all processes and materials.