Add script for short single assembly ExaSMR coupled runs

This commit is contained in:
Paul Romano 2019-06-05 14:08:49 -05:00
parent 9ea73dbce3
commit ec21e23981
4 changed files with 221 additions and 27 deletions

174
smr/build-assembly-short.py Normal file
View file

@ -0,0 +1,174 @@
#!/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, mats
from smr.surfaces import surfs, lattice_pitch, pin_pitch, bottom_fuel_stack, \
top_active_core
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=196,
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 ]
xy_bounds = openmc.model.get_rectangular_prism(
lattice_pitch, lattice_pitch, boundary_type='reflective')
z_bounds = +surfs['bot active core'] & -surfs['top active core']
surfs['bot active core'].boundary_type = 'reflective'
surfs['top active core'].boundary_type = 'reflective'
# 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
cell.region = surfs['lat grid box inner'] & z_bounds
root_universe.add_cell(cell)
# Add outer water cell
cell = openmc.Cell(name='outer water')
cell.fill = mats['H2O']
cell.region = ~surfs['lat grid box inner'] & xy_bounds & z_bounds
root_universe.add_cell(cell)
# 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
if args.tallies == 'mat':
# Count the number of instances for each cell and material
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}
for cell in tqdm(geometry.get_all_material_cells().values(),
desc='Differentiating materials'):
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)]
#### 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'))

View file

@ -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 = {}

View file

@ -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,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 3.1 {}'.format(fuel)])
else:
fuel_fill = mats['UO2 3.1 {}'.format(fuel)]
univs['Fuel pin (3.1%) 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 3.1% enriched fuel pin cell
univs['Fuel pin (3.1%) stack'] = make_pin_stack(
@ -907,7 +924,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'],

View file

@ -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