Merge branch 'fission_q' into poly_and_fission_q

This commit is contained in:
Sterling Harper 2016-08-08 10:07:00 -05:00
commit ea9d4efdc9
21 changed files with 1022 additions and 39 deletions

Binary file not shown.

View file

@ -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)

View file

@ -0,0 +1,52 @@
.. _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
openmc/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

View file

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

View file

@ -55,6 +55,27 @@ Incident Neutron Data
from fission. It is formatted as a reaction product, described in
:ref:`product`.
**/<nuclide name>/fission_energy_release/**
:Attributes: - **format** (*char[]*) -- The energy-dependence format. Either
'Madland' or 'Sher-Beck'
:Datasets: - **fragments** (*double[]*) -- Polynomial coefficients for energy
released in the form of fragments
- **prompt_neutrons** (*double[]* or :ref:`tabulated <1d_tabulated>`)
-- Energy released in the form of prompt neutrons. Polynomial if
the format is Madland or a table if Sher-Beck.
- **delayed_neutrons** (*double[]*) -- Polynomial coefficients for
energy released in the form of delayed neutrons
- **prompt_photons** (*double[]*) -- Polynomial coefficients for
energy released in the form of prompt photons
- **delayed_photons** (*double[]*) -- Polynomial coefficients for
energy released in the form of delayed photons
- **betas** (*double[]*) -- Polynomial coefficients for
energy released in the form of betas
- **neutrinos** (*double[]*) -- Polynomial coefficients for
energy released in the form of neutrinos
-------------------------------
Thermal Neutron Scattering Data
-------------------------------

View file

@ -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/

View file

@ -1809,6 +1809,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``

View file

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

44
openmc/data/endf_utils.py Normal file
View file

@ -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}

View file

@ -0,0 +1,620 @@
from collections import Callable
from copy import deepcopy
import sys
import h5py
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from .data import ATOMIC_SYMBOL
from .endf_utils import read_float, read_CONT_line, identify_nuclide
from .function import Tabulated1D, 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).
form : str
Format used to compute the energy-dependence of the data. Either
'Sher-Beck' or 'Madland'.
"""
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
self._form = 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])
@property
def form(self):
return self._form
@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
@form.setter
def form(self, form):
cv.check_value('format', form, ('Madland', 'Sher-Beck'))
self._form = form
@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.form = 'Madland'
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:
out.form = 'Sher-Beck'
# EFR and ENP are energy independent. Polynomial is used because it
# has a __call__ attribute that handles Iterable inputs. 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 = Polynomial(group['fragments'].value)
obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value)
obj.prompt_photons = Polynomial(group['prompt_photons'].value)
obj.delayed_photons = Polynomial(group['delayed_photons'].value)
obj.betas = Polynomial(group['betas'].value)
obj.neutrinos = Polynomial(group['neutrinos'].value)
if group.attrs['format'].decode() == 'Madland':
obj.form = 'Madland'
obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value)
elif group.attrs['format'].decode() == 'Sher-Beck':
obj.form = 'Sher-Beck'
obj.prompt_neutrons = Tabulated1D.from_hdf5(
group['prompt_neutrons'])
else:
raise ValueError('Unrecognized energy release format')
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
"""
group.create_dataset('fragments', data=self.fragments.coef)
group.create_dataset('delayed_neutrons',
data=self.delayed_neutrons.coef)
group.create_dataset('prompt_photons',
data=self.prompt_photons.coef)
group.create_dataset('delayed_photons',
data=self.delayed_photons.coef)
group.create_dataset('betas', data=self.betas.coef)
group.create_dataset('neutrinos', data=self.neutrinos.coef)
if self.form == 'Madland':
group.attrs['format'] = np.string_('Madland')
group.create_dataset('prompt_neutrons',
data=self.prompt_neutrons.coef)
q_prompt = (self.fragments + self.prompt_neutrons +
self.prompt_photons + Polynomial((0.0, -1.0)))
group.create_dataset('q_prompt', data=q_prompt.coef)
q_recoverable = (self.fragments + self.prompt_neutrons +
self.delayed_neutrons + self.prompt_photons +
self.delayed_photons + self.betas +
Polynomial((0.0, -1.0)))
group.create_dataset('q_recoverable', data=q_recoverable.coef)
elif self.form == 'Sher-Beck':
group.attrs['format'] = np.string_('Sher-Beck')
self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons')
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.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')

View file

@ -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

View file

@ -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):
@ -111,6 +120,14 @@ for filename in ace_libraries:
# Continuous-energy neutron data
neutron = openmc.data.IncidentNeutron.from_ace(
table, args.metastable)
# 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(neutron.name)
# Determine filename

View file

@ -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

View file

@ -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)

View file

@ -3651,6 +3651,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

View file

@ -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
@ -294,6 +300,47 @@ 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')
call read_attribute(temp, fer_group, 'format')
if (temp == 'Madland') then
! The data uses the Madland format, i.e. polynomials
! Read the prompt Q-value
allocate(Polynomial :: this % fission_q_prompt)
fer_dset = open_dataset(fer_group, '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 == 'Sher-Beck') then
! The data uses the Sher-Beck format. Python has handily converted this
! format to Tabulated1Ds.
! Read the prompt Q-value
allocate(Tabulated1D :: this % fission_q_prompt)
fer_dset = open_dataset(fer_group, '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()

View file

@ -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"

View file

@ -681,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
@ -697,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
@ -710,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
@ -750,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

View file

@ -1 +1 @@
930af242a043f2676a000dbc5a2db6b148edcb31ed8c87dbaa35a8efb37a3be8cff30cdf4dc03f9c5c7eb4021f7e4c3327e64681cdd8fd8722c95c69db850227
1bef757d276362fdcd9405096b4cdcbd894f9215ed406493486a45193729be446c9a12242c887f89b6e209ec5beaaacb04dee2fd61e72b4f5c6a8712b776ed6e

View file

@ -1 +1 @@
a51db2a4efc681805f85968e04411dc33beee0532c202f5179b9a82880ab60a75e53fa9141c81045ea1d2842372f2d8da900326f09382ea61dd80a3c9b43bba1
a6e5480c66e6510687bf281983b3f387da45005bb08d8cd0aa629e53c5a42a1b2fa8d76345ad2df497d85f2b273fbc5edc36105a71a73c1107479ca0af8329f6

View file

@ -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'