Added Python scripts to analyze state points.

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 "<Mesh: {0}>".format(tuple(self.dimension))
+        else:
+            return "<Mesh>"
+
+class Filter(object):
+    def __init__(self):
+        self.type = 0
+        self.bins = []
+
+    def __repr__(self):
+        return "<Filter: {0}>".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)]

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()

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()