From 0c02e58b7854e8df82b73f261ff1916cadb51779 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 13 Sep 2019 14:06:42 -0500 Subject: [PATCH 01/10] Documentation for modifications to total heating --- docs/source/methods/energy_deposition.rst | 134 ++++++++++++++++++++++ docs/source/methods/index.rst | 1 + docs/source/usersguide/tallies.rst | 5 +- 3 files changed, 138 insertions(+), 2 deletions(-) create mode 100644 docs/source/methods/energy_deposition.rst diff --git a/docs/source/methods/energy_deposition.rst b/docs/source/methods/energy_deposition.rst new file mode 100644 index 0000000000..5bd16802db --- /dev/null +++ b/docs/source/methods/energy_deposition.rst @@ -0,0 +1,134 @@ +.. _methods_heating: + +============================= +Heating and Energy Deposition +============================= + +As particles traverse a problem, some portion of their energy is deposited at +collision sites. There are a variety of mechanisms that contribute to the +energy deposition, including down-scattering from higher to lower energies, +fission events, and reactions with positive Q-values. 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 numbers are the product of reaction-specific coefficients and +a reaction cross section + +.. math:: + + H(E) = \phi(E)\sum_i\rho_i\sum_rk_{i, r} + +and has units energy per time, typically eV / s. +Here, :math:`k_{i, r}` are the KERMA [Kinetic Energy Release in Materials] +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 secondary neutrons and photons from the sum of +incident neutron energy, :math:`E`, and the reaction :math:`Q` value. + +--------------------------- +The Special Case of 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 MT=458 data in a standard ENDF/B-VII file. +All these quantities have some energy dependence, 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 are 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 number into +a fission and non-fission section, one can write + +.. math:: + + H_i(E) = H_{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 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 remove the assumption of local photon energy. +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 ``301`` is included in :attr:`openmc.Tally.scores`. 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/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 | From 016fc0e43e8de91827f5409411839b7473efacb5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 10 Sep 2019 15:22:23 -0500 Subject: [PATCH 02/10] 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 d566f3080f05bf9765ebc4b7210c134c0c7bb8aa) --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/helpers.py | 42 +++++++++++++++++++++++++++++++ openmc/deplete/operator.py | 31 ++++++++++++++++++----- 3 files changed, 68 insertions(+), 6 deletions(-) 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/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index c50ffd6ab6..5fa9f0654d 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,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 # ------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 86dbeb7550..110d7a7a3f 100644 --- a/openmc/deplete/operator.py +++ b/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] From ca410a82fe88c4103e1b6a3715881d3fe5eefae3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 16 Sep 2019 15:12:42 -0500 Subject: [PATCH 03/10] Ensure strings when setting openmc.lib.Tally.scores In using lists of containing integer reactions as the scores, e.g. ``` t = openmc.lib.Tally() t.scores = [901] ``` The setter would fail as the integers do not have an encode method. This converts all incoming scores to strings prior to passing values to the shared library. --- openmc/lib/tally.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index 97b7437ea3..e6dd3ddbff 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -319,7 +319,7 @@ class Tally(_FortranObjectWithID): @scores.setter def scores(self, scores): scores_ = (c_char_p * len(scores))() - scores_[:] = [x.encode() for x in scores] + scores_[:] = [str(x).encode() for x in scores] _dll.openmc_tally_set_scores(self._index, len(scores), scores_) @property From 74e48521a8cfddf921b680a5b322c911521a0129 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 16 Sep 2019 15:24:06 -0500 Subject: [PATCH 04/10] Pass MT301 or MT901 to EnergyHelper from Operator Depending on settings.photon_transport, obtain the total system energy from the 301 [heating] or 901 [heating-local] scores. The former is used in coupled neutron-photon transport, as this does not include any energy from neutrons or photons, assuming these particles deposit their energy along their life. MT901 is used for neutron transport and is the default, e.g. if scores is None. This score includes energy from prompt and delayed photons taken from MT458 data. Related PR: #1344 - Ability to generate KERMAs assuming local photon energy deposition --- openmc/deplete/helpers.py | 26 ++++++++++++++++++++++++-- openmc/deplete/operator.py | 8 +++++--- 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 5fa9f0654d..d994ab68c5 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -10,6 +10,7 @@ from collections import defaultdict from numpy import dot, zeros, newaxis from . import comm +import openmc.data from openmc.checkvalue import check_type, check_greater_than from openmc.lib import ( Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) @@ -161,6 +162,12 @@ class ChainFissionHelper(EnergyHelper): class EnergyScoreHelper(EnergyHelper): """Class responsible for obtaining system energy via a tally score + Parameters + ---------- + reaction_mt : int or None + Valid score to use when obtaining system energy from openmc. + Defaults to 901 [heating assuming local photons] + Attributes ---------- nuclides : list of str @@ -170,13 +177,28 @@ class EnergyScoreHelper(EnergyHelper): 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 : int + MT reaction number that is scored. """ - def __init__(self): + def __init__(self, score=None): super().__init__() + self.score = score self._tally = None + @property + def score(self): + return self._score + + @score.setter + def score(self, value): + if value is None: + self._score = 901 + else: + check_type("score", value, int) + self._score = value + def prepare(self, *args, **kwargs): """Create a tally for system energy production @@ -185,7 +207,7 @@ class EnergyScoreHelper(EnergyHelper): """ self._tally = Tally() - self._tally.scores = ["energy-deposition"] + self._tally.scores = [self.score] def reset(self): """Obtain system energy from tally diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 110d7a7a3f..eef4acdedf 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -221,9 +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" - else EnergyScoreHelper()) + if energy_mode == "fission-q": + self._energy_helper = ChainFissionHelper() + else: + score = 301 if settings.photon_transport else 901 + self._energy_helper = EnergyScoreHelper(score) # Select and create fission yield helper fission_helper = self._fission_helpers[fission_yield_mode] From 452ad5ac7dd231c5f77b662c7b7e08c0a35822e3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 16 Sep 2019 16:32:29 -0500 Subject: [PATCH 05/10] Avoid zero division when normalizing reaction rates If the system energy is found to be zero, then Operator.__call__ will exit the simulation using comm.Abort / sys.exit if no mpi4py. Without this, the reaction rates will be scaled by power / 0, causing errors downstream in setting material compositions. --- openmc/deplete/dummy_comm.py | 5 +++++ openmc/deplete/operator.py | 20 +++++++++++++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) 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/operator.py b/openmc/deplete/operator.py index eef4acdedf..c9b991d1d1 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -7,6 +7,7 @@ densities is all done in-memory instead of through the filesystem. """ +import sys import copy from collections import OrderedDict from itertools import chain @@ -237,6 +238,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 @@ -272,7 +277,16 @@ class Operator(TransportOperator): time_openmc = time.time() # Extract results - op_result = self._unpack_tallies_and_normalize(power) + try: + op_result = self._unpack_tallies_and_normalize(power) + except ZeroDivisionError: + if comm.rank == 0: + sys.stderr.flush() + print(" No energy reported from openmc tallies. Do you have " + "MT901 data?\n", file=sys.stderr, flush=True) + comm.barrier() + comm.Abort(1) + return copy.deepcopy(op_result) @@ -655,6 +669,10 @@ class Operator(TransportOperator): # J / s / source neutron energy = comm.allreduce(self._energy_helper.energy) + # Guard against divide by zero + if energy == 0: + raise ZeroDivisionError + # Scale reaction rates to obtain units of reactions/sec rates *= power / energy From 952cd2a3cf590c762198161bbec4dfc4b12d842f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 16 Sep 2019 16:49:23 -0500 Subject: [PATCH 06/10] Doc fix for ResultsList.get_depletion_time --- openmc/deplete/results_list.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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. From 6b386f8447071212561be030f7bb9b7ed237406f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 17 Sep 2019 11:34:20 -0500 Subject: [PATCH 07/10] Use string score names for EnergyScoreHelper Operator defaults to using "heating-local" if not running in coupled photon transport. --- openmc/deplete/helpers.py | 25 ++++++------------------- openmc/deplete/operator.py | 2 +- openmc/lib/tally.py | 2 +- 3 files changed, 8 insertions(+), 21 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index d994ab68c5..8e87c532ab 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -10,7 +10,6 @@ from collections import defaultdict from numpy import dot, zeros, newaxis from . import comm -import openmc.data from openmc.checkvalue import check_type, check_greater_than from openmc.lib import ( Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) @@ -164,9 +163,9 @@ class EnergyScoreHelper(EnergyHelper): Parameters ---------- - reaction_mt : int or None - Valid score to use when obtaining system energy from openmc. - Defaults to 901 [heating assuming local photons] + score : string + Valid score to use when obtaining system energy from OpenMC. + Defaults to "heating-local" Attributes ---------- @@ -177,28 +176,16 @@ class EnergyScoreHelper(EnergyHelper): 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 : int - MT reaction number that is scored. + score : str + Score used to obtain system energy """ - def __init__(self, score=None): + def __init__(self, score="heating-local"): super().__init__() self.score = score self._tally = None - @property - def score(self): - return self._score - - @score.setter - def score(self, value): - if value is None: - self._score = 901 - else: - check_type("score", value, int) - self._score = value - def prepare(self, *args, **kwargs): """Create a tally for system energy production diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index c9b991d1d1..3558759b4e 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -225,7 +225,7 @@ class Operator(TransportOperator): if energy_mode == "fission-q": self._energy_helper = ChainFissionHelper() else: - score = 301 if settings.photon_transport else 901 + score = "heating" if settings.photon_transport else "heating-local" self._energy_helper = EnergyScoreHelper(score) # Select and create fission yield helper diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index e6dd3ddbff..97b7437ea3 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -319,7 +319,7 @@ class Tally(_FortranObjectWithID): @scores.setter def scores(self, scores): scores_ = (c_char_p * len(scores))() - scores_[:] = [str(x).encode() for x in scores] + scores_[:] = [x.encode() for x in scores] _dll.openmc_tally_set_scores(self._index, len(scores), scores_) @property From cfbd4f2fd698ce92ef0bc0ca8590f18c86aa6b9e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 17 Sep 2019 12:17:49 -0500 Subject: [PATCH 08/10] Operator aborts in unpacking, not call, if no energy produced --- openmc/deplete/operator.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 3558759b4e..5ea1f120de 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -277,16 +277,7 @@ class Operator(TransportOperator): time_openmc = time.time() # Extract results - try: - op_result = self._unpack_tallies_and_normalize(power) - except ZeroDivisionError: - if comm.rank == 0: - sys.stderr.flush() - print(" No energy reported from openmc tallies. Do you have " - "MT901 data?\n", file=sys.stderr, flush=True) - comm.barrier() - comm.Abort(1) - + op_result = self._unpack_tallies_and_normalize(power) return copy.deepcopy(op_result) @@ -671,7 +662,12 @@ class Operator(TransportOperator): # Guard against divide by zero if energy == 0: - raise ZeroDivisionError + 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 From 2ab3ad7b789b0e1b773eb0a1ad56a4fc430c63a7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 17 Sep 2019 12:18:51 -0500 Subject: [PATCH 09/10] Minor tweaks to Operator, EnergyScoreHelper Per reviewer comments --- openmc/deplete/helpers.py | 2 +- openmc/deplete/operator.py | 7 +------ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 8e87c532ab..9e3e379d58 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -205,7 +205,7 @@ class EnergyScoreHelper(EnergyHelper): the number of MPI processes """ super().reset() - if not comm.rank: + if comm.rank == 0: self._energy = self._tally.results[0, 0, 1] # ------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 5ea1f120de..743615bbe9 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -84,8 +84,7 @@ class Operator(TransportOperator): 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. + 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``. Only applicable @@ -261,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() @@ -274,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) From 17ded504432043e3c4fc6443403b99a552762137 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 17 Sep 2019 12:19:30 -0500 Subject: [PATCH 10/10] Improve energy deposition methodology --- docs/source/methods/energy_deposition.rst | 95 +++++++++++++---------- 1 file changed, 55 insertions(+), 40 deletions(-) diff --git a/docs/source/methods/energy_deposition.rst b/docs/source/methods/energy_deposition.rst index 5bd16802db..af4962c553 100644 --- a/docs/source/methods/energy_deposition.rst +++ b/docs/source/methods/energy_deposition.rst @@ -5,23 +5,23 @@ Heating and Energy Deposition ============================= As particles traverse a problem, some portion of their energy is deposited at -collision sites. There are a variety of mechanisms that contribute to the -energy deposition, including down-scattering from higher to lower energies, -fission events, and reactions with positive Q-values. The information describing -how much energy is deposited for a specific reaction is referred to as +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 numbers are the product of reaction-specific coefficients and +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} + 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] -coefficients for reaction :math:`r` of isotope :math:`i`. +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. @@ -33,14 +33,16 @@ 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) + 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 secondary neutrons and photons from the sum of -incident neutron energy, :math:`E`, and the reaction :math:`Q` value. +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. ---------------------------- -The Special Case of Fission ---------------------------- +------- +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 @@ -54,14 +56,16 @@ broken up into the following categories: - :math:`E_{\beta}` - energy of released :math:`\beta` particles - :math:`E_{\nu}` - energy of neutrinos -These components are defined in MT=458 data in a standard ENDF/B-VII file. -All these quantities have some energy dependence, 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 +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} + 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. @@ -75,34 +79,35 @@ NJOY computes the fission KERMA coefficient using this energy-balance method to .. note:: - The energy from delayed neutrons and photons and beta particles are intentionally - left out from the NJOY calculations + 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 number into -a fission and non-fission section, one can write +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:: - H_i(E) = H_{i, nf}(E) + \left[E_{fr}(E) + E_{\gamma, p}\right]\sigma_{i, f}(E) + 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 MT=458 data. +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. +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 @@ -119,11 +124,12 @@ if ``901`` is included in :attr:`openmc.Tally.scores`. Coupled neutron-photon transport -------------------------------- -Here, OpenMC instructs ``heatr`` to remove the assumption of local photon energy. -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 +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:: @@ -131,4 +137,13 @@ heating number returned from this ``heatr`` run and :math:`M301` be + \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 ``301`` is included in :attr:`openmc.Tally.scores`. +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