Assign volumes to all materials

smr/build-assembly-short.py

@@ -0,0 +1,192 @@
+#!/usr/bin/env python3
+
+import argparse
+import copy
+from math import pi, isclose
+from pathlib import Path
+
+import numpy as np
+from tqdm import tqdm
+import openmc
+
+from smr.materials import materials, mats
+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('-m', '--multipole', action='store_true',
+                    help='Whether to 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)
+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)
+
+
+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
+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}
+
+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/smr/assemblies.py

@@ -120,7 +120,7 @@ def make_assembly(name, universes):
     return universe
 
 
-def assembly_universes(num_rings, num_axial, depleted):
+def assembly_universes(rings, num_axial, depleted):
     """Generate universes for SMR fuel assemblies.
 
     Parameters
@@ -138,7 +138,7 @@ 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(rings, num_axial, depleted)
 
     # Create dictionary to store assembly universes
     univs = {}

smr/smr/pins.py

@@ -142,7 +142,7 @@ 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
@@ -230,8 +230,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'],
@@ -706,12 +706,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 +718,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 +736,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 +761,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 = [
@@ -815,7 +820,7 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
                    univs['SS pin'],
                    univs['SS pin'],
                    univs['end plug'],
-                   univs['Fuel pin (1.6%) stack'],
+                   univs['Fuel pin (1.6%) no grid'],
                    univs['pin plenum'],
                    univs['pin plenum grid (intermediate)'],
                    univs['pin plenum'],
@@ -829,11 +834,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(
@@ -861,7 +872,7 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
                    univs['SS pin'],
                    univs['SS pin'],
                    univs['end plug'],
-                   univs['Fuel pin (2.4%) stack'],
+                   univs['Fuel pin (2.4%) no grid'],
                    univs['pin plenum'],
                    univs['pin plenum grid (intermediate)'],
                    univs['pin plenum'],
@@ -875,11 +886,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(
@@ -907,7 +933,7 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False):
                    univs['SS pin'],
                    univs['SS pin'],
                    univs['end plug'],
-                   univs['Fuel pin (3.1%) stack'],
+                   univs['Fuel pin (3.1%) no grid'],
                    univs['pin plenum'],
                    univs['pin plenum grid (intermediate)'],
                    univs['pin plenum'],

smr/smr/surfaces.py

@@ -37,7 +37,7 @@ pellet_OR          = 0.3195*INCHES/2  # ML17013A274, Table 4.1-2
 pellet_length      = 0.4*INCHES       # ML17013A274, Table 4.1-2
 clad_IR            = 0.326*INCHES/2   # ML17013A274, Table 4.1-2
 clad_OR            = 0.374*INCHES/2   # ML17013A274, Table 4.1-2
-active_fuel_length = 78.74*INCHES     # ML17013A274, Figure 4.2-10
+#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
@@ -75,6 +75,8 @@ grid_strap_side   = 21.47270
 top_nozzle_height = 3.551*INCHES  # ML17013A274, Figure 4.2-2
 top_nozzle_width  = 8.406*INCHES  # ML17013A274, Figure 4.2-2
 
+active_fuel_length  = 10.0*pin_pitch
+
 # core radial parameters
 core_barrel_IR     = 74*INCHES/2  # ML17013A274, Table 4.1-2
 core_barrel_OR     = 78*INCHES/2  # ML17013A274, Table 4.1-2
@@ -83,15 +85,16 @@ rpv_IR             = 120.0  # Estimate?
 rpv_OR             = 135.0  # Estimate?
 
 # 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.007
+lowest_extent        =      0.000 - _reference_z
+bottom_support_plate =     20.000 - _reference_z
+top_support_plate    =     25.000 - _reference_z
+bottom_lower_nozzle  =     25.000 - _reference_z
+top_lower_nozzle     =     35.160 - _reference_z
+bottom_fuel_rod      =     35.160 - _reference_z
+top_lower_thimble    =     36.007 - _reference_z
+bottom_fuel_stack    =     36.007 - _reference_z
+bot_burn_abs         =     41.087 - _reference_z
 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