Fix subdivided fuel (position of lower z-plane, add outer universe)

smr/build-core-fresh.py

@@ -15,6 +15,13 @@ from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core
 from smr.core import core_geometry
 
 
+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',
@@ -53,7 +60,7 @@ if args.tallies == 'mat':
     for cell in geometry.get_all_cells().values():
         if cell.fill in fuel_mats:
             # Fill cell with list of "differentiated" materials
-            cell.fill = [cell.fill.clone() for i in range(cell.num_instances)]
+            cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
 
 #### Create OpenMC "materials.xml" file
 all_materials = geometry.get_all_materials()

smr/smr/pins.py

@@ -7,7 +7,7 @@ import openmc
 from openmc.model import subdivide
 
 from .materials import mats
-from .surfaces import surfs, pellet_OR, bottom_fuel_rod, top_active_core
+from .surfaces import surfs, pellet_OR, bottom_fuel_stack, top_active_core
 
 
 def make_pin(name, surfaces, materials, grid=None):
@@ -700,12 +700,6 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
 
     #### 1.6% ENRICHED FUEL PIN CELL
 
-    if num_axial > 1:
-        # Determine z position between each fuel pellet, omitting the surfaces
-        # corresponding to the very bottom and top of the active fuel length
-        axial_splits = np.linspace(bottom_fuel_rod, top_active_core, num_axial + 1)[1:-1]
-        axial_surfs = [openmc.ZPlane(z0=z) for z in axial_splits]
-
     if num_rings > 1:
         # Get z-cylinder surfaces for each ring
         rings = []
@@ -715,25 +709,28 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
             rings.append(cyl)
 
     def subdivided_fuel(fill):
-        # Create universe for UO2 alone with axial/radial subdivision
-        uo2_cells = []
-        if num_axial > 1:
-            for axial_region in subdivide(axial_surfs):
-                if num_rings > 1:
-                    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))
+        """Create universe for UO2 alone with axial/radial subdivision"""
+        # Create universe with radial rings
+        univ = openmc.Universe()
+        if num_rings > 1:
+            for ring_region in subdivide(rings):
+                cell = openmc.Cell(fill=fill, region=ring_region)
+                univ.add_cell(cell)
         else:
-            if num_rings > 1:
-                for ring_region in subdivide(rings):
-                    cell = openmc.Cell(fill=fill, region=ring_region)
-                    uo2_cells.append(cell)
-            else:
-                raise RuntimeError("Shouldn't call with 1 ring and 1 axial segment")
+            univ.add_cell(openmc.Cell(fill=fill))
 
-        return openmc.Universe(cells=uo2_cells)
+        # If multiple axial segments are needed, use a lattice to avoid
+        # performance cost of neighbor lookups
+        if num_axial > 1:
+            zlattice = openmc.RectLattice()
+            zlattice.lower_left = (-100., -100., bottom_fuel_stack)
+            zlattice.pitch = (200., 200., (top_active_core - bottom_fuel_stack)/num_axial)
+            zlattice.dimension = (1, 1, num_axial)
+            zlattice.universes = np.full((num_axial, 1, 1), univ)
+            zlattice.outer = univ
+            return zlattice
+        else:
+            return univ
 
     # If rings/axial segments are present, create a universe for the subdivided
     # fuel. Otherwise just use a plain material.

smr/smr/surfaces.py

@@ -192,7 +192,7 @@ surfs['bot active core'] = openmc.ZPlane(
 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['top lower thimble'] = surfs['bot active core']
 surfs['BA bot'] = openmc.ZPlane(
     z0=bot_burn_abs, name='bottom of BA')