Merge pull request #698 from smharper/poly_and_fission_q

Function1D's and fission Q-values

.travis.yml

@@ -42,7 +42,7 @@ install: true
 
 before_script:
   - if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
-      wget https://anl.box.com/shared/static/6pwyfjnufam0sb96kqwwrve6vdn8m7u4.xz -O - | tar -C $HOME -xvJ;
+      wget https://anl.box.com/shared/static/dqkwdl7o4lauo91h3mgrn9qno6a3c8mp.xz -O - | tar -C $HOME -xvJ;
     fi
   - export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
 

data/get_nndc_data.py

@@ -151,5 +151,6 @@ if not response or response.lower().startswith('y'):
     env = os.environ.copy()
     env['PYTHONPATH'] = os.path.join(cwd, '..')
 
-    subprocess.call(['../scripts/openmc-ace-to-hdf5', '-d', 'nndc_hdf5']
+    subprocess.call(['../scripts/openmc-ace-to-hdf5', '-d', 'nndc_hdf5',
+                     '--fission_energy_release', 'fission_Q_data_endfb71.h5']
                     + ace_files, env=env)

docs/source/conf.py

@@ -28,6 +28,9 @@ MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
                 'h5py', 'pandas', 'opencg']
 sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
 
+import numpy as np
+np.polynomial.Polynomial = MagicMock
+
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the

docs/source/io_formats/fission_energy.rst

@@ -0,0 +1,53 @@
+.. _usersguide_fission_energy:
+
+==================================
+Fission Energy Release File Format
+==================================
+
+This file is a compact HDF5 representation of the ENDF MT=1, MF=458 data (see
+ENDF-102_ for details).  It gives the information needed to compute the energy
+carried away from fission reactions by each reaction product (e.g. fragment
+nuclei, neutrons) which depends on the incident neutron energy.  OpenMC is
+distributed with one of these files under
+data/fission_Q_data_endfb71.h5.  More files of this format can be created from
+ENDF files with the
+``openmc.data.write_compact_458_library`` function.  They can be read with the
+``openmc.data.FissionEnergyRelease.from_compact_hdf5`` class method.
+
+:Attributes: - **comment** (*char[]*) -- An optional text comment
+             - **component order** (*char[][]*) -- An array of strings
+               specifying the order each reaction product occurs in the data
+               arrays.  The components use the 2-3 letter abbreviations
+               specified in ENDF-102 e.g. EFR for fission fragments and ENP for
+               prompt neutrons.
+
+**/<nuclide name>/**
+    Nuclides are named by concatenating their atomic symbol and mass number. For
+    example, 'U235' or 'Pu239'.  Metastable nuclides are appended with an
+    '_m' and their metastable number.  For example, 'Am242_m1'
+
+:Datasets:
+          - **data** (*double[][][]*) -- The energy release coefficients.  The
+             first axis indexes the component type.  The second axis specifies
+             values or uncertainties.  The third axis indexes the polynomial
+             order.  If the data uses the Sher-Beck format, then the last axis
+             will have a length of one and ENDF-102 should be consulted for
+             energy dependence.  Otherwise, the data uses the Madland format
+             which is a polynomial of incident energy.
+
+             For example, if 'EFR' is given first in the **component order**
+             attribute and the data uses the Madland format, then the energy
+             released in the form of fission fragments at an incident energy
+             :math:`E` is given by
+
+             .. math::
+                \text{data}[0, 0, 0] + \text{data}[0, 0, 1] \cdot E
+                + \text{data}[0, 0, 2] \cdot E^2 + \ldots
+              
+             And its uncertainty is
+
+             .. math::
+                \text{data}[0, 1, 0] + \text{data}[0, 1, 1] \cdot E
+                + \text{data}[0, 1, 2] \cdot E^2 + \ldots
+
+.. _ENDF-102: http://www.nndc.bnl.gov/endfdocs/ENDF-102-2012.pdf

docs/source/io_formats/index.rst

@@ -15,6 +15,7 @@ Data Files
    nuclear_data
    mgxs_library
    data_wmp
+   fission_energy
 
 ------------
 Output Files

docs/source/io_formats/nuclear_data.rst

@@ -55,6 +55,36 @@ Incident Neutron Data
    from fission. It is formatted as a reaction product, described in
    :ref:`product`.
 
+**/<nuclide name>/fission_energy_release/**
+
+:Datasets: - **fragments** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of fragments as a function of incident
+             neutron energy.
+           - **prompt_neutrons** (:ref:`polynomial <1d_polynomial>` or
+             :ref:`tabulated <1d_tabulated>`) -- Energy released in the form of
+             prompt neutrons as a function of incident neutron energy.
+           - **delayed_neutrons** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of delayed neutrons as a function of incident
+             neutron energy.
+           - **prompt_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of prompt photons as a function of incident
+             neutron energy.
+           - **delayed_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of delayed photons as a function of incident
+             neutron energy.
+           - **betas** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of betas as a function of incident
+             neutron energy.
+           - **neutrinos** (:ref:`polynomial <1d_polynomial>`) -- Energy
+             released in the form of neutrinos as a function of incident
+             neutron energy.
+           - **q_prompt** (:ref:`polynomial <1d_polynomial>` or
+             :ref:`tabulated <1d_tabulated>`) -- The prompt fission Q-value
+             (fragments + prompt neutrons + prompt photons - incident energy)
+           - **q_recoverable** (:ref:`polynomial <1d_polynomial>` or
+             :ref:`tabulated <1d_tabulated>`) -- The recoverable fission Q-value
+             (Q_prompt + delayed neutrons + delayed photons + betas)
+
 -------------------------------
 Thermal Neutron Scattering Data
 -------------------------------
@@ -142,17 +172,19 @@ Tabulated
 :Object type: Dataset
 :Datatype: *double[2][]*
 :Description: x-values are listed first followed by corresponding y-values
-:Attributes: - **type** (*char[]*) -- 'tabulated'
+:Attributes: - **type** (*char[]*) -- 'Tabulated1D'
              - **breakpoints** (*int[]*) -- Region breakpoints
              - **interpolation** (*int[]*) -- Region interpolation codes
 
+.. _1d_polynomial:
+
 Polynomial
 ----------
 
 :Object type: Dataset
 :Datatype: *double[]*
 :Description: Polynomial coefficients listed in order of increasing power
-:Attributes: - **type** (*char[]*) -- 'polynomial'
+:Attributes: - **type** (*char[]*) -- 'Polynomial'
 
 Coherent elastic scattering
 ---------------------------

docs/source/pythonapi/index.rst

@@ -348,6 +348,7 @@ Core Classes
     openmc.data.Tabulated1D
     openmc.data.ThermalScattering
     openmc.data.CoherentElastic
+    openmc.data.FissionEnergyRelease
 
 Angle-Energy Distributions
 --------------------------
@@ -381,21 +382,22 @@ Classes
 +++++++
 
 .. autosummary::
-   :toctree: generated
-   :nosignatures:
-   :template: myclass.rst
+    :toctree: generated
+    :nosignatures:
+    :template: myclass.rst
 
-   openmc.data.ace.Library
-   openmc.data.ace.Table
+    openmc.data.ace.Library
+    openmc.data.ace.Table
 
 Functions
 +++++++++
 
 .. autosummary::
-   :toctree: generated
-   :nosignatures:
+    :toctree: generated
+    :nosignatures:
 
-   openmc.data.ace.ascii_to_binary
+    openmc.data.ace.ascii_to_binary
+    openmc.data.write_compact_458_library
 
 .. _Jupyter: https://jupyter.org/
 .. _NumPy: http://www.numpy.org/

docs/source/usersguide/input.rst

@@ -1838,6 +1838,27 @@ The ``<tally>`` element accepts the following sub-elements:
         |                      |:math:`\gamma`-rays are assumed to deposit their   |
         |                      |energy locally. Units are MeV per source particle. |
         +----------------------+---------------------------------------------------+
+        |fission-q-prompt      |The prompt fission energy production rate. This    |
+        |                      |energy comes in the form of fission fragment       |
+        |                      |nuclei, prompt neutrons, and prompt                |
+        |                      |:math:`\gamma`-rays. This value depends on the     |
+        |                      |incident energy and it requires that the nuclear   |
+        |                      |data library contains the optional fission energy  |
+        |                      |release data. Energy is assumed to be deposited    |
+        |                      |locally. Units are MeV per source particle.        |
+        +----------------------+---------------------------------------------------+
+        |fission-q-recoverable |The recoverable fission energy production rate.    |
+        |                      |This energy comes in the form of fission fragment  |
+        |                      |nuclei, prompt and delayed neutrons, prompt and    |
+        |                      |delayed :math:`\gamma`-rays, and delayed           |
+        |                      |:math:`\beta`-rays. This tally differs from the    |
+        |                      |kappa-fission tally in that it is dependent on     |
+        |                      |incident neutron energy and it requires that the   |
+        |                      |nuclear data library contains the optional fission |
+        |                      |energy release data. Energy is assumed to be       |
+        |                      |deposited locally. Units are MeV per source        |
+        |                      |paticle.                                           |
+        +----------------------+---------------------------------------------------+
 
     .. note::
        The ``analog`` estimator is actually identical to the ``collision``

openmc/data/__init__.py

@@ -14,3 +14,4 @@ from .nbody import *
 from .thermal import *
 from .urr import *
 from .library import *
+from .fission_energy import *

openmc/data/endf_utils.py

@@ -0,0 +1,44 @@
+"""This module contains a few utility functions for reading ENDF_ data. It is by
+no means enough to read an entire ENDF file.  For a more complete ENDF reader,
+see Pyne_.
+
+.. _ENDF: http://www.nndc.bnl.gov/endf
+.. _Pyne: http://www.pyne.io
+
+"""
+
+import re
+
+def read_float(float_string):
+    """Parse ENDF 6E11.0 formatted string into a float."""
+    assert len(float_string) == 11
+    pattern = r'([\s\-]\d+\.\d+)([\+\-]\d+)'
+    return float(re.sub(pattern, r'\1e\2', float_string))
+
+
+def read_CONT_line(line):
+    """Parse 80-column line from ENDF CONT record into floats and ints."""
+    return (read_float(line[0:11]), read_float(line[11:22]), int(line[22:33]),
+            int(line[33:44]), int(line[44:55]), int(line[55:66]),
+            int(line[66:70]), int(line[70:72]), int(line[72:75]),
+            int(line[75:80]))
+
+
+def identify_nuclide(fname):
+    """Read the header of an ENDF file and extract identifying information."""
+    with open(fname, 'r') as fh:
+        # Skip the tape id (TPID).
+        line = fh.readline()
+
+        # Read the first HEAD and CONT info.
+        line = fh.readline()
+        ZA, AW, LRP, LFI, NLIB, NMOD, MAT, MF, MT, NS = read_CONT_line(line)
+        line = fh.readline()
+        ELIS, STA, LIS, LISO, junk, NFOR, MAT, MF, MT, NS = read_CONT_line(line)
+
+    # Return dictionary of the most important identifying information.
+    return {'Z': int(ZA) // 1000,
+            'A': int(ZA) % 1000,
+            'LFI': bool(LFI),
+            'LIS': LIS,
+            'LISO': LISO}

openmc/data/fission_energy.py

@@ -0,0 +1,592 @@
+from collections import Callable
+from copy import deepcopy
+import sys
+
+import h5py
+import numpy as np
+
+from .data import ATOMIC_SYMBOL
+from .endf_utils import read_float, read_CONT_line, identify_nuclide
+from .function import Function1D, Tabulated1D, Polynomial, Sum
+import openmc.checkvalue as cv
+
+if sys.version_info[0] >= 3:
+    basestring = str
+
+
+def _extract_458_data(filename):
+    """Read an ENDF file and extract the MF=1, MT=458 values.
+
+    Parameters
+    ----------
+    filename : str
+        Path to and ENDF file
+
+    Returns
+    -------
+    value : dict of str to list of float
+        Dictionary that gives lists of coefficients for each energy component.
+        The keys are the 2-3 letter strings used in ENDF-102, e.g. 'EFR' and
+        'ET'.  The list will have a length of 1 for Sher-Beck data, more for
+        polynomial data.
+    uncertainty : dict of str to list of float
+        A dictionary with the same format as above.  This is probably a
+        one-standard deviation value, but that is not specified explicitly in
+        ENDF-102.  Also, some evaluations will give zero uncertainty.  Use with
+        caution.
+
+    """
+    ident = identify_nuclide(filename)
+
+    if not ident['LFI']:
+        # This nuclide isn't fissionable.
+        return None
+
+    # Extract the MF=1, MT=458 section.
+    lines = []
+    with open(filename, 'r') as fh:
+        line = fh.readline()
+        while line != '':
+            if line[70:75] == ' 1458':
+                lines.append(line)
+            line = fh.readline()
+
+    if len(lines) == 0:
+        # No 458 data here.
+        return None
+
+    # Read the number of coefficients in this LIST record.
+    NPL = read_CONT_line(lines[1])[4]
+
+    # Parse the ENDF LIST into an array.
+    data = []
+    for i in range(NPL):
+        row, column = divmod(i, 6)
+        data.append(read_float(lines[2 + row][11*column:11*(column+1)]))
+
+    # Declare the coefficient names and the order they are given in.  The LIST
+    # contains a value followed immediately by an uncertainty for each of these
+    # components, times the polynomial order + 1.
+    labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET')
+
+    # Associate each set of values and uncertainties with its label.
+    value = {}
+    uncertainty = {}
+    for i, label in enumerate(labels):
+        value[label] = data[2*i::18]
+        uncertainty[label] = data[2*i + 1::18]
+
+    # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not
+    # converted from MeV to eV.  Check for this error and fix it if present.
+    n_coeffs = len(value['EFR'])
+    if n_coeffs == 3:  # Only check 2nd-order data.
+        # Check each energy component for the error.  If a 1 MeV neutron
+        # causes a change of more than 100 MeV, we know something is wrong.
+        error_present = False
+        for coeffs in value.values():
+            second_order = coeffs[2]
+            if abs(second_order) * 1e12 > 1e8:
+                error_present = True
+                break
+
+        # If we found the error, reduce all 2nd-order coeffs by 10**6.
+        if error_present:
+            for coeffs in value.values(): coeffs[2] *= 1e-6
+            for coeffs in uncertainty.values(): coeffs[2] *= 1e-6
+
+    # Convert eV to MeV.
+    for coeffs in value.values():
+        for i in range(len(coeffs)):
+            coeffs[i] *= 10**(-6 + 6*i)
+    for coeffs in uncertainty.values():
+        for i in range(len(coeffs)):
+            coeffs[i] *= 10**(-6 + 6*i)
+
+    return value, uncertainty
+
+
+def write_compact_458_library(endf_files, output_name='fission_Q_data.h5',
+                              comment=None, verbose=False):
+    """Read ENDF files, strip the MF=1 MT=458 data and write to small HDF5.
+
+    Parameters
+    ----------
+    endf_files : Collection of str
+        Strings giving the paths to the ENDF files that will be parsed for data.
+    output_name : str
+        Name of the output HDF5 file.  Default is 'fission_Q_data.h5'.
+    comment : str
+        Comment to write in the output HDF5 file.  Defaults to no comment.
+    verbose : bool
+        If True, print the name of each isomer as it is read.  Defaults to
+        False.
+
+    """
+    # Open the output file.
+    out = h5py.File(output_name, 'w', libver='latest')
+
+    # Write comments, if given.  This commented out comment is the one used for
+    # the library distributed with OpenMC.
+    #comment = ('This data is extracted from ENDF/B-VII.1 library.  Thanks '
+    #           'evaluators, for all your hard work :)  Citation:  '
+    #           'M. B. Chadwick, M. Herman, P. Oblozinsky, '
+    #           'M. E. Dunn, Y. Danon, A. C. Kahler, D. L. Smith, '
+    #           'B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, '
+    #           'D. A. Brown, R. Capote, A. D. Carlson, Y. S. Cho, H. Derrien, '
+    #           'K. Guber, G. M. Hale, S. Hoblit, S. Holloway, T. D. Johnson, '
+    #           'T. Kawano, B. C. Kiedrowski, H. Kim, S. Kunieda, '
+    #           'N. M. Larson, L. Leal, J. P. Lestone, R. C. Little, '
+    #           'E. A. McCutchan, R. E. MacFarlane, M. MacInnes, '
+    #           'C. M. Mattoon, R. D. McKnight, S. F. Mughabghab, '
+    #           'G. P. A. Nobre, G. Palmiotti, A. Palumbo, M. T. Pigni, '
+    #           'V. G. Pronyaev, R. O. Sayer, A. A. Sonzogni, N. C. Summers, '
+    #           'P. Talou, I. J. Thompson, A. Trkov, R. L. Vogt, '
+    #           'S. C. van der Marck, A. Wallner, M. C. White, D. Wiarda, '
+    #           'and P. G. Young. ENDF/B-VII.1 nuclear data for science and '
+    #           'technology: Cross sections, covariances, fission product '
+    #           'yields and decay data", Nuclear Data Sheets, '
+    #           '112(12):2887-2996 (2011).')
+    if comment is not None:
+        out.attrs['comment'] = np.string_(comment)
+
+    # Declare the order of the components.  Use fixed-length numpy strings
+    # because they work well with h5py.
+    labels = np.array(('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER',
+                       'ET'), dtype='S3')
+    out.attrs['component order'] = labels
+
+    # Iterate over the given files.
+    if verbose: print('Reading ENDF files:')
+    for fname in endf_files:
+        if verbose: print(fname)
+
+        ident = identify_nuclide(fname)
+
+        # Skip non-fissionable nuclides.
+        if not ident['LFI']: continue
+
+        # Get the important bits.
+        data = _extract_458_data(fname)
+        if data is None: continue
+        value, uncertainty = data
+
+        # Make a group for this isomer.
+        name = ATOMIC_SYMBOL[ident['Z']] + str(ident['A'])
+        if ident['LISO'] != 0:
+            name += '_m' + str(ident['LISO'])
+        nuclide_group = out.create_group(name)
+
+        # Write all the coefficients into one array.  The first dimension gives
+        # the component (e.g. fragments or prompt neutrons); the second switches
+        # between value and uncertainty; the third gives the polynomial order.
+        n_coeffs = len(value['EFR'])
+        data_out = np.zeros((len(labels), 2, n_coeffs))
+        for i, label in enumerate(labels):
+            data_out[i, 0, :] = value[label.decode()]
+            data_out[i, 1, :] = uncertainty[label.decode()]
+        nuclide_group.create_dataset('data', data=data_out)
+
+    out.close()
+
+
+class FissionEnergyRelease(object):
+    """Energy relased by fission reactions.
+
+    Energy is carried away from fission reactions by many different particles.
+    The attributes of this class specify how much energy is released in the form
+    of fission fragments, neutrons, photons, etc.  Each component is also (in
+    general) a function of the incident neutron energy.
+
+    Following a fission reaction, most of the energy release is carried by the
+    daughter nuclei fragments.  These fragments accelerate apart from the
+    Coulomb force on the time scale of ~10^-20 s [1].  Those fragments emit
+    prompt neutrons between ~10^-18 and ~10^-13 s after scission (although some
+    prompt neutrons may come directly from the scission point) [1].  Prompt
+    photons follow with a time scale of ~10^-14 to ~10^-7 s [1].  The fission
+    products then emit delayed neutrons with half lives between 0.1 and 100 s.
+    The remaining fission energy comes from beta decays of the fission products
+    which release beta particles, photons, and neutrinos (that escape the
+    reactor and do not produce usable heat).
+
+    Use the class methods to instantiate this class from an HDF5 or ENDF
+    dataset.  The :meth:`FissionEnergyRelease.from_hdf5` method builds this
+    class from the usual OpenMC HDF5 data files.
+    :meth:`FissionEnergyRelease.from_endf` uses ENDF-formatted data.
+    :meth:`FissionEnergyRelease.from_compact_hdf5` uses a different HDF5 format
+    that is meant to be compact and store the exact same data as the ENDF
+    format.  Files with this format can be generated with the
+    :func:`openmc.data.write_compact_458_library` function.
+
+    References
+    ----------
+    [1] D. G. Madland, "Total prompt energy release in the neutron-induced
+    fission of ^235U, ^238U, and ^239Pu", Nuclear Physics A 772:113--137 (2006).
+    <http://dx.doi.org/10.1016/j.nuclphysa.2006.03.013>
+
+    Attributes
+    ----------
+    fragments : Callable
+        Function that accepts incident neutron energy value(s) and returns the
+        kinetic energy of the fission daughter nuclides (after prompt neutron
+        emission).
+    prompt_neutrons : Callable
+        Function of energy that returns the kinetic energy of prompt fission
+        neutrons.
+    delayed_neutrons : Callable
+        Function of energy that returns the kinetic energy of delayed neutrons
+        emitted from fission products.
+    prompt_photons : Callable
+        Function of energy that returns the kinetic energy of prompt fission
+        photons.
+    delayed_photons : Callable
+        Function of energy that returns the kinetic energy of delayed photons.
+    betas : Callable
+        Function of energy that returns the kinetic energy of delayed beta
+        particles.
+    neutrinos : Callable
+        Function of energy that returns the kinetic energy of neutrinos.
+    recoverable : Callable
+        Function of energy that returns the kinetic energy of all products that
+        can be absorbed in the reactor (all of the energy except for the
+        neutrinos).
+    total : Callable
+        Function of energy that returns the kinetic energy of all products.
+    q_prompt : Callable
+        Function of energy that returns the prompt fission Q-value (fragments +
+        prompt neutrons + prompt photons - incident neutron energy).
+    q_recoverable : Callable
+        Function of energy that returns the recoverable fission Q-value
+        (total release - neutrinos - incident neutron energy).  This value is
+        sometimes referred to as the pseudo-Q-value.
+    q_total : Callable
+        Function of energy that returns the total fission Q-value (total release
+        - incident neutron energy).
+
+    """
+    def __init__(self):
+        self._fragments = None
+        self._prompt_neutrons = None
+        self._delayed_neutrons = None
+        self._prompt_photons = None
+        self._delayed_photons = None
+        self._betas = None
+        self._neutrinos = None
+
+    @property
+    def fragments(self):
+        return self._fragments
+
+    @property
+    def prompt_neutrons(self):
+        return self._prompt_neutrons
+
+    @property
+    def delayed_neutrons(self):
+        return self._delayed_neutrons
+
+    @property
+    def prompt_photons(self):
+        return self._prompt_photons
+
+    @property
+    def delayed_photons(self):
+        return self._delayed_photons
+
+    @property
+    def betas(self):
+        return self._betas
+
+    @property
+    def neutrinos(self):
+        return self._neutrinos
+
+    @property
+    def recoverable(self):
+        return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons,
+                    self.prompt_photons, self.delayed_photons, self.betas])
+
+    @property
+    def total(self):
+        return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons,
+                    self.prompt_photons, self.delayed_photons, self.betas,
+                    self.neutrinos])
+
+    @property
+    def q_prompt(self):
+        return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons,
+                    lambda E: -E])
+
+    @property
+    def q_recoverable(self):
+        return Sum([self.recoverable, lambda E: -E])
+
+    @property
+    def q_total(self):
+        return Sum([self.total, lambda E: -E])
+
+    @fragments.setter
+    def fragments(self, energy_release):
+        cv.check_type('fragments', energy_release, Callable)
+        self._fragments = energy_release
+
+    @prompt_neutrons.setter
+    def prompt_neutrons(self, energy_release):
+        cv.check_type('prompt_neutrons', energy_release, Callable)
+        self._prompt_neutrons = energy_release
+
+    @delayed_neutrons.setter
+    def delayed_neutrons(self, energy_release):
+        cv.check_type('delayed_neutrons', energy_release, Callable)
+        self._delayed_neutrons = energy_release
+
+    @prompt_photons.setter
+    def prompt_photons(self, energy_release):
+        cv.check_type('prompt_photons', energy_release, Callable)
+        self._prompt_photons = energy_release
+
+    @delayed_photons.setter
+    def delayed_photons(self, energy_release):
+        cv.check_type('delayed_photons', energy_release, Callable)
+        self._delayed_photons = energy_release
+
+    @betas.setter
+    def betas(self, energy_release):
+        cv.check_type('betas', energy_release, Callable)
+        self._betas = energy_release
+
+    @neutrinos.setter
+    def neutrinos(self, energy_release):
+        cv.check_type('neutrinos', energy_release, Callable)
+        self._neutrinos = energy_release
+
+    @classmethod
+    def _from_dictionary(cls, energy_release, incident_neutron):
+        """Generate fission energy release data from a dictionary.
+
+        Parameters
+        ----------
+        energy_release : dict of str to list of float
+            Dictionary that gives lists of coefficients for each energy
+            component.  The keys are the 2-3 letter strings used in ENDF-102,
+            e.g. 'EFR' and 'ET'.  The list will have a length of 1 for Sher-Beck
+            data, more for polynomial data.
+
+        incident_neutron : openmc.data.IncidentNeutron
+            Corresponding incident neutron dataset
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease
+            Fission energy release data
+
+        """
+        out = cls()
+
+        # How many coefficients are given for each component?  If we only find
+        # one value for each, then we need to use the Sher-Beck formula for
+        # energy dependence.  Otherwise, it is a polynomial.
+        n_coeffs = len(energy_release['EFR'])
+        if n_coeffs > 1:
+            out.fragments = Polynomial(energy_release['EFR'])
+            out.prompt_neutrons = Polynomial(energy_release['ENP'])
+            out.delayed_neutrons = Polynomial(energy_release['END'])
+            out.prompt_photons = Polynomial(energy_release['EGP'])
+            out.delayed_photons = Polynomial(energy_release['EGD'])
+            out.betas = Polynomial(energy_release['EB'])
+            out.neutrinos = Polynomial(energy_release['ENU'])
+        else:
+            # EFR and ENP are energy independent.  Use 0-order polynomials to
+            # make a constant function.  The energy-dependence of END is
+            # unspecified in ENDF-102 so assume it is independent.
+            out.fragments = Polynomial((energy_release['EFR'][0]))
+            out.prompt_photons = Polynomial((energy_release['EGP'][0]))
+            out.delayed_neutrons = Polynomial((energy_release['END'][0]))
+
+            # EDP, EB, and ENU are linear.
+            out.delayed_photons = Polynomial((energy_release['EGD'][0], -0.075))
+            out.betas = Polynomial((energy_release['EB'][0], -0.075))
+            out.neutrinos = Polynomial((energy_release['ENU'][0], -0.105))
+
+            # Prompt neutrons require nu-data.  It is not clear from ENDF-102
+            # whether prompt or total nu value should be used, but the delayed
+            # neutron fraction is so small that the difference is negligible.
+            # MT=18 (n, fission) might not be available so try MT=19 (n, f) as
+            # well.
+            if 18 in incident_neutron.reactions:
+                nu_prompt = [p for p in incident_neutron[18].products
+                             if p.particle == 'neutron'
+                             and p.emission_mode == 'prompt']
+            elif 19 in incident_neutron.reactions:
+                nu_prompt = [p for p in incident_neutron[19].products
+                             if p.particle == 'neutron'
+                             and p.emission_mode == 'prompt']
+            else:
+                raise ValueError('IncidentNeutron data has no fission '
+                                 'reaction.')
+            if len(nu_prompt) == 0:
+                raise ValueError('Nu data is needed to compute fission energy '
+                                 'release with the Sher-Beck format.')
+            if len(nu_prompt) > 1:
+                raise ValueError('Ambiguous prompt value.')
+            if not isinstance(nu_prompt[0].yield_, Tabulated1D):
+                raise TypeError('Sher-Beck fission energy release currently '
+                                'only supports Tabulated1D nu data.')
+            ENP = deepcopy(nu_prompt[0].yield_)
+            ENP.y = (energy_release['ENP'] + 1.307 * ENP.x
+                     - 8.07 * (ENP.y - ENP.y[0]))
+            out.prompt_neutrons = ENP
+
+        return out
+
+    @classmethod
+    def from_endf(cls, filename, incident_neutron):
+        """Generate fission energy release data from an ENDF file.
+
+        Parameters
+        ----------
+        filename : str
+            Name of the ENDF file containing fission energy release data
+
+        incident_neutron : openmc.data.IncidentNeutron
+            Corresponding incident neutron dataset
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease
+            Fission energy release data
+
+        """
+
+        # Check to make sure this ENDF file matches the expected isomer.
+        ident = identify_nuclide(filename)
+        if ident['Z'] != incident_neutron.atomic_number:
+            raise ValueError('The atomic number of the ENDF evaluation does '
+                             'not match the given IncidentNeutron.')
+        if ident['A'] != incident_neutron.mass_number:
+            raise ValueError('The atomic mass of the ENDF evaluation does '
+                             'not match the given IncidentNeutron.')
+        if ident['LISO'] != incident_neutron.metastable:
+            raise ValueError('The metastable state of the ENDF evaluation does '
+                             'not match the given IncidentNeutron.')
+        if not ident['LFI']:
+            raise ValueError('The ENDF evaluation is not fissionable.')
+
+        # Read the 458 data from the ENDF file.
+        value, uncertainty = _extract_458_data(filename)
+
+        # Build the object.
+        return cls._from_dictionary(value, incident_neutron)
+
+    @classmethod
+    def from_hdf5(cls, group):
+        """Generate fission energy release data from an HDF5 group.
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to read from
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease
+            Fission energy release data
+
+        """
+
+        obj = cls()
+
+        obj.fragments = Function1D.from_hdf5(group['fragments'])
+        obj.prompt_neutrons = Function1D.from_hdf5(group['prompt_neutrons'])
+        obj.delayed_neutrons = Function1D.from_hdf5(group['delayed_neutrons'])
+        obj.prompt_photons = Function1D.from_hdf5(group['prompt_photons'])
+        obj.delayed_photons = Function1D.from_hdf5(group['delayed_photons'])
+        obj.betas = Function1D.from_hdf5(group['betas'])
+        obj.neutrinos = Function1D.from_hdf5(group['neutrinos'])
+
+        return obj
+
+    @classmethod
+    def from_compact_hdf5(cls, fname, incident_neutron):
+        """Generate fission energy release data from a small HDF5 library.
+
+        Parameters
+        ----------
+        fname : str
+            Path to an HDF5 file containing fission energy release data.  This
+            file should have been generated form the
+            :func:`openmc.data.write_compact_458_library` function.
+
+        incident_neutron : openmc.data.IncidentNeutron
+            Corresponding incident neutron dataset
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease or None
+            Fission energy release data for the given nuclide if it is present
+            in the data file
+
+        """
+
+        fin = h5py.File(fname, 'r')
+
+        components = [s.decode() for s in fin.attrs['component order']]
+
+        nuclide_name = ATOMIC_SYMBOL[incident_neutron.atomic_number]
+        nuclide_name += str(incident_neutron.mass_number)
+        if incident_neutron.metastable != 0:
+            nuclide_name += '_m' + str(incident_neutron.metastable)
+
+        if nuclide_name not in fin: return None
+
+        data = {c: fin[nuclide_name + '/data'][i, 0, :]
+                for i, c in enumerate(components)}
+
+        return cls._from_dictionary(data, incident_neutron)
+
+    def to_hdf5(self, group):
+        """Write energy release data to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
+
+        self.fragments.to_hdf5(group, 'fragments')
+        self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons')
+        self.delayed_neutrons.to_hdf5(group, 'delayed_neutrons')
+        self.prompt_photons.to_hdf5(group, 'prompt_photons')
+        self.delayed_photons.to_hdf5(group, 'delayed_photons')
+        self.betas.to_hdf5(group, 'betas')
+        self.neutrinos.to_hdf5(group, 'neutrinos')
+        
+        if isinstance(self.prompt_neutrons, Polynomial):
+            # Add the polynomials for the relevant components together. Use a
+            # Polynomial((0.0, -1.0)) to subtract incident energy.
+            q_prompt = (self.fragments + self.prompt_neutrons +
+                        self.prompt_photons + Polynomial((0.0, -1.0)))
+            q_prompt.to_hdf5(group, 'q_prompt')
+            q_recoverable = (self.fragments + self.prompt_neutrons +
+                             self.delayed_neutrons + self.prompt_photons +
+                             self.delayed_photons + self.betas +
+                             Polynomial((0.0, -1.0)))
+            q_recoverable.to_hdf5(group, 'q_recoverable')
+
+        elif isinstance(self.prompt_neutrons, Tabulated1D):
+            # Make a Tabulated1D and evaluate the polynomial components at the
+            # table x points to get new y points.  Subtract x from y to remove
+            # incident energy.
+            q_prompt = deepcopy(self.prompt_neutrons)
+            q_prompt.y += self.fragments(q_prompt.x)
+            q_prompt.y += self.prompt_photons(q_prompt.x)
+            q_prompt.y -= q_prompt.x
+            q_prompt.to_hdf5(group, 'q_prompt')
+            q_recoverable = q_prompt
+            q_recoverable.y += self.delayed_neutrons(q_recoverable.x)
+            q_recoverable.y += self.delayed_photons(q_recoverable.x)
+            q_recoverable.y += self.betas(q_recoverable.x)
+            q_recoverable.to_hdf5(group, 'q_recoverable')
+
+        else:
+            raise ValueError('Unrecognized energy release format')

openmc/data/function.py

@@ -1,3 +1,4 @@
+from abc import ABCMeta, abstractmethod
 from collections import Iterable, Callable
 from numbers import Real, Integral
 
@@ -9,7 +10,51 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log',
                         4: 'log-linear', 5: 'log-log'}
 
 
-class Tabulated1D(object):
+class Function1D(object):
+    """A function of one independent variable with HDF5 support."""
+
+    __metaclass__ = ABCMeta
+
+    @abstractmethod
+    def __call__(self): pass
+
+    @abstractmethod
+    def to_hdf5(self, group, name='xy'):
+        """Write function to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+        name : str
+            Name of the dataset to create
+
+        """
+        pass
+
+    @classmethod
+    def from_hdf5(cls, dataset):
+        """Generate function from an HDF5 dataset
+
+        Parameters
+        ----------
+        dataset : h5py.Dataset
+            Dataset to read from
+
+        Returns
+        -------
+        openmc.data.Function1D
+            Function read from dataset
+
+        """
+        for subclass in cls.__subclasses__():
+            if dataset.attrs['type'].decode() == subclass.__name__:
+                return subclass.from_hdf5(dataset)
+        raise ValueError("Unrecognized Function1D class: '"
+                         + dataset.attrs['type'].decode() + "'")
+
+
+class Tabulated1D(Function1D):
     """A one-dimensional tabulated function.
 
     This class mirrors the TAB1 type from the ENDF-6 format. A tabulated
@@ -239,7 +284,7 @@ class Tabulated1D(object):
         """
         dataset = group.create_dataset(name, data=np.vstack(
             [self.x, self.y]))
-        dataset.attrs['type'] = np.string_('tab1')
+        dataset.attrs['type'] = np.string_(type(self).__name__)
         dataset.attrs['breakpoints'] = self.breakpoints
         dataset.attrs['interpolation'] = self.interpolation
 
@@ -258,6 +303,10 @@ class Tabulated1D(object):
             Function read from dataset
 
         """
+        if dataset.attrs['type'].decode() != cls.__name__:
+            raise ValueError("Expected an HDF5 attribute 'type' equal to '"
+                             + cls.__name__ + "'")
+
         x = dataset.value[0, :]
         y = dataset.value[1, :]
         breakpoints = dataset.attrs['breakpoints']
@@ -304,6 +353,42 @@ class Tabulated1D(object):
         return Tabulated1D(x, y, breakpoints, interpolation)
 
 
+class Polynomial(np.polynomial.Polynomial, Function1D):
+    def to_hdf5(self, group, name='xy'):
+        """Write polynomial function to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+        name : str
+            Name of the dataset to create
+
+        """
+        dataset = group.create_dataset(name, data=self.coef)
+        dataset.attrs['type'] = np.string_(type(self).__name__)
+
+    @classmethod
+    def from_hdf5(cls, dataset):
+        """Generate function from an HDF5 dataset
+
+        Parameters
+        ----------
+        dataset : h5py.Dataset
+            Dataset to read from
+
+        Returns
+        -------
+        openmc.data.Function1D
+            Function read from dataset
+
+        """
+        if dataset.attrs['type'].decode() != cls.__name__:
+            raise ValueError("Expected an HDF5 attribute 'type' equal to '"
+                             + cls.__name__ + "'")
+        return cls(dataset.value)
+
+
 class Sum(object):
     """Sum of multiple functions.
 

openmc/data/neutron.py

@@ -9,6 +9,7 @@ import h5py
 
 from .data import ATOMIC_SYMBOL, SUM_RULES
 from .ace import Table, get_table
+from .fission_energy import FissionEnergyRelease
 from .function import Tabulated1D, Sum
 from .product import Product
 from .reaction import Reaction, _get_photon_products
@@ -51,6 +52,9 @@ class IncidentNeutron(object):
         Atomic weight ratio of the target nuclide.
     energy : numpy.ndarray
         The energy values (MeV) at which reaction cross-sections are tabulated.
+    fission_energy : None or openmc.data.FissionEnergyRelease
+        The energy released by fission, tabulated by component (e.g. prompt
+        neutrons or beta particles) and dependent on incident neutron energy
     mass_number : int
         Number of nucleons in the nucleus
     metastable : int
@@ -81,6 +85,7 @@ class IncidentNeutron(object):
         self.temperature = temperature
 
         self._energy = None
+        self._fission_energy = None
         self.reactions = OrderedDict()
         self.summed_reactions = OrderedDict()
         self.urr = None
@@ -126,6 +131,10 @@ class IncidentNeutron(object):
     def energy(self):
         return self._energy
 
+    @property
+    def fission_energy(self):
+        return self._fission_energy
+
     @property
     def temperature(self):
         return self._temperature
@@ -186,6 +195,12 @@ class IncidentNeutron(object):
         cv.check_type('energy grid', energy, Iterable, Real)
         self._energy = energy
 
+    @fission_energy.setter
+    def fission_energy(self, fission_energy):
+        cv.check_type('fission energy release', fission_energy,
+                      FissionEnergyRelease)
+        self._fission_energy = fission_energy
+
     @reactions.setter
     def reactions(self, reactions):
         cv.check_type('reactions', reactions, Mapping)
@@ -276,6 +291,11 @@ class IncidentNeutron(object):
             urr_group = g.create_group('urr')
             self.urr.to_hdf5(urr_group)
 
+        # Write fission energy release data
+        if self.fission_energy is not None:
+            fer_group = g.create_group('fission_energy_release')
+            self.fission_energy.to_hdf5(fer_group)
+
         f.close()
 
     @classmethod
@@ -342,6 +362,11 @@ class IncidentNeutron(object):
             urr_group = group['urr']
             data.urr = ProbabilityTables.from_hdf5(urr_group)
 
+        # Read fission energy release data
+        if 'fission_energy_release' in group:
+            fer_group = group['fission_energy_release']
+            data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
+
         return data
 
     @classmethod

openmc/data/product.py

@@ -3,10 +3,9 @@ from numbers import Real
 import sys
 
 import numpy as np
-from numpy.polynomial.polynomial import Polynomial
 
 import openmc.checkvalue as cv
-from .function import Tabulated1D
+from .function import Tabulated1D, Polynomial, Function1D
 from .angle_energy import AngleEnergy
 
 if sys.version_info[0] >= 3:
@@ -36,7 +35,7 @@ class Product(object):
         yield represents particles from prompt and delayed sources.
     particle : str
         What particle the reaction product is.
-    yield_ : float or openmc.data.Tabulated1D or numpy.polynomial.Polynomial
+    yield_ : openmc.data.Function1D
         Yield of secondary particle in the reaction.
 
     """
@@ -47,7 +46,7 @@ class Product(object):
         self.emission_mode = 'prompt'
         self.distribution = []
         self.applicability = []
-        self.yield_ = 1
+        self.yield_ = Polynomial((1,))  # 0-order polynomial i.e. a constant
 
     def __repr__(self):
         if isinstance(self.yield_, Real):
@@ -119,8 +118,7 @@ class Product(object):
 
     @yield_.setter
     def yield_(self, yield_):
-        cv.check_type('product yield', yield_,
-                      (Real, Tabulated1D, Polynomial))
+        cv.check_type('product yield', yield_, Function1D)
         self._yield = yield_
 
     def to_hdf5(self, group):
@@ -138,16 +136,7 @@ class Product(object):
             group.attrs['decay_rate'] = self.decay_rate
 
         # Write yield
-        if isinstance(self.yield_, Tabulated1D):
-            self.yield_.to_hdf5(group, 'yield')
-            dset = group['yield']
-            dset.attrs['type'] = np.string_('tabulated')
-        elif isinstance(self.yield_, Polynomial):
-            dset = group.create_dataset('yield', data=self.yield_.coef)
-            dset.attrs['type'] = np.string_('polynomial')
-        else:
-            dset = group.create_dataset('yield', data=float(self.yield_))
-            dset.attrs['type'] = np.string_('constant')
+        self.yield_.to_hdf5(group, 'yield')
 
         # Write applicability/distribution
         group.attrs['n_distribution'] = len(self.distribution)
@@ -180,13 +169,7 @@ class Product(object):
             p.decay_rate = group.attrs['decay_rate']
 
         # Read yield
-        yield_type = group['yield'].attrs['type'].decode()
-        if yield_type == 'constant':
-            p.yield_ = group['yield'].value
-        elif yield_type == 'polynomial':
-            p.yield_ = Polynomial(group['yield'].value)
-        elif yield_type == 'tabulated':
-            p.yield_ = Tabulated1D.from_hdf5(group['yield'])
+        p.yield_ = Function1D.from_hdf5(group['yield'])
 
         # Read applicability/distribution
         n_distribution = group.attrs['n_distribution']

openmc/data/reaction.py

@@ -5,13 +5,12 @@ from numbers import Real
 from warnings import warn
 
 import numpy as np
-from numpy.polynomial import Polynomial
 
 import openmc.checkvalue as cv
 from openmc.stats import Uniform
 from .angle_distribution import AngleDistribution
 from .angle_energy import AngleEnergy
-from .function import Tabulated1D
+from .function import Tabulated1D, Polynomial
 from .data import REACTION_NAME
 from .product import Product
 from .uncorrelated import UncorrelatedAngleEnergy
@@ -465,7 +464,8 @@ class Reaction(object):
                         idx = ace.jxs[11] + abs(ty) - 101
                         yield_ = Tabulated1D.from_ace(ace, idx)
                     else:
-                        yield_ = abs(ty)
+                        # 0-order polynomial i.e. a constant
+                        yield_ = Polynomial((abs(ty),))
 
                     neutron = Product('neutron')
                     neutron.yield_ = yield_

scripts/openmc-ace-to-hdf5

@@ -25,6 +25,13 @@ follows the NNDC data convention (1000*Z + A + 300 + 100*m), or the MCNP data
 convention (essentially the same as NNDC, except that the first metastable state
 of Am242 is 95242 and the ground state is 95642).
 
+The optional --fission_energy_release argument will accept an HDF5 file
+containing a library of fission energy release (ENDF MF=1 MT=458) data. A
+library built from ENDF/B-VII.1 data is released with OpenMC and can be found at
+openmc/data/fission_Q_data_endb71.h5. This data is necessary for
+'fission-q-prompt' and 'fission-q-recoverable' tallies, but is not needed
+otherwise.
+
 """
 
 class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
@@ -47,6 +54,8 @@ parser.add_argument('--xsdir', help='MCNP xsdir file that lists '
                     'ACE libraries')
 parser.add_argument('--xsdata', help='Serpent xsdata file that lists '
                     'ACE libraries')
+parser.add_argument('--fission_energy_release', help='HDF5 file containing '
+                    'fission energy release data')
 args = parser.parse_args()
 
 if not os.path.isdir(args.destination):
@@ -124,7 +133,16 @@ for filename in ace_libraries:
             except Exception as e:
                 print('Failed to convert {}: {}'.format(table.name, e))
                 continue
-            print('Converting {} (ACE) to {} (HDF5)'.format(table.name, neutron.name))
+
+            # Fission energy release data, if available
+            if args.fission_energy_release is not None:
+                fer = openmc.data.FissionEnergyRelease.from_compact_hdf5(
+                    args.fission_energy_release, neutron)
+                if fer is not None:
+                    neutron.fission_energy = fer
+
+            print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
+                                                            neutron.name))
 
             # Determine filename
             outfile = os.path.join(args.destination,
@@ -137,7 +155,8 @@ for filename in ace_libraries:
         elif table.name.endswith('t'):
             # Thermal scattering data
             thermal = openmc.data.ThermalScattering.from_ace(table)
-            print('Converting {} (ACE) to {} (HDF5)'.format(table.name, thermal.name))
+            print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
+                                                            thermal.name))
 
             # Determine filename
             outfile = os.path.join(args.destination,

src/constants.F90

@@ -289,7 +289,7 @@ module constants
        EVENT_ABSORB  =  2
 
   ! Tally score type
-  integer, parameter :: N_SCORE_TYPES = 21
+  integer, parameter :: N_SCORE_TYPES = 23
   integer, parameter :: &
        SCORE_FLUX               = -1,  & ! flux
        SCORE_TOTAL              = -2,  & ! total reaction rate
@@ -311,7 +311,9 @@ module constants
        SCORE_EVENTS             = -18, & ! number of events
        SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate
        SCORE_PROMPT_NU_FISSION  = -20, & ! prompt neutron production rate
-       SCORE_INVERSE_VELOCITY   = -21    ! flux-weighted inverse velocity
+       SCORE_INVERSE_VELOCITY   = -21, & ! flux-weighted inverse velocity
+       SCORE_FISS_Q_PROMPT      = -22, & ! prompt fission Q-value
+       SCORE_FISS_Q_RECOV       = -23    ! recoverable fission Q-value
 
   ! Maximum scattering order supported
   integer, parameter :: MAX_ANG_ORDER = 10

src/endf.F90

@@ -14,7 +14,7 @@ contains
   pure function reaction_name(MT) result(string)
 
     integer, intent(in) :: MT
-    character(20)       :: string
+    character(MAX_WORD_LEN) :: string
 
     select case (MT)
     ! Special reactions for tallies
@@ -60,6 +60,10 @@ contains
       string = "events"
     case (SCORE_INVERSE_VELOCITY)
       string = "inverse-velocity"
+    case (SCORE_FISS_Q_PROMPT)
+      string = "fission-q-prompt"
+    case (SCORE_FISS_Q_RECOV)
+      string = "fission-q-recoverable"
 
     ! Normal ENDF-based reactions
     case (TOTAL_XS)

src/endf_header.F90

@@ -30,17 +30,6 @@ module endf_header
     end subroutine function1d_from_hdf5_
   end interface
 
-!===============================================================================
-! CONSTANT1D represents a constant one-dimensional function
-!===============================================================================
-
-  type, extends(Function1D) :: Constant1D
-    real(8) :: y
-  contains
-    procedure :: from_hdf5 => constant1d_from_hdf5
-    procedure :: evaluate => constant1d_evaluate
-  end type Constant1D
-
 !===============================================================================
 ! POLYNOMIAL represents a one-dimensional function expressed as a polynomial
 !===============================================================================
@@ -72,25 +61,6 @@ module endf_header
 
 contains
 
-!===============================================================================
-! Constant1D implementation
-!===============================================================================
-
-  subroutine constant1d_from_hdf5(this, dset_id)
-    class(Constant1D), intent(inout) :: this
-    integer(HID_T),    intent(in)    :: dset_id
-
-    call read_dataset(this % y, dset_id)
-  end subroutine constant1d_from_hdf5
-
-  pure function constant1d_evaluate(this, x) result(y)
-    class(Constant1D), intent(in) :: this
-    real(8),           intent(in) :: x
-    real(8)                       :: y
-
-    y = this % y
-  end function constant1d_evaluate
-
 !===============================================================================
 ! Polynomial implementation
 !===============================================================================

src/input_xml.F90

@@ -3661,6 +3661,10 @@ contains
             t % score_bins(j) = SCORE_KAPPA_FISSION
           case ('inverse-velocity')
             t % score_bins(j) = SCORE_INVERSE_VELOCITY
+          case ('fission-q-prompt')
+            t % score_bins(j) = SCORE_FISS_Q_PROMPT
+          case ('fission-q-recoverable')
+            t % score_bins(j) = SCORE_FISS_Q_RECOV
           case ('current')
             t % score_bins(j) = SCORE_CURRENT
             t % type = TALLY_SURFACE_CURRENT

src/nuclide_header.F90

@@ -10,7 +10,7 @@ module nuclide_header
   use constants
   use dict_header, only: DictIntInt
   use endf,        only: reaction_name, is_fission, is_disappearance
-  use endf_header, only: Function1D, Constant1D, Polynomial, Tabulated1D
+  use endf_header, only: Function1D, Polynomial, Tabulated1D
   use error,       only: fatal_error, warning
   use hdf5_interface, only: read_attribute, open_group, close_group, &
        open_dataset, read_dataset, close_dataset, get_shape
@@ -88,6 +88,10 @@ module nuclide_header
     type(DictIntInt) :: reaction_index ! map MT values to index in reactions
                                        ! array; used at tally-time
 
+    ! Fission energy release
+    class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas
+    class(Function1D), allocatable :: fission_q_recov  ! neutrons, gammas, betas
+
   contains
     procedure :: clear => nuclide_clear
     procedure :: print => nuclide_print
@@ -192,6 +196,8 @@ module nuclide_header
     integer(HID_T) :: rxs_group
     integer(HID_T) :: rx_group
     integer(HID_T) :: total_nu
+    integer(HID_T) :: fer_group                 ! fission_energy_release group
+    integer(HID_T) :: fer_dset
     integer(SIZE_T) :: name_len, name_file_len
     integer(HSIZE_T) :: j
     integer(HSIZE_T) :: dims(1)
@@ -251,8 +257,8 @@ module nuclide_header
       call this % urr_data % from_hdf5(urr_group)
 
       ! if the inelastic competition flag indicates that the inelastic cross
-      ! section should be determined from a normal reaction cross section, we need
-      ! to get the index of the reaction
+      ! section should be determined from a normal reaction cross section, we
+      ! need to get the index of the reaction
       if (this % urr_data % inelastic_flag > 0) then
         do i = 1, size(this % reactions)
           if (this % reactions(i) % MT == this % urr_data % inelastic_flag) then
@@ -283,11 +289,9 @@ module nuclide_header
       total_nu = open_dataset(nu_group, 'yield')
       call read_attribute(temp, total_nu, 'type')
       select case (temp)
-      case ('constant')
-        allocate(Constant1D :: this % total_nu)
-      case ('tabulated')
+      case ('Tabulated1D')
         allocate(Tabulated1D :: this % total_nu)
-      case ('polynomial')
+      case ('Polynomial')
         allocate(Polynomial :: this % total_nu)
       end select
       call this % total_nu % from_hdf5(total_nu)
@@ -296,6 +300,43 @@ module nuclide_header
       call close_group(nu_group)
     end if
 
+    ! Read fission energy release data if present
+    call h5ltpath_valid_f(group_id, 'fission_energy_release', .true., exists, &
+                          hdf5_err)
+    if (exists) then
+      fer_group = open_group(group_id, 'fission_energy_release')
+
+      ! Check to see if this is polynomial or tabulated data
+      fer_dset = open_dataset(fer_group, 'q_prompt')
+      call read_attribute(temp, fer_dset, 'type')
+      if (temp == 'Polynomial') then
+        ! Read the prompt Q-value
+        allocate(Polynomial :: this % fission_q_prompt)
+        call this % fission_q_prompt % from_hdf5(fer_dset)
+        call close_dataset(fer_dset)
+
+        ! Read the recoverable energy Q-value
+        allocate(Polynomial :: this % fission_q_recov)
+        fer_dset = open_dataset(fer_group, 'q_recoverable')
+        call this % fission_q_recov % from_hdf5(fer_dset)
+        call close_dataset(fer_dset)
+      else if (temp == 'Tabulated1D') then
+        ! Read the prompt Q-value
+        allocate(Tabulated1D :: this % fission_q_prompt)
+        call this % fission_q_prompt % from_hdf5(fer_dset)
+        call close_dataset(fer_dset)
+
+        ! Read the recoverable energy Q-value
+        allocate(Tabulated1D :: this % fission_q_recov)
+        fer_dset = open_dataset(fer_group, 'q_recoverable')
+        call this % fission_q_recov % from_hdf5(fer_dset)
+        call close_dataset(fer_dset)
+      else
+        call fatal_error('Unrecognized fission energy release format.')
+      end if
+      call close_group(fer_group)
+    end if
+
     ! Create derived cross section data
     call this % create_derived()
 

src/output.F90

@@ -776,6 +776,8 @@ contains
     score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
     score_names(abs(SCORE_PROMPT_NU_FISSION))  = "Prompt-Nu-Fission Rate"
     score_names(abs(SCORE_INVERSE_VELOCITY))   = "Flux-Weighted Inverse Velocity"
+    score_names(abs(SCORE_FISS_Q_PROMPT))      = "Prompt fission power"
+    score_names(abs(SCORE_FISS_Q_RECOV))       = "Recoverable fission power"
 
     ! Create filename for tally output
     filename = trim(path_output) // "tallies.out"

src/product_header.F90

@@ -5,7 +5,7 @@ module product_header
   use angleenergy_header, only: AngleEnergyContainer
   use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, &
        EMISSION_TOTAL, NEUTRON, PHOTON
-  use endf_header, only: Tabulated1D, Function1D, Constant1D, Polynomial
+  use endf_header, only: Tabulated1D, Function1D, Polynomial
   use hdf5_interface, only: read_attribute, open_group, close_group, &
        open_dataset, close_dataset, read_dataset
   use random_lcg, only: prn
@@ -109,11 +109,9 @@ contains
     yield = open_dataset(group_id, 'yield')
     call read_attribute(temp, yield, 'type')
     select case (temp)
-    case ('constant')
-      allocate(Constant1D :: this % yield)
-    case ('tabulated')
+    case ('Tabulated1D')
       allocate(Tabulated1D :: this % yield)
-    case ('polynomial')
+    case ('Polynomial')
       allocate(Polynomial :: this % yield)
     end select
     call this % yield % from_hdf5(yield)

src/state_point.F90

@@ -51,7 +51,7 @@ contains
     integer(HID_T) :: file_id
     integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
                       mesh_group, filter_group, runtime_group
-    character(20), allocatable :: str_array(:)
+    character(MAX_WORD_LEN), allocatable :: str_array(:)
     character(MAX_FILE_LEN)    :: filename
     type(RegularMesh), pointer :: meshp
     type(TallyObject), pointer    :: tally

src/tally.F90

@@ -1,7 +1,6 @@
 module tally
 
   use constants
-  use endf_header,      only: Constant1D
   use error,            only: fatal_error
   use geometry_header
   use global
@@ -247,24 +246,17 @@ contains
         ! reaction with neutrons in the exit channel
         if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
              (p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
-          ! Don't waste time on very common reactions we know have multiplicities
-          ! of one.
+          ! Don't waste time on very common reactions we know have
+          ! multiplicities of one.
           score = p % last_wgt * flux
         else
-          m = nuclides(p%event_nuclide)%reaction_index% &
+          m = nuclides(p % event_nuclide) % reaction_index % &
                get_key(p % event_MT)
 
           ! Get yield and apply to score
-          associate (rxn => nuclides(p%event_nuclide)%reactions(m))
-            select type (yield => rxn % products(1) % yield)
-            type is (Constant1D)
-              ! Grab the yield from the reaction
-              score = p % last_wgt * yield % y * flux
-            class default
-              ! the yield was already incorporated in to p % wgt per the
-              ! scattering routine
-              score = p % wgt * flux
-            end select
+          associate (rxn => nuclides(p % event_nuclide) % reactions(m))
+            score = p % last_wgt * flux &
+                 * rxn % products(1) % yield % evaluate(p % last_E)
           end associate
         end if
 
@@ -289,16 +281,9 @@ contains
                get_key(p % event_MT)
 
           ! Get yield and apply to score
-          associate (rxn => nuclides(p%event_nuclide)%reactions(m))
-            select type (yield => rxn % products(1) % yield)
-            type is (Constant1D)
-              ! Grab the yield from the reaction
-              score = p % last_wgt * yield % y * flux
-            class default
-              ! the yield was already incorporated in to p % wgt per the
-              ! scattering routine
-              score = p % wgt * flux
-            end select
+          associate (rxn => nuclides(p % event_nuclide) % reactions(m))
+            score = p % last_wgt * flux &
+                 * rxn % products(1) % yield % evaluate(p % last_E)
           end associate
         end if
 
@@ -324,15 +309,8 @@ contains
 
           ! Get yield and apply to score
           associate (rxn => nuclides(p%event_nuclide)%reactions(m))
-            select type (yield => rxn % products(1) % yield)
-            type is (Constant1D)
-              ! Grab the yield from the reaction
-              score = p % last_wgt * yield % y * flux
-            class default
-              ! the yield was already incorporated in to p % wgt per the
-              ! scattering routine
-              score = p % wgt * flux
-            end select
+            score = p % last_wgt * flux &
+                 * rxn % products(1) % yield % evaluate(p % last_E)
           end associate
         end if
 
@@ -703,14 +681,14 @@ contains
           if (survival_biasing) then
             ! No fission events occur if survival biasing is on -- need to
             ! calculate fraction of absorptions that would have resulted in
-            ! fission scale by kappa-fission
-            associate (nuc => nuclides(p%event_nuclide))
-              if (micro_xs(p%event_nuclide)%absorption > ZERO .and. &
-                   nuc%fissionable) then
-                score = p%absorb_wgt * &
-                     nuc%reactions(nuc%index_fission(1))%Q_value * &
-                     micro_xs(p%event_nuclide)%fission / &
-                     micro_xs(p%event_nuclide)%absorption * flux
+            ! fission scaled by kappa-fission
+            associate (nuc => nuclides(p % event_nuclide))
+              if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
+                   nuc % fissionable) then
+                score = p % absorb_wgt * &
+                     nuc % reactions(nuc % index_fission(1)) % Q_value * &
+                     micro_xs(p % event_nuclide) % fission / &
+                     micro_xs(p % event_nuclide) % absorption * flux
               end if
             end associate
           else
@@ -719,12 +697,12 @@ contains
             ! All fission events will contribute, so again we can use
             ! particle's weight entering the collision as the estimate for
             ! the fission energy production rate
-            associate (nuc => nuclides(p%event_nuclide))
-              if (nuc%fissionable) then
-                score = p%last_wgt * &
-                     nuc%reactions(nuc%index_fission(1))%Q_value * &
-                     micro_xs(p%event_nuclide)%fission / &
-                     micro_xs(p%event_nuclide)%absorption * flux
+            associate (nuc => nuclides(p % event_nuclide))
+              if (nuc % fissionable) then
+                score = p % last_wgt * &
+                     nuc % reactions(nuc % index_fission(1)) % Q_value * &
+                     micro_xs(p % event_nuclide) % fission / &
+                     micro_xs(p % event_nuclide) % absorption * flux
               end if
             end associate
           end if
@@ -732,22 +710,23 @@ contains
         else
           if (i_nuclide > 0) then
             associate (nuc => nuclides(i_nuclide))
-              if (nuc%fissionable) then
-                score = nuc%reactions(nuc%index_fission(1))%Q_value * &
-                     micro_xs(i_nuclide)%fission * atom_density * flux
+              if (nuc % fissionable) then
+                score = nuc % reactions(nuc % index_fission(1)) % Q_value * &
+                     micro_xs(i_nuclide) % fission * atom_density * flux
               end if
             end associate
           else
-            do l = 1, materials(p%material)%n_nuclides
+            do l = 1, materials(p % material) % n_nuclides
               ! Determine atom density and index of nuclide
-              atom_density_ = materials(p%material)%atom_density(l)
-              i_nuc = materials(p%material)%nuclide(l)
+              atom_density_ = materials(p % material) % atom_density(l)
+              i_nuc = materials(p % material) % nuclide(l)
 
               ! If nuclide is fissionable, accumulate kappa fission
               associate(nuc => nuclides(i_nuc))
                 if (nuc % fissionable) then
-                  score = score + nuc%reactions(nuc%index_fission(1))%Q_value * &
-                       micro_xs(i_nuc)%fission * atom_density_ * flux
+                  score = score + &
+                       nuc % reactions(nuc % index_fission(1)) % Q_value * &
+                       micro_xs(i_nuc) % fission * atom_density_ * flux
                 end if
               end associate
             end do
@@ -772,6 +751,123 @@ contains
           end if
         end if
 
+      case (SCORE_FISS_Q_PROMPT)
+        if (t % estimator == ESTIMATOR_ANALOG) then
+          if (survival_biasing) then
+            ! No fission events occur if survival biasing is on -- need to
+            ! calculate fraction of absorptions that would have resulted in
+            ! fission scaled by Q-value
+            associate (nuc => nuclides(p % event_nuclide))
+              if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
+                   allocated(nuc % fission_q_prompt)) then
+                score = p % absorb_wgt &
+                     * nuc % fission_q_prompt % evaluate(p % last_E) &
+                     * micro_xs(p % event_nuclide) % fission &
+                     / micro_xs(p % event_nuclide) % absorption * flux
+              end if
+            end associate
+          else
+            ! Skip any non-absorption events
+            if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
+            ! All fission events will contribute, so again we can use
+            ! particle's weight entering the collision as the estimate for
+            ! the fission energy production rate
+            associate (nuc => nuclides(p % event_nuclide))
+              if (allocated(nuc % fission_q_prompt)) then
+                score = p % last_wgt &
+                     * nuc % fission_q_prompt % evaluate(p % last_E) &
+                     * micro_xs(p % event_nuclide) % fission &
+                     / micro_xs(p % event_nuclide) % absorption * flux
+              end if
+            end associate
+          end if
+
+        else
+          if (t % estimator == ESTIMATOR_COLLISION) then
+            E = p % last_E
+          else
+            E = p % E
+          end if
+
+          if (i_nuclide > 0) then
+            if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then
+              score = micro_xs(i_nuclide) % fission * atom_density * flux &
+                      * nuclides(i_nuclide) % fission_q_prompt % evaluate(E)
+            else
+              score = ZERO
+            end if
+          else
+            score = ZERO
+            do l = 1, materials(p % material) % n_nuclides
+              atom_density_ = materials(p % material) % atom_density(l)
+              i_nuc = materials(p % material) % nuclide(l)
+              if (allocated(nuclides(i_nuc) % fission_q_prompt)) then
+                score = score + micro_xs(i_nuc) % fission * atom_density_ &
+                        * flux &
+                        * nuclides(i_nuc) % fission_q_prompt % evaluate(E)
+              end if
+            end do
+          end if
+        end if
+
+      case (SCORE_FISS_Q_RECOV)
+        if (t % estimator == ESTIMATOR_ANALOG) then
+          if (survival_biasing) then
+            ! No fission events occur if survival biasing is on -- need to
+            ! calculate fraction of absorptions that would have resulted in
+            ! fission scaled by Q-value
+            associate (nuc => nuclides(p % event_nuclide))
+              if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
+                   allocated(nuc % fission_q_recov)) then
+                score = p % absorb_wgt &
+                     * nuc % fission_q_recov % evaluate(p % last_E) &
+                     * micro_xs(p % event_nuclide) % fission &
+                     / micro_xs(p % event_nuclide) % absorption * flux
+              end if
+            end associate
+          else
+            ! Skip any non-absorption events
+            if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
+            ! All fission events will contribute, so again we can use
+            ! particle's weight entering the collision as the estimate for
+            ! the fission energy production rate
+            associate (nuc => nuclides(p % event_nuclide))
+              if (allocated(nuc % fission_q_recov)) then
+                score = p % last_wgt &
+                     * nuc % fission_q_recov % evaluate(p % last_E) &
+                     * micro_xs(p % event_nuclide) % fission &
+                     / micro_xs(p % event_nuclide) % absorption * flux
+              end if
+            end associate
+          end if
+
+        else
+          if (t % estimator == ESTIMATOR_COLLISION) then
+            E = p % last_E
+          else
+            E = p % E
+          end if
+
+          if (i_nuclide > 0) then
+            if (allocated(nuclides(i_nuclide) % fission_q_recov)) then
+              score = micro_xs(i_nuclide) % fission * atom_density * flux &
+                      * nuclides(i_nuclide) % fission_q_recov % evaluate(E)
+            else
+              score = ZERO
+            end if
+          else
+            score = ZERO
+            do l = 1, materials(p % material) % n_nuclides
+              atom_density_ = materials(p % material) % atom_density(l)
+              i_nuc = materials(p % material) % nuclide(l)
+              if (allocated(nuclides(i_nuc) % fission_q_recov)) then
+                score = score + micro_xs(i_nuc) % fission * atom_density_ &
+                        * flux * nuclides(i_nuc) % fission_q_recov % evaluate(E)
+              end if
+            end do
+          end if
+        end if
+
       case default
         if (t % estimator == ESTIMATOR_ANALOG) then
           ! Any other score is assumed to be a MT number. Thus, we just need

tests/test_tallies/inputs_true.dat

@@ -1 +1 @@
-930af242a043f2676a000dbc5a2db6b148edcb31ed8c87dbaa35a8efb37a3be8cff30cdf4dc03f9c5c7eb4021f7e4c3327e64681cdd8fd8722c95c69db850227
+1bef757d276362fdcd9405096b4cdcbd894f9215ed406493486a45193729be446c9a12242c887f89b6e209ec5beaaacb04dee2fd61e72b4f5c6a8712b776ed6e

tests/test_tallies/results_true.dat

@@ -1 +1 @@
-a51db2a4efc681805f85968e04411dc33beee0532c202f5179b9a82880ab60a75e53fa9141c81045ea1d2842372f2d8da900326f09382ea61dd80a3c9b43bba1
+5a0f3f1ae244ada7d8c9f444d7a98c2589a9720e78174a276cf96162df6728bab6d534f136f65cb0386c3235eb7d5db47a0dd6504636ca77e7eb375e6978a84d

tests/test_tallies/test_tallies.py

@@ -123,8 +123,9 @@ class TalliesTestHarness(PyAPITestHarness):
             t.filters = [cell_filter]
             t.scores = ['absorption', 'delayed-nu-fission', 'events', 'fission',
                         'inverse-velocity', 'kappa-fission', '(n,2n)', '(n,n1)',
-                        '(n,gamma)', 'nu-fission', 'scatter', 'elastic', 'total',
-                        'prompt-nu-fission']
+                        '(n,gamma)', 'nu-fission', 'scatter', 'elastic', 
+                        'total', 'prompt-nu-fission', 'fission-q-prompt',
+                        'fission-q-recoverable']
         score_tallies[0].estimator = 'tracklength'
         score_tallies[1].estimator = 'analog'
         score_tallies[2].estimator = 'collision'