Fix sourcepoint write in core scripts

Fix the build-assembly.py script

smr/build-assembly-long.py

@@ -0,0 +1,186 @@
+#!/usr/bin/env python3
+
+import argparse
+from math import pi, isclose
+from pathlib import Path
+
+import numpy as np
+from tqdm import tqdm
+import openmc
+
+from smr.materials import materials, clone
+from smr.surfaces import surfs, lattice_pitch, pin_pitch, bottom_fuel_stack, \
+    top_active_core, pellet_OR, active_fuel_length
+from smr.pins import pin_universes, make_stack
+
+
+# Define command-line options
+parser = argparse.ArgumentParser()
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', dest='multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('--clone', action='store_true',
+                    help='Clone materials for each cell instance')
+parser.add_argument('--no-clone', dest='clone', action='store_false',
+                    help='Do not clone materials for each cell instance')
+parser.add_argument('-a', '--axial', type=int, default=100,
+                    help='Number of axial subdivisions in fuel')
+parser.add_argument('-d', '--depleted', action='store_true',
+                    help='Whether UO2 compositions should represent depleted fuel')
+parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(clone=False, multipole=True)
+args = parser.parse_args()
+
+# Make directory for inputs
+if args.output_dir is None:
+    if args.depleted:
+        directory = Path('assembly-long-depleted')
+    else:
+        directory = Path('assembly-long-fresh')
+else:
+    directory = args.output_dir
+directory.mkdir(exist_ok=True)
+
+rings = [0.1*pin_pitch, 0.2*pin_pitch]
+
+assembly_long_surfs = [
+    surfs['bottom FR'],
+    surfs['bot active core'],
+    surfs['top active core'],
+    surfs['top pin plenum'],
+    surfs['top FR'],
+    surfs['bot upper nozzle'],
+    surfs['top upper nozzle']
+]
+
+univs = pin_universes(rings, args.axial, args.depleted)
+fuel_univ = make_stack(
+    'Fuel (3.1%) stack no grid',
+    surfaces=assembly_long_surfs,
+    universes=[
+        univs['water pin'],
+        univs['end plug'],
+        univs['Fuel pin (3.1%) no grid'],
+        univs['pin plenum'],
+        univs['end plug'],
+        univs['water pin']
+    ]
+)
+
+# Define the NumPy array indices for assembly locations where there
+# may be CR guide tubes, instrument tubes and burnable absorbers
+nonfuel_y = np.array([2,2,2,3,3,5,5,5,5,5,8,8,8,8,8,11,11,11,11,11,13,13,14,14,14])
+nonfuel_x = np.array([5,8,11,3,13,2,5,8,11,14,2,5,8,11,14,2,5,8,11,14,3,13,5,8,11])
+
+universes = np.full((17,17), fuel_univ)
+universes[nonfuel_y, nonfuel_x] = univs['GT empty']
+
+# Instantiate the lattice
+lattice = openmc.RectLattice(name='Pin lattice')
+lattice.lower_left = (-17.*pin_pitch/2., -17.*pin_pitch/2.)
+lattice.pitch = (pin_pitch, pin_pitch)
+lattice.universes = universes
+
+# Add lattice to bounding cell
+root_universe = openmc.Universe(name='Root universe')
+cell = openmc.Cell(name='Lattice cell')
+cell.fill = lattice
+z_bounds = +surfs['bottom FR'] & -surfs['top FR']
+cell.region = surfs['lat grid box inner'] & z_bounds
+root_universe.add_cell(cell)
+
+# Apply reflective boundaries on sides and vacuum on bottom/top
+surfs['bottom FR'].boundary_type = 'vacuum'
+surfs['top FR'].boundary_type = 'vacuum'
+for halfspace in surfs['lat grid box inner']:
+    halfspace.surface.boundary_type = 'reflective'
+
+# Define geometry with a single assembly
+geometry = openmc.Geometry(root_universe)
+
+h = active_fuel_length / args.axial
+
+fuel_mats = {}
+
+# Count the number of instances for each cell and material
+if args.clone:
+    geometry.determine_paths(instances_only=True)
+
+for cell in tqdm(geometry.get_all_material_cells().values(),
+                 desc='Differentiating materials / assigning volume'):
+    if cell.fill in materials:
+        # Determine if this material is fuel
+        is_fuel = 'UO2 Fuel' in cell.fill.name
+
+        # Fill cell with list of "differentiated" materials if requested
+        if args.clone:
+            cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
+
+        # Determine volume of each fuel material
+        if is_fuel:
+            upper_right = cell.region.bounding_box[1]
+            if isclose(upper_right[0], rings[0]):
+                ri, ro = 0.0, rings[0]
+            elif isclose(upper_right[0], rings[1]):
+                ri, ro = rings[0], rings[1]
+            else:
+                ri, ro = rings[1], pellet_OR
+
+            if args.clone:
+                for mat in cell.fill:
+                    mat.volume = pi * (ro*ro - ri*ri) * h
+            else:
+                # In non-clone mode, we still need to create a copy of the
+                # material for each ring since they get different volumes
+                if ri not in fuel_mats:
+                    cell.fill = cell.fill.clone()
+                    cell.fill.volume = pi * (ro*ro - ri*ri) * h
+                    fuel_mats[ri] = cell.fill
+                else:
+                    cell.fill = fuel_mats[ri]
+        else:
+            if args.clone:
+                for mat in cell.fill:
+                    mat.volume = 1.0
+            else:
+                cell.fill.volume = 1.0
+
+#### Create OpenMC "materials.xml" file
+print('Getting materials...')
+all_materials = geometry.get_all_materials()
+print('Creating materials collection...')
+materials = openmc.Materials(all_materials.values())
+print('Exporting materials to XML...')
+materials.export_to_xml(str(directory / 'materials.xml'))
+
+
+#### Create OpenMC "geometry.xml" file
+geometry.export_to_xml(str(directory / 'geometry.xml'))
+
+
+#### Create OpenMC "settings.xml" file
+
+# Construct uniform initial source distribution over fissionable zones
+lower_left = (-lattice_pitch/2, -lattice_pitch/2, bottom_fuel_stack)
+upper_right = (lattice_pitch/2, lattice_pitch/2, top_active_core)
+source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
+source.space.only_fissionable = True
+
+settings = openmc.Settings()
+settings.batches = 200
+settings.inactive = 100
+settings.particles = 10000
+settings.output = {'tallies': False, 'summary': False}
+settings.source = source
+settings.sourcepoint = {'write': False}
+settings.temperature = {
+    'default': 531.5,
+    'method': 'interpolation',
+        'range': (500.0, 1300.0)
+}
+if args.multipole:
+    settings.temperature['multipole'] = True
+    settings.temperature['tolerance'] = 1000
+
+settings.export_to_xml(str(directory / 'settings.xml'))

smr/build-assembly-short.py

@@ -0,0 +1,192 @@
+#!/usr/bin/env python3
+
+import argparse
+from math import pi, isclose
+from pathlib import Path
+
+import numpy as np
+from tqdm import tqdm
+import openmc
+
+from smr.materials import materials, clone
+from smr.surfaces import surfs, lattice_pitch, pin_pitch, bottom_fuel_stack, \
+    top_active_core, pellet_OR, clad_OR, clad_IR, guide_tube_IR, guide_tube_OR
+from smr.pins import pin_universes
+
+
+# Define command-line options
+parser = argparse.ArgumentParser()
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='mat',
+                    help='Whether to use distribmats or distribcells for tallies')
+parser.add_argument('-a', '--axial', type=int, default=10,
+                    help='Number of axial subdivisions in fuel')
+parser.add_argument('-d', '--depleted', action='store_true',
+                    help='Whether UO2 compositions should represent depleted fuel')
+parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(multipole=True)
+args = parser.parse_args()
+
+# Make directory for inputs
+if args.output_dir is None:
+    if args.depleted:
+        directory = Path('assembly-short-depleted')
+    else:
+        directory = Path('assembly-short-fresh')
+else:
+    directory = args.output_dir
+directory.mkdir(exist_ok=True)
+
+rings = [0.1*pin_pitch, 0.2*pin_pitch]
+
+# Define the NumPy array indices for assembly locations where there
+# may be CR guide tubes, instrument tubes and burnable absorbers
+nonfuel_y = np.array([2,2,2,3,3,5,5,5,5,5,8,8,8,8,8,11,11,11,11,11,13,13,14,14,14])
+nonfuel_x = np.array([5,8,11,3,13,2,5,8,11,14,2,5,8,11,14,2,5,8,11,14,3,13,5,8,11])
+
+# NO BURNABLE ABSORBERS
+pins = pin_universes(rings, args.axial, args.depleted)
+gtu = pins['GT empty']
+#gti = pins['GT empty instr']
+universes = np.empty((17,17), dtype=openmc.Universe)
+universes[:,:] = pins['Fuel pin (3.1%) no grid']
+universes[nonfuel_y, nonfuel_x] = [    gtu,   gtu,   gtu,
+                                     gtu,              gtu,
+                                   gtu, gtu,  gtu,  gtu, gtu,
+                                   gtu, gtu,  gtu,  gtu, gtu,
+                                   gtu, gtu,  gtu,  gtu, gtu,
+                                     gtu,              gtu,
+                                       gtu,   gtu,   gtu     ]
+
+# Instantiate the lattice
+lattice = openmc.RectLattice(name='Pin lattice')
+lattice.lower_left = (-17.*pin_pitch/2., -17.*pin_pitch/2.)
+lattice.pitch = (pin_pitch, pin_pitch)
+lattice.universes = universes
+
+# Add lattice to bounding cell
+root_universe = openmc.Universe(name='Root universe')
+cell = openmc.Cell(name='Lattice cell')
+cell.fill = lattice
+z_bounds = +surfs['bot active core'] & -surfs['top active core']
+cell.region = surfs['lat grid box inner'] & z_bounds
+root_universe.add_cell(cell)
+
+# Apply reflective boundaries
+surfs['bot active core'].boundary_type = 'reflective'
+surfs['top active core'].boundary_type = 'reflective'
+for halfspace in surfs['lat grid box inner']:
+    halfspace.surface.boundary_type = 'reflective'
+
+# Define geometry with a single assembly
+geometry = openmc.Geometry(root_universe)
+
+#### "Differentiate" the geometry if using distribmats
+h = 10.0*pin_pitch / args.axial
+if args.tallies == 'mat':
+    # Count the number of instances for each cell and material
+    geometry.determine_paths(instances_only=True)
+
+    for cell in tqdm(geometry.get_all_material_cells().values(),
+                     desc='Differentiating materials'):
+        if cell.fill in materials:
+            # Fill cell with list of "differentiated" materials
+            cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
+
+            # Determine volume of each fuel material
+            if 'UO2 Fuel' in cell.fill[0].name:
+                lower_left, _ = cell.region.bounding_box
+                if isclose(lower_left[0], rings[0]):
+                    ri, ro = 0.0, rings[0]
+                elif isclose(lower_left[0], rings[1]):
+                    ri, ro = rings[0], rings[1]
+                else:
+                    ri, ro = rings[1], pellet_OR
+                for mat in cell.fill:
+                    mat.volume = pi * (ro*ro - ri*ri) * h
+            elif cell.fill[0].name == 'Borated Water':
+                for mat in cell.fill:
+                    mat.volume = pin_pitch**2 - pi*clad_OR**2 * h
+            elif cell.fill[0].name == 'Helium':
+                for mat in cell.fill:
+                    mat.volume = pi * (clad_IR**2 - pellet_OR**2) * h
+            elif cell.fill[0].name == 'M5':
+                for mat in cell.fill:
+                    mat.volume = pi * (clad_OR**2 - clad_IR**2) * h
+            elif cell.fill[0].name == 'Zircaloy-4':
+                for mat in cell.fill:
+                    mat.volume = pi * (guide_tube_OR**2 - guide_tube_IR**2) * h
+
+#### Create OpenMC "materials.xml" file
+print('Getting materials...')
+all_materials = geometry.get_all_materials()
+print('Creating materials collection...')
+materials = openmc.Materials(all_materials.values())
+print('Exporting materials to XML...')
+materials.export_to_xml(str(directory / 'materials.xml'))
+
+
+#### Create OpenMC "geometry.xml" file
+geometry.export_to_xml(str(directory / 'geometry.xml'))
+
+
+#### Create OpenMC "settings.xml" file
+
+# Construct uniform initial source distribution over fissionable zones
+lower_left = (-lattice_pitch/2, -lattice_pitch/2, bottom_fuel_stack)
+upper_right = (lattice_pitch/2, lattice_pitch/2, top_active_core)
+source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
+source.space.only_fissionable = True
+
+settings = openmc.Settings()
+settings.batches = 200
+settings.inactive = 100
+settings.particles = 10000
+settings.output = {'tallies': False, 'summary': False}
+settings.source = source
+settings.sourcepoint = {'write': False}
+
+if args.multipole:
+    settings.temperature = {
+        'multipole': True,
+        'tolerance': 1000,
+        'default': 531.5,
+        'method': 'interpolation',
+        'range': (500.0, 1300.0)
+    }
+
+settings.export_to_xml(str(directory / 'settings.xml'))
+
+
+####  Create OpenMC "tallies.xml" file
+tallies = openmc.Tallies()
+
+# Extract all fuel materials
+materials = geometry.get_materials_by_name(name='Fuel', matching=False)
+
+# If using distribcells, create distribcell tally needed for depletion
+if args.tallies == 'cell':
+    # Extract all cells filled by a fuel material
+    fuel_cells = []
+    for cell in geometry.get_all_cells().values():
+        if cell.fill in materials:
+            tally = openmc.Tally(name='depletion tally')
+            tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
+                            'fission', '(n,2n)', '(n,3n)', '(n,4n)']
+            tally.nuclides = cell.fill.get_nuclides()
+            tally.filters.append(openmc.DistribcellFilter([cell]))
+            tallies.append(tally)
+
+# If using distribmats, create material tally needed for depletion
+elif args.tallies == 'mat':
+    tally = openmc.Tally(name='depletion tally')
+    tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
+                    'fission', '(n,2n)', '(n,3n)', '(n,4n)']
+    tally.nuclides = materials[0].get_nuclides()
+    tally.filters = [openmc.MaterialFilter(materials)]
+    tallies.append(tally)
+
+tallies.export_to_xml(str(directory / 'tallies.xml'))

smr/build-assembly.py

@@ -1,24 +1,28 @@
 #!/usr/bin/env python3
 
 import argparse
-import copy
 from pathlib import Path
 
 import numpy as np
 from tqdm import tqdm
 
 import openmc
-from smr.materials import materials
-from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core
+from smr.materials import materials, clone
+from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core, pellet_OR
 from smr.assemblies import assembly_universes
-from smr.plots import assembly_plots
 from smr import inlet_temperature
 
 
 # Define command-line options
 parser = argparse.ArgumentParser()
-parser.add_argument('-m', '--multipole', action='store_true',
-                    help='Whether to use multipole cross sections')
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', dest='multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('--clone', action='store_true',
+                    help='Clone materials for each cell instance')
+parser.add_argument('--no-clone', dest='clone', action='store_false',
+                    help='Do not clone materials for each cell instance')
 parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='mat',
                     help='Whether to use distribmats or distribcells for tallies')
 parser.add_argument('-r', '--rings', type=int, default=10,
@@ -28,6 +32,7 @@ parser.add_argument('-a', '--axial', type=int, default=196,
 parser.add_argument('-d', '--depleted', action='store_true',
                     help='Whether UO2 compositions should represent depleted fuel')
 parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(clone=False, multipole=True)
 args = parser.parse_args()
 
 # Make directory for inputs
@@ -41,26 +46,22 @@ else:
 directory.mkdir(exist_ok=True)
 
 # Define geometry with a single assembly
-assembly = assembly_universes(args.rings, args.axial, args.depleted)
-lattice_sides = openmc.model.get_rectangular_prism(lattice_pitch, lattice_pitch,
-                                                   boundary_type='reflective')
+if args.rings > 1:
+    ring_radii = np.sqrt(np.arange(1, args.rings)*pellet_OR**2 / args.rings)
+else:
+    ring_radii = None
+assembly = assembly_universes(ring_radii, args.axial, args.depleted)
+lattice_sides = openmc.model.rectangular_prism(lattice_pitch, lattice_pitch,
+                                               boundary_type='reflective')
 main_cell = openmc.Cell(
-    fill=assembly['Assembly (3.1%) 16BA'],
+    fill=assembly['Assembly (3.1%)'],
     region=lattice_sides & +surfs['lower bound'] & -surfs['upper bound']
 )
 root_univ = openmc.Universe(cells=[main_cell])
 geometry = openmc.Geometry(root_univ)
 
-
-def clone(material):
-    """Perform copy of material but share nuclide densities"""
-    shared_mat = copy.copy(material)
-    shared_mat.id = None
-    return shared_mat
-
-
 #### "Differentiate" the geometry if using distribmats
-if args.tallies == 'mat':
+if args.clone:
     # Count the number of instances for each cell and material
     geometry.determine_paths(instances_only=True)
 
@@ -100,7 +101,7 @@ settings.inactive = 100
 settings.particles = 10000
 settings.output = {'tallies': False, 'summary': False}
 settings.source = source
-settings.sourcepoint_write = False
+settings.sourcepoint = {'write': False}
 settings.temperature = {
     'default': inlet_temperature,
     'method': 'interpolation',
@@ -142,7 +143,3 @@ elif args.tallies == 'mat':
     tallies.append(tally)
 
 tallies.export_to_xml(str(directory / 'tallies.xml'))
-
-# Create plots
-plots = assembly_plots(main_cell.fill)
-plots.export_to_xml(str(directory / 'plots.xml'))

smr/build-core-fresh.py

@@ -1,34 +1,30 @@
 #!/usr/bin/env python3
 
-import os
-import shutil
-import copy
 import argparse
+from math import pi
 from pathlib import Path
 
 import numpy as np
-
 import openmc
-from smr.materials import materials
-from smr.plots import core_plots
-from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core
+from tqdm import tqdm
+
+from smr.materials import materials, clone
+from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core, \
+    pellet_OR, active_fuel_length
 from smr.core import core_geometry
 from smr import inlet_temperature
 
 
-def clone(mat):
-    """Make a shallow copy of a material (sharing compositions)."""
-    mat_copy = copy.copy(mat)
-    mat_copy.id = None
-    return mat_copy
-
-
 # Define command-line options
 parser = argparse.ArgumentParser()
-parser.add_argument('-m', '--multipole', action='store_true',
-                    help='Whether to use multipole cross sections')
-parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='cell',
-                    help='Whether to use distribmats or distribcells for tallies')
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', dest='multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('--clone', action='store_true',
+                    help='Clone materials for each cell instance')
+parser.add_argument('--no-clone', dest='clone', action='store_false',
+                    help='Do not clone materials for each cell instance')
 parser.add_argument('-r', '--rings', type=int, default=10,
                     help='Number of annular regions in fuel')
 parser.add_argument('-a', '--axial', type=int, default=196,
@@ -36,6 +32,7 @@ parser.add_argument('-a', '--axial', type=int, default=196,
 parser.add_argument('-d', '--depleted', action='store_true',
                     help='Whether UO2 compositions should represent depleted fuel')
 parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(clone=False, multipole=True)
 args = parser.parse_args()
 
 # Make directory for inputs
@@ -48,20 +45,44 @@ else:
     directory = args.output_dir
 directory.mkdir(exist_ok=True)
 
-geometry = core_geometry(args.rings, args.axial, args.depleted)
+if args.rings > 1:
+    ring_radii = np.sqrt(np.arange(1, args.rings)*pellet_OR**2 / args.rings)
+else:
+    ring_radii = None
+geometry = core_geometry(ring_radii, args.axial, args.depleted)
 
-#### "Differentiate" the geometry if using distribmats
-if args.tallies == 'mat':
-    # Count the number of instances for each cell and material
+h = active_fuel_length / args.axial
+fuel_mats = {}
+
+# Count the number of instances for each cell and material
+if args.clone:
     geometry.determine_paths(instances_only=True)
 
-    # Extract all cells filled by a fuel material
-    fuel_mats = {m for m in materials if 'UO2 Fuel' in m.name}
+fuel_volume = pi * pellet_OR**2 * h / args.rings
+for cell in tqdm(geometry.get_all_cells().values(),
+                 desc='Differentiating materials / assigning volume'):
+    if cell.fill in materials:
+        # Determine if this material is fuel
+        name = cell.fill.name
+        is_fuel = 'UO2 Fuel' in name
 
-    for cell in geometry.get_all_cells().values():
-        if cell.fill in fuel_mats:
-            # Fill cell with list of "differentiated" materials
-            cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
+        # Determine volume of each fuel material
+        if is_fuel:
+            if args.clone:
+                # Fill cell with list of "differentiated" materials if requested
+                cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
+                for mat in cell.fill:
+                    mat.volume = fuel_volume
+            else:
+                r_o = cell.region.bounding_box[1][0]
+                if (name, r_o) not in fuel_mats:
+                    cell.fill = cell.fill.clone()
+                    cell.fill.volume = fuel_volume
+                    fuel_mats[name, r_o] = cell.fill
+                else:
+                    cell.fill = fuel_mats[name, r_o]
+        else:
+            cell.fill.volume = 1.0
 
 #### Create OpenMC "materials.xml" file
 all_materials = geometry.get_all_materials()
@@ -87,8 +108,7 @@ settings.inactive = 100
 settings.particles = 10000
 settings.output = {'tallies': False, 'summary': False}
 settings.source = source
-settings.sourcepoint_write = False
-
+settings.sourcepoint = {'write': False}
 settings.temperature = {
     'default': inlet_temperature,
     'method': 'interpolation',
@@ -99,39 +119,3 @@ if args.multipole:
     settings.temperature['tolerance'] = 1000
 
 settings.export_to_xml(str(directory / 'settings.xml'))
-
-
-#### Create OpenMC "plots.xml" file
-plots = core_plots()
-plots.export_to_xml(str(directory / 'plots.xml'))
-
-
-####  Create OpenMC "tallies.xml" file
-tallies = openmc.Tallies()
-
-# Extract all fuel materials
-materials = geometry.get_materials_by_name(name='Fuel', matching=False)
-
-# If using distribcells, create distribcell tally needed for depletion
-if args.tallies == 'cell':
-    # Extract all cells filled by a fuel material
-    fuel_cells = []
-    for cell in geometry.get_all_cells().values():
-        if cell.fill in materials:
-            tally = openmc.Tally(name='depletion tally')
-            tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
-                            'fission', '(n,2n)', '(n,3n)', '(n,4n)']
-            tally.nuclides = cell.fill.get_nuclides()
-            tally.filters.append(openmc.DistribcellFilter([cell]))
-            tallies.append(tally)
-
-# If using distribmats, create material tally needed for depletion
-elif args.tallies == 'mat':
-    tally = openmc.Tally(name='depletion tally')
-    tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
-                    'fission', '(n,2n)', '(n,3n)', '(n,4n)']
-    tally.nuclides = materials[0].get_nuclides()
-    tally.filters = [openmc.MaterialFilter(materials)]
-    tallies.append(tally)
-
-tallies.export_to_xml(str(directory / 'tallies.xml'))

smr/build-core-long.py

@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+
+import argparse
+from math import pi, isclose
+from pathlib import Path
+
+import openmc
+from smr.materials import materials
+from smr.surfaces import bottom_fuel_stack, top_active_core, \
+    pellet_OR, pin_pitch, clad_IR, clad_OR, active_fuel_length
+from smr.core import core_geometry
+from smr import inlet_temperature
+import smr.surfaces
+
+# Define command-line options
+parser = argparse.ArgumentParser()
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', dest='multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('-a', '--axial', type=int, default=100,
+                    help='Number of axial subdivisions in fuel')
+parser.add_argument('-d', '--depleted', action='store_true',
+                    help='Whether UO2 compositions should represent depleted fuel')
+parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(multipole=True)
+args = parser.parse_args()
+
+# Make directory for inputs
+if args.output_dir is None:
+    if args.depleted:
+        directory = Path('core-long-depleted')
+    else:
+        directory = Path('core-long-fresh')
+else:
+    directory = args.output_dir
+directory.mkdir(exist_ok=True)
+
+# Modify lattice pitch
+smr.surfaces.lattice_pitch = lattice_pitch = 17*smr.surfaces.pin_pitch
+
+ring_radii = [0.1*pin_pitch, 0.2*pin_pitch]
+
+geometry = core_geometry(ring_radii, args.axial, args.depleted)
+
+h = active_fuel_length / args.axial
+fuel_mats = {}
+
+for cell in geometry.get_all_cells().values():
+    if cell.fill in materials:
+        # Determine volume of each fuel material
+        name = cell.fill.name
+        if 'UO2 Fuel' in name:
+            upper_right = cell.region.bounding_box[1][0]
+            if isclose(upper_right, ring_radii[0]):
+                ri, ro = 0.0, ring_radii[0]
+            elif isclose(upper_right, ring_radii[1]):
+                ri, ro = ring_radii[0], ring_radii[1]
+            else:
+                ri, ro = ring_radii[1], pellet_OR
+            if (name, ri) not in fuel_mats:
+                cell.fill = cell.fill.clone()
+                cell.fill.volume = pi * (ro*ro - ri*ri) * h
+                fuel_mats[name, ri] = cell.fill
+            else:
+                cell.fill = fuel_mats[name, ri]
+        elif name == 'Helium':
+            cell.fill.volume = pi * (clad_IR**2 - pellet_OR**2) * h
+        elif name == 'M5':
+            # Clad is not subdivided
+            cell.fill.volume = pi * (clad_OR**2 - clad_IR**2) * active_fuel_length
+        else:
+            cell.fill.volume = 1.0
+
+
+#### Create OpenMC "materials.xml" file
+all_materials = geometry.get_all_materials()
+materials = openmc.Materials(all_materials.values())
+materials.export_to_xml(str(directory / 'materials.xml'))
+
+#### Create OpenMC "geometry.xml" file
+geometry.export_to_xml(str(directory / 'geometry.xml'))
+
+
+#### Create OpenMC "settings.xml" file
+
+# Construct uniform initial source distribution over fissionable zones
+lower_left = [-7.*lattice_pitch/2., -7.*lattice_pitch/2., bottom_fuel_stack]
+upper_right = [+7.*lattice_pitch/2., +7.*lattice_pitch/2., top_active_core]
+source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
+source.space.only_fissionable = True
+
+settings = openmc.Settings()
+settings.batches = 200
+settings.inactive = 100
+settings.particles = 20_000_000
+settings.output = {'tallies': False, 'summary': False}
+settings.source = source
+settings.sourcepoint = {'write': False}
+settings.temperature = {
+    'default': inlet_temperature,
+    'method': 'interpolation',
+    'range': (300.0, 1500.0),
+}
+if args.multipole:
+    settings.temperature['multipole'] = True
+    settings.temperature['tolerance'] = 1000
+
+settings.export_to_xml(str(directory / 'settings.xml'))
+
+# Check assembly power distribution
+core_lattice = geometry.get_cells_by_fill_name('Main core')[0].fill
+mesh = openmc.RegularMesh.from_rect_lattice(core_lattice)
+assembly_power = openmc.Tally()
+assembly_power.filters = [openmc.MeshFilter(mesh)]
+assembly_power.scores = ['nu-fission']
+tallies = openmc.Tallies([assembly_power])
+tallies.export_to_xml(directory / 'tallies.xml')

smr/build-core-short.py

@@ -0,0 +1,131 @@
+#!/usr/bin/env python3
+
+import argparse
+from math import pi, isclose
+from pathlib import Path
+
+import openmc
+from smr.materials import materials
+from smr.surfaces import bottom_fuel_stack, top_active_core, \
+    pellet_OR, surfs, pin_pitch, clad_IR, clad_OR
+import smr.surfaces
+import smr.pins
+from smr.core import core_geometry
+from smr import inlet_temperature
+
+# Define command-line options
+parser = argparse.ArgumentParser()
+parser.add_argument('--multipole', action='store_true',
+                    help='Use multipole cross sections')
+parser.add_argument('--no-multipole', action='store_false',
+                    help='Do not use multipole cross sections')
+parser.add_argument('-a', '--axial', type=int, default=3,
+                    help='Number of axial subdivisions in fuel')
+parser.add_argument('-d', '--depleted', action='store_true',
+                    help='Whether UO2 compositions should represent depleted fuel')
+parser.add_argument('-o', '--output-dir', type=Path, default=None)
+parser.set_defaults(multipole=True)
+args = parser.parse_args()
+
+# Make directory for inputs
+if args.output_dir is None:
+    if args.depleted:
+        directory = Path('core-short-depleted')
+    else:
+        directory = Path('core-short-fresh')
+else:
+    directory = args.output_dir
+directory.mkdir(exist_ok=True)
+
+# Modify lattice pitch
+smr.surfaces.lattice_pitch = lattice_pitch = 17*smr.surfaces.pin_pitch
+
+# Modify fuel length
+length = 3. * pin_pitch
+smr.surfaces.active_fuel_length = length
+smr.pins.top_active_core = length
+surfs['top active core'].z0 = length
+
+# Change top and bottom of model to contain only fuel
+surfs['lower bound'].z0 = 0.0
+surfs['lower bound'].boundary_type = 'reflective'
+surfs['upper bound'].z0 = length
+surfs['upper bound'].boundary_type = 'reflective'
+
+ring_radii = [0.1*pin_pitch, 0.2*pin_pitch]
+
+geometry = core_geometry(ring_radii, args.axial, args.depleted)
+
+h = length / args.axial
+fuel_mats = {}
+
+for cell in geometry.get_all_cells().values():
+    if cell.fill in materials:
+        # Determine volume of each fuel material
+        name = cell.fill.name
+        if 'UO2 Fuel' in name:
+            upper_right = cell.region.bounding_box[1][0]
+            if isclose(upper_right, ring_radii[0]):
+                ri, ro = 0.0, ring_radii[0]
+            elif isclose(upper_right, ring_radii[1]):
+                ri, ro = ring_radii[0], ring_radii[1]
+            else:
+                ri, ro = ring_radii[1], pellet_OR
+            if (name, ri) not in fuel_mats:
+                cell.fill = cell.fill.clone()
+                cell.fill.volume = pi * (ro*ro - ri*ri) * h
+                fuel_mats[name, ri] = cell.fill
+            else:
+                cell.fill = fuel_mats[name, ri]
+        elif name == 'Helium':
+            cell.fill.volume = pi * (clad_IR**2 - pellet_OR**2) * h
+        elif name == 'M5':
+            # Clad is not subdivided
+            cell.fill.volume = pi * (clad_OR**2 - clad_IR**2) * length
+        else:
+            cell.fill.volume = 1.0
+
+
+#### Create OpenMC "materials.xml" file
+all_materials = geometry.get_all_materials()
+materials = openmc.Materials(all_materials.values())
+materials.export_to_xml(str(directory / 'materials.xml'))
+
+#### Create OpenMC "geometry.xml" file
+geometry.export_to_xml(str(directory / 'geometry.xml'))
+
+
+#### Create OpenMC "settings.xml" file
+
+# Construct uniform initial source distribution over fissionable zones
+lower_left = [-7.*lattice_pitch/2., -7.*lattice_pitch/2., bottom_fuel_stack]
+upper_right = [+7.*lattice_pitch/2., +7.*lattice_pitch/2., top_active_core]
+source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
+source.space.only_fissionable = True
+
+settings = openmc.Settings()
+settings.batches = 200
+settings.inactive = 100
+settings.particles = 10000
+settings.output = {'tallies': False, 'summary': False}
+settings.source = source
+settings.sourcepoint = {'write': False}
+settings.temperature = {
+    'default': inlet_temperature,
+    'method': 'interpolation',
+    'range': (300.0, 1500.0),
+}
+if args.multipole:
+    settings.temperature['multipole'] = True
+    settings.temperature['tolerance'] = 1000
+
+settings.export_to_xml(str(directory / 'settings.xml'))
+
+# Check assembly power distribution
+core_lattice = geometry.get_cells_by_fill_name('Main core')[0].fill
+mesh = openmc.RegularMesh.from_rect_lattice(core_lattice)
+assembly_power = openmc.Tally()
+assembly_power.filters = [openmc.MeshFilter(mesh)]
+assembly_power.scores = ['nu-fission']
+tallies = openmc.Tallies([assembly_power])
+tallies.export_to_xml(directory / 'tallies.xml')

smr/milestones.md

@@ -0,0 +1,20 @@
+# Milestone Models
+
+- **AD-SE-08-61, Coupled Multiphysics Driver Implementation** --- This milestone
+  used the singlerod short/long problems. Generating the model was done with the
+  script `tests/singlerod/make_openmc_model.py` from the ENRICO repository
+  (there is a `--short` command line option to generate the short version)
+
+- **AD-SE-08-66, Coupled Assembly Analysis** --- This milestone used the
+  assembly short and long (v2) problems. Generating the model was done with
+  `smr/build-assembly-long.py -a 100 --clone` on the ecp-benchmarks repository
+  (git commit `631fefe`, after pull request #14). In these models, materials are
+  fully differentiated across each fuel ring/axial segment.
+
+- **AD-SE-08-73, Full core coupled-physics simulation** --- This milestone used
+  the core-short and core-long (90 layer) models. Generating the models was done
+  with `smr/build-core-short.py` and `smr/build-assembly-long.py -a 90` on the
+  ecp-benchmarks repository (git commit `c7b89db`, after pull request #16). In
+  these models, materials are not differentiated and no grid spacers are
+  present. The lattice pitch is modified to be exactly 17 times the pin pitch
+  (slightly different than NuScale specification).

smr/smr/assemblies.py

@@ -120,13 +120,14 @@ def make_assembly(name, universes):
     return universe
 
 
-def assembly_universes(num_rings, num_axial, depleted):
+def assembly_universes(ring_radii, num_axial, depleted):
     """Generate universes for SMR fuel assemblies.
 
     Parameters
     ----------
-    num_rings : int
-        Number of annual regions in fuel
+    ring_radii : iterable of float
+        Radii of rings in fuel (note that this doesn't need to include the
+        full fuel pin radius)
     num_axial : int
         Number of axial subdivisions in fuel
     depleted : bool
@@ -138,15 +139,14 @@ def assembly_universes(num_rings, num_axial, depleted):
         Dictionary mapping a universe name to a openmc.Universe object
 
     """
-    pins = pin_universes(num_rings, num_axial, depleted)
+    pins = pin_universes(ring_radii, num_axial, depleted)
 
     # Create dictionary to store assembly universes
     univs = {}
 
     # commonly needed universes
-    gtu = pins['GT empty']
+    gtu = pins['GT empty stack']
     gti = pins['GT empty instr']
-    bas = pins['BA stack']
     ins = pins['IT stack']
     crA = pins['GT CR bank A']
     crB = pins['GT CR bank B']
@@ -215,7 +215,7 @@ def assembly_universes(num_rings, num_axial, depleted):
                                            gtu, gtu,  gtu,  gtu, gtu,
                                              gtu,              gtu,
                                                gtu,   gtu,   gtu     ]
-        univs['Assembly (2.4%) no BAs' + comment] = \
+        univs['Assembly (2.4%)' + comment] = \
             make_assembly('Assembly (2.4%) no BAs' + comment, universes)
 
         # WITH CONTROL ROD D BANK
@@ -231,33 +231,6 @@ def assembly_universes(num_rings, num_axial, depleted):
         univs['Assembly (2.4%) CR D' + comment] = \
             make_assembly('Assembly (2.4%) CR D' + comment, universes)
 
-        # WITH 12 BURNABLE ABSORBERS
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (2.4%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   gtu,   bas,
-                                             bas,              bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                           gtu, gtu,  cent, gtu, gtu,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                             bas,              bas,
-                                               bas,   gtu,   bas     ]
-        univs['Assembly (2.4%) 12BA' + comment] = \
-            make_assembly('Assembly (2.4%) 12BA' + comment, universes)
-
-        # WITH 16 BURNABLE ABSORBERS
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (2.4%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   bas,   bas,
-                                             bas,              bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                           bas, gtu,  cent, gtu, bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                             bas,              bas,
-                                               bas,   bas,   bas     ]
-        univs['Assembly (2.4%) 16BA' + comment] = \
-            make_assembly('Assembly (2.4%) 16BA' + comment, universes)
-
-
     #### 3.1% ENRICHED ASSEMBLIES
 
     for cent, comment in [(gti, ''), (ins, ' instr')]:
@@ -288,134 +261,5 @@ def assembly_universes(num_rings, num_axial, depleted):
         univs['Assembly (3.1%) CR SA' + comment] = \
             make_assembly('Assembly (3.1%) CR SA' + comment, universes)
 
-        # WITH 20 BURNABLE ABSORBERS
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   bas,   bas,
-                                             bas,              bas,
-                                           bas, bas,  gtu,  bas, bas,
-                                           bas, gtu,  cent, gtu, bas,
-                                           bas, bas,  gtu,  bas, bas,
-                                             bas,              bas,
-                                               bas,   bas,   bas     ]
-        univs['Assembly (3.1%) 20BA' + comment] = \
-            make_assembly('Assembly (3.1%) 20BA' + comment, universes)
-
-        # WITH 16 BURNABLE ABSORBERS
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   bas,   bas,
-                                             bas,              bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                           bas, gtu,  cent, gtu, bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                             bas,              bas,
-                                               bas,   bas,   bas     ]
-        univs['Assembly (3.1%) 16BA' + comment] = \
-            make_assembly('Assembly (3.1%) 16BA' + comment, universes)
-
-        # WITH 15 BURNABLE ABSORBERS NW
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    gtu,   gtu,   gtu,
-                                             gtu,              gtu,
-                                           gtu, bas,  bas,  bas, bas,
-                                           gtu, bas,  cent, bas, bas,
-                                           gtu, bas,  bas,  bas, bas,
-                                             gtu,              bas,
-                                               bas,   bas,   bas     ]
-        univs['Assembly (3.1%) 15BANW' + comment] = \
-            make_assembly('Assembly (3.1%) 15BANW' + comment, universes)
-
-        # WITH 15 BURNABLE ABSORBERS NE
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    gtu,   gtu,   gtu,
-                                             gtu,              gtu,
-                                           bas, bas,  bas,  bas, gtu,
-                                           bas, bas,  cent, bas, gtu,
-                                           bas, bas,  bas,  bas, gtu,
-                                             bas,              gtu,
-                                               bas,   bas,   bas     ]
-        univs['Assembly (3.1%) 15BANE' + comment] = \
-            make_assembly('Assembly (3.1%) 15BANE' + comment, universes)
-
-        # WITH 15 BURNABLE ABSORBERS SW
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   bas,   bas,
-                                             gtu,              bas,
-                                           gtu, bas,  bas,  bas, bas,
-                                           gtu, bas,  cent, bas, bas,
-                                           gtu, bas,  bas,  bas, bas,
-                                             gtu,              gtu,
-                                               gtu,   gtu,   gtu     ]
-        univs['Assembly (3.1%) 15BASW' + comment] = \
-            make_assembly('Assembly (3.1%) 15BASW' + comment, universes)
-
-        # WITH 15 BURNABLE ABSORBERS SE
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   bas,   bas,
-                                             bas,              gtu,
-                                           bas, bas,  bas,  bas, gtu,
-                                           bas, bas,  cent, bas, gtu,
-                                           bas, bas,  bas,  bas, gtu,
-                                             gtu,              gtu,
-                                               gtu,   gtu,   gtu     ]
-        univs['Assembly (3.1%) 15BASE' + comment] = \
-            make_assembly('Assembly (3.1%) 15BASE' + comment, universes)
-
-        # WITH 6 BURNABLE ABSORBERS N
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    gtu,   gtu,   gtu,
-                                             gtu,              gtu,
-                                           gtu, gtu,  gtu,  gtu, gtu,
-                                           gtu, gtu,  cent, gtu, gtu,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                             bas,              bas,
-                                               bas,   gtu,   bas     ]
-        univs['Assembly (3.1%) 6BAN' + comment] = \
-            make_assembly('Assembly (3.1%) 6BAN' + comment, universes)
-
-        # WITH 6 BURNABLE ABSORBERS S
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   gtu,   bas,
-                                             bas,              bas,
-                                           bas, gtu,  gtu,  gtu, bas,
-                                           gtu, gtu,  cent, gtu, gtu,
-                                           gtu, gtu,  gtu,  gtu, gtu,
-                                             gtu,              gtu,
-                                               gtu,   gtu,   gtu     ]
-        univs['Assembly (3.1%) 6BAS' + comment] = \
-            make_assembly('Assembly (3.1%) 6BAS' + comment, universes)
-
-        # WITH 6 BURNABLE ABSORBERS W
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    gtu,   gtu,   bas,
-                                             gtu,              bas,
-                                           gtu, gtu,  gtu,  gtu, bas,
-                                           gtu, gtu,  cent, gtu, gtu,
-                                           gtu, gtu,  gtu,  gtu, bas,
-                                             gtu,              bas,
-                                               gtu,   gtu,   bas     ]
-        univs['Assembly (3.1%) 6BAW' + comment] = \
-            make_assembly('Assembly (3.1%) 6BAW' + comment, universes)
-
-        # WITH 6 BURNABLE ABSORBERS E
-        universes = np.empty((17,17), dtype=openmc.Universe)
-        universes[:,:] = pins['Fuel (3.1%) stack']
-        universes[nonfuel_y, nonfuel_x] = [    bas,   gtu,   gtu,
-                                             bas,              gtu,
-                                           bas, gtu,  gtu,  gtu, gtu,
-                                           gtu, gtu,  cent, gtu, gtu,
-                                           bas, gtu,  gtu,  gtu, gtu,
-                                             bas,              gtu,
-                                               bas,   gtu,   gtu     ]
-        univs['Assembly (3.1%) 6BAE' + comment] = \
-            make_assembly('Assembly (3.1%) 6BAE' + comment, universes)
 
     return univs

smr/smr/core.py

@@ -5,18 +5,19 @@ import numpy as np
 import openmc
 
 from .materials import mats
-from .surfaces import surfs, lattice_pitch
 from .reflector import reflector_universes
 from .assemblies import assembly_universes
+from smr import surfaces
 
 
-def core_geometry(num_rings, num_axial, depleted):
+def core_geometry(ring_radii, num_axial, depleted):
     """Generate full core SMR geometry.
 
     Parameters
     ----------
-    num_rings : int
-        Number of annual regions in fuel
+    ring_radii : iterable of float
+        Radii of rings in fuel (note that this doesn't need to include the
+        full fuel pin radius)
     num_axial : int
         Number of axial subdivisions in fuel
     depleted : bool
@@ -28,11 +29,12 @@ def core_geometry(num_rings, num_axial, depleted):
         SMR full core geometry
 
     """
-    assembly = assembly_universes(num_rings, num_axial, depleted)
+    assembly = assembly_universes(ring_radii, num_axial, depleted)
     reflector = reflector_universes()
 
     # Construct main core lattice
     core = openmc.RectLattice(name='Main core')
+    lattice_pitch = surfaces.lattice_pitch
     core.lower_left = (-9*lattice_pitch/2, -9*lattice_pitch/2)
     core.pitch = (lattice_pitch, lattice_pitch)
     universes = np.tile(reflector['solid'], (9, 9))
@@ -45,67 +47,67 @@ def core_geometry(num_rings, num_axial, depleted):
 
     universes[1, 1] = reflector['1,1']
     universes[1, 2] = reflector['NW']
-    universes[1, 3] = assembly['Assembly (3.1%) instr']
-    universes[1, 4] = assembly['Assembly (2.4%) CR D']
-    universes[1, 5] = assembly['Assembly (3.1%) instr']
+    universes[1, 3] = assembly['Assembly (3.1%)']
+    universes[1, 4] = assembly['Assembly (3.1%)']
+    universes[1, 5] = assembly['Assembly (3.1%)']
     universes[1, 6] = reflector['NE']
     universes[1, 7] = reflector['1,7']
 
     universes[2, 0] = reflector['2,0']
     universes[2, 1] = reflector['NW']
-    universes[2, 2] = assembly['Assembly (3.1%) instr']
-    universes[2, 3] = assembly['Assembly (2.4%) CR D']
-    universes[2, 4] = assembly['Assembly (3.1%) 16BA']
-    universes[2, 5] = assembly['Assembly (2.4%) CR D']
-    universes[2, 6] = assembly['Assembly (3.1%) instr']
+    universes[2, 2] = assembly['Assembly (3.1%)']
+    universes[2, 3] = assembly['Assembly (2.4%)']
+    universes[2, 4] = assembly['Assembly (1.6%)']
+    universes[2, 5] = assembly['Assembly (2.4%)']
+    universes[2, 6] = assembly['Assembly (3.1%)']
     universes[2, 7] = reflector['NE']
     universes[2, 8] = reflector['2,8']
 
     universes[3, 0] = reflector['3,0']
-    universes[3, 1] = assembly['Assembly (3.1%) instr']
-    universes[3, 2] = assembly['Assembly (2.4%) CR D']
-    universes[3, 3] = assembly['Assembly (3.1%) 16BA']
-    universes[3, 4] = assembly['Assembly (2.4%) CR D']
-    universes[3, 5] = assembly['Assembly (3.1%) 16BA']
-    universes[3, 6] = assembly['Assembly (2.4%) CR D']
-    universes[3, 7] = assembly['Assembly (3.1%) instr']
+    universes[3, 1] = assembly['Assembly (3.1%)']
+    universes[3, 2] = assembly['Assembly (2.4%)']
+    universes[3, 3] = assembly['Assembly (1.6%)']
+    universes[3, 4] = assembly['Assembly (1.6%)']
+    universes[3, 5] = assembly['Assembly (1.6%)']
+    universes[3, 6] = assembly['Assembly (2.4%)']
+    universes[3, 7] = assembly['Assembly (3.1%)']
     universes[3, 8] = reflector['3,8']
 
     universes[4, 0] = reflector['4,0']
-    universes[4, 1] = assembly['Assembly (2.4%) CR D']
-    universes[4, 2] = assembly['Assembly (3.1%) 16BA']
-    universes[4, 3] = assembly['Assembly (2.4%) CR D']
-    universes[4, 4] = assembly['Assembly (1.6%) instr']
-    universes[4, 5] = assembly['Assembly (2.4%) CR D']
-    universes[4, 6] = assembly['Assembly (3.1%) 16BA']
-    universes[4, 7] = assembly['Assembly (2.4%) CR D']
+    universes[4, 1] = assembly['Assembly (3.1%)']
+    universes[4, 2] = assembly['Assembly (1.6%)']
+    universes[4, 3] = assembly['Assembly (1.6%)']
+    universes[4, 4] = assembly['Assembly (2.4%)']
+    universes[4, 5] = assembly['Assembly (1.6%)']
+    universes[4, 6] = assembly['Assembly (1.6%)']
+    universes[4, 7] = assembly['Assembly (3.1%)']
     universes[4, 8] = reflector['4,8']
 
     universes[5, 0] = reflector['5,0']
-    universes[5, 1] = assembly['Assembly (3.1%) instr']
-    universes[5, 2] = assembly['Assembly (2.4%) CR D']
-    universes[5, 3] = assembly['Assembly (3.1%) 16BA']
-    universes[5, 4] = assembly['Assembly (2.4%) CR D']
-    universes[5, 5] = assembly['Assembly (3.1%) 16BA']
-    universes[5, 6] = assembly['Assembly (2.4%) CR D']
-    universes[5, 7] = assembly['Assembly (3.1%) instr']
+    universes[5, 1] = assembly['Assembly (3.1%)']
+    universes[5, 2] = assembly['Assembly (2.4%)']
+    universes[5, 3] = assembly['Assembly (1.6%)']
+    universes[5, 4] = assembly['Assembly (1.6%)']
+    universes[5, 5] = assembly['Assembly (1.6%)']
+    universes[5, 6] = assembly['Assembly (2.4%)']
+    universes[5, 7] = assembly['Assembly (3.1%)']
     universes[5, 8] = reflector['5,8']
 
     universes[6, 0] = reflector['6,0']
     universes[6, 1] = reflector['SW']
-    universes[6, 2] = assembly['Assembly (3.1%) instr']
-    universes[6, 3] = assembly['Assembly (2.4%) CR D']
-    universes[6, 4] = assembly['Assembly (3.1%) 16BA']
-    universes[6, 5] = assembly['Assembly (2.4%) CR D']
-    universes[6, 6] = assembly['Assembly (3.1%) instr']
+    universes[6, 2] = assembly['Assembly (3.1%)']
+    universes[6, 3] = assembly['Assembly (2.4%)']
+    universes[6, 4] = assembly['Assembly (1.6%)']
+    universes[6, 5] = assembly['Assembly (2.4%)']
+    universes[6, 6] = assembly['Assembly (3.1%)']
     universes[6, 7] = reflector['SE']
     universes[6, 8] = reflector['6,8']
 
     universes[7, 1] = reflector['7,1']
     universes[7, 2] = reflector['SW']
-    universes[7, 3] = assembly['Assembly (3.1%) instr']
-    universes[7, 4] = assembly['Assembly (2.4%) CR D']
-    universes[7, 5] = assembly['Assembly (3.1%) instr']
+    universes[7, 3] = assembly['Assembly (3.1%)']
+    universes[7, 4] = assembly['Assembly (3.1%)']
+    universes[7, 5] = assembly['Assembly (3.1%)']
     universes[7, 6] = reflector['SE']
     universes[7, 7] = reflector['7,7']
 
@@ -118,6 +120,7 @@ def core_geometry(num_rings, num_axial, depleted):
     core.universes = universes
 
     root_univ = openmc.Universe(universe_id=0, name='root universe')
+    surfs = surfaces.surfs
 
     # Cylinder filled with core lattice
     cell = openmc.Cell(name='Main core')
@@ -133,75 +136,10 @@ def core_geometry(num_rings, num_axial, depleted):
                    +surfs['lower bound'] & -surfs['upper bound'])
     root_univ.add_cell(cell)
 
-    # Neutron shield panels
-    cell = openmc.Cell(name='neutron shield panel NW')
-    cell.fill = mats['SS']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   +surfs['neutron shield NWbot SEtop'] &
-                   -surfs['neutron shield NWtop SEbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel N')
-    cell.fill = mats['H2O']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   +surfs['neutron shield NWtop SEbot'] &
-                   -surfs['neutron shield NEtop SWbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel SE')
-    cell.fill = mats['SS']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   -surfs['neutron shield NWbot SEtop'] &
-                   +surfs['neutron shield NWtop SEbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel E')
-    cell.fill = mats['H2O']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   +surfs['neutron shield NWbot SEtop'] &
-                   +surfs['neutron shield NEbot SWtop'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel NE')
-    cell.fill = mats['SS']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   +surfs['neutron shield NEbot SWtop'] &
-                   -surfs['neutron shield NEtop SWbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel S')
-    cell.fill = mats['H2O']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   -surfs['neutron shield NWtop SEbot'] &
-                   +surfs['neutron shield NEtop SWbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel SW')
-    cell.fill = mats['SS']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   -surfs['neutron shield NEbot SWtop'] &
-                   +surfs['neutron shield NEtop SWbot'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
-    cell = openmc.Cell(name='neutron shield panel W')
-    cell.fill = mats['H2O']
-    cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-                   -surfs['neutron shield NWbot SEtop'] &
-                   -surfs['neutron shield NEbot SWtop'] &
-                   +surfs['lower bound'] & -surfs['upper bound'])
-    root_univ.add_cell(cell)
-
     # Downcomer
     cell = openmc.Cell(name='downcomer')
     cell.fill = mats['H2O']
-    cell.region = (+surfs['neutron shield OR'] & -surfs['RPV IR'] &
+    cell.region = (+surfs['core barrel OR'] & -surfs['RPV IR'] &
                    +surfs['lower bound'] & -surfs['upper bound'])
     root_univ.add_cell(cell)
 

smr/smr/materials.py

@@ -1,5 +1,7 @@
 """Instantiate the OpenMC Materials needed by the core model."""
 
+import copy
+
 import openmc
 from openmc.data import atomic_weight, atomic_mass, water_density
 
@@ -63,6 +65,15 @@ mats['In'].add_element('Mn', 0.0087, 'wo')
 mats['In'].add_element('Fe', 0.2863, 'wo')
 mats['In'].add_element('Ni', 0.5119, 'wo')
 
+# Create stainless steel 302
+mats['SS302'] = openmc.Material(name='SS302')
+mats['SS302'].set_density('g/cm3', 7.86)
+mats['SS302'].add_element('Si', 0.01, 'wo')
+mats['SS302'].add_element('Cr', 0.18, 'wo')
+mats['SS302'].add_element('Mn', 0.02, 'wo')
+mats['SS302'].add_element('Fe', 0.70, 'wo')
+mats['SS302'].add_element('Ni', 0.09, 'wo')
+
 # Create stainless steel material
 mats['SS'] = openmc.Material(name='SS304')
 mats['SS'].set_density('g/cc', 8.03)
@@ -183,7 +194,6 @@ aB_bsg = aB10_bsg + aB11_bsg
 
 # Create borosilicate glass material
 mats['BSG'] = openmc.Material(name='Borosilicate Glass')
-mats['BSG'].temperature = 300
 mats['BSG'].set_density('g/cc', 2.26)
 mats['BSG'].add_element('O', aO_bsg, 'ao')
 mats['BSG'].add_element('Si', aSi_bsg, 'ao')
@@ -196,7 +206,6 @@ mats['BSG'].add_nuclide('B11', aB11_bsg, 'ao')
 
 # Create 1.6% enriched UO2 fuel material
 mat = openmc.Material(name='1.6% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.31341)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=1.61006)
@@ -204,7 +213,6 @@ mats['UO2 1.6 fresh'] = mat
 
 # Create 2.4% enriched UO2 fuel material
 mat = openmc.Material(name='2.4% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.29748)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=2.39993)
@@ -212,7 +220,6 @@ mats['UO2 2.4 fresh'] = mat
 
 # Create 3.1% enriched UO2 fuel material
 mat = openmc.Material(name='3.1% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.30166)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=3.10221)
@@ -220,7 +227,6 @@ mats['UO2 3.1 fresh'] = mat
 
 # Depleted versions of 1.6%, 2.4%, 3.1% fuel
 mat = openmc.Material(name='2.4% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.29748)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=2.39993)
@@ -229,7 +235,6 @@ for nuc in _DEPLETION_NUCLIDES:
 mats['UO2 2.4 depleted'] = mat
 
 mat = openmc.Material(name='1.6% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.31341)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=1.61006)
@@ -238,7 +243,6 @@ for nuc in _DEPLETION_NUCLIDES:
 mats['UO2 1.6 depleted'] = mat
 
 mat = openmc.Material(name='3.1% Enr. UO2 Fuel')
-mat.temperature = 300
 mat.set_density('g/cc', 10.30166)
 mat.add_element('O', 2., 'ao')
 mat.add_element('U', 1., 'ao', enrichment=3.10221)
@@ -250,3 +254,10 @@ mats['UO2 3.1 depleted'] = mat
 # Construct a collection of Materials to export to XML
 
 materials = openmc.Materials(mats.values())
+
+
+def clone(material):
+    """Perform copy of material but share nuclide densities"""
+    shared_mat = copy.copy(material)
+    shared_mat.id = None
+    return shared_mat

smr/smr/pins.py

@@ -142,13 +142,14 @@ def make_pin_stack(name, zsurfaces, universes, boundary, fuel_fill):
     return universe
 
 
-def pin_universes(num_rings=10, num_axial=196, depleted=False):
+def pin_universes(ring_radii=None, num_axial=196, depleted=False):
     """Generate universes for SMR fuel pins.
 
     Parameters
     ----------
-    num_rings : int
-        Number of annual regions in fuel
+    ring_radii : iterable of float
+        Radii of rings in fuel (note that this doesn't need to include the
+        full fuel pin radius)
     num_axial : int
         Number of axial subdivisions in fuel
     depleted : bool
@@ -230,8 +231,8 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
         surfs['top upper nozzle']
     ]
 
-    univs['GT empty'] = make_stack(
-        'GT empty', surfaces=stack_surfs,
+    univs['GT empty stack'] = make_stack(
+        'GT empty stack', surfaces=stack_surfs,
         universes=[univs['water pin'],
                    univs['water pin'],
                    univs['water pin'],
@@ -472,199 +473,6 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
                        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['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']])
-
-
     # Fuel pin cells
     univs['SS pin'] = make_pin(
         'SS pin',
@@ -674,16 +482,16 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
     univs['end plug'] = make_pin(
         'end plug',
         [surfs['clad OR']],
-        [mats['Zr'], mats['H2O']])
+        [mats['M5'], mats['H2O']])
 
     univs['pin plenum'] = make_pin(
         'pin plenum',
         surfaces=[surfs['plenum spring OR'],
                   surfs['clad IR'],
                   surfs['clad OR']],
-        materials=[mats['In'],
+        materials=[mats['SS302'],
                    mats['He'],
-                   mats['Zr'],
+                   mats['M5'],
                    mats['H2O']])
 
     univs['pin plenum grid (intermediate)'] = make_pin(
@@ -706,12 +514,11 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
         axial_splits = np.linspace(bottom_fuel_stack, top_active_core, num_axial + 1)[1:-1]
         axial_surfs = [openmc.ZPlane(z0=z) for z in axial_splits]
 
-    if num_rings > 1:
+    if ring_radii is not None:
         # Get z-cylinder surfaces for each ring
         rings = []
-        for i in range(1, num_rings):
-            R = sqrt(i*pellet_OR**2/num_rings)
-            cyl = openmc.ZCylinder(r=R, name='fuel ring {}'.format(i))
+        for i, r in enumerate(ring_radii):
+            cyl = openmc.ZCylinder(r=r, name='fuel ring {}'.format(i))
             rings.append(cyl)
 
     def subdivided_fuel(fill):
@@ -719,14 +526,14 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
         uo2_cells = []
         if num_axial > 1:
             for axial_region in subdivide(axial_surfs):
-                if num_rings > 1:
+                if ring_radii is not None:
                     for ring_region in subdivide(rings):
                         cell = openmc.Cell(fill=fill, region=axial_region & ring_region)
                         uo2_cells.append(cell)
                 else:
                     uo2_cells.append(openmc.Cell(fill=fill, region=axial_region))
         else:
-            if num_rings > 1:
+            if ring_radii is not None:
                 for ring_region in subdivide(rings):
                     cell = openmc.Cell(fill=fill, region=ring_region)
                     uo2_cells.append(cell)
@@ -737,7 +544,7 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
 
     # If rings/axial segments are present, create a universe for the subdivided
     # fuel. Otherwise just use a plain material.
-    if num_rings > 1 or num_axial > 1:
+    if ring_radii is not None or num_axial > 1:
         fuel_fill = subdivided_fuel(mats['UO2 1.6 {}'.format(fuel)])
     else:
         fuel_fill = mats['UO2 1.6 {}'.format(fuel)]
@@ -762,6 +569,12 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
         materials=outside_pin_mats,
         grid='intermediate')
 
+    univs['Fuel pin (1.6%) no grid'] = make_pin(
+        'Pin no grid',
+        surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
+        materials=[fuel_fill] + outside_pin_mats
+    )
+
     # Stack all axial pieces of 1.6% enriched fuel pin cell
 
     within_fuel_surfs = [
@@ -829,11 +642,17 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
 
     # If rings/axial segments are present, create a universe for the subdivided
     # fuel. Otherwise just use a plain material.
-    if num_rings > 1 or num_axial > 1:
+    if ring_radii is not None or num_axial > 1:
         fuel_fill = subdivided_fuel(mats['UO2 2.4 {}'.format(fuel)])
     else:
         fuel_fill = mats['UO2 2.4 {}'.format(fuel)]
 
+    univs['Fuel pin (2.4%) no grid'] = make_pin(
+        'Pin no grid',
+        surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
+        materials=[fuel_fill] + outside_pin_mats
+    )
+
     # Stack all axial pieces of 2.4% enriched fuel pin cell
 
     univs['Fuel pin (2.4%) stack'] = make_pin_stack(
@@ -875,11 +694,26 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
 
     # If rings/axial segments are present, create a universe for the subdivided
     # fuel. Otherwise just use a plain material.
-    if num_rings > 1 or num_axial > 1:
+    if ring_radii is not None or num_axial > 1:
         fuel_fill = subdivided_fuel(mats['UO2 3.1 {}'.format(fuel)])
     else:
         fuel_fill = mats['UO2 3.1 {}'.format(fuel)]
 
+    if num_axial > 1:
+        water_cells = []
+        for i, r in enumerate(subdivide(axial_surfs)):
+            cell = openmc.Cell(fill=mats['H2O'], region=r, name=f'Water ({i})')
+            water_cells.append(cell)
+        water_fill = openmc.Universe(cells=water_cells)
+    else:
+        water_fill = mats['H2O']
+
+    univs['Fuel pin (3.1%) no grid'] = make_pin(
+        'Pin no grid',
+        surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
+        materials=[fuel_fill, mats['He'], mats['M5'], water_fill]
+    )
+
     # Stack all axial pieces of 3.1% enriched fuel pin cell
 
     univs['Fuel pin (3.1%) stack'] = make_pin_stack(

smr/smr/reflector.py

@@ -12,7 +12,7 @@ converted to actual dimensions by scaling according to the width of an assembly.
 import openmc
 
 from .materials import mats
-from .surfaces import lattice_pitch
+from smr import surfaces
 
 
 def make_reflector(name, parameters):
@@ -91,6 +91,7 @@ def reflector_universes():
 
     # All pixel widths are scaled according to the actual width of an assembly
     # divided by the width of an assembly in pixels
+    lattice_pitch = surfaces.lattice_pitch
     scale = lattice_pitch/width
 
     # Physical positions

smr/smr/surfaces.py

@@ -13,7 +13,6 @@ NuScale DC application, chapter 4: https://www.nrc.gov/docs/ML1701/ML17013A274.p
 
 """
 
-import copy
 from math import tan, pi
 
 import numpy as np
@@ -40,7 +39,7 @@ clad_OR            = 0.374*INCHES/2   # ML17013A274, Table 4.1-2
 active_fuel_length = 78.74*INCHES     # ML17013A274, Figure 4.2-10
 plenum_length      = 5.311*INCHES     # ML17013A274, Figure 4.2-10
 fuel_rod_length    = 85.00*INCHES     # ML17013A274, Table 4.1-2
-lower_end_cap      = 0.575*INCHES     # ML17007A001, Table 3-2
+lower_end_cap_length = 0.575*INCHES   # ML17007A001, Table 3-2
 
 # pin cell parameters
 guide_tube_IR      = 0.450*INCHES/2  # ML17013A274, Table 4.1-2
@@ -79,19 +78,20 @@ top_nozzle_width  = 8.406*INCHES  # ML17013A274, Figure 4.2-2
 core_barrel_IR     = 74*INCHES/2  # ML17013A274, Table 4.1-2
 core_barrel_OR     = 78*INCHES/2  # ML17013A274, Table 4.1-2
 neutron_shield_OR  = core_barrel_OR + 2.0
-rpv_IR             = 120.0  # Estimate?
-rpv_OR             = 135.0  # Estimate?
+rpv_IR             = 96.5*INCHES/2  # ML17013A274, Table 5.3-1
+rpv_OR             = 105*INCHES/2   # ML17013A274, Table 5.3-1
 
 # axial parameters
-lowest_extent        =      0.000
-bottom_support_plate =     20.000
-top_support_plate    =     25.000
-bottom_lower_nozzle  =     25.000
-top_lower_nozzle     =     35.160
-bottom_fuel_rod      =     35.160
-top_lower_thimble    =     36.007
-bottom_fuel_stack    =     36.007
-bot_burn_abs         =     41.087
+reference_z = -36.6205
+lowest_extent        = reference_z
+bottom_support_plate = lowest_extent + 20.000
+top_support_plate    = bottom_support_plate + 5.000
+bottom_lower_nozzle  = bottom_support_plate + 5.000
+top_lower_nozzle     = bottom_lower_nozzle + 4.0*INCHES
+bottom_fuel_rod      = bottom_lower_nozzle + 4.0*INCHES
+top_lower_thimble    = bottom_fuel_rod + lower_end_cap_length
+bottom_fuel_stack    = bottom_fuel_rod + lower_end_cap_length
+bot_burn_abs         = bottom_fuel_stack + 2.0*INCHES
 top_active_core = bottom_fuel_stack + active_fuel_length
 top_plenum = top_active_core + plenum_length
 top_fuel_rod = bottom_fuel_rod + fuel_rod_length
@@ -127,74 +127,47 @@ neutron_shield_NEtop_SWbot = tan(-pi/6)
 
 surfs = {}
 
-surfs['pellet OR'] = openmc.ZCylinder(
-    r=pellet_OR, name='Pellet OR')
-surfs['plenum spring OR'] = openmc.ZCylinder(
-    r=plenum_spring_OR, name='FR Plenum Spring OR')
-surfs['clad IR'] = openmc.ZCylinder(
-    r=clad_IR, name='Clad IR')
-surfs['clad OR'] = openmc.ZCylinder(
-    r=clad_OR, name='Clad OR')
-surfs['GT IR'] = openmc.ZCylinder(
-    r=guide_tube_IR, name='GT IR (above dashpot)')
-surfs['GT OR'] = openmc.ZCylinder(
-    r=guide_tube_OR, name='GT OR (above dashpot)')
-surfs['GT dashpot IR'] = openmc.ZCylinder(
-    r=guide_tube_dash_IR, name='GT IR (at dashpot)')
-surfs['GT dashpot OR'] = openmc.ZCylinder(
-    r=guide_tube_dash_OR, name='GT OR (at dashpot)')
-surfs['CP OR'] = openmc.ZCylinder(
-    r=boron_carbide_OR, name='Control Poison OR')
-surfs['CR IR'] = openmc.ZCylinder(
-    r=control_rod_IR, name='CR Clad IR')
-surfs['CR OR'] = openmc.ZCylinder(
-    r=control_rod_OR, name='CR Clad OR')
-surfs['BA IR 1'] = openmc.ZCylinder(
-    r=burn_abs_r1, name='BA IR 1')
-surfs['BA IR 2'] = openmc.ZCylinder(
-    r=burn_abs_r2, name='BA IR 2')
-surfs['BA IR 3'] = openmc.ZCylinder(
-    r=burn_abs_r3, name='BA IR 3')
-surfs['BA IR 4'] = openmc.ZCylinder(
-    r=burn_abs_r4, name='BA IR 4')
-surfs['BA IR 5'] = openmc.ZCylinder(
-    r=burn_abs_r5, name='BA IR 5')
-surfs['BA IR 6'] = openmc.ZCylinder(
-    r=burn_abs_r6, name='BA IR 6')
-surfs['BA IR 7'] = openmc.ZCylinder(
-    r=burn_abs_r7, name='BA IR 7')
-surfs['BA IR 8'] = openmc.ZCylinder(
-    r=burn_abs_r8, name='BA IR 8')
-surfs['IT IR'] = copy.deepcopy(surfs['BA IR 5'])
-surfs['IT OR'] = copy.deepcopy(surfs['BA IR 6'])
+surfs['pellet OR'] = openmc.ZCylinder(r=pellet_OR, name='Pellet OR')
+surfs['plenum spring OR'] = openmc.ZCylinder(r=plenum_spring_OR, name='FR Plenum Spring OR')
+surfs['clad IR'] = openmc.ZCylinder(r=clad_IR, name='Clad IR')
+surfs['clad OR'] = openmc.ZCylinder(r=clad_OR, name='Clad OR')
+surfs['GT IR'] = openmc.ZCylinder(r=guide_tube_IR, name='GT IR (above dashpot)')
+surfs['GT OR'] = openmc.ZCylinder(r=guide_tube_OR, name='GT OR (above dashpot)')
+surfs['GT dashpot IR'] = openmc.ZCylinder(r=guide_tube_dash_IR, name='GT IR (at dashpot)')
+surfs['GT dashpot OR'] = openmc.ZCylinder(r=guide_tube_dash_OR, name='GT OR (at dashpot)')
+surfs['CP OR'] = openmc.ZCylinder(r=boron_carbide_OR, name='Control Poison OR')
+surfs['CR IR'] = openmc.ZCylinder(r=control_rod_IR, name='CR Clad IR')
+surfs['CR OR'] = openmc.ZCylinder(r=control_rod_OR, name='CR Clad OR')
+surfs['BA IR 1'] = openmc.ZCylinder(r=burn_abs_r1, name='BA IR 1')
+surfs['BA IR 2'] = openmc.ZCylinder(r=burn_abs_r2, name='BA IR 2')
+surfs['BA IR 3'] = openmc.ZCylinder(r=burn_abs_r3, name='BA IR 3')
+surfs['BA IR 4'] = openmc.ZCylinder(r=burn_abs_r4, name='BA IR 4')
+surfs['BA IR 5'] = openmc.ZCylinder(r=burn_abs_r5, name='BA IR 5')
+surfs['BA IR 6'] = openmc.ZCylinder(r=burn_abs_r6, name='BA IR 6')
+surfs['BA IR 7'] = openmc.ZCylinder(r=burn_abs_r7, name='BA IR 7')
+surfs['BA IR 8'] = openmc.ZCylinder(r=burn_abs_r8, name='BA IR 8')
+surfs['IT IR'] = surfs['BA IR 5']
+surfs['IT OR'] = surfs['BA IR 6']
 
 # Rectangular prisms for grid spacers
-surfs['rod grid box'] = \
-    openmc.get_rectangular_prism(rod_grid_side, rod_grid_side)
+surfs['rod grid box'] = openmc.rectangular_prism(rod_grid_side, rod_grid_side)
 
 # Rectangular prisms for lattice grid sleeves
-surfs['lat grid box inner'] = \
-    openmc.get_rectangular_prism(17.*pin_pitch, 17.*pin_pitch)
-surfs['lat grid box outer'] = \
-    openmc.get_rectangular_prism(grid_strap_side, grid_strap_side)
+surfs['lat grid box inner'] = openmc.rectangular_prism(17.*pin_pitch, 17.*pin_pitch)
+surfs['lat grid box outer'] = openmc.rectangular_prism(grid_strap_side, grid_strap_side)
 
-surfs['bot support plate'] = openmc.ZPlane(
-    z0=bottom_support_plate, name='bot support plate')
-surfs['top support plate'] = openmc.ZPlane(
-    z0=top_support_plate, name='top support plate')
+surfs['bot support plate'] = openmc.ZPlane(z0=bottom_support_plate, name='bot support plate')
+surfs['top support plate'] = openmc.ZPlane(z0=top_support_plate, name='top support plate')
 surfs['bottom FR'] = openmc.ZPlane(z0=bottom_fuel_rod, name='bottom FR')
-surfs['top lower nozzle'] = copy.deepcopy(surfs['bottom FR'])
-surfs['bot lower nozzle'] = copy.deepcopy(surfs['top support plate'])
+surfs['top lower nozzle'] = surfs['bottom FR']
+surfs['bot lower nozzle'] = surfs['top support plate']
 
 # axial surfaces
-surfs['bot active core'] = openmc.ZPlane(
-    z0=bottom_fuel_stack, name='bot active core')
-surfs['top active core'] = openmc.ZPlane(
-    z0=top_active_core, name='top active core')
+surfs['bot active core'] = openmc.ZPlane(z0=bottom_fuel_stack, name='bot active core')
+surfs['top active core'] = openmc.ZPlane(z0=top_active_core, name='top active core')
 
-surfs['top lower thimble'] = copy.deepcopy(surfs['bot active core'])
-surfs['BA bot'] = openmc.ZPlane(
-    z0=bot_burn_abs, name='bottom of BA')
+surfs['top lower thimble'] = surfs['bot active core']
+surfs['BA bot'] = openmc.ZPlane(z0=bot_burn_abs, name='bottom of BA')
 
 for i, (bottom, top) in enumerate(zip(grid_bottom, grid_top)):
     # Create plane for bottom of spacer grid
@@ -207,17 +180,12 @@ for i, (bottom, top) in enumerate(zip(grid_bottom, grid_top)):
     name = 'top of grid {}'.format(i + 1)
     surfs[key] = openmc.ZPlane(z0=top, name=name)
 
-surfs['dashpot top'] = openmc.ZPlane(
-    z0=step0H, name='top dashpot')
+surfs['dashpot top'] = openmc.ZPlane(z0=step0H, name='top dashpot')
 
-surfs['top pin plenum'] = openmc.ZPlane(
-    z0=top_plenum, name='top pin plenum')
-surfs['top FR'] = openmc.ZPlane(
-    z0=top_fuel_rod, name='top FR')
-surfs['bot upper nozzle'] = openmc.ZPlane(
-    z0=bottom_upper_nozzle, name='bottom upper nozzle')
-surfs['top upper nozzle'] = openmc.ZPlane(
-    z0=top_upper_nozzle, name='top upper nozzle')
+surfs['top pin plenum'] = openmc.ZPlane(z0=top_plenum, name='top pin plenum')
+surfs['top FR'] = openmc.ZPlane(z0=top_fuel_rod, name='top FR')
+surfs['bot upper nozzle'] = openmc.ZPlane(z0=bottom_upper_nozzle, name='bottom upper nozzle')
+surfs['top upper nozzle'] = openmc.ZPlane(z0=top_upper_nozzle, name='top upper nozzle')
 
 # Control rod bank surfaces for ARO configuration
 for bank in ['A','B','C','D','E',]:
@@ -226,30 +194,19 @@ for bank in ['A','B','C','D','E',]:
     surfs['bankS{} bot'.format(bank)] = openmc.ZPlane(
         z0=step248H, name='CR bankS{} bottom'.format(bank))
 
-surfs['bankA top'] = openmc.ZPlane(
-    z0=bank_top, name='CR bank A top')
-surfs['bankA bot'] = openmc.ZPlane(
-    z0=bank_bot, name='CR bank A bottom')
-surfs['bankB top'] = openmc.ZPlane(
-    z0=bank_top, name='CR bank B top')
-surfs['bankB bot'] = openmc.ZPlane(
-    z0=bank_bot, name='CR bank B bottom')
-surfs['bankC top'] = openmc.ZPlane(
-    z0=bank_top, name='CR bank C top')
-surfs['bankC bot'] = openmc.ZPlane(
-    z0=bank_bot, name='CR bank C bottom')
-surfs['bankD top'] = openmc.ZPlane(
-    z0=bank_top, name='CR bank D top')
-surfs['bankD bot'] = openmc.ZPlane(
-    z0=bank_bot, name='CR bank D bottom')
+surfs['bankA top'] = openmc.ZPlane(z0=bank_top, name='CR bank A top')
+surfs['bankA bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank A bottom')
+surfs['bankB top'] = openmc.ZPlane(z0=bank_top, name='CR bank B top')
+surfs['bankB bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank B bottom')
+surfs['bankC top'] = openmc.ZPlane(z0=bank_top, name='CR bank C top')
+surfs['bankC bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank C bottom')
+surfs['bankD top'] = openmc.ZPlane(z0=bank_top, name='CR bank D top')
+surfs['bankD bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank D bottom')
 
 # outer radial surfaces
-surfs['core barrel IR'] = openmc.ZCylinder(
-    r=core_barrel_IR, name='core barrel IR')
-surfs['core barrel OR'] = openmc.ZCylinder(
-    r=core_barrel_OR, name='core barrel OR')
-surfs['neutron shield OR'] = openmc.ZCylinder(
-    r=neutron_shield_OR, name='neutron shield OR')
+surfs['core barrel IR'] = openmc.ZCylinder(r=core_barrel_IR, name='core barrel IR')
+surfs['core barrel OR'] = openmc.ZCylinder(r=core_barrel_OR, name='core barrel OR')
+surfs['neutron shield OR'] = openmc.ZCylinder(r=neutron_shield_OR, name='neutron shield OR')
 
 # neutron shield planes
 surfs['neutron shield NWbot SEtop'] = openmc.Plane(
@@ -266,10 +223,8 @@ surfs['neutron shield NEtop SWbot'] = openmc.Plane(
     name='neutron shield NEtop SWbot')
 
 # outer radial surfaces
-surfs['RPV IR'] = openmc.ZCylinder(
-    r=rpv_IR, name='RPV IR')
-surfs['RPV OR'] = openmc.ZCylinder(
-    r=rpv_OR, name='RPV OR', boundary_type='vacuum')
+surfs['RPV IR'] = openmc.ZCylinder(r=rpv_IR, name='RPV IR')
+surfs['RPV OR'] = openmc.ZCylinder(r=rpv_OR, name='RPV OR', boundary_type='vacuum')
 
 # outer axial surfaces
 surfs['upper bound'] = openmc.ZPlane(