Merge branch 'develop' into more-updates

This commit is contained in:
Paul Romano 2019-03-28 14:11:46 -05:00
commit 61445ea60f
7 changed files with 341 additions and 147 deletions

View file

@ -369,7 +369,7 @@ class WindowedMultipole(EqualityMixin):
out.windows = group['windows'][()]
out.broaden_poly = group['broaden_poly'][()].astype(np.bool)
out.broaden_poly = group['broaden_poly'][...].astype(np.bool)
if out.broaden_poly.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('broaden_poly', 'windows'))

View file

@ -19,71 +19,62 @@ from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
from .function import Tabulated1D
_SUBSHELLS = ['K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2',
'O3', 'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2',
'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11',
'Q1', 'Q2', 'Q3']
# Helper function to map designator to subshell string or None
def _subshell(i):
if i == 0:
return None
else:
return _SUBSHELLS[i - 1]
# Electron subshell labels
_SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3',
'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', 'P3', 'P4',
'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11','Q1', 'Q2', 'Q3']
_REACTION_NAME = {
501: 'Total photon interaction',
502: 'Photon coherent scattering',
504: 'Photon incoherent scattering',
515: 'Pair production, electron field',
516: 'Total pair production',
517: 'Pair production, nuclear field',
522: 'Photoelectric absorption',
526: 'Electro-atomic scattering',
527: 'Electro-atomic bremsstrahlung',
528: 'Electro-atomic excitation',
534: 'K (1s1/2) subshell photoelectric',
535: 'L1 (2s1/2) subshell photoelectric',
536: 'L2 (2p1/2) subshell photoelectric',
537: 'L3 (2p3/2) subshell photoelectric',
538: 'M1 (3s1/2) subshell photoelectric',
539: 'M2 (3p1/2) subshell photoelectric',
540: 'M3 (3p3/2) subshell photoelectric',
541: 'M4 (3d3/2) subshell photoelectric',
542: 'M5 (3d5/2) subshell photoelectric',
543: 'N1 (4s1/2) subshell photoelectric',
544: 'N2 (4p1/2) subshell photoelectric',
545: 'N3 (4p3/2) subshell photoelectric',
546: 'N4 (4d3/2) subshell photoelectric',
547: 'N5 (4d5/2) subshell photoelectric',
548: 'N6 (4f5/2) subshell photoelectric',
549: 'N7 (4f7/2) subshell photoelectric',
550: 'O1 (5s1/2) subshell photoelectric',
551: 'O2 (5p1/2) subshell photoelectric',
552: 'O3 (5p3/2) subshell photoelectric',
553: 'O4 (5d3/2) subshell photoelectric',
554: 'O5 (5d5/2) subshell photoelectric',
555: 'O6 (5f5/2) subshell photoelectric',
556: 'O7 (5f7/2) subshell photoelectric',
557: 'O8 (5g7/2) subshell photoelectric',
558: 'O9 (5g9/2) subshell photoelectric',
559: 'P1 (6s1/2) subshell photoelectric',
560: 'P2 (6p1/2) subshell photoelectric',
561: 'P3 (6p3/2) subshell photoelectric',
562: 'P4 (6d3/2) subshell photoelectric',
563: 'P5 (6d5/2) subshell photoelectric',
564: 'P6 (6f5/2) subshell photoelectric',
565: 'P7 (6f7/2) subshell photoelectric',
566: 'P8 (6g7/2) subshell photoelectric',
567: 'P9 (6g9/2) subshell photoelectric',
568: 'P10 (6h9/2) subshell photoelectric',
569: 'P11 (6h11/2) subshell photoelectric',
570: 'Q1 (7s1/2) subshell photoelectric',
571: 'Q2 (7p1/2) subshell photoelectric',
572: 'Q3 (7p3/2) subshell photoelectric'
501: ('Total photon interaction', 'total'),
502: ('Photon coherent scattering', 'coherent'),
504: ('Photon incoherent scattering', 'incoherent'),
515: ('Pair production, electron field', 'pair_production_electron'),
516: ('Total pair production', 'pair_production_total'),
517: ('Pair production, nuclear field', 'pair_production_nuclear'),
522: ('Photoelectric absorption', 'photoelectric'),
526: ('Electro-atomic scattering', 'electro_atomic_scat'),
527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'),
528: ('Electro-atomic excitation', 'electro_atomic_excit'),
534: ('K (1s1/2) subshell photoelectric', 'K'),
535: ('L1 (2s1/2) subshell photoelectric', 'L1'),
536: ('L2 (2p1/2) subshell photoelectric', 'L2'),
537: ('L3 (2p3/2) subshell photoelectric', 'L3'),
538: ('M1 (3s1/2) subshell photoelectric', 'M1'),
539: ('M2 (3p1/2) subshell photoelectric', 'M2'),
540: ('M3 (3p3/2) subshell photoelectric', 'M3'),
541: ('M4 (3d3/2) subshell photoelectric', 'M4'),
542: ('M5 (3d5/2) subshell photoelectric', 'M5'),
543: ('N1 (4s1/2) subshell photoelectric', 'N1'),
544: ('N2 (4p1/2) subshell photoelectric', 'N2'),
545: ('N3 (4p3/2) subshell photoelectric', 'N3'),
546: ('N4 (4d3/2) subshell photoelectric', 'N4'),
547: ('N5 (4d5/2) subshell photoelectric', 'N5'),
548: ('N6 (4f5/2) subshell photoelectric', 'N6'),
549: ('N7 (4f7/2) subshell photoelectric', 'N7'),
550: ('O1 (5s1/2) subshell photoelectric', 'O1'),
551: ('O2 (5p1/2) subshell photoelectric', 'O2'),
552: ('O3 (5p3/2) subshell photoelectric', 'O3'),
553: ('O4 (5d3/2) subshell photoelectric', 'O4'),
554: ('O5 (5d5/2) subshell photoelectric', 'O5'),
555: ('O6 (5f5/2) subshell photoelectric', 'O6'),
556: ('O7 (5f7/2) subshell photoelectric', 'O7'),
557: ('O8 (5g7/2) subshell photoelectric', 'O8'),
558: ('O9 (5g9/2) subshell photoelectric', 'O9'),
559: ('P1 (6s1/2) subshell photoelectric', 'P1'),
560: ('P2 (6p1/2) subshell photoelectric', 'P2'),
561: ('P3 (6p3/2) subshell photoelectric', 'P3'),
562: ('P4 (6d3/2) subshell photoelectric', 'P4'),
563: ('P5 (6d5/2) subshell photoelectric', 'P5'),
564: ('P6 (6f5/2) subshell photoelectric', 'P6'),
565: ('P7 (6f7/2) subshell photoelectric', 'P7'),
566: ('P8 (6g7/2) subshell photoelectric', 'P8'),
567: ('P9 (6g9/2) subshell photoelectric', 'P9'),
568: ('P10 (6h9/2) subshell photoelectric', 'P10'),
569: ('P11 (6h11/2) subshell photoelectric', 'P11'),
570: ('Q1 (7s1/2) subshell photoelectric', 'Q1'),
571: ('Q2 (7p1/2) subshell photoelectric', 'Q2'),
572: ('Q3 (7p3/2) subshell photoelectric', 'Q3')
}
# Compton profiles are read from a pre-generated HDF5 file when they are first
@ -225,7 +216,7 @@ class AtomicRelaxation(EqualityMixin):
# Get shell designators
n = ace.nxs[7]
idx = ace.jxs[11]
shells = [_subshell(int(i)) for i in ace.xss[idx : idx+n]]
shells = [_SUBSHELLS[int(i)] for i in ace.xss[idx : idx+n]]
# Get number of electrons for each shell
idx = ace.jxs[12]
@ -245,8 +236,8 @@ class AtomicRelaxation(EqualityMixin):
if n_transitions > 0:
records = []
for j in range(n_transitions):
subj = _subshell(int(ace.xss[idx]))
subk = _subshell(int(ace.xss[idx + 1]))
subj = _SUBSHELLS[int(ace.xss[idx])]
subk = _SUBSHELLS[int(ace.xss[idx + 1])]
etr = ace.xss[idx + 2]*EV_PER_MEV
if j == 0:
ftr = ace.xss[idx + 3]
@ -301,7 +292,7 @@ class AtomicRelaxation(EqualityMixin):
# Read data for each subshell
for i in range(n_subshells):
params, list_items = get_list_record(file_obj)
subi = _subshell(int(params[0]))
subi = _SUBSHELLS[int(params[0])]
n_transitions = int(params[5])
binding_energy[subi] = list_items[0]
num_electrons[subi] = list_items[1]
@ -310,8 +301,8 @@ class AtomicRelaxation(EqualityMixin):
# Read transition data
records = []
for j in range(n_transitions):
subj = _subshell(int(list_items[6*(j+1)]))
subk = _subshell(int(list_items[6*(j+1) + 1]))
subj = _SUBSHELLS[int(list_items[6*(j+1)])]
subk = _SUBSHELLS[int(list_items[6*(j+1) + 1])]
etr = list_items[6*(j+1) + 2]
ftr = list_items[6*(j+1) + 3]
records.append((subj, subk, etr, ftr))
@ -323,8 +314,70 @@ class AtomicRelaxation(EqualityMixin):
# Return instance of class
return cls(binding_energy, num_electrons, transitions)
def to_hdf5(self, group):
raise NotImplementedError
@classmethod
def from_hdf5(cls, group):
"""Generate atomic relaxation data from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
Returns
-------
openmc.data.AtomicRelaxation
Atomic relaxation data
"""
# Create data dictionaries
binding_energy = {}
num_electrons = {}
transitions = {}
designators = [s.decode() for s in group.attrs['designators']]
columns = ['secondary', 'tertiary', 'energy (eV)', 'probability']
for shell in designators:
# Shell group
sub_group = group[shell]
# Read subshell binding energy and number of electrons
if 'binding_energy' in sub_group.attrs:
binding_energy[shell] = sub_group.attrs['binding_energy']
if 'num_electrons' in sub_group.attrs:
num_electrons[shell] = sub_group.attrs['num_electrons']
# Read transition data
if 'transitions' in sub_group:
df = pd.DataFrame(sub_group['transitions'][()],
columns=columns)
# Replace float indexes back to subshell strings
df[columns[:2]] = df[columns[:2]].replace(
np.arange(float(len(_SUBSHELLS))), _SUBSHELLS)
transitions[shell] = df
return cls(binding_energy, num_electrons, transitions)
def to_hdf5(self, group, shell):
"""Write atomic relaxation data to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
shell : str
The subshell to write data for
"""
# Write subshell binding energy and number of electrons
group.attrs['binding_energy'] = self.binding_energy[shell]
group.attrs['num_electrons'] = self.num_electrons[shell]
# Write transition data with replacements
if shell in self.transitions:
df = self.transitions[shell].replace(
_SUBSHELLS, range(len(_SUBSHELLS)))
group.create_dataset('transitions', data=df.values.astype(float))
class IncidentPhoton(EqualityMixin):
@ -509,7 +562,7 @@ class IncidentPhoton(EqualityMixin):
idx += n_energy
# Copy binding energy
shell = _subshell(d)
shell = _SUBSHELLS[d]
e = data.atomic_relaxation.binding_energy[shell]
rx.subshell_binding_energy = e
@ -580,6 +633,89 @@ class IncidentPhoton(EqualityMixin):
return data
@classmethod
def from_hdf5(cls, group_or_filename):
"""Generate photon reaction from an HDF5 group
Parameters
----------
group_or_filename : h5py.Group or str
HDF5 group containing interaction data. If given as a string, it is
assumed to be the filename for the HDF5 file, and the first group is
used to read from.
Returns
-------
openmc.data.IncidentPhoton
Photon interaction data
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(str(group_or_filename), 'r')
# Make sure version matches
if 'version' in h5file.attrs:
major, minor = h5file.attrs['version']
# For now all versions of HDF5 data can be read
else:
raise IOError(
'HDF5 data does not indicate a version. Your installation '
'of the OpenMC Python API expects version {}.x data.'
.format(HDF5_VERSION_MAJOR))
group = list(h5file.values())[0]
Z = group.attrs['Z']
data = cls(Z)
# Read energy grid
energy = group['energy'][()]
# Read cross section data
for mt, (name, key) in _REACTION_NAME.items():
if key in group:
rgroup = group[key]
elif key in group['subshells']:
rgroup = group['subshells'][key]
else:
continue
data.reactions[mt] = PhotonReaction.from_hdf5(rgroup, mt, energy)
# Check for necessary reactions
for mt in (502, 504, 522):
assert mt in data, "Reaction {} not found".format(mt)
# Read atomic relaxation
data.atomic_relaxation = AtomicRelaxation.from_hdf5(group['subshells'])
# Read Compton profiles
if 'compton_profiles' in group:
rgroup = group['compton_profiles']
profile = data.compton_profiles
profile['num_electrons'] = rgroup['num_electrons'][()]
profile['binding_energy'] = rgroup['binding_energy'][()]
# Get electron momentum values
pz = rgroup['pz'][()]
J = rgroup['J'][()]
if pz.size != J.shape[1]:
raise ValueError("'J' array shape is not consistent with the "
"'pz' array shape")
profile['J'] = [Tabulated1D(pz, Jk) for Jk in J]
# Read bremsstrahlung
if 'bremsstrahlung' in group:
rgroup = group['bremsstrahlung']
data.bremsstrahlung['I'] = rgroup.attrs['I']
for key in ('dcs', 'electron_energy', 'ionization_energy',
'num_electrons', 'photon_energy'):
data.bremsstrahlung[key] = rgroup[key][()]
return data
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export incident photon data to an HDF5 file.
@ -590,6 +726,9 @@ class IncidentPhoton(EqualityMixin):
mode : {'r', r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
that are less backwards compatible but have performance benefits.
"""
# Open file and write version
@ -607,77 +746,25 @@ class IncidentPhoton(EqualityMixin):
union_grid = np.union1d(union_grid, rx.xs.x)
group.create_dataset('energy', data=union_grid)
# Write coherent scattering cross section
rx = self.reactions[502]
coh_group = group.create_group('coherent')
coh_group.create_dataset('xs', data=rx.xs(union_grid))
if rx.scattering_factor is not None:
# Create integrated form factor
ff = deepcopy(rx.scattering_factor)
ff.x *= ff.x
ff.y *= ff.y/Z**2
int_ff = Tabulated1D(ff.x, ff.integral())
int_ff.to_hdf5(coh_group, 'integrated_scattering_factor')
if rx.anomalous_real is not None:
rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real')
if rx.anomalous_imag is not None:
rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag')
# Write incoherent scattering cross section
rx = self[504]
incoh_group = group.create_group('incoherent')
incoh_group.create_dataset('xs', data=rx.xs(union_grid))
if rx.scattering_factor is not None:
rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor')
# Write electron-field pair production cross section
if 515 in self:
pair_group = group.create_group('pair_production_electron')
pair_group.create_dataset('xs', data=self[515].xs(union_grid))
# Write nuclear-field pair production cross section
if 517 in self:
pair_group = group.create_group('pair_production_nuclear')
pair_group.create_dataset('xs', data=self[517].xs(union_grid))
# Write photoelectric cross section
photoelec_group = group.create_group('photoelectric')
photoelec_group.create_dataset('xs', data=self[522].xs(union_grid))
# Write photoionization cross sections
# Write cross sections
shell_group = group.create_group('subshells')
designators = []
for mt, rx in self.reactions.items():
if mt >= 534 and mt <= 572:
# Get name of subshell
shell = _SUBSHELLS[mt - 534]
designators.append(shell)
sub_group = shell_group.create_group(shell)
name, key = _REACTION_NAME[mt]
if mt in [502, 504, 515, 517, 522]:
sub_group = group.create_group(key)
elif mt >= 534 and mt <= 572:
# Subshell
designators.append(key)
sub_group = shell_group.create_group(key)
if self.atomic_relaxation is not None:
relax = self.atomic_relaxation
# Write subshell binding energy and number of electrons
sub_group.attrs['binding_energy'] = relax.binding_energy[shell]
sub_group.attrs['num_electrons'] = relax.num_electrons[shell]
# Write atomic relaxation
if key in self.atomic_relaxation.subshells:
self.atomic_relaxation.to_hdf5(sub_group, key)
else:
continue
# Write transition data with replacements
if shell in relax.transitions:
shell_values = _SUBSHELLS.copy()
shell_values.insert(0, None)
df = relax.transitions[shell].replace(
shell_values, range(len(shell_values)))
sub_group.create_dataset(
'transitions', data=df.values.astype(float))
# Determine threshold
threshold = rx.xs.x[0]
idx = np.searchsorted(union_grid, threshold, side='right') - 1
# Interpolate cross section onto union grid and write
photoionization = rx.xs(union_grid[idx:])
sub_group.create_dataset('xs', data=photoionization)
assert len(union_grid) == len(photoionization) + idx
sub_group['xs'].attrs['threshold_idx'] = idx
rx.to_hdf5(sub_group, union_grid, Z)
shell_group.attrs['designators'] = np.array(designators, dtype='S')
@ -805,7 +892,7 @@ class PhotonReaction(EqualityMixin):
def __repr__(self):
if self.mt in _REACTION_NAME:
return "<Photon Reaction: MT={} {}>".format(
self.mt, _REACTION_NAME[self.mt])
self.mt, _REACTION_NAME[self.mt][0])
else:
return "<Photon Reaction: MT={}>".format(self.mt)
@ -982,3 +1069,93 @@ class PhotonReaction(EqualityMixin):
params, rx.anomalous_imag = get_tab1_record(file_obj)
return rx
@classmethod
def from_hdf5(cls, group, mt, energy):
"""Generate photon reaction from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
mt : int
The MT value of the reaction to get data for
energy : Iterable of float
arrays of energies at which cross sections are tabulated at
Returns
-------
openmc.data.PhotonReaction
Photon reaction data
"""
# Create instance
rx = cls(mt)
# Cross sections
xs = group['xs'][()]
# Replace zero elements to small non-zero to enable log-log
xs[xs == 0.0] = np.exp(-500.0)
# Threshold
threshold_idx = 0
if 'threshold_idx' in group['xs'].attrs:
threshold_idx = group['xs'].attrs['threshold_idx']
# Store cross section
rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5])
# Check for anomalous scattering factor
if 'anomalous_real' in group:
rx.anomalous_real = Tabulated1D.from_hdf5(group['anomalous_real'])
if 'anomalous_imag' in group:
rx.anomalous_imag = Tabulated1D.from_hdf5(group['anomalous_imag'])
# Check for factors / scattering functions
if 'scattering_factor' in group:
rx.scattering_factor = Tabulated1D.from_hdf5(group['scattering_factor'])
return rx
def to_hdf5(self, group, energy, Z):
"""Write photon reaction to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
energy : Iterable of float
arrays of energies at which cross sections are tabulated at
Z : int
atomic number
"""
# Write cross sections
if self.mt >= 534 and self.mt <= 572:
# Determine threshold
threshold = self.xs.x[0]
idx = np.searchsorted(energy, threshold, side='right') - 1
# Interpolate cross section onto union grid and write
photoionization = self.xs(energy[idx:])
group.create_dataset('xs', data=photoionization)
assert len(energy) == len(photoionization) + idx
group['xs'].attrs['threshold_idx'] = idx
else:
group.create_dataset('xs', data=self.xs(energy))
# Write scattering factor
if self.scattering_factor is not None:
if self.mt == 502:
# Create integrated form factor
ff = deepcopy(self.scattering_factor)
ff.x *= ff.x
ff.y *= ff.y/Z**2
int_ff = Tabulated1D(ff.x, ff.integral())
int_ff.to_hdf5(group, 'integrated_scattering_factor')
self.scattering_factor.to_hdf5(group, 'scattering_factor')
if self.anomalous_real is not None:
self.anomalous_real.to_hdf5(group, 'anomalous_real')
if self.anomalous_imag is not None:
self.anomalous_imag.to_hdf5(group, 'anomalous_imag')

View file

@ -20,7 +20,7 @@ class ResultsList(list):
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
# Get number of results stored
n = fh["number"].value.shape[0]
n = fh["number"][...].shape[0]
for i in range(n):
self.append(Results.from_hdf5(fh, i))

View file

@ -104,11 +104,11 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
lattice_type = group['type'][()].decode()
if lattice_type == 'rectangular':
dimension = group['dimension'][()]
lower_left = group['lower_left'][()]
pitch = group['pitch'][()]
dimension = group['dimension'][...]
lower_left = group['lower_left'][...]
pitch = group['pitch'][...]
outer = group['outer'][()]
universe_ids = group['universes'][()]
universe_ids = group['universes'][...]
# Create the Lattice
lattice = openmc.RectLattice(lattice_id, name)

View file

@ -2183,7 +2183,7 @@ class XSdata(object):
kTs_group = group['kTs']
float_temperatures = []
for temperature in temperatures:
kT = kTs_group[temperature].value
kT = kTs_group[temperature][()]
float_temperatures.append(kT / openmc.data.K_BOLTZMANN)
attrs = group.attrs.keys()

View file

@ -272,7 +272,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta):
name = group['name'][()].decode() if 'name' in group else ''
surf_type = group['type'][()].decode()
bc = group['boundary_type'][()].decode()
coeffs = group['coefficients'][()]
coeffs = group['coefficients'][...]
# Create the Surface based on its type
if surf_type == 'x-plane':

View file

@ -129,3 +129,20 @@ def test_export_to_hdf5(tmpdir, element):
filename = str(tmpdir.join('tmp.h5'))
element.export_to_hdf5(filename)
assert os.path.exists(filename)
# Read in data from hdf5
element2 = openmc.data.IncidentPhoton.from_hdf5(filename)
# Check for some cross section and datasets of element and element2
energy = np.logspace(np.log10(1.0), np.log10(1.0e10), num=100)
for mt in (502, 504, 515, 517, 522, 541, 570):
xs = element[mt].xs(energy)
xs2 = element2[mt].xs(energy)
assert np.allclose(xs, xs2)
assert element[502].scattering_factor == element2[502].scattering_factor
assert element.atomic_relaxation.transitions['O3'].equals(
element2.atomic_relaxation.transitions['O3'])
assert (element.compton_profiles['binding_energy'] ==
element2.compton_profiles['binding_energy']).all()
assert (element.bremsstrahlung['electron_energy'] ==
element2.bremsstrahlung['electron_energy']).all()
# Export to hdf5 again
element2.export_to_hdf5(filename, 'w')