diff --git a/docs/source/methods/energy_deposition.rst b/docs/source/methods/energy_deposition.rst new file mode 100644 index 0000000000..af4962c553 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index d8e6ee8190..5eff6c50e6 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -18,3 +18,4 @@ Theory and Methodology eigenvalue parallelization cmfd + energy_deposition diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 5b26eb38aa..985dc05d61 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/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 diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 1676417c66..c147858eb1 100644 --- a/docs/source/usersguide/tallies.rst +++ b/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 | diff --git a/openmc/deplete/dummy_comm.py b/openmc/deplete/dummy_comm.py index b3fa272648..2648fdedce 100644 --- a/openmc/deplete/dummy_comm.py +++ b/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) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index c50ffd6ab6..9e3e379d58 100644 --- a/openmc/deplete/helpers.py +++ b/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 # ------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 86dbeb7550..743615bbe9 100644 --- a/openmc/deplete/operator.py +++ b/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 diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 3b5166ed8a..1530d3955c 100644 --- a/openmc/deplete/results_list.py +++ b/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.