From 0d874cc3b3e7be7683aae9cfa6b6a6cea5afeeda Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Mar 2013 17:08:43 -0400 Subject: [PATCH] Added HDF5 support in statepoint.py. Closes #152 on github. --- src/hdf5_interface.F90 | 18 ++-- src/utils/statepoint.py | 205 ++++++++++++++++++++++-------------- src/utils/statepoint_cmp.py | 20 ++-- 3 files changed, 146 insertions(+), 97 deletions(-) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 4340e30aaa..29c36aaae0 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -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", & diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 1227d89cbe..a717a17270 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -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 "".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] diff --git a/src/utils/statepoint_cmp.py b/src/utils/statepoint_cmp.py index c4ef18a810..0e3d8a94f6 100755 --- a/src/utils/statepoint_cmp.py +++ b/src/utils/statepoint_cmp.py @@ -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