diff --git a/smr/fresh/geometry.xml b/smr/fresh/geometry.xml
index a3e377f..2434ba1 100644
--- a/smr/fresh/geometry.xml
+++ b/smr/fresh/geometry.xml
@@ -6251,7 +6251,7 @@
   <cell fill="91" id="6642" name="Fuel (3.1%) stack (10)" region="60 -61" universe="106" />
   <cell fill="1" id="6643" name="Fuel (3.1%) stack (11)" region="61" universe="106" />
   <cell fill="125" id="6761" name="Assembly (1.6%) no BAs instr lattice" region="33 -34 35 -36" universe="126" />
-  <cell fill="1" id="6762" name="Assembly (1.6%) no BAs instr lattice outer water" region="~(37 -38 39 -40)" universe="126" />
+  <cell id="6762" material="8" name="Assembly (1.6%) no BAs instr lattice outer water" region="~(37 -38 39 -40)" universe="126" />
   <cell id="6763" material="8" name="Assembly (1.6%) no BAs instr lattice axial (0)" region="37 -38 39 -40 ~(33 -34 35 -36) -47" universe="126" />
   <cell id="6764" material="3" name="Assembly (1.6%) no BAs instr lattice axial (1)" region="37 -38 39 -40 ~(33 -34 35 -36) 47 -48" universe="126" />
   <cell id="6765" material="8" name="Assembly (1.6%) no BAs instr lattice axial (2)" region="37 -38 39 -40 ~(33 -34 35 -36) 48 -49" universe="126" />
@@ -6264,7 +6264,7 @@
   <cell id="6772" material="6" name="Assembly (1.6%) no BAs instr lattice axial (9)" region="37 -38 39 -40 ~(33 -34 35 -36) 55 -56" universe="126" />
   <cell id="6773" material="8" name="Assembly (1.6%) no BAs instr lattice axial (last)" region="37 -38 39 -40 ~(33 -34 35 -36) 56" universe="126" />
   <cell fill="145" id="6891" name="Assembly (2.4%) CR D lattice" region="33 -34 35 -36" universe="146" />
-  <cell fill="1" id="6892" name="Assembly (2.4%) CR D lattice outer water" region="~(37 -38 39 -40)" universe="146" />
+  <cell id="6892" material="8" name="Assembly (2.4%) CR D lattice outer water" region="~(37 -38 39 -40)" universe="146" />
   <cell id="6893" material="8" name="Assembly (2.4%) CR D lattice axial (0)" region="37 -38 39 -40 ~(33 -34 35 -36) -47" universe="146" />
   <cell id="6894" material="3" name="Assembly (2.4%) CR D lattice axial (1)" region="37 -38 39 -40 ~(33 -34 35 -36) 47 -48" universe="146" />
   <cell id="6895" material="8" name="Assembly (2.4%) CR D lattice axial (2)" region="37 -38 39 -40 ~(33 -34 35 -36) 48 -49" universe="146" />
@@ -6277,7 +6277,7 @@
   <cell id="6902" material="6" name="Assembly (2.4%) CR D lattice axial (9)" region="37 -38 39 -40 ~(33 -34 35 -36) 55 -56" universe="146" />
   <cell id="6903" material="8" name="Assembly (2.4%) CR D lattice axial (last)" region="37 -38 39 -40 ~(33 -34 35 -36) 56" universe="146" />
   <cell fill="165" id="7021" name="Assembly (3.1%) 16BA lattice" region="33 -34 35 -36" universe="166" />
-  <cell fill="1" id="7022" name="Assembly (3.1%) 16BA lattice outer water" region="~(37 -38 39 -40)" universe="166" />
+  <cell id="7022" material="8" name="Assembly (3.1%) 16BA lattice outer water" region="~(37 -38 39 -40)" universe="166" />
   <cell id="7023" material="8" name="Assembly (3.1%) 16BA lattice axial (0)" region="37 -38 39 -40 ~(33 -34 35 -36) -47" universe="166" />
   <cell id="7024" material="3" name="Assembly (3.1%) 16BA lattice axial (1)" region="37 -38 39 -40 ~(33 -34 35 -36) 47 -48" universe="166" />
   <cell id="7025" material="8" name="Assembly (3.1%) 16BA lattice axial (2)" region="37 -38 39 -40 ~(33 -34 35 -36) 48 -49" universe="166" />
@@ -6290,7 +6290,7 @@
   <cell id="7032" material="6" name="Assembly (3.1%) 16BA lattice axial (9)" region="37 -38 39 -40 ~(33 -34 35 -36) 55 -56" universe="166" />
   <cell id="7033" material="8" name="Assembly (3.1%) 16BA lattice axial (last)" region="37 -38 39 -40 ~(33 -34 35 -36) 56" universe="166" />
   <cell fill="183" id="7138" name="Assembly (3.1%) no BAs instr lattice" region="33 -34 35 -36" universe="184" />
-  <cell fill="1" id="7139" name="Assembly (3.1%) no BAs instr lattice outer water" region="~(37 -38 39 -40)" universe="184" />
+  <cell id="7139" material="8" name="Assembly (3.1%) no BAs instr lattice outer water" region="~(37 -38 39 -40)" universe="184" />
   <cell id="7140" material="8" name="Assembly (3.1%) no BAs instr lattice axial (0)" region="37 -38 39 -40 ~(33 -34 35 -36) -47" universe="184" />
   <cell id="7141" material="3" name="Assembly (3.1%) no BAs instr lattice axial (1)" region="37 -38 39 -40 ~(33 -34 35 -36) 47 -48" universe="184" />
   <cell id="7142" material="8" name="Assembly (3.1%) no BAs instr lattice axial (2)" region="37 -38 39 -40 ~(33 -34 35 -36) 48 -49" universe="184" />
diff --git a/smr/smr/assemblies.py b/smr/smr/assemblies.py
index 3d9c996..d673127 100644
--- a/smr/smr/assemblies.py
+++ b/smr/smr/assemblies.py
@@ -6,7 +6,7 @@ import openmc
 
 from .materials import mats
 from .surfaces import surfs, pin_pitch
-from .pins import univs as pins
+from .pins import pin_universes
 
 
 def make_assembly(name, universes):
@@ -47,7 +47,7 @@ def make_assembly(name, universes):
 
     # Add outer water cell
     cell = openmc.Cell(name=univ_name + ' outer water')
-    cell.fill = pins['water pin']
+    cell.fill = mats['H2O']
     cell.region = ~surfs['lat grid box outer']
     universe.add_cell(cell)
 
@@ -121,6 +121,8 @@ def make_assembly(name, universes):
 
 
 def assembly_universes():
+    pins = pin_universes()
+
     # Create dictionary to store assembly universes
     univs = {}
 
diff --git a/smr/smr/pins.py b/smr/smr/pins.py
index f7c6119..01ae73d 100644
--- a/smr/smr/pins.py
+++ b/smr/smr/pins.py
@@ -140,737 +140,742 @@ def make_pin_stack(name, zsurfaces, universes, boundary, pin_universe):
     return universe
 
 
-univs = {}
+def pin_universes():
+    # Create dictionary to store pin universes
+    univs = {}
 
+    # Dummy water cell
+    cell = openmc.Cell(name='water pin', fill=mats['H2O'])
+    univs['water pin'] = openmc.Universe(name='Empty water pin cell universe')
+    univs['water pin'].add_cell(cell)
 
-#### DUMMY WATER CELL
+    # Guide tube pin cells
+    univs['GT empty'] = make_pin(
+        'GT empty',
+        [surfs['GT IR'], surfs['GT OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']])
+    univs['GT empty grid (bottom)'] = make_pin(
+        'GT empty grid (bottom)',
+        [surfs['GT IR'], surfs['GT OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']], grid='bottom')
+    univs['GT empty grid (intermediate)'] = make_pin(
+        'GT empty grid (intermediate)',
+        [surfs['GT IR'], surfs['GT OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']], grid='intermediate')
+    univs['GT empty nozzle'] = make_pin(
+        'GT empty nozzle',
+        [surfs['GT IR'], surfs['GT OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']])
 
-cell = openmc.Cell(name='water pin')
-cell.fill = mats['H2O']
-univs['water pin'] = openmc.Universe(name='Empty water pin cell universe')
-univs['water pin'].add_cell(cell)
+    univs['GTd empty'] = make_pin(
+        'GT empty at dashpot',
+        [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']])
+    univs['GTd empty grid (bottom)'] = make_pin(
+        'GT empty at dashpot grid (bottom)',
+        [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']], grid='bottom')
+    univs['GTd empty grid (intermediate)'] = make_pin(
+        'GT empty at dashpot grid (intermediate)',
+        [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']], grid='intermediate')
+    univs['GTd empty nozzle'] = make_pin(
+        'GT empty nozzle',
+        [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
+        [mats['H2O'], mats['Zr'], mats['H2O']])
 
+    # Stack all axial pieces of guide tube together
 
-#### GUIDE TUBE PIN CELLS
+    stack_surfs = [
+        surfs['bot support plate'],
+        surfs['top support plate'],
+        surfs['top lower nozzle'],
+        surfs['top lower thimble'],
+        surfs['grid1bot'],
+        surfs['grid1top'],
+        surfs['dashpot top'],
+        surfs['grid2bot'],
+        surfs['grid2top'],
+        surfs['grid3bot'],
+        surfs['grid3top'],
+        surfs['grid4bot'],
+        surfs['grid4top'],
+        surfs['top active core'],
+        surfs['grid5bot'],
+        surfs['grid5top'],
+        surfs['top pin plenum'],
+        surfs['top FR'],
+        surfs['bot upper nozzle'],
+        surfs['top upper nozzle']
+    ]
 
-univs['GT empty'] = make_pin(
-    'GT empty', [surfs['GT IR'], surfs['GT OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']])
-univs['GT empty grid (bottom)'] = make_pin(
-    'GT empty grid (bottom)', [surfs['GT IR'], surfs['GT OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']], grid='bottom')
-univs['GT empty grid (intermediate)'] = make_pin(
-    'GT empty grid (intermediate)', [surfs['GT IR'], surfs['GT OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']], grid='intermediate')
-univs['GT empty nozzle'] = make_pin(
-    'GT empty nozzle', [surfs['GT IR'], surfs['GT OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']])
-
-univs['GTd empty'] = make_pin(
-    'GT empty at dashpot', [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']])
-univs['GTd empty grid (bottom)'] = make_pin(
-    'GT empty at dashpot grid (bottom)',
-    [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']], grid='bottom')
-univs['GTd empty grid (intermediate)'] = make_pin(
-    'GT empty at dashpot grid (intermediate)',
-    [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']], grid='intermediate')
-univs['GTd empty nozzle'] = make_pin(
-    'GT empty nozzle', [surfs['GT dashpot IR'], surfs['GT dashpot OR']],
-    [mats['H2O'], mats['Zr'], mats['H2O']])
-
-# Stack all axial pieces of guide tube together
-
-stack_surfs = [
-    surfs['bot support plate'],
-    surfs['top support plate'],
-    surfs['top lower nozzle'],
-    surfs['top lower thimble'],
-    surfs['grid1bot'],
-    surfs['grid1top'],
-    surfs['dashpot top'],
-    surfs['grid2bot'],
-    surfs['grid2top'],
-    surfs['grid3bot'],
-    surfs['grid3top'],
-    surfs['grid4bot'],
-    surfs['grid4top'],
-    surfs['top active core'],
-    surfs['grid5bot'],
-    surfs['grid5top'],
-    surfs['top pin plenum'],
-    surfs['top FR'],
-    surfs['bot upper nozzle'],
-    surfs['top upper nozzle']
-]
-
-univs['GT empty'] = make_stack(
-    'GT empty', surfaces=stack_surfs,
-    universes=[univs['water pin'],
-               univs['water pin'],
-               univs['water pin'],
-               univs['GTd empty'],
-               univs['GTd empty'],
-               univs['GTd empty grid (bottom)'],
-               univs['GTd empty'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['water pin'],
-               univs['water pin']])
-
-univs['GT empty instr'] = make_stack(
-    'GT empty instr', surfaces=stack_surfs,
-    universes=[univs['water pin'],
-               univs['water pin'],
-               univs['water pin'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty grid (bottom)'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty grid (intermediate)'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['GT empty'],
-               univs['water pin'],
-               univs['water pin']])
-
-
-#### INSTRUMENT TUBE PIN CELL
-
-univs['IT'] = make_pin(
-    'IT', [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']])
-univs['IT grid (bottom)'] = make_pin(
-    'IT grid (bottom)',
-    [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='bottom')
-univs['IT grid (intermediate)'] = make_pin(
-    'IT grid (intermediate)',
-    [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='intermediate')
-
-univs['IT nozzle'] = make_pin(
-    'IT nozzle',
-    [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['SS']])
-univs['IT dashpot'] = make_pin(
-    'IT dashpot', [surfs['IT IR'], surfs['IT OR']],
-    [mats['Air'], mats['Zr'], mats['H2O']])
-
-# Stack all axial pieces of instrument tube together
-
-univs['IT stack'] = make_stack(
-    'GT instr', surfaces=stack_surfs,
-    universes=[univs['IT dashpot'],
-               univs['IT dashpot'],
-               univs['IT dashpot'],
-               univs['IT'],
-               univs['IT'],
-               univs['IT grid (bottom)'],
-               univs['IT'],
-               univs['IT'],
-               univs['IT grid (intermediate)'],
-               univs['IT'],
-               univs['IT grid (intermediate)'],
-               univs['IT'],
-               univs['IT grid (intermediate)'],
-               univs['IT'],
-               univs['IT'],
-               univs['IT grid (intermediate)'],
-               univs['IT'],
-               univs['IT'],
-               univs['IT'],
-               univs['IT dashpot'],
-               univs['water pin']])
-
-
-#### CONTROL ROD PIN CELLS
-
-univs['CR'] = make_pin(
-    'CR', [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']])
-univs['CR grid (bottom)'] = make_pin(
-    'CR grid (bottom)', [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='bottom')
-univs['CR grid (intermediate)'] = make_pin(
-    'CR grid (intermediate)',
-    [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='intermediate')
-univs['CR nozzle'] = make_pin(
-    'CR nozzle', [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
-    [mats['AIC'], mats['Air'], mats['SS'], mats['H2O']])
-
-univs['CR blank'] = make_pin(
-    'CR blank',
-    [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']])
-univs['CR blank grid (bottom)'] = make_pin(
-    'CR blank grid (bottom)',
-    [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='bottom')
-univs['CR blank grid (intermediate)'] = make_pin(
-    'CR blank grid (intermediate)',
-    [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
-    [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
-    grid='intermediate')
-univs['CR blank nozzle'] = make_pin(
-    'CR blank nozzle', [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
-    [mats['SS'], mats['Air'], mats['SS'], mats['H2O']])
-univs['CR blank bare'] = make_pin(
-    'CR blank bare', [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
-    [mats['SS'], mats['Air'], mats['SS'], mats['H2O']])
-univs['CR bare'] = make_pin(
-    'CR bare', [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
-    [mats['AIC'], mats['Air'], mats['SS'], mats['H2O']])
-
-# Stack all axial pieces of control rod tubes together for each bank
-
-banks = ['A', 'B', 'C', 'D', 'SA', 'SB', 'SC', 'SD', 'SE']
-
-for b in banks:
-
-    # no grid, no nozzle
-    univs['GT CR bank {} dummy'.format(b)] = make_stack(
-        'GT CR bank {} dummy'.format(b),
-        surfaces=[surfs['bottom FR'],
-                  surfs['dashpot top'],
-                  surfs['bank{} bot'.format(b)],
-                  surfs['bank{} top'.format(b)]],
+    univs['GT empty'] = make_stack(
+        'GT empty', surfaces=stack_surfs,
         universes=[univs['water pin'],
+                   univs['water pin'],
+                   univs['water pin'],
+                   univs['GTd empty'],
+                   univs['GTd empty'],
+                   univs['GTd empty grid (bottom)'],
                    univs['GTd empty'],
                    univs['GT empty'],
-                   univs['CR'],
-                   univs['CR blank']])
-
-    # bottom grid
-    univs['GT CR bank {} dummy grid (bottom)'.format(b)] = make_stack(
-        'GT CR bank {} dummy grid (bottom)'.format(b),
-        surfaces=[surfs['bottom FR'],
-                  surfs['dashpot top'],
-                  surfs['bank{} bot'.format(b)],
-                  surfs['bank{} top'.format(b)]],
-        universes=[univs['water pin'],
-                   univs['GTd empty grid (bottom)'],
-                   univs['GT empty grid (bottom)'],
-                   univs['CR grid (bottom)'],
-                   univs['CR blank grid (bottom)']])
-
-    # intermediate grid
-    univs['GT CR bank {} dummy grid (intermediate)'.format(b)] = make_stack(
-        'GT CR bank {} dummy grid (intermediate)'.format(b),
-        surfaces=[surfs['bottom FR'],
-                  surfs['dashpot top'],
-                  surfs['bank{} bot'.format(b)],
-                  surfs['bank{} top'.format(b)]],
-        universes=[univs['water pin'],
-                   univs['GTd empty grid (intermediate)'],
                    univs['GT empty grid (intermediate)'],
-                   univs['CR grid (intermediate)'],
-                   univs['CR blank grid (intermediate)']])
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['water pin'],
+                   univs['water pin']])
 
-    # nozzle
-    univs['GT CR bank {} dummy nozzle'.format(b)] = make_stack(
-        'GT CR bank {} dummy nozzle'.format(b),
-        surfaces=[surfs['bottom FR'],
-                  surfs['dashpot top'],
-                  surfs['bank{} bot'.format(b)],
-                  surfs['bank{} top'.format(b)]],
+    univs['GT empty instr'] = make_stack(
+        'GT empty instr', surfaces=stack_surfs,
         universes=[univs['water pin'],
-                   univs['GTd empty nozzle'],
-                   univs['GT empty nozzle'],
-                   univs['CR nozzle'],
-                   univs['CR blank nozzle']])
+                   univs['water pin'],
+                   univs['water pin'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty grid (bottom)'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty grid (intermediate)'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['GT empty'],
+                   univs['water pin'],
+                   univs['water pin']])
 
-    # bare
-    univs['GT CR bank {} dummy bare'.format(b)] = make_stack(
-        'GT CR bank {} dummy bare'.format(b),
-        surfaces=[surfs['bottom FR'],
-                  surfs['dashpot top'],
-                  surfs['bank{} bot'.format(b)],
-                  surfs['bank{} top'.format(b)]],
+
+    # Instrument tube pin cell
+    univs['IT'] = make_pin(
+        'IT',
+        [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']])
+    univs['IT grid (bottom)'] = make_pin(
+        'IT grid (bottom)',
+        [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='bottom')
+    univs['IT grid (intermediate)'] = make_pin(
+        'IT grid (intermediate)',
+        [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='intermediate')
+
+    univs['IT nozzle'] = make_pin(
+        'IT nozzle',
+        [surfs['IT IR'], surfs['IT OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['Air'], mats['Zr'], mats['H2O'], mats['Zr'], mats['SS']])
+    univs['IT dashpot'] = make_pin(
+        'IT dashpot',
+        [surfs['IT IR'], surfs['IT OR']],
+        [mats['Air'], mats['Zr'], mats['H2O']])
+
+    # Stack all axial pieces of instrument tube together
+
+    univs['IT stack'] = make_stack(
+        'GT instr',
+        surfaces=stack_surfs,
+        universes=[univs['IT dashpot'],
+                   univs['IT dashpot'],
+                   univs['IT dashpot'],
+                   univs['IT'],
+                   univs['IT'],
+                   univs['IT grid (bottom)'],
+                   univs['IT'],
+                   univs['IT'],
+                   univs['IT grid (intermediate)'],
+                   univs['IT'],
+                   univs['IT grid (intermediate)'],
+                   univs['IT'],
+                   univs['IT grid (intermediate)'],
+                   univs['IT'],
+                   univs['IT'],
+                   univs['IT grid (intermediate)'],
+                   univs['IT'],
+                   univs['IT'],
+                   univs['IT'],
+                   univs['IT dashpot'],
+                   univs['water pin']])
+
+    # Control rod pin cells
+    univs['CR'] = make_pin(
+        'CR',
+        [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']])
+    univs['CR grid (bottom)'] = make_pin(
+        'CR grid (bottom)',
+        [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='bottom')
+    univs['CR grid (intermediate)'] = make_pin(
+        'CR grid (intermediate)',
+        [surfs['CP OR'], surfs['CR IR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['AIC'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='intermediate')
+    univs['CR nozzle'] = make_pin(
+        'CR nozzle',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
+        [mats['AIC'], mats['Air'], mats['SS'], mats['H2O']])
+
+    univs['CR blank'] = make_pin(
+        'CR blank',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']])
+    univs['CR blank grid (bottom)'] = make_pin(
+        'CR blank grid (bottom)',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='bottom')
+    univs['CR blank grid (intermediate)'] = make_pin(
+        'CR blank grid (intermediate)',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR'], surfs['GT IR'], surfs['GT OR']],
+        [mats['SS'], mats['Air'], mats['SS'], mats['H2O'], mats['Zr'], mats['H2O']],
+        grid='intermediate')
+    univs['CR blank nozzle'] = make_pin(
+        'CR blank nozzle',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
+        [mats['SS'], mats['Air'], mats['SS'], mats['H2O']])
+    univs['CR blank bare'] = make_pin(
+        'CR blank bare',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
+        [mats['SS'], mats['Air'], mats['SS'], mats['H2O']])
+    univs['CR bare'] = make_pin(
+        'CR bare',
+        [surfs['CP OR'], surfs['CR IR'], surfs['CR OR']],
+        [mats['AIC'], mats['Air'], mats['SS'], mats['H2O']])
+
+    # Stack all axial pieces of control rod tubes together for each bank
+
+    for b in ['A', 'B', 'C', 'D', 'SA', 'SB', 'SC', 'SD', 'SE']:
+        # no grid, no nozzle
+        univs['GT CR bank {} dummy'.format(b)] = make_stack(
+            'GT CR bank {} dummy'.format(b),
+            surfaces=[surfs['bottom FR'],
+                      surfs['dashpot top'],
+                      surfs['bank{} bot'.format(b)],
+                      surfs['bank{} top'.format(b)]],
+            universes=[univs['water pin'],
+                       univs['GTd empty'],
+                       univs['GT empty'],
+                       univs['CR'],
+                       univs['CR blank']])
+
+        # bottom grid
+        univs['GT CR bank {} dummy grid (bottom)'.format(b)] = make_stack(
+            'GT CR bank {} dummy grid (bottom)'.format(b),
+            surfaces=[surfs['bottom FR'],
+                      surfs['dashpot top'],
+                      surfs['bank{} bot'.format(b)],
+                      surfs['bank{} top'.format(b)]],
+            universes=[univs['water pin'],
+                       univs['GTd empty grid (bottom)'],
+                       univs['GT empty grid (bottom)'],
+                       univs['CR grid (bottom)'],
+                       univs['CR blank grid (bottom)']])
+
+        # intermediate grid
+        univs['GT CR bank {} dummy grid (intermediate)'.format(b)] = make_stack(
+            'GT CR bank {} dummy grid (intermediate)'.format(b),
+            surfaces=[surfs['bottom FR'],
+                      surfs['dashpot top'],
+                      surfs['bank{} bot'.format(b)],
+                      surfs['bank{} top'.format(b)]],
+            universes=[univs['water pin'],
+                       univs['GTd empty grid (intermediate)'],
+                       univs['GT empty grid (intermediate)'],
+                       univs['CR grid (intermediate)'],
+                       univs['CR blank grid (intermediate)']])
+
+        # nozzle
+        univs['GT CR bank {} dummy nozzle'.format(b)] = make_stack(
+            'GT CR bank {} dummy nozzle'.format(b),
+            surfaces=[surfs['bottom FR'],
+                      surfs['dashpot top'],
+                      surfs['bank{} bot'.format(b)],
+                      surfs['bank{} top'.format(b)]],
+            universes=[univs['water pin'],
+                       univs['GTd empty nozzle'],
+                       univs['GT empty nozzle'],
+                       univs['CR nozzle'],
+                       univs['CR blank nozzle']])
+
+        # bare
+        univs['GT CR bank {} dummy bare'.format(b)] = make_stack(
+            'GT CR bank {} dummy bare'.format(b),
+            surfaces=[surfs['bottom FR'],
+                      surfs['dashpot top'],
+                      surfs['bank{} bot'.format(b)],
+                      surfs['bank{} top'.format(b)]],
+            universes=[univs['water pin'],
+                       univs['GTd empty nozzle'],
+                       univs['GT empty nozzle'],
+                       univs['CR bare'],
+                       univs['CR blank bare']])
+
+        # final combination of all axial pieces for control rod bank "b"
+        univs['GT CR bank {}'.format(b)] = make_stack(
+            'GT CR bank {}'.format(b), stack_surfs,
+            universes=[univs['water pin'],
+                       univs['GT CR bank {} dummy nozzle'.format(b)],
+                       univs['GT CR bank {} dummy nozzle'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy grid (bottom)'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy grid (bottom)'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy'.format(b)],
+                       univs['GT CR bank {} dummy bare'.format(b)],
+                       univs['GT CR bank {} dummy bare'.format(b)]])
+
+
+    #### BURNABLE ABSORBER PIN CELLS
+
+    univs['BA'] = make_pin(
+        'BA',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['BA IR 7'],
+                  surfs['BA IR 8']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']])
+
+    univs['BA grid (bottom)'] = make_pin(
+        'BA grid (bottom)',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['BA IR 7'],
+                  surfs['BA IR 8']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']],
+        grid='bottom')
+
+    univs['BA grid (intermediate)'] = make_pin(
+        'BA grid (intermediate)',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['BA IR 7'],
+                  surfs['BA IR 8']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']],
+        grid='intermediate')
+
+    univs['BA dashpot'] = make_pin(
+        'BA dashpot',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['GT dashpot IR'],
+                  surfs['GT dashpot OR']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']])
+
+    univs['BA dashpot grid (bottom)'] = make_pin(
+        'BA dashpot grid (bottom)',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['GT dashpot IR'],
+                  surfs['GT dashpot OR']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']],
+        grid='bottom')
+
+    univs['BA dashpot grid (intermediate)'] = make_pin(
+        'BA dashpot grid (intermediate)',
+        surfaces=[surfs['BA IR 1'],
+                  surfs['BA IR 2'],
+                  surfs['BA IR 3'],
+                  surfs['BA IR 4'],
+                  surfs['BA IR 5'],
+                  surfs['BA IR 6'],
+                  surfs['GT dashpot IR'],
+                  surfs['GT dashpot OR']],
+        materials=[mats['Air'],
+                   mats['SS'],
+                   mats['Air'],
+                   mats['BSG'],
+                   mats['Air'],
+                   mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']],
+        grid='intermediate')
+
+    univs['BA blank SS'] = make_pin(
+        'BA blank SS',
+        surfaces=[surfs['BA IR 6'],
+                  surfs['BA IR 7'],
+                  surfs['BA IR 8']],
+        materials=[mats['SS'],
+                   mats['H2O'],
+                   mats['Zr'],
+                   mats['H2O']])
+
+    univs['BA blank SS bare'] = make_pin(
+        'BA blank SS bare',
+        surfaces=[surfs['BA IR 6']],
+        materials=[mats['SS'],
+                   mats['H2O']])
+
+    stack_surfs_BA = [
+        surfs['bot support plate'],
+        surfs['top support plate'],
+        surfs['top lower nozzle'],
+        surfs['top lower thimble'],
+        surfs['grid1bot'],
+        surfs['BA bot'],
+        surfs['grid1top'],
+        surfs['dashpot top'],
+        surfs['grid2bot'],
+        surfs['grid2top'],
+        surfs['grid3bot'],
+        surfs['grid3top'],
+        surfs['grid4bot'],
+        surfs['grid4top'],
+        surfs['top active core'],
+        surfs['grid5bot'],
+        surfs['grid5top'],
+        surfs['top pin plenum'],
+        surfs['top FR'],
+        surfs['bot upper nozzle'],
+        surfs['top upper nozzle']]
+
+    # Stack all axial pieces of control rod tubes together for each bank
+
+    univs['BA stack'] = make_stack(
+        'BA stack', stack_surfs_BA,
         universes=[univs['water pin'],
-                   univs['GTd empty nozzle'],
-                   univs['GT empty nozzle'],
-                   univs['CR bare'],
-                   univs['CR blank bare']])
+                   univs['water pin'],
+                   univs['water pin'],
+                   univs['GTd empty'],
+                   univs['GTd empty'],
+                   univs['GTd empty grid (bottom)'],
+                   univs['BA dashpot grid (bottom)'],
+                   univs['BA dashpot'],
+                   univs['BA'],
+                   univs['BA grid (intermediate)'],
+                   univs['BA'],
+                   univs['BA grid (intermediate)'],
+                   univs['BA'],
+                   univs['BA grid (intermediate)'],
+                   univs['BA'],
+                   univs['BA blank SS'],
+                   univs['BA blank SS'],
+                   univs['BA blank SS'],
+                   univs['BA blank SS'],
+                   univs['BA blank SS'],
+                   univs['BA blank SS bare'],
+                   univs['water pin']])
 
-    # final combination of all axial pieces for control rod bank "b"
-    univs['GT CR bank {}'.format(b)] = make_stack(
-        'GT CR bank {}'.format(b), stack_surfs,
+
+    # Fuel pin cells
+    univs['SS pin'] = make_pin(
+        'SS pin',
+        [surfs['clad OR']],
+        [mats['SS'], mats['H2O']])
+
+    univs['end plug'] = make_pin(
+        'end plug',
+        [surfs['clad OR']],
+        [mats['Zr'], mats['H2O']])
+
+    univs['pin plenum'] = make_pin(
+        'pin plenum',
+        surfaces=[surfs['plenum spring OR'],
+                  surfs['clad IR'],
+                  surfs['clad OR']],
+        materials=[mats['In'],
+                   mats['He'],
+                   mats['Zr'],
+                   mats['H2O']])
+
+    univs['pin plenum grid (intermediate)'] = make_pin(
+        'pin plenum grid (intermediate)',
+        surfaces=[surfs['plenum spring OR'],
+                  surfs['clad IR'],
+                  surfs['clad OR']],
+        materials=[mats['In'],
+                   mats['He'],
+                   mats['Zr'],
+                   mats['H2O']],
+        grid='intermediate')
+
+
+    #### 1.6% ENRICHED FUEL PIN CELL
+
+    # Determine z position between each fuel pellet, omitting the surfaces
+    # corresponding to the very bottom and top of the active fuel length
+    n_pellets = 196
+    axial_splits = np.linspace(bottom_fuel_rod, top_active_core, n_pellets + 1)[1:-1]
+    axial_surfs = [openmc.ZPlane(z0=z) for z in axial_splits]
+
+    # Get z-cylinder surfaces for each ring
+    rings = [surfs['fuel ring {}'.format(i)] for i in range(1, n_rings)]
+
+    # Create universe for UO2 alone with axial/radial subdivision
+    uo2_cells = []
+    for axial_region in subdivide(axial_surfs):
+        for ring_region in subdivide(rings):
+            uo2_cells.append(openmc.Cell(
+                fill=mats['UO2 1.6'],
+                region=axial_region & ring_region
+            ))
+    univs['Subdivided Fuel (1.6%)'] = openmc.Universe(cells=uo2_cells)
+
+    outside_pin_surfaces = [surfs['clad IR'], surfs['clad OR']]
+    outside_pin_mats = [mats['He'], mats['Zr'], mats['H2O']]
+
+    univs['Outside pin'] = make_pin(
+        'Outside pin',
+        surfaces=outside_pin_surfaces,
+        materials=outside_pin_mats)
+
+    univs['Outside pin grid (bottom)'] = make_pin(
+        'Outside pin grid (bottom)',
+        surfaces=outside_pin_surfaces,
+        materials=outside_pin_mats,
+        grid='bottom')
+
+    univs['Outside pin grid (intermediate)'] = make_pin(
+        'Outside pin grid (intermediate)',
+        surfaces=outside_pin_surfaces,
+        materials=outside_pin_mats,
+        grid='intermediate')
+
+    # Stack all axial pieces of 1.6% enriched fuel pin cell
+
+    within_fuel_surfs = [
+        surfs['grid1bot'],
+        surfs['grid1top'],
+        surfs['dashpot top'],
+        surfs['grid2bot'],
+        surfs['grid2top'],
+        surfs['grid3bot'],
+        surfs['grid3top'],
+        surfs['grid4bot'],
+        surfs['grid4top']
+    ]
+
+    univs['Fuel pin (1.6%) stack'] = make_pin_stack(
+        'Fuel pin (1.6%) stack',
+        zsurfaces=within_fuel_surfs,
+        universes=[
+            univs['Outside pin'],
+            univs['Outside pin grid (bottom)'],
+            univs['Outside pin'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin']
+        ],
+        boundary=surfs['pellet OR'],
+        pin_universe=univs['Subdivided Fuel (1.6%)'])
+
+    fuel_stack_surfs = [
+        surfs['bot support plate'],
+        surfs['top support plate'],
+        surfs['top lower nozzle'],
+        surfs['top lower thimble'],
+        surfs['top active core'],
+        surfs['grid5bot'],
+        surfs['grid5top'],
+        surfs['top pin plenum'],
+        surfs['top FR'],
+        surfs['bot upper nozzle'],
+        surfs['top upper nozzle']
+    ]
+
+    univs['Fuel (1.6%) stack'] = make_stack(
+        'Fuel (1.6%) stack',
+        surfaces=fuel_stack_surfs,
         universes=[univs['water pin'],
-                   univs['GT CR bank {} dummy nozzle'.format(b)],
-                   univs['GT CR bank {} dummy nozzle'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy grid (bottom)'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy grid (intermediate)'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy grid (bottom)'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy'.format(b)],
-                   univs['GT CR bank {} dummy bare'.format(b)],
-                   univs['GT CR bank {} dummy bare'.format(b)]])
+                   univs['SS pin'],
+                   univs['SS pin'],
+                   univs['end plug'],
+                   univs['Fuel pin (1.6%) stack'],
+                   univs['pin plenum'],
+                   univs['pin plenum grid (intermediate)'],
+                   univs['pin plenum'],
+                   univs['end plug'],
+                   univs['water pin'],
+                   univs['SS pin'],
+                   univs['water pin']])
 
 
-#### BURNABLE ABSORBER PIN CELLS
+    #### 2.4% ENRICHED FUEL PIN CELL
 
-univs['BA'] = make_pin(
-    'BA',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['BA IR 7'],
-              surfs['BA IR 8']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']])
+    # Create universe for UO2 alone with axial/radial subdivision
+    uo2_cells = []
+    for axial_region in subdivide(axial_surfs):
+        for ring_region in subdivide(rings):
+            uo2_cells.append(openmc.Cell(
+                fill=mats['UO2 2.4'],
+                region=axial_region & ring_region
+            ))
+    univs['Subdivided Fuel (2.4%)'] = openmc.Universe(cells=uo2_cells)
 
-univs['BA grid (bottom)'] = make_pin(
-    'BA grid (bottom)',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['BA IR 7'],
-              surfs['BA IR 8']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']],
-    grid='bottom')
+    # Stack all axial pieces of 2.4% enriched fuel pin cell
 
-univs['BA grid (intermediate)'] = make_pin(
-    'BA grid (intermediate)',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['BA IR 7'],
-              surfs['BA IR 8']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']],
-    grid='intermediate')
+    univs['Fuel pin (2.4%) stack'] = make_pin_stack(
+        'Fuel pin (2.4%) stack',
+        zsurfaces=within_fuel_surfs,
+        universes=[
+            univs['Outside pin'],
+            univs['Outside pin grid (bottom)'],
+            univs['Outside pin'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin']
+        ],
+        boundary=surfs['pellet OR'],
+        pin_universe=univs['Subdivided Fuel (2.4%)'])
 
-univs['BA dashpot'] = make_pin(
-    'BA dashpot',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['GT dashpot IR'],
-              surfs['GT dashpot OR']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']])
-
-univs['BA dashpot grid (bottom)'] = make_pin(
-    'BA dashpot grid (bottom)',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['GT dashpot IR'],
-              surfs['GT dashpot OR']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']],
-    grid='bottom')
-
-univs['BA dashpot grid (intermediate)'] = make_pin(
-    'BA dashpot grid (intermediate)',
-    surfaces=[surfs['BA IR 1'],
-              surfs['BA IR 2'],
-              surfs['BA IR 3'],
-              surfs['BA IR 4'],
-              surfs['BA IR 5'],
-              surfs['BA IR 6'],
-              surfs['GT dashpot IR'],
-              surfs['GT dashpot OR']],
-    materials=[mats['Air'],
-               mats['SS'],
-               mats['Air'],
-               mats['BSG'],
-               mats['Air'],
-               mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']],
-    grid='intermediate')
-
-univs['BA blank SS'] = make_pin(
-    'BA blank SS',
-    surfaces=[surfs['BA IR 6'],
-              surfs['BA IR 7'],
-              surfs['BA IR 8']],
-    materials=[mats['SS'],
-               mats['H2O'],
-               mats['Zr'],
-               mats['H2O']])
-
-univs['BA blank SS bare'] = make_pin(
-    'BA blank SS bare',
-    surfaces=[surfs['BA IR 6']],
-    materials=[mats['SS'],
-               mats['H2O']])
-
-stack_surfs_BA = [
-    surfs['bot support plate'],
-    surfs['top support plate'],
-    surfs['top lower nozzle'],
-    surfs['top lower thimble'],
-    surfs['grid1bot'],
-    surfs['BA bot'],
-    surfs['grid1top'],
-    surfs['dashpot top'],
-    surfs['grid2bot'],
-    surfs['grid2top'],
-    surfs['grid3bot'],
-    surfs['grid3top'],
-    surfs['grid4bot'],
-    surfs['grid4top'],
-    surfs['top active core'],
-    surfs['grid5bot'],
-    surfs['grid5top'],
-    surfs['top pin plenum'],
-    surfs['top FR'],
-    surfs['bot upper nozzle'],
-    surfs['top upper nozzle']]
-
-# Stack all axial pieces of control rod tubes together for each bank
-
-univs['BA stack'] = make_stack(
-    'BA stack', stack_surfs_BA,
-    universes=[univs['water pin'],
-               univs['water pin'],
-               univs['water pin'],
-               univs['GTd empty'],
-               univs['GTd empty'],
-               univs['GTd empty grid (bottom)'],
-               univs['BA dashpot grid (bottom)'],
-               univs['BA dashpot'],
-               univs['BA'],
-               univs['BA grid (intermediate)'],
-               univs['BA'],
-               univs['BA grid (intermediate)'],
-               univs['BA'],
-               univs['BA grid (intermediate)'],
-               univs['BA'],
-               univs['BA blank SS'],
-               univs['BA blank SS'],
-               univs['BA blank SS'],
-               univs['BA blank SS'],
-               univs['BA blank SS'],
-               univs['BA blank SS bare'],
-               univs['water pin']])
+    univs['Fuel (2.4%) stack'] = make_stack(
+        'Fuel (2.4%) stack',
+        surfaces=fuel_stack_surfs,
+        universes=[univs['water pin'],
+                   univs['SS pin'],
+                   univs['SS pin'],
+                   univs['end plug'],
+                   univs['Fuel pin (2.4%) stack'],
+                   univs['pin plenum'],
+                   univs['pin plenum grid (intermediate)'],
+                   univs['pin plenum'],
+                   univs['end plug'],
+                   univs['water pin'],
+                   univs['SS pin'],
+                   univs['water pin']])
 
 
-#### FUEL PIN CELLS
+    #### 3.1% ENRICHED FUEL PIN CELL
 
-univs['SS pin'] = make_pin(
-    'SS pin',
-    surfaces=[surfs['clad OR']],
-    materials=[mats['SS'],
-               mats['H2O']])
+    # Create universe for UO2 alone with axial/radial subdivision
+    uo2_cells = []
+    for axial_region in subdivide(axial_surfs):
+        for ring_region in subdivide(rings):
+            uo2_cells.append(openmc.Cell(
+                fill=mats['UO2 3.1'],
+                region=axial_region & ring_region
+            ))
+    univs['Subdivided Fuel (3.1%)'] = openmc.Universe(cells=uo2_cells)
 
-univs['end plug'] = make_pin(
-    'end plug',
-    surfaces=[surfs['clad OR']],
-    materials=[mats['Zr'],
-               mats['H2O']])
+    # Stack all axial pieces of 3.1% enriched fuel pin cell
 
-univs['pin plenum'] = make_pin(
-    'pin plenum',
-    surfaces=[surfs['plenum spring OR'],
-              surfs['clad IR'],
-              surfs['clad OR']],
-    materials=[mats['In'],
-               mats['He'],
-               mats['Zr'],
-               mats['H2O']])
+    univs['Fuel pin (3.1%) stack'] = make_pin_stack(
+        'Fuel pin (3.1%) stack',
+        zsurfaces=within_fuel_surfs,
+        universes=[
+            univs['Outside pin'],
+            univs['Outside pin grid (bottom)'],
+            univs['Outside pin'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin'],
+            univs['Outside pin grid (intermediate)'],
+            univs['Outside pin']
+        ],
+        boundary=surfs['pellet OR'],
+        pin_universe=univs['Subdivided Fuel (3.1%)'])
 
-univs['pin plenum grid (intermediate)'] = make_pin(
-    'pin plenum grid (intermediate)',
-    surfaces=[surfs['plenum spring OR'],
-              surfs['clad IR'],
-              surfs['clad OR']],
-    materials=[mats['In'],
-               mats['He'],
-               mats['Zr'],
-               mats['H2O']],
-    grid='intermediate')
+    univs['Fuel (3.1%) stack'] = make_stack(
+        'Fuel (3.1%) stack',
+        surfaces=fuel_stack_surfs,
+        universes=[univs['water pin'],
+                   univs['SS pin'],
+                   univs['SS pin'],
+                   univs['end plug'],
+                   univs['Fuel pin (3.1%) stack'],
+                   univs['pin plenum'],
+                   univs['pin plenum grid (intermediate)'],
+                   univs['pin plenum'],
+                   univs['end plug'],
+                   univs['water pin'],
+                   univs['SS pin'],
+                   univs['water pin']])
 
-
-#### 1.6% ENRICHED FUEL PIN CELL
-
-# Determine z position between each fuel pellet, omitting the surfaces
-# corresponding to the very bottom and top of the active fuel length
-n_pellets = 196
-axial_splits = np.linspace(bottom_fuel_rod, top_active_core, n_pellets + 1)[1:-1]
-axial_surfs = [openmc.ZPlane(z0=z) for z in axial_splits]
-
-# Get z-cylinder surfaces for each ring
-rings = [surfs['fuel ring {}'.format(i)] for i in range(1, n_rings)]
-
-# Create universe for UO2 alone with axial/radial subdivision
-uo2_cells = []
-for axial_region in subdivide(axial_surfs):
-    for ring_region in subdivide(rings):
-        uo2_cells.append(openmc.Cell(
-            fill=mats['UO2 1.6'],
-            region=axial_region & ring_region
-        ))
-univs['Subdivided Fuel (1.6%)'] = openmc.Universe(cells=uo2_cells)
-
-outside_pin_surfaces = [surfs['clad IR'], surfs['clad OR']]
-outside_pin_mats = [mats['He'], mats['Zr'], mats['H2O']]
-
-univs['Outside pin'] = make_pin(
-    'Outside pin',
-    surfaces=outside_pin_surfaces,
-    materials=outside_pin_mats)
-
-univs['Outside pin grid (bottom)'] = make_pin(
-    'Outside pin grid (bottom)',
-    surfaces=outside_pin_surfaces,
-    materials=outside_pin_mats,
-    grid='bottom')
-
-univs['Outside pin grid (intermediate)'] = make_pin(
-    'Outside pin grid (intermediate)',
-    surfaces=outside_pin_surfaces,
-    materials=outside_pin_mats,
-    grid='intermediate')
-
-# Stack all axial pieces of 1.6% enriched fuel pin cell
-
-within_fuel_surfs = [
-    surfs['grid1bot'],
-    surfs['grid1top'],
-    surfs['dashpot top'],
-    surfs['grid2bot'],
-    surfs['grid2top'],
-    surfs['grid3bot'],
-    surfs['grid3top'],
-    surfs['grid4bot'],
-    surfs['grid4top']
-]
-
-univs['Fuel pin (1.6%) stack'] = make_pin_stack(
-    'Fuel pin (1.6%) stack',
-    zsurfaces=within_fuel_surfs,
-    universes=[
-        univs['Outside pin'],
-        univs['Outside pin grid (bottom)'],
-        univs['Outside pin'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin']
-    ],
-    boundary=surfs['pellet OR'],
-    pin_universe=univs['Subdivided Fuel (1.6%)'])
-
-fuel_stack_surfs = [
-    surfs['bot support plate'],
-    surfs['top support plate'],
-    surfs['top lower nozzle'],
-    surfs['top lower thimble'],
-    surfs['top active core'],
-    surfs['grid5bot'],
-    surfs['grid5top'],
-    surfs['top pin plenum'],
-    surfs['top FR'],
-    surfs['bot upper nozzle'],
-    surfs['top upper nozzle']
-]
-
-univs['Fuel (1.6%) stack'] = make_stack(
-    'Fuel (1.6%) stack',
-    surfaces=fuel_stack_surfs,
-    universes=[univs['water pin'],
-               univs['SS pin'],
-               univs['SS pin'],
-               univs['end plug'],
-               univs['Fuel pin (1.6%) stack'],
-               univs['pin plenum'],
-               univs['pin plenum grid (intermediate)'],
-               univs['pin plenum'],
-               univs['end plug'],
-               univs['water pin'],
-               univs['SS pin'],
-               univs['water pin']])
-
-
-#### 2.4% ENRICHED FUEL PIN CELL
-
-# Create universe for UO2 alone with axial/radial subdivision
-uo2_cells = []
-for axial_region in subdivide(axial_surfs):
-    for ring_region in subdivide(rings):
-        uo2_cells.append(openmc.Cell(
-            fill=mats['UO2 2.4'],
-            region=axial_region & ring_region
-        ))
-univs['Subdivided Fuel (2.4%)'] = openmc.Universe(cells=uo2_cells)
-
-# Stack all axial pieces of 2.4% enriched fuel pin cell
-
-univs['Fuel pin (2.4%) stack'] = make_pin_stack(
-    'Fuel pin (2.4%) stack',
-    zsurfaces=within_fuel_surfs,
-    universes=[
-        univs['Outside pin'],
-        univs['Outside pin grid (bottom)'],
-        univs['Outside pin'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin']
-    ],
-    boundary=surfs['pellet OR'],
-    pin_universe=univs['Subdivided Fuel (2.4%)'])
-
-univs['Fuel (2.4%) stack'] = make_stack(
-    'Fuel (2.4%) stack',
-    surfaces=fuel_stack_surfs,
-    universes=[univs['water pin'],
-               univs['SS pin'],
-               univs['SS pin'],
-               univs['end plug'],
-               univs['Fuel pin (2.4%) stack'],
-               univs['pin plenum'],
-               univs['pin plenum grid (intermediate)'],
-               univs['pin plenum'],
-               univs['end plug'],
-               univs['water pin'],
-               univs['SS pin'],
-               univs['water pin']])
-
-
-#### 3.1% ENRICHED FUEL PIN CELL
-
-# Create universe for UO2 alone with axial/radial subdivision
-uo2_cells = []
-for axial_region in subdivide(axial_surfs):
-    for ring_region in subdivide(rings):
-        uo2_cells.append(openmc.Cell(
-            fill=mats['UO2 3.1'],
-            region=axial_region & ring_region
-        ))
-univs['Subdivided Fuel (3.1%)'] = openmc.Universe(cells=uo2_cells)
-
-# Stack all axial pieces of 3.1% enriched fuel pin cell
-
-univs['Fuel pin (3.1%) stack'] = make_pin_stack(
-    'Fuel pin (3.1%) stack',
-    zsurfaces=within_fuel_surfs,
-    universes=[
-        univs['Outside pin'],
-        univs['Outside pin grid (bottom)'],
-        univs['Outside pin'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin'],
-        univs['Outside pin grid (intermediate)'],
-        univs['Outside pin']
-    ],
-    boundary=surfs['pellet OR'],
-    pin_universe=univs['Subdivided Fuel (3.1%)'])
-
-univs['Fuel (3.1%) stack'] = make_stack(
-    'Fuel (3.1%) stack',
-    surfaces=fuel_stack_surfs,
-    universes=[univs['water pin'],
-               univs['SS pin'],
-               univs['SS pin'],
-               univs['end plug'],
-               univs['Fuel pin (3.1%) stack'],
-               univs['pin plenum'],
-               univs['pin plenum grid (intermediate)'],
-               univs['pin plenum'],
-               univs['end plug'],
-               univs['water pin'],
-               univs['SS pin'],
-               univs['water pin']])
+    return univs