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

Folded multiple temperatures in to each h5 data file; now im getting 50% less storage required of endf70. yay! Next I think I can get rid of some s(a,b) information, will investigate that and then on to revising openMC so it can read and use this data

Browse source
This commit is contained in:
Adam Nelson 2016-08-21 06:05:19 -04:00
parent fd6b22b454
commit 3b65dde0fd
5 changed files with 657 additions and 222 deletions
Download patch file Download diff file Expand all files Collapse all files

4
openmc/data/data.py
View file

@ -218,3 +218,7 @@ def atomic_mass(isotope):
isotope = isotope[:isotope.find('_')]
return _ATOMIC_MASS.get(isotope.lower())
def kT_to_K(kT):
K = kT / 8.6173324e-11
return K

288
openmc/data/neutron.py
View file

@ -7,7 +7,7 @@ from warnings import warn
import numpy as np
import h5py
from .data import ATOMIC_SYMBOL, SUM_RULES
from .data import ATOMIC_SYMBOL, SUM_RULES, kT_to_K
from .ace import Table, get_table
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum
@ -21,6 +21,64 @@ if sys.version_info[0] >= 3:
basestring = str
def _get_metadata(zaid, metastable_scheme='nndc'):
"""Method to obtain the complete element name, element, Z, mass_number,
and metastable state
Parameters
----------
zaid : int
ZAID (1000*Z + A) obtained from a library
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
Returns
-------
"""
cv.check_type('zaid', zaid, int)
cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp'])
Z = zaid // 1000
mass_number = zaid % 1000
if metastable_scheme == 'mcnp':
if zaid > 1000000:
# New SZA format
Z = Z % 1000
if zaid == 1095242:
metastable = 0
else:
metastable = zaid // 1000000
else:
if zaid == 95242:
metastable = 1
elif zaid == 95642:
metastable = 0
else:
metastable = 1 if mass_number > 300 else 0
elif metastable_scheme == 'nndc':
metastable = 1 if mass_number > 300 else 0
while mass_number > 3 * Z:
mass_number -= 100
# Determine name
element = ATOMIC_SYMBOL[Z]
name = '{}{}'.format(element, mass_number)
if metastable > 0:
name += '_m{}'.format(metastable)
return (name, element, Z, mass_number, metastable)
class IncidentNeutron(EqualityMixin):
"""Continuous-energy neutron interaction data.
@ -40,8 +98,9 @@ class IncidentNeutron(EqualityMixin):
Metastable state of the nucleus. A value of zero indicates ground state.
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
temperature : float
Temperature of the target nuclide in MeV.
kTs : Iterable float
List of temperatures of the target nuclide in the data set.
The temperatures have units of MeV.
Attributes
----------
@ -51,8 +110,10 @@ class IncidentNeutron(EqualityMixin):
Atomic symbol of the nuclide, e.g., 'Zr'
atomic_weight_ratio : float
Atomic weight ratio of the target nuclide.
energy : numpy.ndarray
energy : dict of numpy.ndarray
The energy values (MeV) at which reaction cross-sections are tabulated.
They keys of the dict are the temperature string ('296.3K') for each
set of energies
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
@ -69,23 +130,28 @@ class IncidentNeutron(EqualityMixin):
summed_reactions : collections.OrderedDict
Contains summed cross sections, e.g., the total cross section. The keys
are the MT values and the values are Reaction objects.
temperature : float
Temperature of the target nuclide in MeV.
temperatures : Iterable of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings with 1 decimal place,
i.e., '293.6K'
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of MeV.
urr : None or openmc.data.ProbabilityTables
Unresolved resonance region probability tables
"""
def __init__(self, name, atomic_number, mass_number, metastable,
atomic_weight_ratio, temperature):
atomic_weight_ratio, kTs):
self.name = name
self.atomic_number = atomic_number
self.mass_number = mass_number
self.metastable = metastable
self.atomic_weight_ratio = atomic_weight_ratio
self.temperature = temperature
self._energy = None
self.kTs = kTs
self.temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs]
self.energy = {}
self._fission_energy = None
self.reactions = OrderedDict()
self.summed_reactions = OrderedDict()
@ -128,18 +194,10 @@ class IncidentNeutron(EqualityMixin):
def atomic_weight_ratio(self):
return self._atomic_weight_ratio
@property
def energy(self):
return self._energy
@property
def fission_energy(self):
return self._fission_energy
@property
def temperature(self):
return self._temperature
@property
def reactions(self):
return self._reactions
@ -159,7 +217,7 @@ class IncidentNeutron(EqualityMixin):
@property
def atomic_symbol(self):
return atomic_symbol[self.atomic_number]
return ATOMIC_SYMBOL[self.atomic_number]
@atomic_number.setter
def atomic_number(self, atomic_number):
@ -185,17 +243,6 @@ class IncidentNeutron(EqualityMixin):
cv.check_greater_than('atomic weight ratio', atomic_weight_ratio, 0.0)
self._atomic_weight_ratio = atomic_weight_ratio
@temperature.setter
def temperature(self, temperature):
cv.check_type('temperature', temperature, Real)
cv.check_greater_than('temperature', temperature, 0.0, True)
self._temperature = temperature
@energy.setter
def energy(self, energy):
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,
@ -218,6 +265,87 @@ class IncidentNeutron(EqualityMixin):
(ProbabilityTables, type(None)))
self._urr = urr
def add_temperature_from_ace(self, ace_or_filename, metastable_scheme='nndc'):
"""Add data to the IncidentNeutron object from an ACE file at a
different temperature.
Parameters
----------
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
"""
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
ace = get_table(ace_or_filename)
# Obtain the information needed to check if this ACE file is for the
# same nuclide.
zaid, xs = ace.name.split('.')
name, element, Z, mass_number, metastable = \
_get_metadata(int(zaid), metastable_scheme)
# If this ACE data matches the data within self then get the data
if ace.temperature not in self.kTs:
if name == self.name:
# Add temperature and kTs
strT = "{0:.1f}K".format(kT_to_K(ace.temperature))
self.temperatures.append(strT)
self.kTs.append(ace.temperature)
# Read energy grid
n_energy = ace.nxs[3]
energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy]
self.energy[strT] = energy
total_xs = \
Tabulated1D(energy, ace.xss[ace.jxs[1] +
n_energy:ace.jxs[1] +
2 * n_energy])
abs_xs = Tabulated1D(energy, ace.xss[ace.jxs[1] + 2 *
n_energy:ace.jxs[1] +
3 * n_energy])
self.summed_reactions[1].add_temperature(strT, 0, total_xs)
if 27 in self.summed_reactions:
self.summed_reactions[27].add_temperature(strT, 0, abs_xs)
# Read each reaction and get the xs data out of it
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
rx = Reaction.from_ace(ace, i)
xsdata = list(rx.T_data.values())[0]
self.reactions[rx.mt].add_temperature(strT,
xsdata.threshold_idx,
xsdata.xs)
# Obtain data for the summed photon reactions
for mt in self.summed_reactions:
if mt not in [1, 27]:
# Create summed appropriate cross section
mts = self.get_reaction_components(mt)
xsdata = Sum([self.reactions[mt_i].T_data[strT].xs
for mt_i in mts])
self.summed_reactions[mt].add_temperature(strT, 0,
xsdata)
else:
raise ValueError('Data provided for an incorrect nuclide')
else:
raise Warning('Temperature data set already within '
'IncidentNeutron object')
def get_reaction_components(self, mt):
"""Determine what reactions make up summed reaction.
@ -271,10 +399,12 @@ class IncidentNeutron(EqualityMixin):
g.attrs['A'] = self.mass_number
g.attrs['metastable'] = self.metastable
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['temperature'] = self.temperature
g.attrs['kTs'] = self.kTs
# Write energy grid
g.create_dataset('energy', data=self.energy)
eg = g.create_group('energy')
for temperature in self.temperatures:
eg.create_dataset(temperature, data=self.energy[temperature])
# Write reaction data
rxs_group = g.create_group('reactions')
@ -327,19 +457,22 @@ class IncidentNeutron(EqualityMixin):
mass_number = group.attrs['A']
metastable = group.attrs['metastable']
atomic_weight_ratio = group.attrs['atomic_weight_ratio']
temperature = group.attrs['temperature']
kTs = group.attrs['kTs'].tolist()
temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs]
data = cls(name, atomic_number, mass_number, metastable,
atomic_weight_ratio, temperature)
atomic_weight_ratio, kTs)
# Read energy grid
data.energy = group['energy'].value
e_group = group['energy']
for temperature in temperatures:
data.energy[temperature] = e_group[temperature].value
# Read reaction data
rxs_group = group['reactions']
for name, obj in sorted(rxs_group.items()):
if name.startswith('reaction_'):
rx = Reaction.from_hdf5(obj, data.energy)
rx = Reaction.from_hdf5(obj, data.energy, temperatures)
data.reactions[rx.mt] = rx
# Read total nu data if available
@ -351,12 +484,15 @@ class IncidentNeutron(EqualityMixin):
# MTs never depend on lower MTs.
for mt_sum in sorted(SUM_RULES, reverse=True):
if mt_sum not in data:
xs_components = [data[mt].xs for mt in SUM_RULES[mt_sum]
if mt in data]
if len(xs_components) > 0:
rxn = Reaction(mt_sum)
rxn.xs = Sum(xs_components)
data.summed_reactions[mt_sum] = rxn
for it, T in enumerate(data.temperatures):
xs_components = \
[data[mt].T_data[T].xs for mt in SUM_RULES[mt_sum]
if mt in data]
if len(xs_components) > 0:
if it == 0:
data.summed_reactions[mt_sum] = Reaction(mt_sum)
data.summed_reactions[mt_sum].add_temperature(
T, 0, Sum(xs_components))
# Read unresolved resonance probability tables
if 'urr' in group:
@ -376,9 +512,9 @@ class IncidentNeutron(EqualityMixin):
Parameters
----------
ace : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If the value is a string, it is assumed to
be the filename for the ACE file.
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
@ -394,6 +530,8 @@ class IncidentNeutron(EqualityMixin):
Incident neutron continuous-energy data
"""
# First obtain the data for the first provided ACE table/file
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
@ -401,54 +539,37 @@ class IncidentNeutron(EqualityMixin):
# If mass number hasn't been specified, make an educated guess
zaid, xs = ace.name.split('.')
zaid = int(zaid)
Z = zaid // 1000
mass_number = zaid % 1000
name, element, Z, mass_number, metastable = \
_get_metadata(int(zaid), metastable_scheme)
if metastable_scheme == 'mcnp':
if zaid > 1000000:
# New SZA format
Z = Z % 1000
if zaid == 1095242:
metastable = 0
else:
metastable = zaid // 1000000
else:
if zaid == 95242:
metastable = 1
elif zaid == 95642:
metastable = 0
else:
metastable = 1 if mass_number > 300 else 0
elif metastable_scheme == 'nndc':
metastable = 1 if mass_number > 300 else 0
# Assign temperature to the running list
kTs = [ace.temperature]
temperatures = ["{0:.1f}K".format(kT_to_K(ace.temperature))]
while mass_number > 3*Z:
mass_number -= 100
# Determine name for group
element = ATOMIC_SYMBOL[Z]
if metastable > 0:
name = '{}{}_m{}.{}'.format(element, mass_number, metastable, xs)
else:
name = '{}{}.{}'.format(element, mass_number, xs)
# If mass number hasn't been specified, make an educated guess
zaid, xs = ace.name.split('.')
name, element, Z, mass_number, metastable = \
_get_metadata(int(zaid), metastable_scheme)
data = cls(name, Z, mass_number, metastable,
ace.atomic_weight_ratio, ace.temperature)
ace.atomic_weight_ratio, kTs)
# Read energy grid
n_energy = ace.nxs[3]
energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy]
data.energy = energy
total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2*n_energy]
absorption_xs = ace.xss[ace.jxs[1] + 2*n_energy:ace.jxs[1] + 3*n_energy]
data.energy[temperatures[0]] = energy
total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2 * n_energy]
absorption_xs = ace.xss[ace.jxs[1] + 2*n_energy:ace.jxs[1] +
3 * n_energy]
# Create summed reactions (total and absorption)
total = Reaction(1)
total.xs = Tabulated1D(energy, total_xs)
total.add_temperature(temperatures[-1], 0,
Tabulated1D(energy, total_xs))
data.summed_reactions[1] = total
absorption = Reaction(27)
absorption.xs = Tabulated1D(energy, absorption_xs)
absorption.add_temperature(temperatures[-1], 0,
Tabulated1D(energy, absorption_xs))
data.summed_reactions[27] = absorption
# Read each reaction
@ -478,7 +599,12 @@ class IncidentNeutron(EqualityMixin):
warn('Photon production is present for MT={} but no '
'reaction components exist.'.format(mt))
continue
rx.xs = Sum([data.reactions[mt_i].xs for mt_i in mts])
threshold_idx = \
np.amin([data.reactions[mt_i].T_data[temperatures[-1]].xs.x[0]
for mt_i in mts])
xsvals = Sum([data.reactions[mt_i].T_data[temperatures[-1]].xs
for mt_i in mts])
rx.add_temperature(temperatures[-1], threshold_idx, xsvals)
# Determine summed cross section
rx.products += _get_photon_products(ace, rx)

141
openmc/data/reaction.py
View file

@ -1,7 +1,7 @@
from __future__ import division, unicode_literals
from collections import Iterable, Callable
from copy import deepcopy
from numbers import Real
from numbers import Real, Integral
from warnings import warn
import numpy as np
@ -12,7 +12,7 @@ from openmc.stats import Uniform
from .angle_distribution import AngleDistribution
from .angle_energy import AngleEnergy
from .function import Tabulated1D, Polynomial
from .data import REACTION_NAME
from .data import REACTION_NAME, kT_to_K
from .product import Product
from .uncorrelated import UncorrelatedAngleEnergy
@ -211,7 +211,7 @@ def _get_photon_products(ace, rx):
# Get photon production cross section
photon_prod_xs = ace.xss[idx + 2:idx + 2 + n_energy]
neutron_xs = rx.xs(energy)
neutron_xs = list(rx.T_data.values())[0].xs(energy)
idx = np.where(neutron_xs > 0.)
# Calculate photon yield
@ -251,6 +251,78 @@ def _get_photon_products(ace, rx):
return photons
class XS(EqualityMixin):
"""An energy-dependent cross-section for a reaction channel
Parameters
----------
threshold_idx : int
The index on the energy grid corresponding to the threshold of this
reaction.
xs : dict of callable
Microscopic cross section for this reaction as a function of incident
energy; these cross sections are provided in a dictionary where the key
is the temperature of the cross section set.
Attributes
----------
threshold : float
Threshold of the reaction in MeV
threshold_idx : int
The index on the energy grid corresponding to the threshold of this
reaction.
xs : dict of callable
Microscopic cross section for this reaction as a function of incident
energy; these cross sections are provided in a dictionary where the key
is the temperature of the cross section set.
"""
def __init__(self, threshold_idx, xs):
self._threshold_idx = threshold_idx
self._xs = xs
@property
def threshold_idx(self):
return self._threshold_idx
@property
def threshold(self):
return self.xs.x[0]
@property
def xs(self):
return self._xs
@threshold_idx.setter
def threshold_idx(self, threshold_idx):
cv.check_type('threshold_idx', threshold_idx, Integral)
cv.check_greater_than('threshold_idx', threshold_idx, 0, equality=True)
self._threshold_idx = threshold_idx
@xs.setter
def xs(self, xs):
cv.check_type('reaction cross section', xs, Callable)
if isinstance(xs, Tabulated1D):
for y in xs.y:
cv.check_greater_than('reaction cross section', y, 0.0, True)
self._xs = xs
def to_hdf5(self, group):
"""Write XS to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['threshold_idx'] = self.threshold_idx + 1
if self.xs is not None:
group.create_dataset('xs', data=self.xs.y)
class Reaction(EqualityMixin):
"""A nuclear reaction
@ -281,9 +353,10 @@ class Reaction(EqualityMixin):
threshold_idx : int
The index on the energy grid corresponding to the threshold of this
reaction.
xs : callable
T_data : dict of openmc.data.XS
Microscopic cross section for this reaction as a function of incident
energy
energy; these cross sections are provided in a dictionary where the key
is the temperature of the cross section set.
products : Iterable of openmc.data.Product
Reaction products
derived_products : Iterable of openmc.data.Product
@ -296,8 +369,7 @@ class Reaction(EqualityMixin):
self.center_of_mass = True
self.mt = mt
self.q_value = 0.
self.threshold_idx = 0
self._xs = None
self.T_data = {}
self.products = []
self.derived_products = []
@ -319,14 +391,6 @@ class Reaction(EqualityMixin):
def products(self):
return self._products
@property
def threshold(self):
return self.xs.x[0]
@property
def xs(self):
return self._xs
@center_of_mass.setter
def center_of_mass(self, center_of_mass):
cv.check_type('center of mass', center_of_mass, (bool, np.bool_))
@ -342,13 +406,9 @@ class Reaction(EqualityMixin):
cv.check_type('reaction products', products, Iterable, Product)
self._products = products
@xs.setter
def xs(self, xs):
cv.check_type('reaction cross section', xs, Callable)
if isinstance(xs, Tabulated1D):
for y in xs.y:
cv.check_greater_than('reaction cross section', y, 0.0, True)
self._xs = xs
def add_temperature(self, temperature, threshold_idx, xs):
cv.check_type('temperature', temperature, str)
self.T_data[temperature] = XS(threshold_idx, xs)
def to_hdf5(self, group):
"""Write reaction to an HDF5 group
@ -366,16 +426,16 @@ class Reaction(EqualityMixin):
else:
group.attrs['label'] = np.string_(self.mt)
group.attrs['Q_value'] = self.q_value
group.attrs['threshold_idx'] = self.threshold_idx + 1
group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0
if self.xs is not None:
group.create_dataset('xs', data=self.xs.y)
for T in self.T_data:
Tgroup = group.create_group(T)
self.T_data[T].to_hdf5(Tgroup)
for i, p in enumerate(self.products):
pgroup = group.create_group('product_{}'.format(i))
p.to_hdf5(pgroup)
@classmethod
def from_hdf5(cls, group, energy):
def from_hdf5(cls, group, energy, temperatures):
"""Generate reaction from an HDF5 group
Parameters
@ -384,6 +444,8 @@ class Reaction(EqualityMixin):
HDF5 group to write to
energy : Iterable of float
Array of energies at which cross sections are tabulated at
temperatures : Iterable of float
Array of temperatures at which to obtain the cross sections
Returns
-------
@ -391,16 +453,19 @@ class Reaction(EqualityMixin):
Reaction data
"""
mt = group.attrs['mt']
rx = cls(mt)
rx.q_value = group.attrs['Q_value']
rx.threshold_idx = group.attrs['threshold_idx'] - 1
rx.center_of_mass = bool(group.attrs['center_of_mass'])
# Read cross section
if 'xs' in group:
xs = group['xs'].value
rx.xs = Tabulated1D(energy[rx.threshold_idx:], xs)
# Read cross section data
for T in temperatures:
Tgroup = group[T]
if 'xs' in Tgroup:
threshold_idx = Tgroup.attrs['threshold_idx'] - 1
xs = Tabulated1D(energy[T][threshold_idx:], Tgroup['xs'].value)
rx.T_data[T] = XS(threshold_idx, xs)
# Determine number of products
n_product = 0
@ -421,6 +486,10 @@ class Reaction(EqualityMixin):
n_grid = ace.nxs[3]
grid = ace.xss[ace.jxs[1]:ace.jxs[1] + n_grid]
# Convert data temperature to a "300.0K" number for indexing
# temperature data
strT = "{0:.1f}K".format(kT_to_K(ace.temperature))
if i_reaction > 0:
mt = int(ace.xss[ace.jxs[3] + i_reaction - 1])
rx = cls(mt)
@ -435,11 +504,11 @@ class Reaction(EqualityMixin):
loc = int(ace.xss[ace.jxs[6] + i_reaction - 1])
# Determine starting index on energy grid
rx.threshold_idx = int(ace.xss[ace.jxs[7] + loc - 1]) - 1
threshold_idx = int(ace.xss[ace.jxs[7] + loc - 1]) - 1
# Determine number of energies in reaction
n_energy = int(ace.xss[ace.jxs[7] + loc])
energy = grid[rx.threshold_idx:rx.threshold_idx + n_energy]
energy = grid[threshold_idx:threshold_idx + n_energy]
# Read reaction cross section
xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy]
@ -450,7 +519,8 @@ class Reaction(EqualityMixin):
"to zero.".format(rx.mt, ace.name))
xs[xs < 0.0] = 0.0
rx.xs = Tabulated1D(energy, xs)
tabulated_xs = Tabulated1D(energy, xs)
rx.T_data[strT] = XS(threshold_idx, tabulated_xs)
# ==================================================================
# YIELD AND ANGLE-ENERGY DISTRIBUTION
@ -509,7 +579,8 @@ class Reaction(EqualityMixin):
"Setting to zero.".format(ace.name))
elastic_xs[elastic_xs < 0.0] = 0.0
rx.xs = Tabulated1D(grid, elastic_xs)
tabulated_xs = Tabulated1D(grid, elastic_xs)
rx.T_data[strT] = XS(0, tabulated_xs)
# No energy distribution for elastic scattering
neutron = Product('neutron')

354
openmc/data/thermal.py
View file

@ -8,6 +8,7 @@ import h5py
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import kT_to_K
from .ace import Table, get_table
from .angle_energy import AngleEnergy
from .function import Tabulated1D
@ -89,7 +90,7 @@ class CoherentElastic(EqualityMixin):
if isinstance(E, Iterable):
E = np.asarray(E)
idx = np.searchsorted(self.bragg_edges, E)
return self.factors[idx]/E
return self.factors[idx] / E
def __len__(self):
@ -159,8 +160,9 @@ class ThermalScattering(EqualityMixin):
ZAID identifier of the table, e.g. lwtr.10t.
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
temperature : float
Temperature of the target nuclide in eV.
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of MeV.
Attributes
----------
@ -174,22 +176,29 @@ class ThermalScattering(EqualityMixin):
approximation
name : str
Name of the table, e.g. lwtr.20t.
temperature : float
Temperature of the target nuclide in eV.
temperatures : Iterable of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings with 1 decimal place,
i.e., '293.6K'
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of MeV.
zaids : Iterable of int
ZAID identifiers that the thermal scattering data applies to
"""
def __init__(self, name, atomic_weight_ratio, temperature):
def __init__(self, name, atomic_weight_ratio, kTs):
self.name = name
self.atomic_weight_ratio = atomic_weight_ratio
self.temperature = temperature
self.elastic_xs = None
self.elastic_mu_out = None
self.inelastic_xs = None
self.inelastic_e_out = None
self.inelastic_mu_out = None
self.kTs = kTs
self.temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs]
self.elastic_xs = {}
self.elastic_mu_out = {}
self.inelastic_xs = {}
self.inelastic_e_out = {}
self.inelastic_mu_out = {}
self.inelastic_dist = {}
self.secondary_mode = None
self.zaids = []
@ -217,26 +226,191 @@ class ThermalScattering(EqualityMixin):
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['temperature'] = self.temperature
g.attrs['kTs'] = self.kTs
g.attrs['zaids'] = self.zaids
g.attrs['secondary_mode'] = np.string_(self.secondary_mode)
# Write thermal elastic scattering
if self.elastic_xs is not None:
elastic_group = g.create_group('elastic')
self.elastic_xs.to_hdf5(elastic_group, 'xs')
if self.elastic_mu_out is not None:
elastic_group.create_dataset('mu_out', data=self.elastic_mu_out)
for T in self.temperatures:
Tg = g.create_group(T)
# Write thermal elastic scattering
if self.elastic_xs:
elastic_group = Tg.create_group('elastic')
# Write thermal inelastic scattering
if self.inelastic_xs is not None:
inelastic_group = g.create_group('inelastic')
self.inelastic_xs.to_hdf5(inelastic_group, 'xs')
inelastic_group.attrs['secondary_mode'] = np.string_(self.secondary_mode)
if self.secondary_mode in ('equal', 'skewed'):
inelastic_group.create_dataset('energy_out', data=self.inelastic_e_out)
inelastic_group.create_dataset('mu_out', data=self.inelastic_mu_out)
elif self.secondary_mode == 'continuous':
self.inelastic_dist.to_hdf5(inelastic_group)
self.elastic_xs[T].to_hdf5(elastic_group, 'xs')
if self.elastic_mu_out:
elastic_group.create_dataset('mu_out',
data=self.elastic_mu_out[T])
# Write thermal inelastic scattering
if self.inelastic_xs:
inelastic_group = Tg.create_group('inelastic')
self.inelastic_xs[T].to_hdf5(inelastic_group, 'xs')
if self.secondary_mode in ('equal', 'skewed'):
inelastic_group.create_dataset('energy_out',
data=self.inelastic_e_out[T])
inelastic_group.create_dataset('mu_out',
data=self.inelastic_mu_out[T])
elif self.secondary_mode == 'continuous':
self.inelastic_dist[T].to_hdf5(inelastic_group)
f.close()
def add_temperature_from_ace(self, ace_or_filename, name=None):
"""Add data to the ThermalScattering object from an ACE file at a
different temperature.
Parameters
----------
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
GND-conforming name of the material, e.g. c_H_in_H2O. If none is
passed, the appropriate name is guessed based on the name of the ACE
table.
Returns
-------
openmc.data.ThermalScattering
Thermal scattering data
"""
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
ace = get_table(ace_or_filename)
# Get new name that is GND-consistent
ace_name, xs = ace.name.split('.')
if name is None:
if ace_name.lower() in _THERMAL_NAMES:
name = _THERMAL_NAMES[ace_name.lower()]
else:
# Make an educated guess? This actually works well for JEFF-3.2
# which stupidly uses names like lw00.32t, lw01.32t, etc. for
# different temperatures
matches = get_close_matches(
ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5)
if len(matches) > 0:
name = _THERMAL_NAMES[matches[0]]
else:
# OK, we give up. Just use the ACE name.
name = 'c_' + ace.name
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(ace.name, name))
# If this ACE data matches the data within self then get the data
if ace.temperature not in self.kTs:
if name == self.name:
# Add temperature and kTs
strT = "{0:.1f}K".format(kT_to_K(ace.temperature))
self.temperatures.append(strT)
self.kTs.append(ace.temperature)
# Incoherent inelastic scattering cross section
idx = ace.jxs[1]
n_energy = int(ace.xss[idx])
energy = ace.xss[idx + 1: idx + 1 + n_energy]
xs = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
self.inelastic_xs[strT] = Tabulated1D(energy, xs)
# Make sure secondary_mode is always equal. This should always
# be the case, but to reduce future debugging should something
# change, this will alert the developers to the issue.
if ace.nxs[7] == 0:
secondary_mode = 'equal'
elif ace.nxs[7] == 1:
secondary_mode = 'skewed'
elif ace.nxs[7] == 2:
secondary_mode = 'continuous'
if secondary_mode != self.secondary_mode:
raise ValueError('Secondary Modes are inconsistent.')
n_energy_out = ace.nxs[4]
if self.secondary_mode in ('equal', 'skewed'):
n_mu = ace.nxs[3]
idx = ace.jxs[3]
self.inelastic_e_out[strT] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2):
n_mu + 2]
self.inelastic_e_out[strT].shape = \
(n_energy, n_energy_out)
self.inelastic_mu_out[strT] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)]
self.inelastic_mu_out[strT].shape = \
(n_energy, n_energy_out, n_mu + 2)
self.inelastic_mu_out[strT] = \
self.inelastic_mu_out[strT][:, :, 1:]
else:
n_mu = ace.nxs[3] - 1
idx = ace.jxs[3]
locc = ace.xss[idx:idx + n_energy].astype(int)
n_energy_out = \
ace.xss[idx + n_energy:idx + 2 * n_energy].astype(int)
energy_out = []
mu_out = []
for i in range(n_energy):
idx = locc[i]
# Outgoing energy distribution for incoming energy i
e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True)
eout_i.c = c
# Outgoing angle distribution for each
# (incoming, outgoing) energy pair
mu_i = []
for j in range(n_energy_out[i]):
mu = ace.xss[idx + 4:idx + 4 + n_mu]
p_mu = 1. / n_mu * np.ones(n_mu)
mu_ij = Discrete(mu, p_mu)
mu_ij.c = np.cumsum(p_mu)
mu_i.append(mu_ij)
idx += 3 + n_mu
energy_out.append(eout_i)
mu_out.append(mu_i)
# Create correlated angle-energy distribution
breakpoints = [n_energy]
interpolation = [2]
energy = self.inelastic_xs[strT].x
self.inelastic_dist[strT] = CorrelatedAngleEnergy(
breakpoints, interpolation, energy, energy_out, mu_out)
# Incoherent/coherent elastic scattering cross section
idx = ace.jxs[4]
if idx != 0:
n_energy = int(ace.xss[idx])
energy = ace.xss[idx + 1: idx + 1 + n_energy]
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
if ace.nxs[5] == 4:
self.elastic_xs[strT] = CoherentElastic(energy, P)
else:
self.elastic_xs[strT] = Tabulated1D(energy, P)
# Angular distribution
n_mu = ace.nxs[6]
if n_mu != -1:
idx = ace.jxs[6]
self.elastic_mu_out[strT] = \
ace.xss[idx:idx + n_energy * n_mu]
self.elastic_mu_out[strT].shape = \
(n_energy, n_mu)
else:
raise ValueError('Data provided for an incorrect library')
else:
raise Warning('Temperature data set already within '
'IncidentNeutron object')
@classmethod
def from_hdf5(cls, group_or_filename):
@ -263,35 +437,46 @@ class ThermalScattering(EqualityMixin):
name = group.name[1:]
atomic_weight_ratio = group.attrs['atomic_weight_ratio']
temperature = group.attrs['temperature']
table = cls(name, atomic_weight_ratio, temperature)
kTs = group.attrs['kTs'].tolist()
temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs]
table = cls(name, atomic_weight_ratio, kTs)
table.zaids = group.attrs['zaids']
table.secondary_mode = group.attrs['secondary_mode'].decode()
# Read thermal elastic scattering
if 'elastic' in group:
elastic_group = group['elastic']
for T in temperatures:
Tgroup = group[T]
if 'elastic' in Tgroup:
elastic_group = Tgroup['elastic']
# Cross section
elastic_xs_type = elastic_group['xs'].attrs['type'].decode()
if elastic_xs_type == 'tab1':
table.elastic_xs = Tabulated1D.from_hdf5(elastic_group['xs'])
elif elastic_xs_type == 'bragg':
table.elastic_xs = CoherentElastic.from_hdf5(elastic_group['xs'])
# Cross section
elastic_xs_type = elastic_group['xs'].attrs['type'].decode()
if elastic_xs_type == 'Tabulated1D':
table.elastic_xs[T] = \
Tabulated1D.from_hdf5(elastic_group['xs'])
elif elastic_xs_type == 'bragg':
table.elastic_xs[T] = \
CoherentElastic.from_hdf5(elastic_group['xs'])
# Angular distribution
if 'mu_out' in elastic_group:
table.elastic_mu_out = elastic_group['mu_out'].value
# Angular distribution
if 'mu_out' in elastic_group:
table.elastic_mu_out[T] = \
elastic_group['mu_out'].value
# Read thermal inelastic scattering
if 'inelastic' in group:
inelastic_group = group['inelastic']
table.secondary_mode = inelastic_group.attrs['secondary_mode'].decode()
table.inelastic_xs = Tabulated1D.from_hdf5(inelastic_group['xs'])
if table.secondary_mode in ('equal', 'skewed'):
table.inelastic_e_out = inelastic_group['energy_out']
table.inelastic_mu_out = inelastic_group['mu_out']
elif table.secondary_mode == 'continuous':
table.inelastic_dist = AngleEnergy.from_hdf5(inelastic_group)
# Read thermal inelastic scattering
if 'inelastic' in Tgroup:
inelastic_group = Tgroup['inelastic']
table.inelastic_xs[T] = \
Tabulated1D.from_hdf5(inelastic_group['xs'])
if table.secondary_mode in ('equal', 'skewed'):
table.inelastic_e_out[T] = \
inelastic_group['energy_out']
table.inelastic_mu_out[T] = \
inelastic_group['mu_out']
elif table.secondary_mode == 'continuous':
table.inelastic_dist[T] = \
AngleEnergy.from_hdf5(inelastic_group)
return table
@ -301,7 +486,7 @@ class ThermalScattering(EqualityMixin):
Parameters
----------
ace : openmc.data.ace.Table or str
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
@ -324,7 +509,7 @@ class ThermalScattering(EqualityMixin):
ace_name, xs = ace.name.split('.')
if name is None:
if ace_name.lower() in _THERMAL_NAMES:
name = _THERMAL_NAMES[ace_name.lower()] + '.' + xs
name = _THERMAL_NAMES[ace_name.lower()]
else:
# Make an educated guess?? This actually works well for JEFF-3.2
# which stupidly uses names like lw00.32t, lw01.32t, etc. for
@ -332,21 +517,25 @@ class ThermalScattering(EqualityMixin):
matches = get_close_matches(
ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5)
if len(matches) > 0:
name = _THERMAL_NAMES[matches[0]] + '.' + xs
name = _THERMAL_NAMES[matches[0]]
else:
# OK, we give up. Just use the ACE name.
name = 'c_' + ace.name
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(ace.name, name))
table = cls(name, ace.atomic_weight_ratio, ace.temperature)
# Assign temperature to the running list
kTs = [ace.temperature]
temperatures = ["{0:.1f}K".format(kT_to_K(ace.temperature))]
table = cls(name, ace.atomic_weight_ratio, kTs)
# Incoherent inelastic scattering cross section
idx = ace.jxs[1]
n_energy = int(ace.xss[idx])
energy = ace.xss[idx+1 : idx+1+n_energy]
xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy]
table.inelastic_xs = Tabulated1D(energy, xs)
table.inelastic_xs[temperatures[0]] = Tabulated1D(energy, xs)
if ace.nxs[7] == 0:
table.secondary_mode = 'equal'
@ -359,34 +548,45 @@ class ThermalScattering(EqualityMixin):
if table.secondary_mode in ('equal', 'skewed'):
n_mu = ace.nxs[3]
idx = ace.jxs[3]
table.inelastic_e_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2):n_mu+2]
table.inelastic_e_out.shape = (n_energy, n_energy_out)
table.inelastic_e_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2):
n_mu + 2]
table.inelastic_e_out[temperatures[0]].shape = \
(n_energy, n_energy_out)
table.inelastic_mu_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2)]
table.inelastic_mu_out.shape = (n_energy, n_energy_out, n_mu+2)
table.inelastic_mu_out = table.inelastic_mu_out[:, :, 1:]
table.inelastic_mu_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)]
table.inelastic_mu_out[temperatures[0]].shape = \
(n_energy, n_energy_out, n_mu+2)
table.inelastic_mu_out[temperatures[0]] = \
table.inelastic_mu_out[temperatures[0]][:, :, 1:]
else:
n_mu = ace.nxs[3] - 1
idx = ace.jxs[3]
locc = ace.xss[idx:idx + n_energy].astype(int)
n_energy_out = ace.xss[idx + n_energy:idx + 2*n_energy].astype(int)
n_energy_out = \
ace.xss[idx + n_energy:idx + 2 * n_energy].astype(int)
energy_out = []
mu_out = []
for i in range(n_energy):
idx = locc[i]
# Outgoing energy distribution for incoming energy i
e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):n_mu + 3]
p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):n_mu + 3]
c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):n_mu + 3]
e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True)
eout_i.c = c
# Outgoing angle distribution for each (incoming, outgoing) energy pair
# Outgoing angle distribution for each
# (incoming, outgoing) energy pair
mu_i = []
for j in range(n_energy_out[i]):
mu = ace.xss[idx + 4:idx + 4 + n_mu]
p_mu = 1./n_mu*np.ones(n_mu)
p_mu = 1. / n_mu * np.ones(n_mu)
mu_ij = Discrete(mu, p_mu)
mu_ij.c = np.cumsum(p_mu)
mu_i.append(mu_ij)
@ -398,28 +598,30 @@ class ThermalScattering(EqualityMixin):
# Create correlated angle-energy distribution
breakpoints = [n_energy]
interpolation = [2]
energy = table.inelastic_xs.x
table.inelastic_dist = CorrelatedAngleEnergy(
energy = table.inelastic_xs[temperatures[0]].x
table.inelastic_dist[temperatures[0]] = CorrelatedAngleEnergy(
breakpoints, interpolation, energy, energy_out, mu_out)
# Incoherent/coherent elastic scattering cross section
idx = ace.jxs[4]
if idx != 0:
n_energy = int(ace.xss[idx])
energy = ace.xss[idx+1 : idx+1+n_energy]
P = ace.xss[idx+1+n_energy : idx+1+2*n_energy]
energy = ace.xss[idx + 1: idx + 1 + n_energy]
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
if ace.nxs[5] == 4:
table.elastic_xs = CoherentElastic(energy, P)
table.elastic_xs[temperatures[0]] = CoherentElastic(energy, P)
else:
table.elastic_xs = Tabulated1D(energy, P)
table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P)
# Angular distribution
n_mu = ace.nxs[6]
if n_mu != -1:
idx = ace.jxs[6]
table.elastic_mu_out = ace.xss[idx:idx + n_energy*n_mu]
table.elastic_mu_out.shape = (n_energy, n_mu)
table.elastic_mu_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_mu]
table.elastic_mu_out[temperatures[0]].shape = \
(n_energy, n_mu)
# Get relevant ZAIDs
pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int)

92
scripts/openmc-ace-to-hdf5
View file

@ -124,47 +124,79 @@ for filename in ace_libraries:
continue
lib = openmc.data.ace.Library(filename)
nuclides = {}
for table in lib.tables:
if table.name.endswith('c'):
name, xs = table.name.split('.')
if xs.endswith('c'):
# Continuous-energy neutron data
try:
if name not in nuclides:
neutron = openmc.data.IncidentNeutron.from_ace(
table, args.metastable)
except Exception as e:
print('Failed to convert {}: {}'.format(table.name, e))
continue
table, args.metastable)
# try:
# neutron = openmc.data.IncidentNeutron.from_ace(
# table, args.metastable)
# except Exception as e:
# print('Failed to convert {}: {}'.format(table.name, e))
# continue
# 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
# 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))
print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
neutron.name))
# Determine filename
outfile = os.path.join(args.destination,
neutron.name.replace('.', '_') + '.h5')
neutron.export_to_hdf5(outfile, 'w')
# Determine filename
outfile = os.path.join(args.destination,
neutron.name.replace('.', '_') + '.h5')
neutron.export_to_hdf5(outfile, 'w')
# Register with library
library.register_file(outfile)
# Register with library
library.register_file(outfile)
elif table.name.endswith('t'):
# Add nuclide to list
nuclides[name] = outfile
else:
# Then we only need to append the data
print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
neutron.name))
neutron = \
openmc.data.IncidentNeutron.from_hdf5(nuclides[name])
neutron.add_temperature_from_ace(table, args.metastable)
neutron.export_to_hdf5(outfile + '_1', 'w')
os.rename(outfile + '_1', outfile)
elif xs.endswith('t'):
# Thermal scattering data
thermal = openmc.data.ThermalScattering.from_ace(table)
print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
thermal.name))
if name not in nuclides:
thermal = openmc.data.ThermalScattering.from_ace(table)
print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
thermal.name))
# Determine filename
outfile = os.path.join(args.destination,
thermal.name.replace('.', '_') + '.h5')
thermal.export_to_hdf5(outfile, 'w')
# Determine filename
outfile = os.path.join(args.destination,
thermal.name.replace('.', '_') + '.h5')
thermal.export_to_hdf5(outfile, 'w')
# Register with library
library.register_file(outfile, 'thermal')
# Register with library
library.register_file(outfile, 'thermal')
# Add data to list
nuclides[name] = outfile
else:
# Then we only need to append the data
print('Converting {} (ACE) to {} (HDF5)'.format(table.name,
thermal.name))
# if table.name == 'poly.11t':
# import pdb; pdb.set_trace()
thermal = openmc.data.ThermalScattering.from_hdf5(nuclides[name])
thermal.add_temperature_from_ace(table)
thermal.export_to_hdf5(outfile + '_1', 'w')
os.rename(outfile + '_1', outfile)
# Write cross_sections.xml
libpath = os.path.join(args.destination, 'cross_sections.xml')