Move Python API into openmc/ directory and move Python scripts into scripts/

directory.

scripts/convert_binary.py

@@ -0,0 +1,77 @@
+#!/usr/bin/env python
+
+from __future__ import division
+
+from struct import pack
+import sys
+
+
+def ascii_to_binary(ascii_file, binary_file):
+    """Convert an ACE file in ASCII format (type 1) to binary format (type 2).
+
+    Parameters
+    ----------
+    ascii_file : str
+        Filename of ASCII ACE file
+    binary_file : str
+        Filename of binary ACE file to be written
+
+    """
+
+    # Open ASCII file
+    ascii = open(ascii_file, 'r')
+
+    # Set default record length
+    record_length = 4096
+
+    # Read data from ASCII file
+    lines = ascii.readlines()
+    ascii.close()
+
+    # Open binary file
+    binary = open(binary_file, 'wb')
+
+    idx = 0
+    while idx < len(lines):
+        # Read/write header block
+        hz = lines[idx][:10].encode('UTF-8')
+        aw0 = float(lines[idx][10:22])
+        tz = float(lines[idx][22:34])
+        hd = lines[idx][35:45].encode('UTF-8')
+        hk = lines[idx + 1][:70].encode('UTF-8')
+        hm = lines[idx + 1][70:80].encode('UTF-8')
+        binary.write(pack('=10sdd10s70s10s', hz, aw0, tz, hd, hk, hm))
+
+        # Read/write IZ/AW pairs
+        data = ' '.join(lines[idx + 2:idx + 6]).split()
+        iz = list(map(int, data[::2]))
+        aw = list(map(float, data[1::2]))
+        izaw = [item for sublist in zip(iz, aw) for item in sublist]
+        binary.write(pack('=' + 16*'id', *izaw))
+
+        # Read/write NXS and JXS arrays. Null bytes are added at the end so
+        # that XSS will start at the second record
+        nxs = list(map(int, ' '.join(lines[idx + 6:idx + 8]).split()))
+        jxs = list(map(int, ' '.join(lines[idx + 8:idx + 12]).split()))
+        binary.write(pack('=16i32i{0}x'.format(record_length - 500), *(nxs + jxs)))
+
+        # Read/write XSS array. Null bytes are added to form a complete record
+        # at the end of the file
+        n_lines = (nxs[0] + 3)//4
+        xss = list(map(float, ' '.join(lines[idx + 12:idx + 12 + n_lines]).split()))
+        extra_bytes = record_length - ((len(xss)*8 - 1) % record_length + 1)
+        binary.write(pack('={0}d{1}x'.format(nxs[0], extra_bytes), *xss))
+
+        # Advance to next table in file
+        idx += 12 + n_lines
+
+    # Close binary file
+    binary.close()
+
+if __name__ == '__main__':
+    # Check for proper number of arguments
+    if len(sys.argv) < 3:
+        sys.exit('Usage: {0} ascii_file binary_file'.format(sys.argv[0]))
+
+    # Convert ASCII file
+    ascii_to_binary(sys.argv[1], sys.argv[2])

scripts/convert_xsdata.py

@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+
+import os
+import sys
+from xml.dom.minidom import getDOMImplementation
+
+types = {1: "neutron", 2: "dosimetry", 3: "thermal"}
+
+
+class Xsdata(object):
+
+    def __init__(self, filename):
+        self._table_dict = {}
+        self.tables = []
+
+        for line in open(filename, 'r'):
+            words = line.split()
+
+            # If this listing is just an alias listing, only assign the alias
+            # attribute
+            name = words[1]
+            alias = words[0]
+            table = self.find_table(name)
+            if table:
+                if name not in table.alias:
+                    table.alias.append(alias)
+                continue
+
+            table = XsdataTable()
+            table.name = name
+            table.type = types[int(words[2])]
+            table.zaid = int(words[3])
+            table.metastable = int(words[4])
+            table.awr = float(words[5])
+            table.temperature = 8.6173423e-11 * float(words[6])
+            table.binary = int(words[7])
+            table.path = words[8]
+
+            self.tables.append(table)
+            self._table_dict[name] = table
+
+        # Check for common directory
+        self.directory = os.path.dirname(self.tables[0].path)
+        for table in self.tables:
+            if not table.path.startswith(self.directory):
+                self.directory = None
+                break
+
+    def to_xml(self):
+        # Create XML document
+        impl = getDOMImplementation()
+        doc = impl.createDocument(None, "cross_sections", None)
+
+        # Get root element
+        root = doc.documentElement
+
+        # Add a directory node
+        if self.directory:
+            directoryNode = doc.createElement("directory")
+            text = doc.createTextNode(self.directory)
+            directoryNode.appendChild(text)
+            root.appendChild(directoryNode)
+
+            for table in self.tables:
+                table.path = os.path.basename(table.path)
+
+        # Add a node for each table
+        for table in self.tables:
+            node = table.to_xml_node(doc)
+            root.appendChild(node)
+
+        return doc
+
+    def find_table(self, name):
+        if name in self._table_dict:
+            return self._table_dict[name]
+        else:
+            return None
+
+
+class XsdataTable(object):
+
+    def __init__(self):
+        self.alias = []
+
+    def to_xml_node(self, doc):
+        node = doc.createElement("ace_table")
+        node.setAttribute("name", self.name)
+        for attribute in ["alias", "zaid", "type", "metastable",
+                          "awr", "temperature", "binary", "path"]:
+            if hasattr(self, attribute):
+                # Join string for alias attribute
+                if attribute == "alias":
+                    if not self.alias:
+                        continue
+                    string = " ".join(self.alias)
+                else:
+                    string = "{0}".format(getattr(self, attribute))
+
+                # Skip metastable and binary if 0
+                if attribute == "metastable" and self.metastable == 0:
+                    continue
+                if attribute == "binary" and self.binary == 0:
+                    continue
+
+                # Create attribute node
+                # nodeAttr = doc.createElement(attribute)
+                # text = doc.createTextNode(string)
+                # nodeAttr.appendChild(text)
+                # node.appendChild(nodeAttr)
+                node.setAttribute(attribute, string)
+        return node
+
+
+if __name__ == '__main__':
+    # Read command line arguments
+    if len(sys.argv) < 3:
+        sys.exit("Usage: convert_xsdata.py  xsdataFile  xmlFile")
+    xsdataFile = sys.argv[1]
+    xmlFile = sys.argv[2]
+
+    # Read xsdata and create XML document object
+    xsdataObject = Xsdata(xsdataFile)
+    doc = xsdataObject.to_xml()
+
+    # Reduce number of lines
+    lines = doc.toprettyxml(indent='  ')
+    lines = lines.replace('<alias>\n      ', '<alias>')
+    lines = lines.replace('\n    </alias>', '</alias>')
+    lines = lines.replace('<zaid>\n      ', '<zaid>')
+    lines = lines.replace('\n    </zaid>', '</zaid>')
+    lines = lines.replace('<type>\n      ', '<type>')
+    lines = lines.replace('\n    </type>', '</type>')
+    lines = lines.replace('<awr>\n      ', '<awr>')
+    lines = lines.replace('\n    </awr>', '</awr>')
+    lines = lines.replace('<temperature>\n      ', '<temperature>')
+    lines = lines.replace('\n    </temperature>', '</temperature>')
+    lines = lines.replace('<path>\n      ', '<path>')
+    lines = lines.replace('\n    </path>', '</path>')
+    lines = lines.replace('<metastable>\n      ', '<metastable>')
+    lines = lines.replace('\n    </metastable>', '</metastable>')
+    lines = lines.replace('<binary>\n      ', '<binary>')
+    lines = lines.replace('\n    </binary>', '</binary>')
+
+    # Write document in pretty XML to specified file
+    f = open(xmlFile, 'w')
+    f.write(lines)
+    f.close()

scripts/convert_xsdir.py

@@ -0,0 +1,288 @@
+#!/usr/bin/env python
+
+import os
+import sys
+from xml.dom.minidom import getDOMImplementation
+
+elements = [None, "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na",
+            "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V",
+            "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se",
+            "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh",
+            "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba",
+            "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho",
+            "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt",
+            "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac",
+            "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
+            "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg",
+            "Cn"]
+
+
+class Xsdir(object):
+
+    def __init__(self, filename):
+        self.f = open(filename, 'r')
+        self.filename = os.path.abspath(filename)
+        self.directory = os.path.dirname(filename)
+        self.awr = {}
+        self.tables = []
+
+        self.filetype = set()
+        self.recordlength = set()
+        self.entries = set()
+
+        # Read first section (DATAPATH)
+        line = self.f.readline()
+        words = line.split()
+        if words:
+            if words[0].lower().startswith('datapath'):
+                if '=' in words[0]:
+                    index = line.index('=')
+                    self.datapath = line[index+1:].strip()
+                else:
+                    if len(line.strip()) > 8:
+                        self.datapath = line[8:].strip()
+            else:
+                self.f.seek(0)
+
+        # Read second section
+        line = self.f.readline()
+        words = line.split()
+        assert len(words) == 3
+        assert words[0].lower() == 'atomic'
+        assert words[1].lower() == 'weight'
+        assert words[2].lower() == 'ratios'
+
+        while True:
+            line = self.f.readline()
+            words = line.split()
+
+            # Check for end of second section
+            if len(words) % 2 != 0 or words[0] == 'directory':
+                break
+
+            for zaid, awr in zip(words[::2], words[1::2]):
+                self.awr[zaid] = awr
+
+        # Read third section
+        while words[0] != 'directory':
+            words = self.f.readline().split()
+
+        while True:
+            words = self.f.readline().split()
+            if not words:
+                break
+
+            # Handle continuation lines
+            while words[-1] == '+':
+                extraWords = self.f.readline().split()
+                words = words[:-1] + extraWords
+            assert len(words) >= 7
+
+            # Create XsdirTable object and add to line
+            table = XsdirTable(self.directory)
+            self.tables.append(table)
+
+            # All tables have at least 7 attributes
+            table.name = words[0]
+            table.awr = float(words[1])
+            table.filename = words[2]
+            table.access = words[3]
+            table.filetype = int(words[4])
+            table.location = int(words[5])
+            table.length = int(words[6])
+
+            self.filetype.add(table.filetype)
+
+            if len(words) > 7:
+                table.recordlength = int(words[7])
+                self.recordlength.add(table.recordlength)
+            if len(words) > 8:
+                table.entries = int(words[8])
+                self.entries.add(table.entries)
+            if len(words) > 9:
+                table.temperature = float(words[9])
+            if len(words) > 10:
+                table.ptable = (words[10] == 'ptable')
+
+        if len(self.filetype) == 1:
+            if 1 in self.filetype:
+                self.filetype = 'ascii'
+            elif 2 in self.filetype:
+                self.filetype = 'binary'
+        else:
+            self.filetype = None
+
+        if len(self.recordlength) == 1:
+            self.recordlength = list(self.recordlength)[0]
+        else:
+            self.recordlength = None
+        if len(self.entries) == 1:
+            self.entries = list(self.entries)[0]
+        else:
+            self.recordlength = None
+
+    def to_xml(self):
+        # Create XML document
+        impl = getDOMImplementation()
+        doc = impl.createDocument(None, "cross_sections", None)
+
+        # Get root element
+        root = doc.documentElement
+
+        # Add a directory node
+        if self.directory:
+            directoryNode = doc.createElement("directory")
+            text = doc.createTextNode(self.directory)
+            directoryNode.appendChild(text)
+            root.appendChild(directoryNode)
+
+            for table in self.tables:
+                table.path = os.path.basename(table.path)
+
+        # Add filetype, record_length and entries nodes
+        if self.filetype:
+            node = doc.createElement("filetype")
+            text = doc.createTextNode(self.filetype)
+            node.appendChild(text)
+            root.appendChild(node)
+        if self.recordlength:
+            node = doc.createElement("record_length")
+            text = doc.createTextNode(str(self.recordlength))
+            node.appendChild(text)
+            root.appendChild(node)
+        if self.entries:
+            node = doc.createElement("entries")
+            text = doc.createTextNode(str(self.entries))
+            node.appendChild(text)
+            root.appendChild(node)
+
+        # Add a node for each table
+        for table in self.tables:
+            if table.name[-1] in ['e', 'p', 'u', 'h', 'g', 'm', 'd']:
+                continue
+            node = table.to_xml_node(doc)
+            root.appendChild(node)
+
+        return doc
+
+
+class XsdirTable(object):
+
+    def __init__(self, directory=None):
+        self.directory = None
+        self.name = None
+        self.awr = None
+        self.filename = None
+        self.access = None
+        self.filetype = None
+        self.location = None
+        self.length = None
+        self.recordlength = None
+        self.entries = None
+        self.temperature = None
+        self.ptable = False
+
+    @property
+    def path(self):
+        if self.directory:
+            return os.path.join(self.directory, self.filename)
+        else:
+            return self.filename
+
+    @path.setter
+    def path(self, value):
+        self.diretory = ''
+        self.filename = value
+
+    @property
+    def metastable(self):
+        # Only valid for neutron cross-sections
+        if not self.name.endswith('c'):
+            return
+
+        # Handle special case of Am-242 and Am-242m
+        if self.zaid == '95242':
+            return 1
+        elif self.zaid == '95642':
+            return 0
+
+        # All other cases
+        A = int(self.zaid) % 1000
+        if A > 300:
+            return 1
+        else:
+            return 0
+
+    @property
+    def alias(self):
+        zaid = self.zaid
+        if zaid:
+            Z = int(zaid[:-3])
+            A = zaid[-3:]
+
+            if A == '000':
+                s = 'Nat'
+            elif zaid == '95242':
+                s = '242m'
+            elif zaid == '95642':
+                s = '242'
+            elif int(A) > 300:
+                s = str(int(A) - 400) + "m"
+            else:
+                s = str(int(A))
+
+            return "{0}-{1}.{2}".format(elements[Z], s, self.xs)
+        else:
+            return None
+
+    @property
+    def zaid(self):
+        if self.name.endswith('c'):
+            return self.name[:self.name.find('.')]
+        else:
+            return 0
+
+    @property
+    def xs(self):
+        return self.name[self.name.find('.')+1:]
+
+    def to_xml_node(self, doc):
+        node = doc.createElement("ace_table")
+        node.setAttribute("name", self.name)
+        for attribute in ["alias", "zaid", "type", "metastable", "awr",
+                          "temperature", "path", "location"]:
+            if hasattr(self, attribute):
+                string = str(getattr(self, attribute))
+
+                # Skip metastable and binary if 0
+                if attribute == "metastable" and self.metastable == 0:
+                    continue
+
+                # Skip any attribute that is none
+                if getattr(self, attribute) is None:
+                    continue
+
+                # Create attribute node
+                node.setAttribute(attribute, string)
+
+        return node
+
+
+if __name__ == '__main__':
+    # Read command line arguments
+    if len(sys.argv) < 3:
+        sys.exit("Usage: convert_xsdir.py  xsdirFile  xmlFile")
+    xsdirFile = sys.argv[1]
+    xmlFile = sys.argv[2]
+
+    # Read xsdata and create XML document object
+    xsdirObject = Xsdir(xsdirFile)
+    doc = xsdirObject.to_xml()
+
+    # Reduce number of lines
+    lines = doc.toprettyxml(indent='  ')
+
+    # Write document in pretty XML to specified file
+    f = open(xmlFile, 'w')
+    f.write(lines)
+    f.close()

scripts/memory_usage.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+"""
+This script reads a cross_sections.out file, adds up the memory usage for each
+nuclide and S(a,b) table, and displays the total memory usage.
+
+"""
+
+from __future__ import print_function
+
+import sys
+import os
+
+if len(sys.argv) > 1:
+    # Get path to cross_sections.out file from command line argument
+    filename = sys.argv[-1]
+else:
+    # Set default path for cross_sections.out
+    filename = 'cross_sections.out'
+    if not os.path.exists(filename):
+        raise OSError('Could not find cross_sections.out file!')
+
+# Open file handle for cross_sections.out file
+f = open(filename, 'r')
+
+# Initialize memory size arrays
+memory_xs = []
+memory_angle = []
+memory_energy = []
+memory_urr = []
+memory_total = []
+memory_sab = []
+
+while True:
+    # Read next line in file
+    line = f.readline()
+
+    # Check for EOF
+    if line == '':
+        break
+
+    # Look for block listing memory usage for a nuclide
+    words = line.split()
+    if len(words) == 2 and words[0] == 'Memory':
+        memory_xs.append(int(f.readline().split()[-2]))
+        memory_angle.append(int(f.readline().split()[-2]))
+        memory_energy.append(int(f.readline().split()[-2]))
+        memory_urr.append(int(f.readline().split()[-2]))
+        memory_total.append(int(f.readline().split()[-2]))
+
+    # Look for memory usage for S(a,b) table
+    if len(words) == 5 and words[1] == 'Used':
+        memory_sab.append(int(words[-2]))
+
+# Write out summary memory usage
+print('Memory Requirements')
+print('  Reaction Cross Sections        = ' + str(sum(memory_xs)))
+print('  Secondary Angle Distributions  = ' + str(sum(memory_angle)))
+print('  Secondary Energy Distributions = ' + str(sum(memory_energy)))
+print('  Probability Tables             = ' + str(sum(memory_urr)))
+print('  S(a,b) Tables                  = ' + str(sum(memory_sab)))
+print('  Total                          = ' + str(sum(memory_total)))

scripts/plot_mesh_tally.py

@@ -0,0 +1,304 @@
+#!/usr/bin/env python
+
+"""Python script to plot tally data generated by OpenMC."""
+
+import os
+import sys
+
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg
+from matplotlib.figure import Figure
+import matplotlib.pyplot as plt
+import numpy as np
+
+from openmc.statepoint import StatePoint
+
+if sys.version_info[0] < 3:
+    import Tkinter as tk
+    import tkFileDialog as filedialog
+    import tkFont as font
+    import tkMessageBox as messagebox
+    import ttk as ttk
+else:
+    import tkinter as tk
+    import tkinter.filedialog as filedialog
+    import tkinter.font as font
+    import tkinter.messagebox as messagebox
+    import tkinter.ttk as ttk
+
+
+class MeshPlotter(tk.Frame):
+    def __init__(self, parent, filename):
+        tk.Frame.__init__(self, parent)
+
+        self.labels = {'cell': 'Cell:', 'cellborn': 'Cell born:',
+                       'surface': 'Surface:', 'material': 'Material:',
+                       'universe': 'Universe:', 'energyin': 'Energy in:',
+                       'energyout': 'Energy out:'}
+
+        self.filterBoxes = {}
+
+        # Read data from source or leakage fraction file
+        self.get_file_data(filename)
+
+        # Set up top-level window
+        top = self.winfo_toplevel()
+        top.title('Mesh Tally Plotter: ' + filename)
+        top.rowconfigure(0, weight=1)
+        top.columnconfigure(0, weight=1)
+        self.grid(sticky=tk.W+tk.N)
+
+        # Create widgets and draw to screen
+        self.create_widgets()
+        self.update()
+
+    def create_widgets(self):
+        figureFrame = tk.Frame(self)
+        figureFrame.grid(row=0, column=0)
+
+        # Create the Figure and Canvas
+        self.dpi = 100
+        self.fig = Figure((5.0, 5.0), dpi=self.dpi)
+        self.canvas = FigureCanvasTkAgg(self.fig, master=figureFrame)
+        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
+
+        # Create the navigation toolbar, tied to the canvas
+        self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame)
+        self.mpl_toolbar.update()
+        self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1)
+
+        # Create frame for comboboxes
+        self.selectFrame = tk.Frame(self)
+        self.selectFrame.grid(row=1, column=0, sticky=tk.W+tk.E)
+
+        # Tally selection
+        labelTally = tk.Label(self.selectFrame, text='Tally:')
+        labelTally.grid(row=0, column=0, sticky=tk.W)
+        self.tallyBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.tallyBox['values'] = [self.datafile.tallies[i].id
+                                   for i in self.meshTallies]
+        self.tallyBox.current(0)
+        self.tallyBox.grid(row=0, column=1, sticky=tk.W+tk.E)
+        self.tallyBox.bind('<<ComboboxSelected>>', self.update)
+
+        # Planar basis selection
+        labelBasis = tk.Label(self.selectFrame, text='Basis:')
+        labelBasis.grid(row=1, column=0, sticky=tk.W)
+        self.basisBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.basisBox['values'] = ('xy', 'yz', 'xz')
+        self.basisBox.current(0)
+        self.basisBox.grid(row=1, column=1, sticky=tk.W+tk.E)
+        self.basisBox.bind('<<ComboboxSelected>>', self.update)
+
+        # Axial level
+        labelAxial = tk.Label(self.selectFrame, text='Axial level:')
+        labelAxial.grid(row=2, column=0, sticky=tk.W)
+        self.axialBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.axialBox.grid(row=2, column=1, sticky=tk.W+tk.E)
+        self.axialBox.bind('<<ComboboxSelected>>', self.redraw)
+
+        # Option for mean/uncertainty
+        labelMean = tk.Label(self.selectFrame, text='Mean/Uncertainty:')
+        labelMean.grid(row=3, column=0, sticky=tk.W)
+        self.meanBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.meanBox['values'] = ('Mean', 'Absolute uncertainty',
+                                  'Relative uncertainty')
+        self.meanBox.current(0)
+        self.meanBox.grid(row=3, column=1, sticky=tk.W+tk.E)
+        self.meanBox.bind('<<ComboboxSelected>>', self.update)
+
+        # Scores
+        labelScore = tk.Label(self.selectFrame, text='Score:')
+        labelScore.grid(row=4, column=0, sticky=tk.W)
+        self.scoreBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.scoreBox.grid(row=4, column=1, sticky=tk.W+tk.E)
+        self.scoreBox.bind('<<ComboboxSelected>>', self.redraw)
+
+        # Filter label
+        boldfont = font.Font(weight='bold')
+        labelFilters = tk.Label(self.selectFrame, text='Filters:',
+                                font=boldfont)
+        labelFilters.grid(row=5, column=0, sticky=tk.W)
+
+    def update(self, event=None):
+        if not event:
+            widget = None
+        else:
+            widget = event.widget
+
+        tally_id = self.meshTallies[self.tallyBox.current()]
+        selectedTally = self.datafile.tallies[tally_id]
+
+        # Get mesh for selected tally
+        self.mesh = selectedTally.filters_by_name['mesh'].mesh
+
+        # Get mesh dimensions
+        self.nx, self.ny, self.nz = self.mesh.dimension
+
+        # Repopulate comboboxes baesd on current basis selection
+        text = self.basisBox.get()
+        if text == 'xy':
+            self.axialBox['values'] = [str(i+1) for i in range(self.nz)]
+        elif text == 'yz':
+            self.axialBox['values'] = [str(i+1) for i in range(self.nx)]
+        else:
+            self.axialBox['values'] = [str(i+1) for i in range(self.ny)]
+        self.axialBox.current(0)
+
+        # If update() was called by a change in the basis combobox, we don't
+        # need to repopulate the filters
+        if widget == self.basisBox:
+            self.redraw()
+            return
+
+        # Update scores
+        self.scoreBox['values'] = selectedTally.scores
+        self.scoreBox.current(0)
+
+        # Remove any filter labels/comboboxes that exist
+        for row in range(6, self.selectFrame.grid_size()[1]):
+            for w in self.selectFrame.grid_slaves(row=row):
+                w.grid_forget()
+                w.destroy()
+
+        # create a label/combobox for each filter in selected tally
+        count = 0
+        for f in selectedTally.filters:
+            filterType = f.type
+            if filterType == 'mesh':
+                continue
+            count += 1
+
+            # Create label and combobox for this filter
+            label = tk.Label(self.selectFrame, text=self.labels[filterType])
+            label.grid(row=count+6, column=0, sticky=tk.W)
+            combobox = ttk.Combobox(self.selectFrame, state='readonly')
+            self.filterBoxes[filterType] = combobox
+
+            # Set combobox items
+            if filterType in ['energyin', 'energyout']:
+                combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2])
+                                      for i in range(f.length)]
+            else:
+                combobox['values'] = [str(i) for i in f.bins]
+
+            combobox.current(0)
+            combobox.grid(row=count+6, column=1, sticky=tk.W+tk.E)
+            combobox.bind('<<ComboboxSelected>>', self.redraw)
+
+        # If There are no filters, leave a 'None available' message
+        if count == 0:
+            count += 1
+            label = tk.Label(self.selectFrame, text="None Available")
+            label.grid(row=count+6, column=0, sticky=tk.W)
+
+        self.redraw()
+
+    def redraw(self, event=None):
+        basis = self.basisBox.current() + 1
+        axial_level = self.axialBox.current() + 1
+        mbvalue = self.meanBox.get()
+
+        # Get selected tally
+        tally_id = self.meshTallies[self.tallyBox.current()]
+        selectedTally = self.datafile.tallies[tally_id]
+
+        # Create spec_list
+        spec_list = []
+        for f in selectedTally.filters:
+            if f.type == 'mesh':
+                mesh_filter = f
+                continue
+            index = self.filterBoxes[f.type].current()
+            spec_list.append((f, index))
+
+        text = self.basisBox.get()
+        if text == 'xy':
+            dims = (self.nx, self.ny)
+        elif text == 'yz':
+            dims = (self.ny, self.nz)
+        else:
+            dims = (self.nx, self.nz)
+        matrix = np.zeros(dims)
+        for i in range(dims[0]):
+            for j in range(dims[1]):
+                if text == 'xy':
+                    meshtuple = (i + 1, j + 1, axial_level)
+                elif text == 'yz':
+                    meshtuple = (axial_level, i + 1, j + 1)
+                else:
+                    meshtuple = (i + 1, axial_level, j + 1)
+                filters, filter_bins = zip(*spec_list + [
+                    (mesh_filter, meshtuple)])
+                mean = selectedTally.get_value(
+                    self.scoreBox.get(), filters, filter_bins)
+                stdev = selectedTally.get_value(
+                    self.scoreBox.get(), filters, filter_bins,
+                    value='std_dev')
+                if mbvalue == 'Mean':
+                    matrix[i, j] = mean
+                elif mbvalue == 'Absolute uncertainty':
+                    matrix[i, j] = stdev
+                else:
+                    if mean > 0.:
+                        matrix[i, j] = stdev/mean
+                    else:
+                        matrix[i, j] = 0.
+
+        # Clear the figure
+        self.fig.clear()
+
+        # Make figure, set up color bar
+        self.axes = self.fig.add_subplot(111)
+        cax = self.axes.imshow(matrix.transpose(), vmin=0.0, vmax=matrix.max(),
+                               interpolation='none', origin='lower')
+        self.fig.colorbar(cax)
+
+        self.axes.set_xticks([])
+        self.axes.set_yticks([])
+        self.axes.set_aspect('equal')
+
+        # Draw canvas
+        self.canvas.draw()
+
+    def get_file_data(self, filename):
+        # Create StatePoint object and read in data
+        self.datafile = StatePoint(filename)
+        self.datafile.read_results()
+        self.datafile.compute_stdev()
+
+        # Find which tallies are mesh tallies
+        self.meshTallies = []
+        for itally, tally in self.datafile.tallies.iteritems():
+            if any([f.type == 'mesh' for f in tally.filters]):
+                self.meshTallies.append(itally)
+                tally.filters_by_name = {f.type: f for f in tally.filters}
+
+        if not self.meshTallies:
+            messagebox.showerror("Invalid StatePoint File",
+                                 "File does not contain mesh tallies!")
+            sys.exit(1)
+
+
+if __name__ == '__main__':
+    # Hide root window
+    root = tk.Tk()
+    root.withdraw()
+
+    # If no filename given as command-line argument, open file dialog
+    if len(sys.argv) < 2:
+        filename = filedialog.askopenfilename(title='Select statepoint file',
+                                              initialdir='.')
+    else:
+        filename = sys.argv[1]
+
+    if filename:
+        # Check to make sure file exists
+        if not os.path.isfile(filename):
+            messagebox.showerror("File not found",
+                                 "Could not find regular file: " + filename)
+            sys.exit(1)
+
+        app = MeshPlotter(root, filename)
+        root.deiconify()
+        root.mainloop()

scripts/statepoint_3d.py

@@ -0,0 +1,393 @@
+#!/usr/bin/env python2
+
+from __future__ import division, print_function
+
+import sys
+import itertools
+import re
+import warnings
+
+from openmc.statepoint import StatePoint
+
+alphanum = re.compile(r"[\W_]+")
+
+err = False
+
+################################################################################
+def parse_options():
+  """Process command line arguments"""
+
+
+
+  def tallies_callback(option, opt, value, parser):
+    """Option parser function for list of tallies"""
+    global err
+    try:
+      setattr(parser.values, option.dest, [int(v) for v in value.split(',')])
+    except:
+      p.print_help()
+      err = True
+
+  def scores_callback(option, opt, value, parser):
+    """Option parser function for list of scores"""
+    global err
+    try:
+      scores = {}
+      entries = value.split(',')
+      for e in entries:
+        tally,score = [int(i) for i in e.split('.')]
+        if not tally in scores: scores[tally] = []
+        scores[tally].append(score)
+      setattr(parser.values, option.dest, scores)
+    except:
+      p.print_help()
+      err = True
+
+  def filters_callback(option, opt, value, parser):
+    """Option parser function for list of filters"""
+    global err
+    try:
+      filters = {}
+      entries = value.split(',')
+      for e in entries:
+        tally,filter_,bin = [i for i in e.split('.')]
+        tally,bin = int(tally),int(bin)
+        if not tally in filters: filters[tally] = {}
+        if not filter_ in filters[tally]: filters[tally][filter_] = []
+        filters[tally][filter_].append(bin)
+      setattr(parser.values, option.dest, filters)
+    except:
+      p.print_help()
+      err = True
+
+  from optparse import OptionParser
+  usage = r"""%prog [options] <statepoint_file>
+
+The default is to process all tallies and all scores into one file. Subsets
+can be chosen using the options.  For example, to only process tallies 2 and 4
+with all scores on tally 2 and only scores 1 and 3 on tally 4:
+
+%prog -t 2,4 -s 4.1,4.3 <statepoint_file>
+
+Likewise if you have additional filters on a tally you can specify a subset of
+bins for each filter for that tally. For example to process all tallies and
+scores, but only energyin bin #1 in tally 2:
+
+%prog -f 2.energyin.1 <statepoint_file>
+
+You can list the available tallies, scores, and filters with the -l option:
+
+%prog -l <statepoint_file>"""
+  p = OptionParser(usage=usage)
+  p.add_option('-t', '--tallies', dest='tallies', type='string', default=None,
+                action='callback', callback=tallies_callback,
+                help='List of tally indices to process, separated by commas.' \
+                     ' Default is to process all tallies.')
+  p.add_option('-s', '--scores', dest='scores', type='string', default=None,
+               action='callback', callback=scores_callback,
+               help='List of score indices to process, separated by commas, ' \
+                 'specified as {tallyid}.{scoreid}.' \
+                 ' Default is to process all scores in each tally.')
+  p.add_option('-f', '--filters', dest='filters', type='string', default=None,
+               action='callback', callback=filters_callback,
+               help='List of filter bins to process, separated by commas, ' \
+                 'specified as {tallyid}.{filter}.{binid}. ' \
+                 'Default is to process all filter combinaiton for each score.')
+  p.add_option('-l', '--list', dest='list', action='store_true',
+               help='List the tally and score indices available in the file.')
+  p.add_option('-o', '--output', action='store', dest='output',
+               default='tally', help='path to output SILO file.')
+  p.add_option('-e', '--error', dest='valerr', default=False,
+               action='store_true',
+               help='Flag to extract errors instead of values.')
+  p.add_option('-v', '--vtk', action='store_true', dest='vtk',
+               default=False, help='Flag to convert to VTK instead of SILO.')
+  parsed = p.parse_args()
+
+  if not parsed[1]:
+    p.print_help()
+    return parsed, err
+
+  if parsed[0].valerr:
+    parsed[0].valerr = 1
+  else:
+    parsed[0].valerr = 0
+
+  return parsed, err
+
+################################################################################
+def main(file_, o):
+  """Main program"""
+
+  sp = StatePoint(file_)
+  sp.read_results()
+
+  validate_options(sp, o)
+
+  if o.list:
+    print_available(sp)
+    return
+
+  if o.vtk:
+    if not o.output[-4:] == ".vtm": o.output += ".vtm"
+  else:
+    if not o.output[-5:] == ".silo": o.output += ".silo"
+
+  if o.vtk:
+    try:
+      import vtk
+    except:
+      print('The vtk python bindings do not appear to be installed properly.\n'
+            'On Ubuntu: sudo apt-get install python-vtk\n'
+            'See: http://www.vtk.org/')
+      return
+  else:
+    try:
+      import silomesh
+    except:
+      print('The silomesh package does not appear to be installed properly.\n'
+            'See: https://github.com/nhorelik/silomesh/')
+      return
+
+  if o.vtk:
+    blocks = vtk.vtkMultiBlockDataSet()
+    blocks.SetNumberOfBlocks(5)
+    block_idx = 0
+  else:
+    silomesh.init_silo(o.output)
+
+  # Tally loop #################################################################
+  for tally in sp.tallies:
+
+    # skip non-mesh tallies or non-user-specified tallies
+    if o.tallies and not tally.id in o.tallies: continue
+    if not 'mesh' in tally.filters: continue
+
+    print("Processing Tally {}...".format(tally.id))
+
+    # extract filter options and mesh parameters for this tally
+    filtercombos = get_filter_combos(tally)
+    meshparms = get_mesh_parms(sp, tally)
+    nx,ny,nz = meshparms[:3]
+    ll = meshparms[3:6]
+    ur = meshparms[6:9]
+
+    if o.vtk:
+      ww = [(u-l)/n for u,l,n in zip(ur,ll,(nx,ny,nz))]
+      grid = grid = vtk.vtkImageData()
+      grid.SetDimensions(nx+1,ny+1,nz+1)
+      grid.SetOrigin(*ll)
+      grid.SetSpacing(*ww)
+    else:
+      silomesh.init_mesh('Tally_{}'.format(tally.id), *meshparms)
+
+    # Score loop ###############################################################
+    for sid,score in enumerate(tally.scores):
+
+      # skip non-user-specified scrores for this tally
+      if o.scores and tally.id in o.scores and not sid in o.scores[tally.id]:
+        continue
+
+      # Filter loop ############################################################
+      for filterspec in filtercombos:
+
+        # skip non-user-specified filter bins
+        skip = False
+        if o.filters and tally.id in o.filters:
+          for filter_,bin in filterspec[1:]:
+            if filter_ in o.filters[tally.id] and \
+               not bin in o.filters[tally.id][filter_]:
+              skip = True
+              break
+        if skip: continue
+
+        # find and sanitize the variable name for this score
+        varname = get_sanitized_filterspec_name(tally, score, filterspec)
+        if o.vtk:
+          vtkdata = vtk.vtkDoubleArray()
+          vtkdata.SetName(varname)
+          dataforvtk = {}
+        else:
+          silomesh.init_var(varname)
+
+        lbl = "\t Score {}.{} {}:\t\t{}".format(tally.id, sid+1, score, varname)
+
+        # Mesh fill loop #######################################################
+        for x in range(1,nx+1):
+          sys.stdout.write(lbl+" {0}%\r".format(int(x/nx*100)))
+          sys.stdout.flush()
+          for y in range(1,ny+1):
+            for z in range(1,nz+1):
+              filterspec[0][1] = (x,y,z)
+              val = sp.get_value(tally.id-1, filterspec, sid)[o.valerr]
+              if o.vtk:
+                # vtk cells go z, y, x, so we store it now and enter it later
+                i = (z-1)*nx*ny + (y-1)*nx + x-1
+                dataforvtk[i] = float(val)
+              else:
+                silomesh.set_value(float(val), x, y, z)
+
+        # end mesh fill loop
+        print()
+        if o.vtk:
+          for i in range(nx*ny*nz):
+            vtkdata.InsertNextValue(dataforvtk[i])
+          grid.GetCellData().AddArray(vtkdata)
+          del vtkdata
+
+        else:
+          silomesh.finalize_var()
+
+      # end filter loop
+
+    # end score loop
+    if o.vtk:
+      blocks.SetBlock(block_idx, grid)
+      block_idx += 1
+    else:
+      silomesh.finalize_mesh()
+
+  # end tally loop
+  if o.vtk:
+    writer = vtk.vtkXMLMultiBlockDataWriter()
+    writer.SetFileName(o.output)
+    writer.SetInput(blocks)
+    writer.Write()
+  else:
+    silomesh.finalize_silo()
+
+################################################################################
+def get_sanitized_filterspec_name(tally, score, filterspec):
+  """Returns a name fit for silo vars for a given filterspec, tally and score"""
+
+  comboname = "_"+" ".join(["{}_{}".format(filter_, bin)
+                        for filter_, bin in filterspec[1:]])
+  if len(filterspec[1:]) == 0: comboname = ''
+  varname = 'Tally_{}_{}{}'.format(tally.id, score, comboname)
+  varname = alphanum.sub('_', varname)
+  return varname
+
+################################################################################
+def get_filter_combos(tally):
+  """Returns a list of all filter spec combinations, excluding meshes
+
+  Each combo has the mesh spec as the first element, to be set later.
+  These filter specs correspond with the second argument to StatePoint.get_value
+  """
+
+  specs = []
+
+  if len(tally.filters) == 1:
+    return [[['mesh', [1, 1, 1]]]]
+
+  filters = list(tally.filters.keys())
+  filters.pop(filters.index('mesh'))
+  nbins = [tally.filters[f].length for f in filters]
+
+  combos = [ [b] for b in range(nbins[0])]
+  for i,b in enumerate(nbins[1:]):
+    prod = list(itertools.product(combos, range(b)))
+    if i == 0:
+      combos = prod
+    else:
+      combos = [[v for v in p[0]] + [p[1]] for p in prod]
+
+  for c in combos:
+    spec = [['mesh', [1, 1, 1]]]
+    for i,bin in enumerate(c):
+      spec.append((filters[i], bin))
+    specs.append(spec)
+
+  return specs
+
+################################################################################
+def get_mesh_parms(sp, tally):
+  meshid = tally.filters['mesh'].bins[0]
+  for i,m in enumerate(sp.meshes):
+    if m.id == meshid:
+      mesh = m
+  return mesh.dimension + mesh.lower_left + mesh.upper_right
+
+################################################################################
+def print_available(sp):
+  """Prints available tallies/scores in a statepoint"""
+
+  print("Available tally and score indices:")
+  for tally in sp.tallies:
+    mesh = ""
+    if not 'mesh' in tally.filters: mesh = "(no mesh)"
+    print("\tTally {} {}".format(tally.id, mesh))
+    scores = ["{}.{}: {}".format(tally.id, sid, score)
+              for sid, score in enumerate(tally.scores)]
+    for score in scores:
+      print("\t\tScore {}".format(score))
+      for filter_ in tally.filters:
+        if filter_ == 'mesh': continue
+        for bin in range(tally.filters[filter_].length):
+          print("\t\t\tFilters: {}.{}.{}".format(tally.id, filter_, bin))
+
+################################################################################
+def validate_options(sp,o):
+  """Validates specified tally/score options for the current statepoint"""
+
+  available_tallies = [t.id for t in sp.tallies]
+  if o.tallies:
+    for otally in o.tallies:
+      if not otally in available_tallies:
+        warnings.warn('Tally {} not in statepoint file'.format(otally))
+        continue
+      else:
+        for tally in sp.tallies:
+          if tally.id == otally: break
+      if not 'mesh' in tally.filters:
+        warnings.warn('Tally {} contains no mesh'.format(otally))
+      if o.scores and otally in o.scores.keys():
+        for oscore in o.scores[otally]:
+          if oscore > len(tally.scores):
+            warnings.warn('No score {} in tally {}'.format(oscore, otally))
+
+  if o.scores:
+    for otally in o.scores.keys():
+      if not otally in available_tallies:
+        warnings.warn('Tally {} not in statepoint file'.format(otally))
+        continue
+      if o.tallies and not otally in o.tallies:
+        warnings.warn(
+          'Skipping scores for tally {}, excluded by tally list'.format(otally))
+        continue
+
+  if o.filters:
+    for otally in o.filters.keys():
+      if not otally in available_tallies:
+        warnings.warn('Tally {} not in statepoint file'.format(otally))
+        continue
+      if o.tallies and not otally in o.tallies:
+        warnings.warn(
+          'Skipping filters for tally {}, excluded by tally list'.format(otally))
+        continue
+      for tally in sp.tallies:
+        if tally.id == otally: break
+      for filter_ in o.filters[otally]:
+        if filter_ == 'mesh':
+          warnings.warn('Cannot specify mesh filter bins')
+          continue
+        if not filter_ in tally.filters.keys():
+          warnings.warn(
+                  'Tally {} does not contain filter {}'.format(otally, filter_))
+          continue
+        for bin in o.filters[otally][filter_]:
+          if bin >= tally.filters[filter_].length:
+            warnings.warn(
+                 'No bin {} in tally {} filter {}'.format(bin, otally, filter_))
+
+################################################################################
+# monkeypatch to suppress the source echo produced by warnings
+def formatwarning(message, category, filename, lineno, line):
+  return "{}:{}: {}: {}\n".format(filename, lineno, category.__name__, message)
+warnings.formatwarning = formatwarning
+
+################################################################################
+if __name__ == '__main__':
+    (options, args), err = parse_options()
+    if args and not err:
+      main(args[0],options)

scripts/statepoint_histogram.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+from sys import argv
+from math import sqrt
+
+import numpy as np
+import scipy.stats
+import matplotlib.pyplot as plt
+
+from openmc.statepoint import StatePoint
+
+# Get filename
+filename = argv[1]
+
+# Create StatePoint object
+sp = StatePoint(filename)
+sp.read_results()
+sp.compute_ci()
+
+# Check if tallies are present
+if not sp.tallies_present:
+    raise Exception("No tally data in state point!")
+
+# Loop over all tallies
+for i, t in sp.tallies.items():
+    # Determine relative error and fraction of bins with less than 1% half-width
+    # of CI
+    n_bins = t.mean.size
+    relative_error = t.std_dev[t.mean > 0.] / t.mean[t.mean > 0.]
+    fraction = float(sum(relative_error < 0.01))/n_bins
+
+    # Display results
+    print("Tally " + str(i))
+    print("  Fraction under 1% = {0}".format(fraction))
+    print("  Min relative error = {0}".format(min(relative_error)))
+    print("  Max relative error = {0}".format(max(relative_error)))
+    print("  Non-scoring bins = {0}".format(
+          1.0 - float(relative_error.size)/n_bins))
+
+    # Plot histogram
+    plt.hist(relative_error, 100)
+    plt.show()

scripts/tally_conv.py

@@ -0,0 +1,375 @@
+#!/usr/bin/env python
+
+# This program takes OpenMC statepoint binary files and creates a variety of
+# outputs from them which should provide the user with an idea of the
+# convergence behavior of all the tallies and filters defined by the user in
+# tallies.xml.  The program can directly plot the value and errors of each
+# tally, filter, score combination; it can save these plots to a file; and 
+# it can also save the data used in these plots to a CSV file for importing in
+# to other plotting packages such as Excel, gnuplot, MathGL, or Veusz.
+
+# To use the program, run this program from the working directory of the openMC
+# problem to analyze.  
+
+# The USER OPTIONS block below provides four options for the user to set:
+# fileType, printxs, showImg, and savetoCSV.  See the options block for more 
+# information.
+
+from __future__ import print_function
+from math import sqrt, pow
+from glob import glob
+
+import numpy as np
+import scipy.stats
+import matplotlib.pyplot as plt
+
+from openmc.statepoint import StatePoint
+
+##################################### USER OPTIONS
+
+# Set filetype (the file extension desired, without the period.)
+# Options are backend dependent, but most backends support png, pdf, ps, eps 
+# and svg.  Write "none" if no saved files are desired.
+fileType = "none"
+
+# Set if cross-sections or reaction rates are desired printxs = True means X/S
+printxs = False
+
+# Set if the figures should be displayed to screen or not (True means show)
+showImg = False
+
+# Save to CSV for use in more advanced plotting programs like GNUPlot, MathGL
+savetoCSV = True
+
+##################################### END USER OPTIONS
+
+## Find if tallies.xml exists.
+#if glob('./tallies.xml') != None:
+#    # It exists
+#    tallyData = talliesXML('tallies.xml')
+#else: 
+#    # It does not exist.
+#    tallyData = None
+
+# Find all statepoints in this directory.
+files = glob('./statepoint.*.binary')
+fileNums = []
+begin = 13
+# Arrange the file list in increasing batch order
+for i in range(len(files)):
+    end = files[i].find(".binary")
+    fileNums.append(int(files[i][begin:end]))
+fileNums.sort()
+# Re-make filenames
+files = []
+for i in range(len(fileNums)):
+    files.append("./statepoint." + str(fileNums[i]) + ".binary")
+
+# Initialize arrays as needed
+mean = [None for x in range(len(files))]
+uncert = [None for x in range(len(files))]
+scoreType = [None for x in range(len(files))]
+active_batches = [None for x in range(len(files))]
+
+for i_batch in range(len(files)):
+    
+    # Get filename    
+    batch_filename = files[i_batch]
+    
+    # Create StatePoint object
+    sp = StatePoint(batch_filename)
+
+    # Read number of realizations for global tallies
+    sp.n_realizations = sp._get_int()[0]
+    
+    # Read global tallies
+    n_global_tallies = sp._get_int()[0]
+    sp.global_tallies = np.array(sp._get_double(2*n_global_tallies))
+    sp.global_tallies.shape = (n_global_tallies, 2)
+
+    # Flag indicating if tallies are present
+    tallies_present = sp._get_int()[0]
+
+    # Check if tallies are present
+    if not tallies_present:
+        raise Exception("No tally data in state point!")
+    
+    # Increase the dimensionality of our main variables
+    mean[i_batch] = [None for x in range(len(sp.tallies))]
+    uncert[i_batch] = [None for x in range(len(sp.tallies))]
+    scoreType[i_batch] = [None for x in range(len(sp.tallies))]     
+    
+    # Loop over all tallies
+    for i_tally, t in enumerate(sp.tallies):
+        # Calculate t-value for 95% two-sided CI
+        n = t.n_realizations
+        t_value = scipy.stats.t.ppf(0.975, n - 1)
+    
+        # Store the batch count    
+        active_batches[i_batch] = n
+    
+        # Resize the 2nd dimension      
+        mean[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
+        uncert[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
+        scoreType[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
+        
+        for i_filter in range(t.total_filter_bins):
+            # Resize the 3rd dimension            
+            mean[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
+            uncert[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
+            scoreType[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
+            print(t.total_filter_bins,t.n_nuclides)            
+            for i_nuclide in range(t.n_nuclides):
+                mean[i_batch][i_tally][i_filter][i_nuclide] = \
+                    [None for x in range(t.n_scores)]
+                uncert[i_batch][i_tally][i_filter][i_nuclide] = \
+                    [None for x in range(t.n_scores)]
+                scoreType[i_batch][i_tally][i_filter][i_nuclide] = \
+                    [None for x in range(t.n_scores)]
+                for i_score in range(t.n_scores):
+                    scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score] = \
+                        t.scores[i_score]                 
+                    s, s2 = sp._get_double(2)
+                    s /= n
+                    mean[i_batch][i_tally][i_filter][i_nuclide][i_score] = s
+                    if s != 0.0:
+                        relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s
+                    else:
+                        relative_error = 0.0
+                    uncert[i_batch][i_tally][i_filter][i_nuclide][i_score] = relative_error
+
+# Reorder the data lists in to a list order more conducive for plotting:
+# The indexing should be: [tally][filter][score][batch]
+meanPlot = [None for x in range(len(mean[0]))] # Set to the number of tallies
+uncertPlot = [None for x in range(len(mean[0]))] # Set to the number of tallies
+absUncertPlot = [None for x in range(len(mean[0]))] # Set to number of tallies
+filterLabel = [None for x in range(len(mean[0]))] #Set to the number of tallies
+fluxLoc = [None for x in range(len(mean[0]))] # Set to the number of tallies
+printxs = [False for x in range(len(mean[0]))] # Set to the number of tallies
+
+# Get and set the correct sizes for the rest of the dimensions
+for i_tally in range(len(meanPlot)):
+    # Set 2nd (score) dimension
+    meanPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
+    uncertPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
+    absUncertPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
+    filterLabel[i_tally] = [None for x in range(len(mean[0][i_tally]))]
+    
+    # Initialize flux location so it will be -1 if not found
+    fluxLoc[i_tally] = -1
+    
+    for i_filter in range(len(meanPlot[i_tally])):
+        # Set 3rd (filter) dimension
+        meanPlot[i_tally][i_filter] = \
+            [None for x in range(len(mean[0][i_tally][i_filter]))]
+        uncertPlot[i_tally][i_filter] = \
+            [None for x in range(len(mean[0][i_tally][i_filter]))]
+        absUncertPlot[i_tally][i_filter] = \
+            [None for x in range(len(mean[0][i_tally][i_filter]))]
+        filterLabel[i_tally][i_filter] = \
+            [None for x in range(len(mean[0][i_tally][i_filter]))]
+            
+        for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
+            # Set 4th (nuclide)) dimension
+            meanPlot[i_tally][i_filter][i_nuclide] = \
+                [None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
+            uncertPlot[i_tally][i_filter][i_nuclide] = \
+                [None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
+            absUncertPlot[i_tally][i_filter][i_nuclide] = \
+                [None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
+        
+            for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
+                # Set 5th (batch) dimension
+                meanPlot[i_tally][i_filter][i_nuclide][i_score] = \
+                    [None for x in range(len(mean))]
+                uncertPlot[i_tally][i_filter][i_nuclide][i_score] = \
+                    [None for x in range(len(mean))]
+                absUncertPlot[i_tally][i_filter][i_nuclide][i_score] = \
+                    [None for x in range(len(mean))]
+                    
+                # Get filterLabel (this should be moved to its own function)
+                #??? How to do?
+            
+        # Set flux location if found
+        # all batches and all tallies will have the same score ordering, hence
+        # the 0's in the 1st, 3rd, and 4th dimensions.
+        if scoreType[0][i_tally][0][0][i_score] == 'flux': 
+            fluxLoc[i_tally] = i_score
+
+# Set printxs array according to the printxs input
+if printxs:
+    for i_tally in range(len(fluxLoc)):
+        if fluxLoc[i_tally] != -1:
+            printxs[i_tally] = True
+        
+# Now rearrange the data as suitable, and perform xs conversion if necessary
+for i_batch in range(len(mean)):
+    for i_tally in range(len(mean[i_batch])):
+        for i_filter in range(len(mean[i_batch][i_tally])):
+            for i_nuclide in range(len(mean[i_batch][i_tally][i_filter])):
+                for i_score in range(len(mean[i_batch][i_tally][i_filter][i_nuclide])):
+                    if (printxs[i_tally] and \
+                        ((scoreType[0][i_tally][i_filter][i_nuclide][i_score] != 'flux') and \
+                        (scoreType[0][i_tally][i_filter][i_nuclide][i_score] != 'current'))):
+                        
+                        # Perform rate to xs conversion
+                        # mean is mean/fluxmean
+                        meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
+                            mean[i_batch][i_tally][i_filter][i_nuclide][i_score] / \
+                            mean[i_batch][i_tally][i_filter][i_nuclide][fluxLoc[i_tally]]
+                        
+                        # Update the relative uncertainty via error propagation
+                        uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
+                            sqrt(pow(uncert[i_batch][i_tally][i_filter][i_nuclide][i_score],2) \
+                            + pow(uncert[i_batch][i_tally][i_filter][i_nuclide][fluxLoc[i_tally]],2))
+                    else: 
+                        
+                        # Do not perform rate to xs conversion
+                        meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
+                            mean[i_batch][i_tally][i_filter][i_nuclide][i_score]
+                        uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
+                            uncert[i_batch][i_tally][i_filter][i_nuclide][i_score]
+
+                    # Both have the same absolute uncertainty calculation
+                    absUncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
+                        uncert[i_batch][i_tally][i_filter][i_nuclide][i_score] * \
+                        mean[i_batch][i_tally][i_filter][i_nuclide][i_score]    
+
+# Set plotting constants
+xLabel = "Batches"
+xLabel = xLabel.title() # not necessary for now, but is left in to handle if 
+# the previous line changes
+
+# Begin plotting
+for i_tally in range(len(meanPlot)):
+    # Set tally string (placeholder until I put tally labels in statePoint)
+    tallyStr = "Tally " + str(i_tally + 1)
+    
+    for i_filter in range(len(meanPlot[i_tally])):
+
+        # Set filter string
+        filterStr = "Filter " + str(i_filter + 1)
+        
+        for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
+            
+            nuclideStr = "Nuclide " + str(i_nuclide + 1)
+        
+            for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
+    
+                # Set score string
+                scoreStr = scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score]
+                scoreStr = scoreStr.title()
+                if (printxs[i_tally] and ((scoreStr != 'Flux') and \
+                    (scoreStr != 'Current'))):
+                    scoreStr = scoreStr + "-XS"
+                
+                # set Title 
+                title = "Convergence of " + scoreStr + " in " + tallyStr + " for "\
+                    + filterStr + " and " + nuclideStr
+                
+                # set yLabel
+                yLabel = scoreStr
+                yLabel = yLabel.title()
+    
+                # Set saving filename
+                fileName = "tally_" + str(i_tally + 1) + "_" + scoreStr + \
+                    "_filter_" + str(i_filter+1) + "_nuclide_" + str(i_nuclide+1) \
+                    + "." + fileType
+                REfileName = "tally_" + str(i_tally + 1) + "_" + scoreStr + \
+                    "RE_filter_" + str(i_filter+1) + "_nuclide_" + str(i_nuclide+1) \
+                    + "." + fileType
+                
+                # Plot mean with absolute error bars
+                plt.errorbar(active_batches, \
+                    meanPlot[i_tally][i_filter][i_nuclide][i_score][:], \
+                    absUncertPlot[i_tally][i_filter][i_nuclide][i_score][:],fmt='o-',aa=True)
+                plt.xlabel(xLabel)
+                plt.ylabel(yLabel)
+                plt.title(title)
+                if (fileType != 'none'):
+                    plt.savefig(fileName)
+                if showImg:            
+                    plt.show()
+                plt.clf()
+                
+                # Plot relative uncertainty
+                plt.plot(active_batches, \
+                    uncertPlot[i_tally][i_filter][i_nuclide][i_score][:],'o-',aa=True)
+                plt.xlabel(xLabel)
+                plt.ylabel("Relative Error of " + yLabel)
+                plt.title("Relative Error of " + title)
+                if (fileType != 'none'):
+                    plt.savefig(REfileName)
+                if showImg:            
+                    plt.show()
+                plt.clf()
+            
+if savetoCSV:
+    # This block loops through each tally, and for each tally:
+        # Creates a new file
+        # Writes the scores and filters for that tally in csv format.
+        # The columns will be: batches,then for each filter: all the scores
+        # The rows, of course, are the data points per batch.
+    
+    for i_tally in range(len(meanPlot)):
+        # Set tally string (placeholder until I put tally labels in statePoint)
+        tallyStr = "Tally " + str(i_tally + 1)
+        CSV_filename = "./tally" + str(i_tally+1)+".csv"
+        # Open the file
+        f = open(CSV_filename, 'w')  
+        
+        # Write the header line        
+        
+        lineText = "Batches" 
+        
+        for i_filter in range(len(meanPlot[i_tally])):
+        
+            # Set filter string
+            filterStr = "Filter " + str(i_filter + 1)
+            
+            for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
+            
+                nuclideStr = "Nuclide " + str(i_nuclide + 1)
+                    
+                for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
+                    
+                    # Set the title
+                    scoreStr = scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score]
+                    scoreStr = scoreStr.title()
+                    if (printxs[i_tally] and ((scoreStr != 'Flux') and \
+                        (scoreStr != 'Current'))):
+                        scoreStr = scoreStr + "-XS"
+                    
+                    # set header 
+                    headerText = scoreStr + " for " + filterStr + " for " + nuclideStr
+                    
+                    lineText = lineText + "," + headerText + \
+                        ",Abs Unc of " + headerText + \
+                        ",Rel Unc of " + headerText
+                    
+        f.write(lineText + "\n")   
+
+        # Write the data lines, each row is a different batch                     
+        
+        for i_batch in range(len(meanPlot[i_tally][0][0][0])):
+        
+            lineText = repr(active_batches[i_batch])        
+        
+            for i_filter in range(len(meanPlot[i_tally])):
+                        
+                for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
+                        
+                    for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
+                        
+                        fieldText = \
+                            repr(meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch]) + \
+                            "," + \
+                            repr(absUncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch]) +\
+                            "," + \
+                            repr(uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch])
+                    
+                        lineText = lineText + "," + fieldText
+            
+            f.write(lineText + "\n")
+    
+   

scripts/track.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python2
+"""Convert binary particle track to VTK poly data.
+
+Usage information can be obtained by running 'track.py --help':
+
+    usage: track.py [-h] [-o OUT] IN [IN ...]
+
+    Convert particle track file to a .pvtp file.
+
+    positional arguments:
+      IN                 Input particle track data filename(s).
+
+    optional arguments:
+      -h, --help         show this help message and exit
+      -o OUT, --out OUT  Output VTK poly data filename.
+
+"""
+
+import os
+import argparse
+import struct
+import vtk
+
+
+def _parse_args():
+    # Create argument parser.
+    parser = argparse.ArgumentParser(
+        description='Convert particle track file to a .pvtp file.')
+    parser.add_argument('input', metavar='IN', type=str, nargs='+',
+                        help='Input particle track data filename(s).')
+    parser.add_argument('-o', '--out', metavar='OUT', type=str, dest='out',
+                        help='Output VTK poly data filename.')
+
+    # Parse and return commandline arguments.
+    return parser.parse_args()
+
+
+def main():
+    # Parse commandline arguments.
+    args = _parse_args()
+
+    # Check input file extensions.
+    for fname in args.input:
+        if not (fname.endswith('.h5') or fname.endswith('.binary')):
+            raise ValueError("Input file names must either end with '.h5' or"
+                             "'.binary'.")
+    
+    # Make sure that the output filename ends with '.pvtp'.
+    if not args.out:
+        args.out = 'tracks.pvtp'
+    elif os.path.splitext(args.out)[1] != '.pvtp':
+        args.out = ''.join([args.out, '.pvtp'])
+
+    # Import HDF library if HDF files are present
+    for fname in args.input:
+        if fname.endswith('.h5'):
+            import h5py
+            break
+
+    # Initialize data arrays and offset.
+    points = vtk.vtkPoints()
+    cells = vtk.vtkCellArray()
+    point_offset = 0
+    for fname in args.input:
+        # Write coordinate values to points array.
+        if fname.endswith('.binary'):
+            track = open(fname, 'rb').read()
+            coords = [struct.unpack("ddd", track[24*i : 24*(i+1)])
+                      for i in range(len(track)/24)]
+            n_points = len(coords)
+            for triplet in coords:
+                points.InsertNextPoint(triplet)
+        else:
+            coords = h5py.File(fname).get('coordinates')
+            n_points = coords.shape[0]
+            for i in range(n_points):
+                points.InsertNextPoint(coords[i,:])
+                
+        # Create VTK line and assign points to line.
+        line = vtk.vtkPolyLine()
+        line.GetPointIds().SetNumberOfIds(n_points)
+        for i in range(n_points):
+            line.GetPointIds().SetId(i, point_offset+i)
+        
+        cells.InsertNextCell(line)
+        point_offset += n_points
+    data = vtk.vtkPolyData()
+    data.SetPoints(points)
+    data.SetLines(cells)
+    
+    writer = vtk.vtkXMLPPolyDataWriter()
+    writer.SetInput(data)
+    writer.SetFileName(args.out)
+    writer.Write()
+
+
+
+if __name__ == '__main__':
+    main()

scripts/update_inputs.py

@@ -0,0 +1,255 @@
+#!/usr/bin/env python
+"""Update OpenMC's input XML files to the latest format.
+
+Usage information can be obtained by running 'update_inputs.py --help':
+
+usage: update_lattices.py [-h] IN [IN ...]
+
+Update lattices in geometry.xml files to the latest format. This will remove
+'outside' attributes/elements and replace them with 'outer' attributes. Note
+that this script will not delete the given files; it will append '.original'
+to the given files and write new ones.
+
+positional arguments:
+  IN          Input geometry.xml file(s).
+
+optional arguments:
+  -h, --help  show this help message and exit
+
+"""
+
+from __future__ import print_function
+
+import argparse
+from itertools import chain
+from random import randint
+from shutil import move
+import xml.etree.ElementTree as ET
+
+
+description = "Update OpenMC's input XML files to the latest format."
+epilog = """\
+If any of the given files do not match the most up-to-date formatting, then they
+will be automatically rewritten.  The old out-of-date files will not be deleted;
+they will be moved to a new file with '.original' appended to their name.
+
+Formatting changes that will be made:
+
+geometry.xml:  Lattices containing 'outside' attributes/tags will be replaced
+  with lattices containing 'outer' attributes, and the appropriate
+  cells/universes will be added.
+"""
+
+
+def parse_args():
+    """Read the input files from the commandline."""
+    # Create argument parser.
+    parser = argparse.ArgumentParser(
+         description=description,
+         epilog=epilog,
+         formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument('input', metavar='IN', type=str, nargs='+',
+                        help='Input geometry.xml file(s).')
+
+    # Parse and return commandline arguments.
+    return parser.parse_args()
+
+
+def get_universe_ids(geometry_root):
+    """Return a set of universe id numbers."""
+    root = geometry_root
+    out = {0}
+
+    # Get the ids of universes defined by cells.
+    for cell in root.iter('cell'):
+        # Get universe attributes.
+        if 'universe' in cell.attrib:
+            uid = cell.attrib['universe']
+            out.add(int(uid))
+
+        # Get universe elements.
+        elif cell.find('universe') is not None:
+            elem = cell.find('universe')
+            uid = elem.text
+            out.add(int(uid))
+
+    # Get the ids of universes defined by lattices.
+    for lat in root.iter('lattice'):
+        # Get id attributes.
+        if 'id' in lat.attrib:
+            uid = lat.attrib['id']
+            out.add(int(uid))
+
+        # Get id elements.
+        elif lat.find('id') is not None:
+            elem = lat.find('id')
+            uid = elem.text
+            out.add(int(uid))
+
+    return out
+
+
+def get_cell_ids(geometry_root):
+    """Return a set of cell id numbers."""
+    root = geometry_root
+    out = set()
+
+    # Get the ids of universes defined by cells.
+    for cell in root.iter('cell'):
+        # Get id attributes.
+        if 'id' in cell.attrib:
+            cid = cell.attrib['id']
+            out.add(int(cid))
+
+        # Get id elements.
+        elif cell.find('id') is not None:
+            elem = cell.find('id')
+            cid = elem.text
+            out.add(int(cid))
+
+    return out
+
+
+def find_new_id(current_ids, preferred=None):
+    """Return a new id that is not already present in current_ids."""
+    distance_from_preferred = 21
+    max_random_attempts = 10000
+
+    # First, try to find an id near the preferred number.
+    if preferred is not None:
+        assert isinstance(preferred, int)
+        for i in range(1, distance_from_preferred):
+            if (preferred - i not in current_ids) and (preferred - i > 0):
+                return preferred - i
+            if (preferred + i not in current_ids) and (preferred + i > 0):
+                return preferred + i
+
+    # If that was unsuccessful, attempt to randomly guess a new id number.
+    for i in range(max_random_attempts):
+        num = randint(1, 2147483647)
+        if num not in current_inds:
+            return num
+
+    # Raise an error if an id was not found.
+    raise RuntimeError('Could not find a unique id number for a new universe.')
+
+
+def get_lat_id(lattice_element):
+    """Return the id integer of the lattice_element."""
+    assert isinstance(lattice_element, ET.Element)
+    if 'id' in lattice_element.attrib:
+        return int(lattice_element.attrib['id'].strip())
+    elif any([child.tag == 'id' for child in lattice_element]):
+        elem = lattice_element.find('id')
+        return int(elem.text.strip())
+    else:
+        raise RuntimeError('Could not find the id for a lattice.')
+
+
+def pop_lat_outside(lattice_element):
+    """Return lattice's outside material and remove from attributes/elements."""
+    assert isinstance(lattice_element, ET.Element)
+
+    # Check attributes.
+    if 'outside' in lattice_element.attrib:
+        material = lattice_element.attrib['outside'].strip()
+        del lattice_element.attrib['outside']
+
+    # Check subelements.
+    elif any([child.tag == 'outside' for child in lattice_element]):
+        elem = lattice_element.find('outside')
+        material = elem.text.strip()
+        lattice_element.remove(elem)
+
+    # No 'outside' specified.  This means the outside is a void.
+    else:
+        material = 'void'
+
+    return material
+
+
+def update_geometry(geometry_root):
+    """Update the given XML geometry tree. Return True if changes were made."""
+    root = geometry_root
+    was_updated = False
+
+    # Ignore files that do not contain lattices.
+    if all([child.tag != 'lattice' for child in root]): return False
+
+    # Get a set of already-used universe and cell ids.
+    uids = get_universe_ids(root)
+    cids = get_cell_ids(root)
+    taken_ids = uids.union(cids)
+
+    # Replace 'outside' with 'outer' in lattices.
+    for lat in chain(root.iter('lattice'), root.iter('hex_lattice')):
+        # Get the lattice's id.
+        lat_id = get_lat_id(lat)
+
+        # Ignore lattices that have 'outer' specified.
+        if any([child.tag == 'outer' for child in lat]): continue
+        if 'outer' in lat.attrib: continue
+
+        # Pop the 'outside' material.
+        material = pop_lat_outside(lat)
+
+        # Get an id number for a new outer universe.  Ideally, the id should
+        # be close to the lattice's id.
+        new_uid = find_new_id(taken_ids, preferred=lat_id)
+        assert new_uid not in taken_ids
+
+        # Add the new universe filled with the old 'outside' material to the
+        # geometry.
+        new_cell = ET.Element('cell')
+        new_cell.attrib['id'] = str(new_uid)
+        new_cell.attrib['universe'] = str(new_uid)
+        new_cell.attrib['material'] = material
+        root.append(new_cell)
+        taken_ids.add(new_uid)
+
+        # Add the new universe to the lattice's 'outer' attribute.
+        lat.attrib['outer'] = str(new_uid)
+
+        was_updated = True
+
+    # Remove 'type' from lattice definitions.
+    for lat in root.iter('lattice'):
+        elem = lat.find('type')
+        if elem is not None:
+            lat.remove(elem)
+            was_updated = True
+        if 'type' in lat.attrib:
+            del lat.attrib['type']
+            was_updated = True
+
+    # Change 'width' to 'pitch' in lattice definitions.
+    for lat in root.iter('lattice'):
+        elem = lat.find('width')
+        if elem is not None:
+            elem.tag = 'pitch'
+            was_updated = True
+        if 'width' in lat.attrib:
+            lat.attrib['pitch'] = lat.attrib['width']
+            del lat.attrib['width']
+            was_updated = True
+
+    return was_updated
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    for fname in args.input:
+        # Parse the XML data.
+        tree = ET.parse(fname)
+        root = tree.getroot()
+        was_updated = False
+
+        if root.tag == 'geometry':
+            was_updated = update_geometry(root)
+
+        if was_updated:
+            # Move the original geometry file to preserve it.
+            move(fname, fname + '.original')
+
+            # Write a new geometry file.
+            tree.write(fname)

scripts/voxel.py

@@ -0,0 +1,120 @@
+#!/usr/bin/env python2
+
+from __future__ import division, print_function
+
+import struct
+import sys
+
+################################################################################
+def parse_options():
+  """Process command line arguments"""
+
+  from optparse import OptionParser
+  usage = r"""%prog [options] <voxel_file>"""
+  p = OptionParser(usage=usage)
+  p.add_option('-o', '--output', action='store', dest='output',
+             default='plot', help='Path to output SILO or VTK file.')
+  p.add_option('-v', '--vtk', action='store_true', dest='vtk',
+           default=False, help='Flag to convert to VTK instead of SILO.')
+  parsed = p.parse_args()
+  if not parsed[1]:
+    p.print_help()
+    return parsed
+  return parsed
+
+################################################################################
+def main(file_, o):
+
+  print(file_)
+  fh = open(file_,'rb')
+  header = get_header(fh)
+  meshparms = header['dimension'] + header['lower_left'] + header['upper_right']
+  nx,ny,nz = meshparms[0], meshparms[1], meshparms[2]
+  ll = header['lower_left']
+
+  if o.vtk:
+    try:
+      import vtk
+    except:
+      print('The vtk python bindings do not appear to be installed properly.\n'
+            'On Ubuntu: sudo apt-get install python-vtk\n'
+            'See: http://www.vtk.org/')
+      return
+
+    origin = [(l+w*n/2.) for n,l,w in zip((nx,ny,nz),ll,header['width'])]
+
+    grid = vtk.vtkImageData()
+    grid.SetDimensions(nx+1,ny+1,nz+1)
+    grid.SetOrigin(*ll)
+    grid.SetSpacing(*header['width'])
+
+    data = vtk.vtkDoubleArray()
+    data.SetName("id")
+    data.SetNumberOfTuples(nx*ny*nz)
+    for x in range(nx):
+      sys.stdout.write(" {0}%\r".format(int(x/nx*100)))
+      sys.stdout.flush()
+      for y in range(ny):
+        for z in range(nz):
+          i = z*nx*ny + y*nx + x
+          id_ = get_int(fh)[0]
+          data.SetValue(i, id_)
+    grid.GetCellData().AddArray(data)
+
+    writer = vtk.vtkXMLImageDataWriter()
+    writer.SetInput(grid)
+    if not o.output[-4:] == ".vti": o.output += ".vti"
+    writer.SetFileName(o.output)
+    writer.Write()
+
+  else:
+
+    try:
+      import silomesh
+    except:
+      print('The silomesh package does not appear to be installed properly.\n'
+            'See: https://github.com/nhorelik/silomesh/')
+      return
+    if not o.output[-5:] == ".silo": o.output += ".silo"
+    silomesh.init_silo(o.output)
+    silomesh.init_mesh('plot', *meshparms)
+    silomesh.init_var("id")
+    for x in range(1,nx+1):
+      sys.stdout.write(" {0}%\r".format(int(x/nx*100)))
+      sys.stdout.flush()
+      for y in range(1,ny+1):
+        for z in range(1,nz+1):
+          id_ = get_int(fh)[0]
+          silomesh.set_value(float(id_), x, y, z)
+    print()
+    silomesh.finalize_var()
+    silomesh.finalize_mesh()
+    silomesh.finalize_silo()
+
+################################################################################
+def get_header(file_):
+  nx,ny,nz = get_int(file_, 3)
+  wx,wy,wz = get_double(file_, 3)
+  lx,ly,lz = get_double(file_, 3)
+  header = {'dimension':[nx,ny,nz], 'width':[wx,wy,wz], 'lower_left':[lx,ly,lz],
+            'upper_right': [lx+wx*nx,ly+wy*ny,lz+wz*nz]}
+  return header
+
+################################################################################
+def get_data(file_, n, typeCode, size):
+  return list(struct.unpack('={0}{1}'.format(n,typeCode),
+                            file_.read(n*size)))
+
+################################################################################
+def get_int(file_, n=1, path=None):
+  return get_data(file_, n, 'i', 4)
+
+################################################################################
+def get_double(file_, n=1, path=None):
+  return get_data(file_, n, 'd', 8)
+
+################################################################################
+if __name__ == '__main__':
+  (options, args) = parse_options()
+  if args:
+    main(args[0],options)

scripts/xml_validate.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+
+from __future__ import print_function
+
+import os
+import sys
+import glob
+import lxml.etree as etree
+from subprocess import call
+from optparse import OptionParser
+
+# Command line parsing
+parser = OptionParser()
+parser.add_option('-r', '--relaxng-path', dest='relaxng',
+                  help="Path to RelaxNG files.")
+parser.add_option('-i', '--input-path', dest='inputs', default=os.getcwd(),
+                  help="Path to OpenMC input files." )
+(options, args) = parser.parse_args()
+
+# Colored output
+if sys.stdout.isatty():
+    OK = '\033[92m'
+    FAIL = '\033[91m'
+    NOT_FOUND = '\033[93m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+else:
+    OK = ''
+    FAIL = ''
+    ENDC = ''
+    BOLD = ''
+    NOT_FOUND = ''
+
+# Get absolute paths
+if options.relaxng is not None:
+    relaxng_path = os.path.abspath(options.relaxng)
+if options.inputs is not None:
+    inputs_path = os.path.abspath(options.inputs)
+
+# Search for relaxng path if not set
+if options.relaxng is None:
+    xml_validate_path = os.path.abspath(os.path.dirname(sys.argv[0]))
+    if "bin" in xml_validate_path:
+        relaxng_path = os.path.join(xml_validate_path, "..", "share", "relaxng")
+    elif os.path.join("src", "utils") in xml_validate_path:
+        relaxng_path = os.path.join(xml_validate_path, "..", "relaxng")
+    else:
+        raise Exception("Set RelaxNG path with -r command line option.")
+if not os.path.exists(relaxng_path):
+    raise Exception("RelaxNG path: {0} does not exist, set with -r "
+                    "command line option.".format(relaxng_path))
+
+# Make sure there are .rng files in RelaxNG path
+rng_files = glob.glob(os.path.join(relaxng_path, "*.rng"))
+if len(rng_files) == 0:
+    raise Exception("No .rng files found in RelaxNG "
+                    "path: {0}.".format(relaxng_path))
+
+# Get list of xml input files
+xml_files = glob.glob(os.path.join(inputs_path, "*.xml"))
+if len(xml_files) == 0:
+    raise Exception("No .xml files found at input path: {0}"
+                    ".".format(inputs_path))
+
+# Begin loop around input files
+for xml_file in xml_files:
+
+    text = "Validating {0}".format(os.path.basename(xml_file))
+    print(text + '.'*(30 - len(text)), end="")
+
+    # Validate the XML file
+    try:
+        xml_tree = etree.parse(xml_file)
+    except etree.XMLSyntaxError as e:
+        print(BOLD + FAIL + '[XML ERROR]' + ENDC)
+        print("    {0}".format(e))
+        continue
+
+    # Get xml_filename prefix
+    xml_prefix = os.path.basename(xml_file)
+    xml_prefix = xml_prefix.split(".")[0]
+
+    # Search for rng file
+    rng_file = os.path.join(relaxng_path, xml_prefix + ".rng")
+    if rng_file in rng_files:
+
+        # read in RelaxNG
+        relaxng_doc = etree.parse(rng_file)
+        relaxng = etree.RelaxNG(relaxng_doc)
+
+        # validate xml file again RelaxNG
+        try:
+            relaxng.assertValid(xml_tree)
+            print(BOLD + OK + '[VALID]' + ENDC)
+        except (etree.DocumentInvalid, TypeError) as e:
+            print(BOLD + FAIL + '[NOT VALID]' + ENDC)
+            print("    {0}".format(e))
+
+    # RNG file does not exist
+    else:
+        print(BOLD + NOT_FOUND + '[NO RELAXNG FOUND]' + ENDC)