mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Added HDF5 support in statepoint.py. Closes #152 on github.
This commit is contained in:
Paul Romano
parent
5ccc7844d7
commit
0d874cc3b3
3 changed files with 146 additions and 97 deletions
|
|
@ -499,7 +499,7 @@ contains
|
|||
! Write information for meshes
|
||||
MESH_LOOP: do i = 1, n_meshes
|
||||
! Create temporary group for each mesh
|
||||
call h5gcreate_f(tallies_group, "mesh " // to_str(meshes(i) % id), &
|
||||
call h5gcreate_f(tallies_group, "mesh" // to_str(i), &
|
||||
temp_group, hdf5_err)
|
||||
|
||||
! Write type and number of dimensions
|
||||
|
|
@ -530,7 +530,7 @@ contains
|
|||
t => tallies(i)
|
||||
|
||||
! Create group for this tally
|
||||
call h5gcreate_f(tallies_group, "tally " // to_str(t % id), &
|
||||
call h5gcreate_f(tallies_group, "tally" // to_str(i), &
|
||||
temp_group, hdf5_err)
|
||||
|
||||
! Write size of each tally
|
||||
|
|
@ -544,7 +544,7 @@ contains
|
|||
|
||||
FILTER_LOOP: do j = 1, t % n_filters
|
||||
! Create filter group
|
||||
call h5gcreate_f(temp_group, "filter " // to_str(j), filter_group, &
|
||||
call h5gcreate_f(temp_group, "filter" // to_str(j), filter_group, &
|
||||
hdf5_err)
|
||||
|
||||
! Write type of filter
|
||||
|
|
@ -904,7 +904,7 @@ contains
|
|||
! Write information for meshes
|
||||
MESH_LOOP: do i = 1, n_meshes
|
||||
! Create temporary group for each mesh
|
||||
call h5gcreate_f(tallies_group, "mesh " // to_str(meshes(i) % id), &
|
||||
call h5gcreate_f(tallies_group, "mesh" // to_str(i), &
|
||||
temp_group, hdf5_err)
|
||||
|
||||
! Write id, type, and number of dimensions
|
||||
|
|
@ -936,7 +936,7 @@ contains
|
|||
t => tallies(i)
|
||||
|
||||
! Create group for this tally
|
||||
call h5gcreate_f(tallies_group, "tally " // to_str(t % id), &
|
||||
call h5gcreate_f(tallies_group, "tally" // to_str(i), &
|
||||
temp_group, hdf5_err)
|
||||
|
||||
! Write id
|
||||
|
|
@ -957,7 +957,7 @@ contains
|
|||
|
||||
FILTER_LOOP: do j = 1, t % n_filters
|
||||
! Create filter group
|
||||
call h5gcreate_f(temp_group, "filter " // to_str(j), filter_group, &
|
||||
call h5gcreate_f(temp_group, "filter" // to_str(j), filter_group, &
|
||||
hdf5_err)
|
||||
|
||||
! Write type of filter
|
||||
|
|
@ -1077,7 +1077,7 @@ contains
|
|||
t => tallies(i)
|
||||
|
||||
! Open group for the i-th tally
|
||||
call h5gopen_f(tallies_group, "tally " // to_str(t % id), &
|
||||
call h5gopen_f(tallies_group, "tally" // to_str(i), &
|
||||
temp_group, hdf5_err)
|
||||
|
||||
! Write sum and sum_sq for each bin
|
||||
|
|
@ -1218,8 +1218,8 @@ contains
|
|||
|
||||
TALLIES_LOOP: do i = 1, n_tallies
|
||||
! Open tally group
|
||||
call h5gopen_f(hdf5_state_point, "tallies/tally " // &
|
||||
to_str(tallies(i) % id), tally_group, hdf5_err)
|
||||
call h5gopen_f(hdf5_state_point, "tallies/tally" // &
|
||||
to_str(i), tally_group, hdf5_err)
|
||||
|
||||
! Read number of realizations
|
||||
call hdf5_read_integer(tally_group, "n_realizations", &
|
||||
|
|
|
|||
|
|
@ -90,31 +90,6 @@ score_types.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
|
|||
score_types.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
|
||||
score_types.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})
|
||||
|
||||
|
||||
class BinaryFile(object):
|
||||
def __init__(self, filename):
|
||||
self._f = open(filename, 'rb')
|
||||
|
||||
def _get_data(self, n, typeCode, size):
|
||||
return list(struct.unpack('={0}{1}'.format(n,typeCode),
|
||||
self._f.read(n*size)))
|
||||
|
||||
def _get_int(self, n=1):
|
||||
return self._get_data(n, 'i', 4)
|
||||
|
||||
def _get_long(self, n=1):
|
||||
return self._get_data(n, 'q', 8)
|
||||
|
||||
def _get_float(self, n=1):
|
||||
return self._get_data(n, 'f', 4)
|
||||
|
||||
def _get_double(self, n=1):
|
||||
return self._get_data(n, 'd', 8)
|
||||
|
||||
def _get_string(self, length, n=1):
|
||||
data = self._get_data(length*n, 's', 1)[0]
|
||||
return [data[i*length:(i+1)*length] for i in range(n)]
|
||||
|
||||
class Mesh(object):
|
||||
def __init__(self):
|
||||
pass
|
||||
|
|
@ -149,9 +124,15 @@ class SourceSite(object):
|
|||
return "<SourceSite: xyz={0} at E={1}>".format(self.xyz, self.E)
|
||||
|
||||
|
||||
class StatePoint(BinaryFile):
|
||||
class StatePoint(object):
|
||||
def __init__(self, filename):
|
||||
super(StatePoint, self).__init__(filename)
|
||||
if filename.endswith('.h5'):
|
||||
import h5py
|
||||
self._f = h5py.File(filename, 'r')
|
||||
self._hdf5 = True
|
||||
else:
|
||||
self._f = open(filename, 'rb')
|
||||
self._hdf5 = False
|
||||
|
||||
# Set flags for what data was read
|
||||
self._metadata = False
|
||||
|
|
@ -168,106 +149,125 @@ class StatePoint(BinaryFile):
|
|||
|
||||
def _read_metadata(self):
|
||||
# Read statepoint revision
|
||||
self.revision = self._get_int()[0]
|
||||
self.revision = self._get_int(path='revision_statepoint')[0]
|
||||
|
||||
# Read OpenMC version
|
||||
self.version = self._get_int(3)
|
||||
if self._hdf5:
|
||||
self.version = [self._get_int(path='version_major')[0],
|
||||
self._get_int(path='version_minor')[0],
|
||||
self._get_int(path='version_release')[0]]
|
||||
else:
|
||||
self.version = self._get_int(3)
|
||||
|
||||
# Read date and time
|
||||
self.date_and_time = self._get_string(19)[0]
|
||||
self.date_and_time = self._get_string(19, path='date_and_time')
|
||||
|
||||
# Read path
|
||||
self.path = self._get_string(255)[0].strip()
|
||||
self.path = self._get_string(255, path='path').strip()
|
||||
|
||||
# Read random number seed
|
||||
self.seed = self._get_long()[0]
|
||||
self.seed = self._get_long(path='seed')[0]
|
||||
|
||||
# Read run information
|
||||
self.run_mode = self._get_int()[0]
|
||||
self.n_particles = self._get_long()[0]
|
||||
self.n_batches = self._get_int()[0]
|
||||
self.run_mode = self._get_int(path='run_mode')[0]
|
||||
self.n_particles = self._get_long(path='n_particles')[0]
|
||||
self.n_batches = self._get_int(path='n_batches')[0]
|
||||
|
||||
# Read current batch
|
||||
self.current_batch = self._get_int()[0]
|
||||
self.current_batch = self._get_int(path='current_batch')[0]
|
||||
|
||||
# Read criticality information
|
||||
if self.run_mode == 2:
|
||||
self.n_inactive, self.gen_per_batch = self._get_int(2)
|
||||
self.k_batch = self._get_double(self.current_batch)
|
||||
self.entropy = self._get_double(self.current_batch*self.gen_per_batch)
|
||||
self.k_col_abs = self._get_double()[0]
|
||||
self.k_col_tra = self._get_double()[0]
|
||||
self.k_abs_tra = self._get_double()[0]
|
||||
self.k_combined = self._get_double(2)
|
||||
self.n_inactive = self._get_int(path='n_inactive')[0]
|
||||
self.gen_per_batch = self._get_int(path='gen_per_batch')[0]
|
||||
self.k_batch = self._get_double(self.current_batch, path='k_batch')
|
||||
self.entropy = self._get_double(
|
||||
self.current_batch*self.gen_per_batch, path='entropy')
|
||||
self.k_col_abs = self._get_double(path='k_col_abs')[0]
|
||||
self.k_col_tra = self._get_double(path='k_col_tra')[0]
|
||||
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
|
||||
self.k_combined = self._get_double(2, path='k_combined')
|
||||
|
||||
# Read number of meshes
|
||||
n_meshes = self._get_int()[0]
|
||||
n_meshes = self._get_int(path='tallies/n_meshes')[0]
|
||||
|
||||
# Read meshes
|
||||
for i in range(n_meshes):
|
||||
m = Mesh()
|
||||
self.meshes.append(m)
|
||||
|
||||
base = 'tallies/mesh' + str(i+1) + '/'
|
||||
|
||||
# Read id, mesh type, and number of dimensions
|
||||
m.id, m.type, n = self._get_int(3)
|
||||
m.id = self._get_int(path=base+'id')[0]
|
||||
m.type = self._get_int(path=base+'type')[0]
|
||||
n = self._get_int(path=base+'n_dimension')[0]
|
||||
|
||||
# Read mesh size, lower-left coordinates, upper-right coordinates,
|
||||
# and width of each mesh cell
|
||||
m.dimension = self._get_int(n)
|
||||
m.lower_left = self._get_double(n)
|
||||
m.upper_right = self._get_double(n)
|
||||
m.width = self._get_double(n)
|
||||
m.dimension = self._get_int(n, path=base+'dimension')
|
||||
m.lower_left = self._get_double(n, path=base+'lower_left')
|
||||
m.upper_right = self._get_double(n, path=base+'upper_right')
|
||||
m.width = self._get_double(n, path=base+'width')
|
||||
|
||||
# Read number of tallies
|
||||
n_tallies = self._get_int()[0]
|
||||
n_tallies = self._get_int(path='tallies/n_tallies')[0]
|
||||
|
||||
for i in range(n_tallies):
|
||||
# Create Tally object and add to list of tallies
|
||||
t = Tally()
|
||||
self.tallies.append(t)
|
||||
|
||||
base = 'tallies/tally' + str(i+1) + '/'
|
||||
|
||||
# Read id and number of realizations
|
||||
t.id, t.n_realizations = self._get_int(2)
|
||||
t.id = self._get_int(path=base+'id')[0]
|
||||
t.n_realizations = self._get_int(path=base+'n_realizations')[0]
|
||||
|
||||
# Read sizes of tallies
|
||||
t.total_score_bins = self._get_int()[0]
|
||||
t.total_filter_bins = self._get_int()[0]
|
||||
t.total_score_bins = self._get_int(path=base+'total_score_bins')[0]
|
||||
t.total_filter_bins = self._get_int(path=base+'total_filter_bins')[0]
|
||||
|
||||
# Read number of filters
|
||||
n_filters = self._get_int()[0]
|
||||
n_filters = self._get_int(path=base+'n_filters')[0]
|
||||
|
||||
for j in range(n_filters):
|
||||
# Create Filter object
|
||||
f = Filter()
|
||||
|
||||
base = 'tallies/tally{0}/filter{1}/'.format(i+1, j+1)
|
||||
|
||||
# Get type of filter
|
||||
f.type = filter_types[self._get_int()[0]]
|
||||
f.type = filter_types[self._get_int(path=base+'type')[0]]
|
||||
|
||||
# Add to filter dictionary
|
||||
t.filters[f.type] = f
|
||||
|
||||
# Determine how many bins are in this filter
|
||||
f.length = self._get_int()[0]
|
||||
f.length = self._get_int(path=base+'n_bins')[0]
|
||||
assert f.length > 0
|
||||
if f.type == 'energyin' or f.type == 'energyout':
|
||||
f.bins = self._get_double(f.length + 1)
|
||||
f.bins = self._get_double(f.length + 1, path=base+'bins')
|
||||
elif f.type == 'mesh':
|
||||
f.bins = self._get_int()
|
||||
f.bins = self._get_int(path=base+'bins')
|
||||
else:
|
||||
f.bins = self._get_int(f.length)
|
||||
f.bins = self._get_int(f.length, path=base+'bins')
|
||||
|
||||
base = 'tallies/tally' + str(i+1) + '/'
|
||||
|
||||
# Read nuclide bins
|
||||
n_nuclides = self._get_int()[0]
|
||||
n_nuclides = self._get_int(path=base+'n_nuclide_bins')[0]
|
||||
t.n_nuclides = n_nuclides
|
||||
t.nuclides = self._get_int(n_nuclides)
|
||||
t.nuclides = self._get_int(n_nuclides, path=base+'nuclide_bins')
|
||||
|
||||
# Read score bins and scattering order
|
||||
t.n_scores = self._get_int()[0]
|
||||
t.scores = [score_types[j] for j in self._get_int(t.n_scores)]
|
||||
t.scatt_order = self._get_int(t.n_scores)
|
||||
t.n_scores = self._get_int(path=base+'n_score_bins')[0]
|
||||
t.scores = [score_types[j] for j in self._get_int(
|
||||
t.n_scores, path=base+'score_bins')]
|
||||
t.scatt_order = self._get_int(t.n_scores, path=base+'scatt_order')
|
||||
|
||||
# Read number of user score bins
|
||||
t.n_user_scores = self._get_int()[0]
|
||||
t.n_user_scores = self._get_int(path=base+'n_user_score_bins')[0]
|
||||
|
||||
# Set up stride
|
||||
stride = 1
|
||||
|
|
@ -284,23 +284,31 @@ class StatePoint(BinaryFile):
|
|||
self._read_metadata()
|
||||
|
||||
# Number of realizations for global tallies
|
||||
self.n_realizations = self._get_int()[0]
|
||||
self.n_realizations = self._get_int(path='n_realizations')[0]
|
||||
|
||||
# Read global tallies
|
||||
n_global_tallies = self._get_int()[0]
|
||||
self.global_tallies = np.array(self._get_double(2*n_global_tallies))
|
||||
self.global_tallies.shape = (n_global_tallies, 2)
|
||||
n_global_tallies = self._get_int(path='n_global_tallies')[0]
|
||||
if self._hdf5:
|
||||
data = self._f['global_tallies'].value
|
||||
self.global_tallies = np.column_stack((data['sum'], data['sum_sq']))
|
||||
else:
|
||||
self.global_tallies = np.array(self._get_double(2*n_global_tallies))
|
||||
self.global_tallies.shape = (n_global_tallies, 2)
|
||||
|
||||
# Flag indicating if tallies are present
|
||||
tallies_present = self._get_int()[0]
|
||||
tallies_present = self._get_int(path='tallies/tallies_present')[0]
|
||||
|
||||
# Read tally results
|
||||
if tallies_present:
|
||||
|
||||
for t in self.tallies:
|
||||
for i, t in enumerate(self.tallies):
|
||||
n = t.total_score_bins * t.total_filter_bins
|
||||
t.results = np.array(self._get_double(2*n))
|
||||
t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
|
||||
if self._hdf5:
|
||||
path = 'tallies/tally{0}/results'.format(i+1)
|
||||
data = self._f[path].value
|
||||
t.results = np.dstack((data['sum'], data['sum_sq']))
|
||||
else:
|
||||
t.results = np.array(self._get_double(2*n))
|
||||
t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
|
||||
|
||||
# Indicate that tally results have been read
|
||||
self._results = True
|
||||
|
|
@ -310,15 +318,22 @@ class StatePoint(BinaryFile):
|
|||
if not self._results:
|
||||
self.read_results()
|
||||
|
||||
# For HDF5 state points, copy entire bank
|
||||
if self._hdf5:
|
||||
source_sites = self._f['source_bank'].value
|
||||
|
||||
for i in range(self.n_particles):
|
||||
s = SourceSite()
|
||||
self.source.append(s)
|
||||
|
||||
# Read position, angle, and energy
|
||||
s.weight = self._get_double()[0]
|
||||
s.xyz = self._get_double(3)
|
||||
s.uvw = self._get_double(3)
|
||||
s.E = self._get_double()[0]
|
||||
if self._hdf5:
|
||||
s.weight, s.xyz, s.uvw, s.E = source_sites[i]
|
||||
else:
|
||||
s.weight = self._get_double()[0]
|
||||
s.xyz = self._get_double(3)
|
||||
s.uvw = self._get_double(3)
|
||||
s.E = self._get_double()[0]
|
||||
|
||||
def generate_ci(self, confidence=0.95):
|
||||
"""Calculates confidence intervals for each tally bin."""
|
||||
|
|
@ -423,3 +438,37 @@ class StatePoint(BinaryFile):
|
|||
# sum of squares, or it could be mean and stdev if self.generate_stdev()
|
||||
# has been called already.
|
||||
return t.results[filter_index, score_index]
|
||||
|
||||
def _get_data(self, n, typeCode, size):
|
||||
return list(struct.unpack('={0}{1}'.format(n,typeCode),
|
||||
self._f.read(n*size)))
|
||||
|
||||
def _get_int(self, n=1, path=None):
|
||||
if self._hdf5:
|
||||
return self._f[path].value
|
||||
else:
|
||||
return self._get_data(n, 'i', 4)
|
||||
|
||||
def _get_long(self, n=1, path=None):
|
||||
if self._hdf5:
|
||||
return self._f[path].value
|
||||
else:
|
||||
return self._get_data(n, 'q', 8)
|
||||
|
||||
def _get_float(self, n=1, path=None):
|
||||
if self._hdf5:
|
||||
return self._f[path].value
|
||||
else:
|
||||
return self._get_data(n, 'f', 4)
|
||||
|
||||
def _get_double(self, n=1, path=None):
|
||||
if self._hdf5:
|
||||
return self._f[path].value
|
||||
else:
|
||||
return self._get_data(n, 'd', 8)
|
||||
|
||||
def _get_string(self, n=1, path=None):
|
||||
if self._hdf5:
|
||||
return self._f[path].value
|
||||
else:
|
||||
return self._get_data(n, 's', 1)[0]
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
|
||||
import sys
|
||||
|
||||
from numpy.testing import assert_allclose
|
||||
from numpy.testing import assert_allclose, assert_equal
|
||||
|
||||
from statepoint import StatePoint
|
||||
|
||||
|
|
@ -47,10 +47,10 @@ assert_allclose(sp1.global_tallies, sp2.global_tallies)
|
|||
assert len(sp1.meshes) == len(sp2.meshes)
|
||||
for m1, m2 in zip(sp1.meshes, sp2.meshes):
|
||||
assert m1.type == m2.type
|
||||
assert m1.dimension == m2.dimension
|
||||
assert m1.lower_left == m2.lower_left
|
||||
assert m1.upper_right == m2.upper_right
|
||||
assert m1.width == m2.width
|
||||
assert_equal(m1.dimension, m2.dimension)
|
||||
assert_equal(m1.lower_left, m2.lower_left)
|
||||
assert_equal(m1.upper_right, m2.upper_right)
|
||||
assert_equal(m1.width, m2.width)
|
||||
|
||||
# Compare tallies
|
||||
assert len(sp1.tallies) == len(sp2.tallies)
|
||||
|
|
@ -64,11 +64,11 @@ for t1, t2 in zip(sp1.tallies, sp2.tallies):
|
|||
for f1, f2 in zip(t1.filters.values(), t2.filters.values()):
|
||||
assert f1.type == f2.type
|
||||
assert f1.length == f2.length
|
||||
assert f1.bins == f2.bins
|
||||
assert_equal(f1.bins, f2.bins)
|
||||
|
||||
# Compare nuclide and score bins
|
||||
assert t1.nuclides == t2.nuclides
|
||||
assert t1.scores == t2.scores
|
||||
assert_equal(t1.nuclides, t2.nuclides)
|
||||
assert_equal(t1.scores, t2.scores)
|
||||
|
||||
# Compare tally results
|
||||
assert_allclose(t1.results, t2.results)
|
||||
|
|
@ -80,7 +80,7 @@ if sp1.run_mode == 2:
|
|||
|
||||
assert len(sp1.source) == len(sp2.source)
|
||||
for s1, s2 in zip(sp1.source, sp2.source):
|
||||
assert_allclose(s1.weight, s2.weight)
|
||||
assert s1.weight == s2.weight
|
||||
assert_allclose(s1.xyz, s2.xyz)
|
||||
assert_allclose(s1.uvw, s2.uvw)
|
||||
assert_allclose(s1.E, s2.E)
|
||||
assert s1.E == s2.E
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue