diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py new file mode 100644 index 0000000000..19b53b9c73 --- /dev/null +++ b/src/utils/statepoint.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python + +import struct + +filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', + 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} + +score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter', + -5: 'scatter-1', -6: 'scatter-2', -7: 'scatter-3', + -8: 'transport', -9: 'diffusion', -10: 'n1n', -11: 'n2n', + -12: 'n3n', -13: 'n4n', -14: 'absorption', -15: 'fission', + -16: 'nu-fission', -17: 'current'} + +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) + +class Mesh(object): + def __init__(self): + pass + + def __repr__(self): + if hasattr(self, "dimension"): + return "".format(tuple(self.dimension)) + else: + return "" + +class Filter(object): + def __init__(self): + self.type = 0 + self.bins = [] + + def __repr__(self): + return "".format(self.type) + +class Tally(object): + def __init__(self): + self.filters = [] + +class StatePoint(BinaryFile): + def __init__(self, filename): + super(StatePoint, self).__init__(filename) + + # Initialize arrays for meshes and tallies + self.meshes = [] + self.tallies = [] + + # Read all metadata + self._read_metadata() + + def _read_metadata(self): + # Read statepoint revision + self.revision = self.get_int()[0] + + # Read OpenMC version + self.version = self.get_int(3) + + # Read random number seed + self.seed = self.get_long()[0] + + # Read run information + self.run_mode = self.get_int()[0] + self.n_particles = self.get_long()[0] + self.n_batches, self.n_inactive, self.gen_per_batch = self.get_int(3) + + # Read current batch + self.current_batch = self.get_int()[0] + + # Read batch keff and entropy + keff = self.get_double(self.current_batch) + entropy = self.get_double(self.current_batch) + + # Read global tallies + self.n_global_tallies = self.get_int()[0] + self.global_tallies = self.get_double(2*self.n_global_tallies) + + # Read number of meshes + n_meshes = self.get_int()[0] + + # Read meshes + for i in range(n_meshes): + m = Mesh() + self.meshes.append(m) + + # Read type of mesh and number of dimensions + m.type = self.get_int()[0] + n = self.get_int()[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) + + + # Read number of tallies + n_tallies = self.get_int()[0] + + for i in range(n_tallies): + # Create Tally object and add to list of tallies + t = Tally() + self.tallies.append(t) + + # Read sizes of tallies + t.n_score_bins = self.get_int()[0] + t.n_filter_bins = self.get_int()[0] + + # Read number of filters + n_filters = self.get_int()[0] + + for j in range(n_filters): + # Create Filter object and add to list of filters + f = Filter() + t.filters.append(f) + + # Get type of filter + f.type = filter_types[self.get_int()[0]] + + # Determine how many bins are in this filter + f.length = self.get_int()[0] + assert f.length > 0 + if f.type == 'energyin' or f.type == 'energyout': + f.bins = self.get_double(f.length + 1) + elif f.type == 'mesh': + f.bins = self.get_int() + else: + f.bins = self.get_int(f.length) + + # Read nuclide bins + n_nuclides = self.get_int()[0] + t.nuclides = self.get_int(n_nuclides) + + # Read score bins + n_scores = self.get_int()[0] + t.scores = [score_types[j] for j in self.get_int(n_scores)] diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py new file mode 100755 index 0000000000..da86e9cecd --- /dev/null +++ b/src/utils/statepoint_histogram.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python + +from sys import argv +from struct import unpack +from math import sqrt + +import scipy.stats +import matplotlib.pyplot as plt + +from statepoint import StatePoint + +# Get filename +filename = argv[1] + +# Determine axial level +axial_level = int(argv[2]) if len(argv) > 2 else 1 +score = int(argv[3]) if len(argv) > 3 else 1 + +# Create StatePoint object +sp = StatePoint(filename) + +# Calculate t-value for 95% two-sided CI +n = sp.current_batch - sp.n_inactive +t_value = scipy.stats.t.ppf(0.975, n) + +# Loop over all tallies +for i, t in enumerate(sp.tallies): + # Create lists for tallies + mean = [] + uncertainties = [] + nonzero = [] + + n_bins = t.n_score_bins * t.n_filter_bins + + # Loop over filter/score bins + for j in range(n_bins): + # Read sum and sum-squared + s, s2 = unpack('=2d', sp._f.read(16)) + s /= n + if s != 0.0: + relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s + uncertainties.append(relative_error) + + # Display number of tallies with less than 1% CIs + fraction = (sum([1.0 if re < 0.01 else 0.0 for re in uncertainties]) + /n_bins) + print("Tally " + str(i+1)) + print(" Fraction under 1% = {0}".format(fraction)) + print(" Min relative error = {0}".format(min(uncertainties))) + print(" Max relative error = {0}".format(max(uncertainties))) + print(" Non-scoring bins = {0}".format(1.0 - float(len(uncertainties))/n_bins)) + + # Plot histogram + plt.hist(uncertainties, 100) + plt.show() diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py new file mode 100755 index 0000000000..0039ff6d73 --- /dev/null +++ b/src/utils/statepoint_meshplot.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +from sys import argv +from struct import unpack +from math import sqrt + +import numpy as np +import scipy.stats +import matplotlib.pyplot as plt + +from statepoint import StatePoint + +# Get filename +filename = argv[1] + +# Determine axial level +axial_level = int(argv[2]) if len(argv) > 2 else 1 +score = int(argv[3]) if len(argv) > 3 else 1 + +# Create StatePoint object +sp = StatePoint(filename) + +# Calculate t-value for 95% two-sided CI +n = sp.current_batch - sp.n_inactive +t_value = scipy.stats.t.ppf(0.975, n) + +# Loop over all tallies +print("Reading data...") +for t in sp.tallies: + n_bins = t.n_score_bins * t.n_filter_bins + i_mesh = [f.type for f in t.filters].index('mesh') + + # Get Mesh object + m = sp.meshes[t.filters[i_mesh].bins[0] - 1] + nx, ny, nz = m.dimension + ns = t.n_score_bins * t.n_filter_bins / (nx*ny*nz) + + assert n_bins == nx*ny*nz*ns + + # Create lists for tallies + mean = np.zeros((nx,ny)) + error = np.zeros((nx,ny)) + criteria = np.zeros((nx,ny)) + + # Determine starting position for data + start = sp._f.tell() + (axial_level-1)*ns*16 + (score - 1)*16 + + for x in range(nx): + for y in range(ny): + # Seek to position of data + sp._f.seek(start + x*ny*nz*ns*16 + y*nz*ns*16) + + # Read sum and sum-squared + s, s2 = unpack('=2d', sp._f.read(16)) + s /= n + mean[x,y] = s + if s != 0.0: + error[x,y] = t_value*sqrt((s2/n - s*s)/(n-1))/s + criteria[x,y] = 1.0 if error[x,y] < 0.05 else 0.0 + + # Make figure + print("Making colorplot...") + plt.imshow(error, interpolation="nearest") + plt.colorbar() + plt.xlim((0,nx)) + plt.ylim((0,ny)) + plt.xticks(np.linspace(0,nx,5), + np.linspace(m.lower_left[0],m.upper_right[0],5)) + plt.yticks(np.linspace(0,ny,5), + np.linspace(m.lower_left[1],m.upper_right[1],5)) + plt.show()