mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Add Madland fission-Q support to openmc.data
This commit is contained in:
Sterling Harper
parent
c8a188ca26
commit
a38f53d26f
4 changed files with 351 additions and 0 deletions
|
|
@ -14,3 +14,4 @@ from .nbody import *
|
|||
from .thermal import *
|
||||
from .urr import *
|
||||
from .library import *
|
||||
from .fission_energy import *
|
||||
|
|
|
|||
43
openmc/data/endf_utils.py
Normal file
43
openmc/data/endf_utils.py
Normal file
|
|
@ -0,0 +1,43 @@
|
|||
"""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 = '([\s\\-]\d+\\.\d+)([\\+\\-]\d+)'
|
||||
mantissa, exponent = re.match(pattern, float_string).groups()
|
||||
return float(mantissa + 'e' + exponent)
|
||||
|
||||
|
||||
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,
|
||||
'LIS': LIS,
|
||||
'LISO': LISO}
|
||||
285
openmc/data/fission_energy.py
Normal file
285
openmc/data/fission_energy.py
Normal file
|
|
@ -0,0 +1,285 @@
|
|||
from collections import Callable
|
||||
import sys
|
||||
#from warnings import warn
|
||||
|
||||
import numpy as np
|
||||
from numpy.polynomial.polynomial import Polynomial
|
||||
|
||||
from .function import Tabulated1D, Sum
|
||||
from .endf_utils import read_float, read_CONT_line, identify_nuclide
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
if sys.version_info[0] >= 3:
|
||||
basestring = str
|
||||
|
||||
|
||||
class FissionEnergyRelease(object):
|
||||
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 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_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:
|
||||
pass
|
||||
if ident['A'] != incident_neutron.mass_number:
|
||||
pass
|
||||
if ident['LISO'] != incident_neutron.metastable:
|
||||
pass
|
||||
|
||||
# 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()
|
||||
|
||||
# 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. If we only find
|
||||
# one value for each of these components, then we need to use the
|
||||
# Sher-Beck formula for energy dependence. Otherwise, it is a
|
||||
# polynomial.
|
||||
labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET')
|
||||
|
||||
# Associate each set of values and uncertainties with its label.
|
||||
value = dict()
|
||||
uncertainty = dict()
|
||||
for i in range(len(labels)):
|
||||
value[labels[i]] = data[2*i::18]
|
||||
uncertainty[labels[i]] = 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
|
||||
|
||||
# Perform the sanity check again... just in case.
|
||||
for coeffs in value.values():
|
||||
second_order = coeffs[2]
|
||||
if abs(second_order) * 1e12 > 1e8:
|
||||
raise ValueError("Encountered a ludicrously large second-"
|
||||
"order polynomial coefficient.")
|
||||
|
||||
# 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)
|
||||
|
||||
out = cls()
|
||||
if n_coeffs > 1:
|
||||
out.form = 'Madland'
|
||||
out.fragments = Polynomial(value['EFR'])
|
||||
out.prompt_neutrons = Polynomial(value['ENP'])
|
||||
out.delayed_neutrons = Polynomial(value['END'])
|
||||
out.prompt_photons = Polynomial(value['EGP'])
|
||||
out.delayed_photons = Polynomial(value['EGD'])
|
||||
out.betas = Polynomial(value['EB'])
|
||||
out.neutrinos = Polynomial(value['ENU'])
|
||||
else:
|
||||
out.form = 'Sher-Beck'
|
||||
raise NotImplemented
|
||||
|
||||
return out
|
||||
|
||||
@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()
|
||||
if group.attrs['format'] == 'Madland':
|
||||
obj.fragments = Polynomial(group['fragments'].value)
|
||||
obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].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)
|
||||
elif group.attrs['format'] == 'Sher-Beck':
|
||||
raise NotImplemented
|
||||
else:
|
||||
raise ValueError('Unrecognized energy release format')
|
||||
|
||||
return obj
|
||||
|
||||
def to_hdf5(self, group):
|
||||
"""Write energy release data to an HDF5 group
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to write to
|
||||
|
||||
"""
|
||||
|
||||
if self.form == 'Madland':
|
||||
group.attrs['format'] = np.string_('Madland')
|
||||
group.create_dataset('fragments', data=self.fragments.coef)
|
||||
group.create_dataset('prompt_neutrons',
|
||||
data=self.prompt_neutrons.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)
|
||||
elif self.form == 'Sher-Beck':
|
||||
group.attrs['format'] = np.string_('Sher-Beck')
|
||||
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')
|
||||
else:
|
||||
raise ValueError('Unrecognized energy release format')
|
||||
|
|
@ -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
|
||||
|
|
@ -81,6 +82,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 +128,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 +192,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 +288,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
|
||||
|
|
@ -331,6 +348,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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue