crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add MT=458 data library and supporting methods

Browse source
This commit is contained in:
Sterling Harper 2016-08-02 14:10:32 -05:00
parent 725ef2d923
commit 84cfd2d6a9
5 changed files with 155 additions and 63 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
data/fission_Q_data_endfb71.h5 Normal file
View file

Binary file not shown.

3
data/get_nndc_data.py
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)

195
openmc/data/fission_energy.py
View file

@ -305,6 +305,89 @@ class FissionEnergyRelease(object):
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.
@ -327,72 +410,22 @@ class FissionEnergyRelease(object):
# Check to make sure this ENDF file matches the expected isomer.
ident = identify_nuclide(filename)
if ident['Z'] != incident_neutron.atomic_number:
pass
raise ValueError('The atomic number of the ENDF evaluation does '
'not match the given IncidentNeutron.')
if ident['A'] != incident_neutron.mass_number:
pass
raise ValueError('The atomic mass of the ENDF evaluation does '
'not match the given IncidentNeutron.')
if ident['LISO'] != incident_neutron.metastable:
pass
if not ident['LIF']:
pass
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)
# Declare the coefficient names. 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')
# How many coefficients are given for each coefficient? 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(value['EFR'])
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'
# 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((value['EFR'][0]))
out.prompt_photons = Polynomial((value['EGP'][0]))
out.delayed_neutrons = Polynomial((value['END'][0]))
# EDP, EB, and ENU are linear.
out.delayed_photons = Polynomial((value['EGD'][0], -0.075))
out.betas = Polynomial((value['EB'][0], -0.075))
out.neutrinos = Polynomial((value['ENU'][0], -0.105))
# Prompt neutrons require nu-data. It is not clear from ENDF-102
# whether prompt or total nu values should be used, but the delayed
# neutron fraction is so small that the difference is negligible.
nu_prompt = [p for p in incident_neutron[18].products
if p.particle == 'neutron'
and p.emission_mode == 'prompt']
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 nu 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 = value['ENP'] + 1.307 * ENP.x - 8.07 * (ENP.y - ENP.y[0])
out.prompt_neutrons = ENP
return out
# Build the object.
return cls._from_dictionary(value, incident_neutron)
@classmethod
def from_hdf5(cls, group):
@ -419,9 +452,9 @@ class FissionEnergyRelease(object):
obj.betas = Polynomial(group['betas'].value)
obj.neutrinos = Polynomial(group['neutrinos'].value)
if group.attrs['format'] == 'Madland':
if group.attrs['format'].decode() == 'Madland':
obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value)
elif group.attrs['format'] == 'Sher-Beck':
elif group.attrs['format'].decode() == 'Sher-Beck':
obj.prompt_neutrons = Tabulated1D.from_hdf5(
group['prompt_neutrons'])
else:
@ -429,6 +462,44 @@ class FissionEnergyRelease(object):
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
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 = str(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

3
openmc/data/neutron.py
View file

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

17
scripts/openmc-ace-to-hdf5
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