Updated depletion scripts to reflect latest OpenMC PyAPI and OpenDeplete

This commit is contained in:
Will Boyd 2017-03-07 18:04:29 -05:00
parent 3a30feea76
commit 377ceb928c
14 changed files with 131 additions and 1431 deletions

View file

@ -1,4 +1,4 @@
from collections import OrderedDict, defaultdict
from collections import OrderedDict
import copy
import os
@ -9,56 +9,7 @@ import opendeplete
from geometry import beavrs, openmc_geometry
#### Create "dummy" inputs to export distribcell paths for burnable cells
# Create OpenMC "materials.xml" file
beavrs.write_openmc_materials()
# Create OpenMC "geometry.xml" file
openmc_geometry.export_to_xml()
# Construct uniform initial source distribution over fissionable zones
lower_left = [-21.41728, -21.41728, +192.5]
upper_right = [+21.41728, +21.41728, +197.5]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
# Create OpenMC "settings.xml" file
settings_file = openmc.Settings()
settings_file.batches = 2
settings_file.inactive = 1
settings_file.particles = 10
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
settings_file.export_to_xml()
# Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Instantiate a "dummy" distribcell tally for each cell we wish to deplete
for cell in fuel_cells:
tally = openmc.Tally(name='dummy distribcell tally')
distribcell_filter = openmc.DistribcellFilter([cell.id])
tally.filters = [distribcell_filter]
tally.scores = ['fission']
tallies.append(tally)
tallies.export_to_xml()
# Run OpenMC to generate summary.h5 file
openmc.run()
# Open "summary.h5" file
su = openmc.Summary('summary.h5')
fuel_cells = su.openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
#### Setup OpenDeplete Materials wrapper
# Setup OpenDeplete Materials wrapper
materials = opendeplete.Materials()
materials.temperature = OrderedDict()
materials.sab = OrderedDict()
@ -66,10 +17,17 @@ materials.initial_density = OrderedDict()
materials.burn = OrderedDict()
materials.cross_sections = os.environ["OPENMC_CROSS_SECTIONS"]
# Count the number of instances for each cell and material
openmc_geometry.determine_paths()
# Extract all cells filled by a fuel material
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Extract cell materials, temperatures and sab
for cell in su.openmc_geometry.get_all_material_cells():
for cell in openmc_geometry.get_all_material_cells().values():
materials.burn[cell.name] = 'fuel' in cell.fill.name.lower()
materials.temperature[cell.name] = cell.temperature[0]
materials.temperature[cell.name] = 300
if len(cell.fill._sab) > 0:
materials.sab[cell.name] = cell.fill._sab[0]
@ -84,9 +42,8 @@ for cell in fuel_cells:
densities[nuclide][1] * 1e24
# Determine the maximum material ID
all_mats = su.openmc_geometry.get_all_materials()
max_material_id = 0
for material in all_mats:
for material in openmc_geometry.get_all_materials().values():
max_material_id = max(max_material_id, material.id)
# FIXME: Automatically extract info needed to calculate burnable cell volumes
@ -94,48 +51,45 @@ for material in all_mats:
radius = 0.39218
height = 5.
# Use defaultdict since OpenDeplete assumes volumes specified for all cells
volumes = defaultdict(lambda: 1)
# Assign distribmats for each material
for cell in fuel_cells:
new_materials = []
num_instances = len(cell.distribcell_paths)
for i in range(num_instances):
for i in range(cell.num_instances):
new_material = copy.deepcopy(cell.fill)
new_material.id = max_material_id + 1
max_material_id += 1
new_materials.append(new_material)
# Store volume of burnable fuel rods cells
volumes[new_material.id] = np.pi * radius**2 * height
new_material.volume = np.pi * radius**2 * height
new_material.depletable = True
new_material.temperature = 300
cell.fill = new_materials
# Create dt vector for 1 month with 15 day timesteps
dt1 = 15*24*60*60 # 15 days
# Create dt vector for 1 month with 5 day timesteps
dt1 = 5*24*60*60 # 5 days
dt2 = 1.*30*24*60*60 # 1 months
N = np.floor(dt2/dt1)
dt = np.repeat([dt1], N)
# Create settings variable
settings = opendeplete.Settings()
settings = opendeplete.OpenMCSettings()
settings.openmc_call = ["mpirun", "openmc"]
settings.particles = 120000
settings.batches = 20
settings.inactive = 10
settings.lower_left = lower_left
settings.upper_right = upper_right
settings.lower_left = [-21.41728, -21.41728, +192.5]
settings.upper_right = [+21.41728, +21.41728, +197.5]
settings.entropy_dimension = [17*2, 17*2, 1]
settings.power = 2.337e15 * ((17.*17.*2.) / 1.5**2) # MeV/second cm from CASMO
# MeV/second cm from CASMO
settings.power = 2.337e15 * ((17.*17.*2.) / 1.5**2) * height
settings.dt_vec = dt
settings.output_dir = 'depleted'
op = opendeplete.Operator()
op.geometry_fill(su.openmc_geometry, volumes, materials, settings)
op = opendeplete.OpenMCOperator(openmc_geometry, settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
opendeplete.integrate(op, opendeplete.ce_cm_c1)

View file

@ -64,9 +64,8 @@ plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
for cell in openmc_geometry.get_all_material_cells().values():
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)

View file

@ -1,287 +0,0 @@
"""Creates a 2D 2x2 assembly colorset with periodic BCs."""
import numpy as np
import opencg
import openmc
import openmc.opencg_compatible as opencg_compatible
from beavrs.builder import BEAVRS
def find_assembly(assembly_name, wrap_geometry=True):
"""Find a fuel assembly with some string name in the BEAVRS OpenCG model.
This method extracts the fuel assembly and wraps it in an OpenCG Geometry.
The returned geometry has reflective boundary conditions along all
boundaries. The z-axis is bounded between z=200 and z=210 cm.
Parameters
----------
assembly_name : str
The name of the fuel assembly lattice
wrap_geometry : bool
If false, the fuel assembly Lattice is returned. If true, the fuel
assembly Lattice is wrapped in an OpenCG Geometry and returned (default).
Returns
-------
fuel_assembly
The OpenCG Lattice or Geometry for the assembly or None if not found
"""
# Get all OpenCG Universes
all_univ = beavrs.main_universe.get_all_universes()
# Iterate over all Universes
fuel_assembly = None
for univ_id, univ in all_univ.items():
if univ.name == assembly_name:
fuel_assembly = univ
# Wrap lattice in a Geometry if requested by the user
if wrap_geometry:
# Create a root Cell
root_cell = opencg.Cell(name='root cell')
root_cell.fill = fuel_assembly
# Make mixed reflective / vacuum boundaries
min_x = opencg.XPlane(x0=root_cell.fill.min_x, boundary='reflective')
max_x = opencg.XPlane(x0=root_cell.fill.max_x, boundary='reflective')
min_y = opencg.YPlane(y0=root_cell.fill.min_y, boundary='reflective')
max_y = opencg.YPlane(y0=root_cell.fill.max_y, boundary='reflective')
max_z = opencg.ZPlane(z0=197.5, boundary='reflective')
min_z = opencg.ZPlane(z0=192.5, boundary='reflective')
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
root_cell.add_surface(surface=min_z, halfspace=+1)
root_cell.add_surface(surface=max_z, halfspace=-1)
# Create a root Universe
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create a Geometry
fuel_assembly = opencg.Geometry()
fuel_assembly.root_universe = root_univ
return fuel_assembly
def build_two_by_two(assembly1_name, assembly2_name):
"""Build a 2x2 fuel assembly geometry.
This routine puts reflective boundary conditions along all boundaries.
Parameters
----------
assembly1_name : str
The BEAVRS fuel assembly to place in the bottom right and top left
assembly2_name : str
The BEAVRS fuel assembly to place in the bottom left and top right
Returns
-------
opencg.Geometry
A 2x2 fuel assembly OpenCG Geometry
"""
fuel_assembly1 = find_assembly(assembly1_name, wrap_geometry=False)
fuel_assembly2 = find_assembly(assembly2_name, wrap_geometry=False)
# Find the water material
all_cells = beavrs.main_universe.get_all_cells()
for cell_uuid, cell in all_cells.items():
if cell.type == 'material' and cell.fill.name == 'water':
water = cell.fill
# Create a Cell/Universe around the first fuel assembly
fuel_cell1 = opencg.Cell(name='assm1 cell')
fuel_cell1.fill = fuel_assembly1
fuel_univ1 = opencg.Universe(name='assm1 universe')
fuel_univ1.add_cell(fuel_cell1)
# Create a Cell/Universe around the second fuel assembly
fuel_cell2 = opencg.Cell(name='assm2 cell')
fuel_cell2.fill = fuel_assembly2
fuel_univ2 = opencg.Universe(name='assm2 universe')
fuel_univ2.add_cell(fuel_cell2)
# Create a 3x3 lattice two fuel assemblies surrounded by a water reflector
two_by_two_lattice = opencg.Lattice(name='reflector')
lat_width = fuel_assembly1.max_x - fuel_assembly1.min_x
two_by_two_lattice.width = [lat_width, lat_width, 1000.]
two_by_two_lattice.offset = [0., 0., 0.]
two_by_two_lattice.dimension = [2, 2, 1]
two_by_two_lattice.universes = [[fuel_univ1, fuel_univ2],
[fuel_univ2, fuel_univ1]]
# Create a Geometry around the reflected lattice
root_cell = opencg.Cell(name='root cell')
root_cell.fill = two_by_two_lattice
# Make mixed reflective / vacuum boundaries
min_x = opencg.XPlane(x0=root_cell.fill.min_x, boundary='periodic')
max_x = opencg.XPlane(x0=root_cell.fill.max_x, boundary='periodic')
min_y = opencg.YPlane(y0=root_cell.fill.min_y, boundary='periodic')
max_y = opencg.YPlane(y0=root_cell.fill.max_y, boundary='periodic')
min_z = opencg.ZPlane(z0=192.5, boundary='reflective')
max_z = opencg.ZPlane(z0=197.5, boundary='reflective')
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
root_cell.add_surface(surface=min_z, halfspace=+1)
root_cell.add_surface(surface=max_z, halfspace=-1)
# Create a root Universe for this fuel assembly
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create an OpenCG Geometry for this fuel assembly
two_by_two = opencg.Geometry()
two_by_two.root_universe = root_univ
return two_by_two
#### Create OpenMC "materials.xml" and "geometry.xml" files
# Instantiate a BEAVRS object
beavrs = BEAVRS(nndc_xs=True)
# Write all BEAVRS materials to materials.xml file
beavrs.write_openmc_materials()
# Extract fuel assemblies of interest from BEAVRS model
two_by_two = build_two_by_two('Fuel 1.6% enr instr no BAs',
'Fuel 3.1% enr instr 20')
openmc_geometry = opencg_compatible.get_openmc_geometry(two_by_two)
openmc_geometry.export_to_xml()
#### Create OpenMC "settings.xml" file
# Query the user on whether to use multipole cross sections
multipole = input('Use multipole cross sections? (y/n): ').lower()
multipole = True if multipole == 'y' else False
# Construct uniform initial source distribution over fissionable zones
lower_left = two_by_two.bounds[:3]
upper_right = two_by_two.bounds[3:]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings_file = openmc.Settings()
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 10000
settings_file.ptables = True
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
if multipole:
settings_file.temperature = {'multipole': True, 'tolerance': 1000}
settings_file.export_to_xml()
#### Create OpenMC "plots.xml" file
# Initialize the BEAVRS color mapping scheme
beavrs.write_openmc_plots()
# Create a plot colored by materials
plot = openmc.Plot()
bounds = two_by_two.bounds
plot.width = [two_by_two.max_x - two_by_two.min_x,
two_by_two.max_y - two_by_two.min_y]
plot.origin = [bounds[0] + (bounds[3] - bounds[0]) / 2.,
bounds[1] + (bounds[4] - bounds[1]) / 2.,
bounds[2] + (bounds[5] - bounds[2]) / 2.]
plot.color = 'mat'
plot.filename = '2x2-periodic'
plot.col_spec = beavrs.plots.colspec_mat
plot.pixels = [1000, 1000]
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)
# CASMO 70-group structure
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.array([
0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.042, 0.05, 0.058, 0.067,
0.08, 0.1, 0.14, 0.18, 0.22, 0.25, 0.28, 0.3, 0.32, 0.35, 0.4, 0.5, 0.625,
0.78, 0.85, 0.91, 0.95, 0.972, 0.996, 1.02, 1.045, 1.071, 1.097, 1.123,
1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4., 9.877, 15.968, 27.7, 48.052,
75.501, 148.73, 367.26001, 906.90002, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 15.03e3, 24.78e3, 40.85e3, 67.34e3, 111.e3, 183e3, 302.5e3, 500e3,
821e3, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2e7])
# Initialize a 70-group "distribcell" MGXS library
cell_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
cell_mgxs_lib.energy_groups = energy_groups
cell_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
cell_mgxs_lib.domain_type = 'distribcell'
cell_mgxs_lib.domains = fuel_cells
cell_mgxs_lib.correction = None
cell_mgxs_lib.build_library()
# Initialize a 70-group "material" MGXS library
mat_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
mat_mgxs_lib.energy_groups = energy_groups
mat_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
mat_mgxs_lib.domain_type = 'material'
mat_mgxs_lib.correction = None
mat_mgxs_lib.build_library()
#### Create mesh tallies for verification of pin-wise reaction rates
# Instantiate a tally Mesh
mesh = openmc.Mesh(name='assembly mesh')
mesh.type = 'regular'
mesh.dimension = [34, 34, 1]
mesh.lower_left = lower_left
mesh.width = (np.array(upper_right) - np.array(lower_left))
mesh.width[:2] /= 34
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate energy-integrated fission rate mesh Tally
fission_rates = openmc.Tally(name='fission rates')
fission_rates.filters = [mesh_filter]
fission_rates.scores = ['fission']
# Instantiate energy-wise U-238 capture rate mesh Tally
capture_rates = openmc.Tally(name='u-238 capture')
capture_rates.filters = [mesh_filter]
capture_rates.nuclides = ['U238']
capture_rates.scores = ['absorption', 'fission']
#### Create OpenMC "tallies.xml" file
# Create a "tallies.xml" file for the mesh tallies
tallies_file = openmc.Tallies([fission_rates, capture_rates])
cell_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
mat_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()

View file

@ -1,4 +1,4 @@
from collections import OrderedDict, defaultdict
from collections import OrderedDict
import copy
import os
@ -9,56 +9,7 @@ import opendeplete
from geometry import beavrs, openmc_geometry
#### Create "dummy" inputs to export distribcell paths for burnable cells
# Create OpenMC "materials.xml" file
beavrs.write_openmc_materials()
# Create OpenMC "geometry.xml" file
openmc_geometry.export_to_xml()
# Construct uniform initial source distribution over fissionable zones
lower_left = [-32.12592, -32.12592, 192.5]
upper_right = [32.12592, 32.12592, 197.5]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
# Create OpenMC "settings.xml" file
settings_file = openmc.Settings()
settings_file.batches = 2
settings_file.inactive = 1
settings_file.particles = 10
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
settings_file.export_to_xml()
# Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Instantiate a "dummy" distribcell tally for each cell we wish to deplete
for cell in fuel_cells:
tally = openmc.Tally(name='dummy distribcell tally')
distribcell_filter = openmc.DistribcellFilter([cell.id])
tally.filters = [distribcell_filter]
tally.scores = ['fission']
tallies.append(tally)
tallies.export_to_xml()
# Run OpenMC to generate summary.h5 file
openmc.run()
# Open "summary.h5" file
su = openmc.Summary('summary.h5')
fuel_cells = su.openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
#### Setup OpenDeplete Materials wrapper
# Setup OpenDeplete Materials wrapper
materials = opendeplete.Materials()
materials.temperature = OrderedDict()
materials.sab = OrderedDict()
@ -66,10 +17,17 @@ materials.initial_density = OrderedDict()
materials.burn = OrderedDict()
materials.cross_sections = os.environ["OPENMC_CROSS_SECTIONS"]
# Count the number of instances for each cell and material
openmc_geometry.determine_paths()
# Extract all cells filled by a fuel material
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Extract cell materials, temperatures and sab
for cell in su.openmc_geometry.get_all_material_cells():
for cell in openmc_geometry.get_all_material_cells().values():
materials.burn[cell.name] = 'fuel' in cell.fill.name.lower()
materials.temperature[cell.name] = cell.temperature[0]
materials.temperature[cell.name] = 300
if len(cell.fill._sab) > 0:
materials.sab[cell.name] = cell.fill._sab[0]
@ -84,9 +42,8 @@ for cell in fuel_cells:
densities[nuclide][1] * 1e24
# Determine the maximum material ID
all_mats = su.openmc_geometry.get_all_materials()
max_material_id = 0
for material in all_mats:
for material in openmc_geometry.get_all_materials().values():
max_material_id = max(max_material_id, material.id)
# FIXME: Automatically extract info needed to calculate burnable cell volumes
@ -94,48 +51,45 @@ for material in all_mats:
radius = 0.39218
height = 5.
# Use defaultdict since OpenDeplete assumes volumes specified for all cells
volumes = defaultdict(lambda: 1)
# Assign distribmats for each material
for cell in fuel_cells:
new_materials = []
num_instances = len(cell.distribcell_paths)
for i in range(num_instances):
for i in range(cell.num_instances):
new_material = copy.deepcopy(cell.fill)
new_material.id = max_material_id + 1
max_material_id += 1
new_materials.append(new_material)
# Store volume of burnable fuel rods cells
volumes[new_material.id] = np.pi * radius**2 * height
new_material.volume = np.pi * radius**2 * height
new_material.depletable = True
new_material.temperature = 300
cell.fill = new_materials
# Create dt vector for 1 month with 15 day timesteps
dt1 = 15*24*60*60 # 15 days
# Create dt vector for 1 month with 5 day timesteps
dt1 = 5*24*60*60 # 5 days
dt2 = 1.*30*24*60*60 # 1 months
N = np.floor(dt2/dt1)
dt = np.repeat([dt1], N)
# Create settings variable
settings = opendeplete.Settings()
settings = opendeplete.OpenMCSettings()
settings.openmc_call = ["mpirun", "openmc"]
settings.particles = 120000
settings.batches = 30
settings.inactive = 20
settings.lower_left = lower_left
settings.upper_right = upper_right
settings.lower_left = [-32.12592, -32.12592, 192.5]
settings.upper_right = [32.12592, 32.12592, 197.5]
settings.entropy_dimension = [17*3, 17*3, 1]
settings.power = 2.337e15 * ((17.*17.*2.) / 1.5**2) # MeV/second cm from CASMO
# MeV/second cm from CASMO
settings.power = 2.337e15 * ((17.*17.*2.) / 1.5**2) * height
settings.dt_vec = dt
settings.output_dir = 'depleted'
op = opendeplete.Operator()
op.geometry_fill(su.openmc_geometry, volumes, materials, settings)
op = opendeplete.OpenMCOperator(openmc_geometry, settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
opendeplete.integrate(op, opendeplete.ce_cm_c1)

View file

@ -67,9 +67,8 @@ plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
for cell in openmc_geometry.get_all_material_cells().values():
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)

View file

@ -1,299 +0,0 @@
"""Creates a 2D 2x2 assembly colorset with a water reflector."""
import numpy as np
import opencg
import openmc
import openmc.opencg_compatible as opencg_compatible
from beavrs.builder import BEAVRS
def find_assembly(assembly_name, wrap_geometry=True):
"""Find a fuel assembly with some string name in the BEAVRS OpenCG model.
This method extracts the fuel assembly and wraps it in an OpenCG Geometry.
The returned geometry has reflective boundary conditions along all
boundaries. The z-axis is bounded between z=200 and z=210 cm.
Parameters
----------
assembly_name : str
The name of the fuel assembly lattice
wrap_geometry : bool
If false, the fuel assembly Lattice is returned. If true, the fuel
assembly Lattice is wrapped in an OpenCG Geometry and returned (default).
Returns
-------
fuel_assembly
The OpenCG Lattice or Geometry for the assembly or None if not found
"""
# Get all OpenCG Universes
all_univ = beavrs.main_universe.get_all_universes()
# Iterate over all Universes
fuel_assembly = None
for univ_id, univ in all_univ.items():
if univ.name == assembly_name:
fuel_assembly = univ
# Wrap lattice in a Geometry if requested by the user
if wrap_geometry:
# Create a root Cell
root_cell = opencg.Cell(name='root cell')
root_cell.fill = fuel_assembly
# Make mixed reflective / vacuum boundaries
min_x = opencg.XPlane(x0=root_cell.fill.min_x, boundary='reflective')
max_x = opencg.XPlane(x0=root_cell.fill.max_x, boundary='reflective')
min_y = opencg.YPlane(y0=root_cell.fill.min_y, boundary='reflective')
max_y = opencg.YPlane(y0=root_cell.fill.max_y, boundary='reflective')
max_z = opencg.ZPlane(z0=197.5, boundary='reflective')
min_z = opencg.ZPlane(z0=192.5, boundary='reflective')
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
root_cell.add_surface(surface=min_z, halfspace=+1)
root_cell.add_surface(surface=max_z, halfspace=-1)
# Create a root Universe
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create a Geometry
fuel_assembly = opencg.Geometry()
fuel_assembly.root_universe = root_univ
return fuel_assembly
def build_reflector(assembly1_name, assembly2_name):
"""Build a 2x2 fuel assembly geometry with a water reflector on the
bottom and right.
This routine puts reflective boundary conditions along min x, max y and z
and vacuum boundary conditions along the max x and min y boundaries.
Parameters
----------
assembly1_name : str
The BEAVRS fuel assembly to place in the bottom right and top left
assembly2_name : str
The BEAVRS fuel assembly to place in the bottom left and top right
Returns
-------
opencg.Geometry
A 2x2 fuel assembly and reflector OpenCG Geometry
"""
fuel_assembly1 = find_assembly(assembly1_name, wrap_geometry=False)
fuel_assembly2 = find_assembly(assembly2_name, wrap_geometry=False)
# Find the water material
all_cells = beavrs.main_universe.get_all_cells()
for cell_uuid, cell in all_cells.items():
if cell.type == 'material' and cell.fill.name == 'water':
water = cell.fill
# Create a Cell/Universe around the first fuel assembly
fuel_cell1 = opencg.Cell(name='assm1 cell')
fuel_cell1.fill = fuel_assembly1
fuel_univ1 = opencg.Universe(name='assm1 universe')
fuel_univ1.add_cell(fuel_cell1)
# Create a Cell/Universe around the second fuel assembly
fuel_cell2 = opencg.Cell(name='assm2 cell')
fuel_cell2.fill = fuel_assembly2
fuel_univ2 = opencg.Universe(name='assm2 universe')
fuel_univ2.add_cell(fuel_cell2)
# Create a Cell/Universe with water
water_cell = opencg.Cell(name='water cell', fill=water)
water_univ = opencg.Universe(name='water universe')
water_univ.add_cell(water_cell)
# Create a 3x3 lattice two fuel assemblies surrounded by a water reflector
reflector_lattice = opencg.Lattice(name='reflector')
lat_width = fuel_assembly1.max_x - fuel_assembly1.min_x
reflector_lattice.width = [lat_width, lat_width, 1000.]
reflector_lattice.dimension = [3, 3, 1]
reflector_lattice.offset = [0., 0., 0.]
reflector_lattice.universes = [[fuel_univ1, fuel_univ2, water_univ],
[fuel_univ2, fuel_univ1, water_univ],
[water_univ, water_univ, water_univ]]
# Create a Geometry around the reflected lattice
root_cell = opencg.Cell(name='root cell')
root_cell.fill = reflector_lattice
# Make mixed reflective / vacuum boundaries
min_x = opencg.XPlane(x0=root_cell.fill.min_x, boundary='reflective')
max_x = opencg.XPlane(x0=root_cell.fill.max_x, boundary='vacuum')
min_y = opencg.YPlane(y0=root_cell.fill.min_y, boundary='vacuum')
max_y = opencg.YPlane(y0=root_cell.fill.max_y, boundary='reflective')
min_z = opencg.ZPlane(z0=192.5, boundary='reflective')
max_z = opencg.ZPlane(z0=197.5, boundary='reflective')
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
root_cell.add_surface(surface=min_z, halfspace=+1)
root_cell.add_surface(surface=max_z, halfspace=-1)
# Create a root Universe for this fuel assembly
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create an OpenCG Geometry for this fuel assembly
reflector = opencg.Geometry()
reflector.root_universe = root_univ
return reflector
#### Create OpenMC "materials.xml" and "geometry.xml" files
# Instantiate a BEAVRS object
beavrs = BEAVRS(nndc_xs=True)
# Write all BEAVRS materials to materials.xml file
beavrs.write_openmc_materials()
# Extract fuel assemblies of interest from BEAVRS model
reflector = build_reflector('Fuel 1.6% enr instr no BAs',
'Fuel 3.1% enr instr 20')
openmc_geometry = opencg_compatible.get_openmc_geometry(reflector)
openmc_geometry.export_to_xml()
#### Create OpenMC "settings.xml" file
# Query the user on whether to use multipole cross sections
multipole = input('Use multipole cross sections? (y/n): ').lower()
multipole = True if multipole == 'y' else False
# Construct uniform initial source distribution over fissionable zones
lower_left = reflector.bounds[:3]
upper_right = reflector.bounds[3:]
lat_width = (np.array(upper_right) - np.array(lower_left))
lat_width[:2] /= 3.
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings_file = openmc.Settings()
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 10000
settings_file.ptables = True
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
if multipole:
settings_file.temperature = {'multipole': True, 'tolerance': 1000}
settings_file.export_to_xml()
#### Create OpenMC "plots.xml" file
# Initialize the BEAVRS color mapping scheme
beavrs.write_openmc_plots()
# Create a plot colored by materials
plot = openmc.Plot()
bounds = reflector.bounds
plot.width = [reflector.max_x - reflector.min_x,
reflector.max_y - reflector.min_y]
plot.origin = [bounds[0] + (bounds[3] - bounds[0]) / 2.,
bounds[1] + (bounds[4] - bounds[1]) / 2.,
bounds[2] + (bounds[5] - bounds[2]) / 2.]
plot.color = 'mat'
plot.filename = '2x2-reflector'
plot.col_spec = beavrs.plots.colspec_mat
plot.pixels = [1000, 1000]
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)
# CASMO 70-group structure
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.array([
0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.042, 0.05, 0.058, 0.067,
0.08, 0.1, 0.14, 0.18, 0.22, 0.25, 0.28, 0.3, 0.32, 0.35, 0.4, 0.5, 0.625,
0.78, 0.85, 0.91, 0.95, 0.972, 0.996, 1.02, 1.045, 1.071, 1.097, 1.123,
1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4., 9.877, 15.968, 27.7, 48.052,
75.501, 148.73, 367.26001, 906.90002, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 15.03e3, 24.78e3, 40.85e3, 67.34e3, 111.e3, 183e3, 302.5e3, 500e3,
821e3, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2e7])
# Initialize a 70-group "distribcell" MGXS library
cell_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
cell_mgxs_lib.energy_groups = energy_groups
cell_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
cell_mgxs_lib.domain_type = 'distribcell'
cell_mgxs_lib.domains = fuel_cells
cell_mgxs_lib.correction = None
cell_mgxs_lib.build_library()
# Initialize a 70-group "material" MGXS library
mat_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
mat_mgxs_lib.energy_groups = energy_groups
mat_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
mat_mgxs_lib.domain_type = 'material'
mat_mgxs_lib.correction = None
mat_mgxs_lib.build_library()
#### Create mesh tallies for verification of pin-wise reaction rates
# Instantiate a tally Mesh
mesh = openmc.Mesh(name='assembly mesh')
mesh.type = 'regular'
mesh.dimension = [34, 34, 1]
mesh.lower_left = [lower_left[0], lower_left[1] + lat_width[1], lower_left[2]]
mesh.width = np.array(lat_width)
mesh.width[:2] /= 17.
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate energy-integrated fission rate mesh Tally
fission_rates = openmc.Tally(name='fission rates')
fission_rates.filters = [mesh_filter]
fission_rates.scores = ['fission']
# Instantiate energy-wise U-238 capture rate mesh Tally
capture_rates = openmc.Tally(name='u-238 capture')
capture_rates.filters = [mesh_filter]
capture_rates.nuclides = ['U238']
capture_rates.scores = ['absorption', 'fission']
#### Create OpenMC "tallies.xml" file
# Create a "tallies.xml" file for the mesh tallies
tallies_file = openmc.Tallies([fission_rates, capture_rates])
cell_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
mat_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()

View file

@ -1,4 +1,4 @@
from collections import OrderedDict, defaultdict
from collections import OrderedDict
import copy
import os
@ -9,55 +9,7 @@ import opendeplete
from geometry import beavrs, openmc_geometry
#### Create "dummy" inputs to export distribcell paths for burnable cells
# Create OpenMC "materials.xml" file
beavrs.write_openmc_materials()
# Create OpenMC "geometry.xml" file
openmc_geometry.export_to_xml()
# Construct uniform initial source distribution over fissionable zones
lower_left = [-10.70864, -10.70864, +192.5]
upper_right = [+10.70864, +10.70864, +197.5]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
# Create OpenMC "settings.xml" file
settings_file = openmc.Settings()
settings_file.batches = 2
settings_file.inactive = 1
settings_file.particles = 10
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
settings_file.export_to_xml()
# Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Instantiate a "dummy" distribcell tally for each cell we wish to deplete
for cell in fuel_cells:
tally = openmc.Tally(name='dummy distribcell tally')
distribcell_filter = openmc.DistribcellFilter([cell.id])
tally.filters = [distribcell_filter]
tally.scores = ['fission']
tallies.append(tally)
tallies.export_to_xml()
# Run OpenMC to generate summary.h5 file
openmc.run()
# Open "summary.h5" file
su = openmc.Summary('summary.h5')
fuel_cells = su.openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
#### Setup OpenDeplete Materials wrapper
# Setup OpenDeplete Materials wrapper
materials = opendeplete.Materials()
materials.temperature = OrderedDict()
materials.sab = OrderedDict()
@ -65,10 +17,17 @@ materials.initial_density = OrderedDict()
materials.burn = OrderedDict()
materials.cross_sections = os.environ["OPENMC_CROSS_SECTIONS"]
# Count the number of instances for each cell and material
openmc_geometry.determine_paths()
# Extract all cells filled by a fuel material
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Extract cell materials, temperatures and sab
for cell in su.openmc_geometry.get_all_material_cells():
for cell in openmc_geometry.get_all_material_cells().values():
materials.burn[cell.name] = 'fuel' in cell.fill.name.lower()
materials.temperature[cell.name] = cell.temperature[0]
materials.temperature[cell.name] = 300
if len(cell.fill._sab) > 0:
materials.sab[cell.name] = cell.fill._sab[0]
@ -83,9 +42,8 @@ for cell in fuel_cells:
densities[nuclide][1] * 1e24
# Determine the maximum material ID
all_mats = su.openmc_geometry.get_all_materials()
max_material_id = 0
for material in all_mats:
for material in openmc_geometry.get_all_materials().values():
max_material_id = max(max_material_id, material.id)
# FIXME: Automatically extract info needed to calculate burnable cell volumes
@ -93,48 +51,45 @@ for material in all_mats:
radius = 0.39218
height = 5.
# Use defaultdict since OpenDeplete assumes volumes specified for all cells
volumes = defaultdict(lambda: 1)
# Assign distribmats for each material
for cell in fuel_cells:
new_materials = []
num_instances = len(cell.distribcell_paths)
for i in range(num_instances):
for i in range(cell.num_instances):
new_material = copy.deepcopy(cell.fill)
new_material.id = max_material_id + 1
max_material_id += 1
new_materials.append(new_material)
# Store volume of burnable fuel rods cells
volumes[new_material.id] = np.pi * radius**2 * height
new_material.volume = np.pi * radius**2 * height
new_material.depletable = True
new_material.temperature = 300
cell.fill = new_materials
# Create dt vector for 1 month with 15 day timesteps
dt1 = 15*24*60*60 # 15 days
# Create dt vector for 1 month with 5 day timesteps
dt1 = 5*24*60*60 # 5 days
dt2 = 1.*30*24*60*60 # 1 months
N = np.floor(dt2/dt1)
dt = np.repeat([dt1], N)
# Create settings variable
settings = opendeplete.Settings()
settings = opendeplete.OpenMCSettings()
settings.openmc_call = ["mpirun", "openmc"]
settings.particles = 30000
settings.batches = 20
settings.inactive = 10
settings.lower_left = lower_left
settings.upper_right = upper_right
settings.lower_left = [-10.70864, -10.70864, +192.5]
settings.upper_right = [+10.70864, +10.70864, +197.5]
settings.entropy_dimension = [17, 17, 1]
settings.power = 2.337e15 * (17.**2 / 1.5**2) # MeV/second cm from CASMO
# MeV/second cm from CASMO
settings.power = 2.337e15 * (17.**2 / 1.5**2) * height
settings.dt_vec = dt
settings.output_dir = 'depleted'
op = opendeplete.Operator()
op.geometry_fill(su.openmc_geometry, volumes, materials, settings)
op = opendeplete.OpenMCOperator(openmc_geometry, settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
opendeplete.integrate(op, opendeplete.ce_cm_c1)

View file

@ -64,9 +64,8 @@ plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
for cell in openmc_geometry.get_all_material_cells().values():
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)

View file

@ -1,207 +0,0 @@
"""Creates a 2D fuel assembly with reflective BCs."""
import numpy as np
import opencg
import openmc
import openmc.opencg_compatible as opencg_compatible
from beavrs.builder import BEAVRS
def find_assembly(assembly_name, wrap_geometry=True):
"""Find a fuel assembly with some string name in the BEAVRS OpenCG model.
This method extracts the fuel assembly and wraps it in an OpenCG Geometry.
The returned geometry has reflective boundary conditions along all
boundaries. The z-axis is bounded between z=200 and z=210 cm.
Parameters
----------
assembly_name : str
The name of the fuel assembly lattice
wrap_geometry : bool
If false, the fuel assembly Lattice is returned. If true, the fuel
assembly Lattice is wrapped in an OpenCG Geometry and returned (default).
Returns
-------
fuel_assembly
The OpenCG Lattice or Geometry for the assembly or None if not found
"""
# Get all OpenCG Universes
all_univ = beavrs.main_universe.get_all_universes()
# Iterate over all Universes
fuel_assembly = None
for univ_id, univ in all_univ.items():
if univ.name == assembly_name:
fuel_assembly = univ
# Wrap lattice in a Geometry if requested by the user
if wrap_geometry:
# Create a root Cell
root_cell = opencg.Cell(name='root cell')
root_cell.fill = fuel_assembly
# Make mixed reflective / vacuum boundaries
min_x = opencg.XPlane(x0=root_cell.fill.min_x, boundary='reflective')
max_x = opencg.XPlane(x0=root_cell.fill.max_x, boundary='reflective')
min_y = opencg.YPlane(y0=root_cell.fill.min_y, boundary='reflective')
max_y = opencg.YPlane(y0=root_cell.fill.max_y, boundary='reflective')
max_z = opencg.ZPlane(z0=197.5, boundary='reflective')
min_z = opencg.ZPlane(z0=192.5, boundary='reflective')
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
root_cell.add_surface(surface=min_z, halfspace=+1)
root_cell.add_surface(surface=max_z, halfspace=-1)
# Create a root Universe
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create a Geometry
fuel_assembly = opencg.Geometry()
fuel_assembly.root_universe = root_univ
return fuel_assembly
#### Create OpenMC "materials.xml" and "geometry.xml" files
# Instantiate a BEAVRS object
beavrs = BEAVRS(nndc_xs=True)
# Write all BEAVRS materials to materials.xml file
beavrs.write_openmc_materials()
# Extract fuel assembly of interest from BEAVRS model
assm_name = 'Fuel 1.6% enr instr no BAs'
fuel_assembly = find_assembly(assm_name)
openmc_geometry = opencg_compatible.get_openmc_geometry(fuel_assembly)
openmc_geometry.export_to_xml()
#### Create OpenMC "settings.xml" file
# Query the user on whether to use multipole cross sections
multipole = input('Use multipole cross sections? (y/n): ').lower()
multipole = True if multipole == 'y' else False
# Construct uniform initial source distribution over fissionable zones
lower_left = fuel_assembly.bounds[:3]
upper_right = fuel_assembly.bounds[3:]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings_file = openmc.Settings()
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 10000
settings_file.ptables = True
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
if multipole:
settings_file.temperature = {'multipole': True, 'tolerance': 1000}
settings_file.export_to_xml()
#### Create OpenMC "plots.xml" file
# Initialize the BEAVRS color mapping scheme
beavrs.write_openmc_plots()
# Create a plot colored by materials
plot = openmc.Plot()
bounds = fuel_assembly.bounds
plot.width = [fuel_assembly.max_x - fuel_assembly.min_x,
fuel_assembly.max_y - fuel_assembly.min_y]
plot.origin = [bounds[0] + (bounds[3] - bounds[0]) / 2.,
bounds[1] + (bounds[4] - bounds[1]) / 2.,
bounds[2] + (bounds[5] - bounds[2]) / 2.]
plot.color = 'mat'
plot.filename = 'assembly'
plot.col_spec = beavrs.plots.colspec_mat
plot.pixels = [1000, 1000]
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = openmc_geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)
# CASMO 70-group structure
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.array([
0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.042, 0.05, 0.058, 0.067,
0.08, 0.1, 0.14, 0.18, 0.22, 0.25, 0.28, 0.3, 0.32, 0.35, 0.4, 0.5, 0.625,
0.78, 0.85, 0.91, 0.95, 0.972, 0.996, 1.02, 1.045, 1.071, 1.097, 1.123,
1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4., 9.877, 15.968, 27.7, 48.052,
75.501, 148.73, 367.26001, 906.90002, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 15.03e3, 24.78e3, 40.85e3, 67.34e3, 111.e3, 183e3, 302.5e3, 500e3,
821e3, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2e7])
# Initialize a 70-group "distribcell" MGXS library
cell_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
cell_mgxs_lib.energy_groups = energy_groups
cell_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
cell_mgxs_lib.domain_type = 'distribcell'
cell_mgxs_lib.domains = fuel_cells
cell_mgxs_lib.correction = None
cell_mgxs_lib.build_library()
# Initialize a 70-group "material" MGXS library
mat_mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
mat_mgxs_lib.energy_groups = energy_groups
mat_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
mat_mgxs_lib.domain_type = 'material'
mat_mgxs_lib.correction = None
mat_mgxs_lib.build_library()
#### Create mesh tallies for verification of pin-wise reaction rates
# Instantiate a tally Mesh
mesh = openmc.Mesh(name='assembly mesh')
mesh.type = 'regular'
mesh.dimension = [17, 17, 1]
mesh.lower_left = lower_left
mesh.width = (np.array(upper_right) - np.array(lower_left))
mesh.width[:2] /= 17
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate energy-integrated fission rate mesh Tally
fission_rates = openmc.Tally(name='fission rates')
fission_rates.filters = [mesh_filter]
fission_rates.scores = ['fission']
# Instantiate energy-wise U-238 capture rate mesh Tally
capture_rates = openmc.Tally(name='u-238 capture')
capture_rates.filters = [mesh_filter]
capture_rates.nuclides = ['U238']
capture_rates.scores = ['absorption', 'fission']
#### Create OpenMC "tallies.xml" file
# Create a "tallies.xml" file for the mesh tallies
tallies_file = openmc.Tallies([fission_rates, capture_rates])
cell_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
mat_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()

View file

@ -1,4 +1,4 @@
from collections import OrderedDict, defaultdict
from collections import OrderedDict
import copy
import os
@ -9,55 +9,7 @@ import opendeplete
from geometry import beavrs, openmc_geometry
#### Create "dummy" inputs to export distribcell paths for burnable cells
# Create OpenMC "materials.xml" file
beavrs.write_openmc_materials()
# Create OpenMC "geometry.xml" file
openmc_geometry.export_to_xml()
# Construct uniform initial source distribution over fissionable zones
lower_left = [-0.62992, -0.62992, -10.0]
upper_right = [+0.62992, +0.62992, +10.0]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
# Create OpenMC "settings.xml" file
settings_file = openmc.Settings()
settings_file.batches = 2
settings_file.inactive = 1
settings_file.particles = 10
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
settings_file.export_to_xml()
# Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Instantiate a "dummy" distribcell tally for each cell we wish to deplete
for cell in fuel_cells:
tally = openmc.Tally(name='dummy distribcell tally')
distribcell_filter = openmc.DistribcellFilter([cell.id])
tally.filters = [distribcell_filter]
tally.scores = ['fission']
tallies.append(tally)
tallies.export_to_xml()
# Run OpenMC to generate summary.h5 file
openmc.run()
# Open "summary.h5" file
su = openmc.Summary('summary.h5')
fuel_cells = su.openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
#### Setup OpenDeplete Materials wrapper
# Setup OpenDeplete Materials wrapper
materials = opendeplete.Materials()
materials.temperature = OrderedDict()
materials.sab = OrderedDict()
@ -65,10 +17,17 @@ materials.initial_density = OrderedDict()
materials.burn = OrderedDict()
materials.cross_sections = os.environ["OPENMC_CROSS_SECTIONS"]
# Count the number of instances for each cell and material
openmc_geometry.determine_paths()
# Extract all cells filled by a fuel material
fuel_cells = openmc_geometry.get_cells_by_name(
name='enr radial 0: Fuel', case_sensitive=True)
# Extract cell materials, temperatures and sab
for cell in su.openmc_geometry.get_all_material_cells():
for cell in openmc_geometry.get_all_material_cells().values():
materials.burn[cell.name] = 'fuel' in cell.fill.name.lower()
materials.temperature[cell.name] = cell.temperature[0]
materials.temperature[cell.name] = 300
if len(cell.fill._sab) > 0:
materials.sab[cell.name] = cell.fill._sab[0]
@ -83,9 +42,8 @@ for cell in fuel_cells:
densities[nuclide][1] * 1e24
# Determine the maximum material ID
all_mats = su.openmc_geometry.get_all_materials()
max_material_id = 0
for material in all_mats:
for material in openmc_geometry.get_all_materials().values():
max_material_id = max(max_material_id, material.id)
# FIXME: Automatically extract info needed to calculate burnable cell volumes
@ -93,48 +51,45 @@ for material in all_mats:
radius = 0.39218
height = 5.
# Use defaultdict since OpenDeplete assumes volumes specified for all cells
volumes = defaultdict(lambda: 1)
# Assign distribmats for each material
for cell in fuel_cells:
new_materials = []
num_instances = len(cell.distribcell_paths)
for i in range(num_instances):
for i in range(cell.num_instances):
new_material = copy.deepcopy(cell.fill)
new_material.id = max_material_id + 1
max_material_id += 1
new_materials.append(new_material)
# Store volume of burnable fuel rods cells
volumes[new_material.id] = np.pi * radius**2 * height
new_material.volume = np.pi * radius**2 * height
new_material.depletable = True
new_material.temperature = 300
cell.fill = new_materials
# Create dt vector for 1 month with 15 day timesteps
dt1 = 15*24*60*60 # 15 days
# Create dt vector for 1 month with 5 day timesteps
dt1 = 5*24*60*60 # 5 days
dt2 = 1.*30*24*60*60 # 1 months
N = np.floor(dt2/dt1)
dt = np.repeat([dt1], N)
# Create settings variable
settings = opendeplete.Settings()
settings = opendeplete.OpenMCSettings()
settings.openmc_call = ["mpirun", "openmc"]
settings.particles = 10000
settings.batches = 20
settings.inactive = 10
settings.lower_left = lower_left
settings.upper_right = upper_right
settings.lower_left = [-0.62992, -0.62992, -10.0]
settings.upper_right = [+0.62992, +0.62992, +10.0]
settings.entropy_dimension = [1, 1, 1]
settings.power = 2.337e15 * (1.**2 / 1.5**2) # MeV/second cm from CASMO
# MeV/second cm from CASMO
settings.power = 2.337e15 * (1.**2 / 1.5**2) * height
settings.dt_vec = dt
settings.output_dir = 'depleted'
op = opendeplete.Operator()
op.geometry_fill(su.openmc_geometry, volumes, materials, settings)
op = opendeplete.OpenMCOperator(openmc_geometry, settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
opendeplete.integrate(op, opendeplete.ce_cm_c1)

View file

@ -1,164 +0,0 @@
"""Creates a 2D fuel pin cell with reflective BCs."""
import numpy as np
import opencg
import openmc
import openmc.opencg_compatible as opencg_compatible
from beavrs.builder import BEAVRS
def find_pin(pin_name, wrap_geometry=True):
"""Find a fuel pin with some string name in the BEAVRS OpenCG model.
This method extracts the pin cell and wraps it in an OpenCG Geometry.
The returned geometry has reflective boundary conditions along the x and y
boundaries. The z-axis left unbounded.
Parameters
----------
pin_name : str
The name of the fuel pin universe
wrap_geometry : bool
If false, the pin cell Universe is returned. If true, the pin cell
Universe is wrapped in an OpenCG Geometry and returned (default).
Returns
-------
opencg.Universe
The OpenCG Universe or Geometry for this fuel pin or None if not found
"""
# Get all OpenCG Universes
all_univ = beavrs.main_universe.get_all_universes()
# Iterate over all Universes
fuel_pin = None
for univ_id, univ in all_univ.items():
if univ._name == pin_name:
fuel_pin = univ
# Wrap pin cell Universe in a Geometry if requested by the user
if wrap_geometry:
# Make reflective boundaries
pin_pitch = 0.62992
min_x = opencg.XPlane(x0=-pin_pitch, boundary='reflective')
max_x = opencg.XPlane(x0=pin_pitch, boundary='reflective')
min_y = opencg.YPlane(y0=-pin_pitch, boundary='reflective')
max_y = opencg.YPlane(y0=pin_pitch, boundary='reflective')
# Create a root Cell
root_cell = opencg.Cell(name='root cell')
root_cell.fill = fuel_pin
# Add boundaries to the root Cell
root_cell.add_surface(surface=min_x, halfspace=+1)
root_cell.add_surface(surface=max_x, halfspace=-1)
root_cell.add_surface(surface=min_y, halfspace=+1)
root_cell.add_surface(surface=max_y, halfspace=-1)
# Create a root Universe
root_univ = opencg.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Create a Geometry
fuel_pin = opencg.Geometry()
fuel_pin.root_universe = root_univ
return fuel_pin
#### Create OpenMC "materials.xml" and "geometry.xml" files
# User-specified enrichment of 1.6, 2.4 or 3.1 percent
enrichment = 1.6
# Instantiate a BEAVRS object
beavrs = BEAVRS(nndc_xs=True)
# Write all BEAVRS materials to materials.xml file
beavrs.write_openmc_materials()
# Extract fuel pin of interest from BEAVRS model
pin_name = 'Fuel rod active region - {}% enr'.format(enrichment)
pin_geometry = find_pin(pin_name)
openmc_geometry = opencg_compatible.get_openmc_geometry(pin_geometry)
openmc_geometry.export_to_xml()
#### Create OpenMC "settings.xml" file
# Query the user on whether to use multipole cross sections
multipole = input('Use multipole cross sections? (y/n): ').lower()
multipole = True if multipole == 'y' else False
# Construct uniform initial source distribution over fissionable zones
lower_left = pin_geometry.bounds[:2] + [-10.]
upper_right = pin_geometry.bounds[3:5] + [10.]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings_file = openmc.Settings()
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 10000
settings_file.ptables = True
settings_file.output = {'tallies': False}
settings_file.source = source
settings_file.sourcepoint_write = False
if multipole:
settings_file.temperature = {'multipole': True, 'tolerance': 1000}
settings_file.export_to_xml()
#### Create OpenMC "plots.xml" file
# Initialize the BEAVRS color mapping scheme
beavrs.write_openmc_plots()
# Create a plot colored by materials
plot = openmc.Plot()
bounds = pin_geometry.bounds
plot.width = [pin_geometry.max_x - pin_geometry.min_x,
pin_geometry.max_y - pin_geometry.min_y]
plot.origin = [bounds[0] + (bounds[3] - bounds[0]) / 2.,
bounds[1] + (bounds[4] - bounds[1]) / 2.,
bounds[2] + (bounds[5] - bounds[2]) / 2.]
plot.color = 'mat'
plot.filename = 'fuel-pin'
plot.col_spec = beavrs.plots.colspec_mat
plot.pixels = [1000, 1000]
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()
#### Create OpenMC MGXS library and "tallies.xml" file
# CASMO 70-group structure
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.array([
0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.042, 0.05, 0.058, 0.067,
0.08, 0.1, 0.14, 0.18, 0.22, 0.25, 0.28, 0.3, 0.32, 0.35, 0.4, 0.5, 0.625,
0.78, 0.85, 0.91, 0.95, 0.972, 0.996, 1.02, 1.045, 1.071, 1.097, 1.123,
1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4., 9.877, 15.968, 27.7, 48.052,
75.501, 148.73, 367.26001, 906.90002, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 15.03e3, 24.78e3, 40.85e3, 67.34e3, 111.e3, 183e3, 302.5e3, 500e3,
821e3, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2e7])
# Initialize a 70-group MGXS library
mgxs_lib = openmc.mgxs.Library(openmc_geometry, by_nuclide=True)
mgxs_lib.energy_groups = energy_groups
mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
mgxs_lib.domain_type = 'material'
mgxs_lib.correction = None
mgxs_lib.build_library()
# Create a "tallies.xml" file for the MGXS Library
tallies_file = openmc.Tallies()
mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()

View file

@ -1,4 +1,4 @@
from collections import OrderedDict, defaultdict
from collections import OrderedDict
import copy
import os
@ -10,64 +10,7 @@ from smr.materials import materials
from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core
from smr.core import geometry
#### Create "dummy" inputs to export distribcell paths for burnable cells
# Create OpenMC "geometry.xml" file
geometry.export_to_xml()
# Create OpenMC "materials.xml" file
materials.export_to_xml()
# 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
# Create OpenMC "settings.xml" file
settings = openmc.Settings()
settings.batches = 2
settings.inactive = 1
settings.particles = 10
settings.output = {'tallies': False}
settings.source = source
settings.sourcepoint_write = False
settings.export_to_xml()
# Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
fuel_cells = geometry.get_cells_by_name(
name='(1.6%) (0)', case_sensitive=True)
fuel_cells.extend(geometry.get_cells_by_name(
name='(2.4%) (0)', case_sensitive=True))
fuel_cells.extend(geometry.get_cells_by_name(
name='(3.1%) (0)', case_sensitive=True))
# Instantiate a "dummy" distribcell tally for each cell we wish to deplete
for cell in fuel_cells:
tally = openmc.Tally(name='dummy distribcell tally')
distribcell_filter = openmc.DistribcellFilter([cell.id])
tally.filters = [distribcell_filter]
tally.scores = ['fission']
tallies.append(tally)
tallies.export_to_xml()
# Run OpenMC to generate summary.h5 file
openmc.run()
# Open "summary.h5" file
su = openmc.Summary('summary.h5')
fuel_cells = su.openmc_geometry.get_cells_by_name(
name='(1.6%) (0)', case_sensitive=True)
fuel_cells.extend(su.openmc_geometry.get_cells_by_name(
name='(2.4%) (0)', case_sensitive=True))
fuel_cells.extend(su.openmc_geometry.get_cells_by_name(
name='(3.1%) (0)', case_sensitive=True))
#### Setup OpenDeplete Materials wrapper
# Setup OpenDeplete Materials wrapper
materials = opendeplete.Materials()
materials.temperature = OrderedDict()
materials.sab = OrderedDict()
@ -75,10 +18,21 @@ materials.initial_density = OrderedDict()
materials.burn = OrderedDict()
materials.cross_sections = os.environ["OPENMC_CROSS_SECTIONS"]
# Count the number of instances for each cell and material
geometry.determine_paths()
# Extract all cells filled by a fuel material
fuel_cells = geometry.get_cells_by_name(
name='(1.6%) (0)', case_sensitive=True)
fuel_cells.extend(geometry.get_cells_by_name(
name='(2.4%) (0)', case_sensitive=True))
fuel_cells.extend(geometry.get_cells_by_name(
name='(3.1%) (0)', case_sensitive=True))
# Extract cell materials, temperatures and sab
for cell in su.openmc_geometry.get_all_material_cells():
for cell in geometry.get_all_material_cells().values():
materials.burn[cell.name] = 'fuel' in cell.fill.name.lower()
materials.temperature[cell.name] = cell.temperature[0]
materials.temperature[cell.name] = 300
if len(cell.fill._sab) > 0:
materials.sab[cell.name] = cell.fill._sab[0]
@ -93,58 +47,56 @@ for cell in fuel_cells:
densities[nuclide][1] * 1e24
# Determine the maximum material ID
all_mats = su.openmc_geometry.get_all_materials()
max_material_id = 0
for material in all_mats:
for material in geometry.get_all_materials().values():
max_material_id = max(max_material_id, material.id)
# FIXME: Automatically extract info needed to calculate burnable cell volumes
# Fuel rod geometric parameters
radius = 0.39218
height = 5.
# Use defaultdict since OpenDeplete assumes volumes specified for all cells
volumes = defaultdict(lambda: 1)
height = 200.
# Assign distribmats for each material
for cell in fuel_cells:
new_materials = []
num_instances = len(cell.distribcell_paths)
for i in range(num_instances):
for i in range(cell.num_instances):
new_material = copy.deepcopy(cell.fill)
new_material.id = max_material_id + 1
max_material_id += 1
new_materials.append(new_material)
# Store volume of burnable fuel rods cells
volumes[new_material.id] = np.pi * radius**2 * height
new_material.volume = np.pi * radius**2 * height
new_material.depletable = True
new_material.temperature = 300
cell.fill = new_materials
# Create dt vector for 1 month with 15 day timesteps
dt1 = 15*24*60*60 # 15 days
# Create dt vector for 1 month with 5 day timesteps
dt1 = 5*24*60*60 # 5 days
dt2 = 1.*30*24*60*60 # 1 months
N = np.floor(dt2/dt1)
dt = np.repeat([dt1], N)
# Create settings variable
settings = opendeplete.Settings()
settings = opendeplete.OpenMCSettings()
settings.openmc_call = ["mpirun", "openmc"]
settings.particles = 1000000
settings.batches = 200
settings.inactive = 100
settings.lower_left = lower_left
settings.upper_right = upper_right
settings.lower_left = \
[-7.*lattice_pitch/2., -7.*lattice_pitch/2., bottom_fuel_stack]
settings.upper_right = \
[+7.*lattice_pitch/2., +7.*lattice_pitch/2., top_active_core]
settings.entropy_dimension = [15, 15, 1]
settings.power = 2.337e15 * ((17.*17.*37.) / 1.5**2) # MeV/second cm from CASMO
# MeV/second cm from CASMO
settings.power = 2.337e15 * ((17.*17.*37.) / 1.5**2) * height
settings.dt_vec = dt
settings.output_dir = 'depleted'
op = opendeplete.Operator()
op.geometry_fill(su.openmc_geometry, volumes, materials, settings)
op = opendeplete.OpenMCOperator(geometry, settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
opendeplete.integrate(op, opendeplete.ce_cm_c1)

View file

@ -51,9 +51,8 @@ plots.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
for cell in geometry.get_all_material_cells().values():
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)

View file

@ -1,109 +0,0 @@
import numpy as np
import openmc
from smr.materials import materials
from smr.plots import plots
from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core
from smr.core import geometry
#### Create OpenMC "geometry.xml" file
geometry.export_to_xml()
#### Create OpenMC "materials.xml" file
materials.export_to_xml()
#### Create OpenMC "settings.xml" file
# Query the user on whether to use multipole cross sections
multipole = input('Use multipole cross sections? (y/n): ').lower()
multipole = (multipole == 'y')
# 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}
settings.source = source
settings.sourcepoint_write = False
if multipole:
settings.temperature = {'multipole': True, 'tolerance': 1000}
settings.export_to_xml()
#### Create OpenMC "plots.xml" file
plots.export_to_xml()
#### Create OpenMC MGXS libraries
# Get all cells filled with a "fuel" material
mat_cells = geometry.get_all_material_cells()
fuel_cells = []
for cell in mat_cells:
if 'fuel' in cell.fill.name.lower():
fuel_cells.append(cell)
# CASMO 8-group structure
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.array([0., 0.058e-6, 0.14e-6, 0.28e-6,
0.625e-6, 4.e-6, 5.53e-3, 821.e-3, 20.])
# Initialize a 70-group "distribcell" MGXS library
cell_mgxs_lib = openmc.mgxs.Library(geometry, by_nuclide=True)
cell_mgxs_lib.energy_groups = energy_groups
cell_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
cell_mgxs_lib.domain_type = 'distribcell'
cell_mgxs_lib.domains = fuel_cells
cell_mgxs_lib.correction = None
cell_mgxs_lib.build_library()
# Initialize a 70-group "material" MGXS library
mat_mgxs_lib = openmc.mgxs.Library(geometry, by_nuclide=True)
mat_mgxs_lib.energy_groups = energy_groups
mat_mgxs_lib.mgxs_types = ['total', 'nu-fission', 'nu-scatter matrix', 'chi']
mat_mgxs_lib.domain_type = 'material'
mat_mgxs_lib.correction = None
mat_mgxs_lib.build_library()
#### Create mesh tallies for verification of pin-wise reaction rates
# Instantiate a tally Mesh
mesh = openmc.Mesh(name='assembly mesh')
mesh.type = 'regular'
mesh.dimension = [7*17, 7*17, 100]
mesh.lower_left = lower_left
mesh.width = (np.array(upper_right) - np.array(lower_left))
mesh.width[:2] /= (7*17)
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate energy-integrated fission rate mesh Tally
fission_rates = openmc.Tally(name='fission rates')
fission_rates.filters = [mesh_filter]
fission_rates.scores = ['fission']
# Instantiate energy-wise U-238 capture rate mesh Tally
capture_rates = openmc.Tally(name='u-238 capture')
capture_rates.filters = [mesh_filter]
capture_rates.nuclides = ['U238']
capture_rates.scores = ['absorption', 'fission']
#### Create OpenMC "tallies.xml" file
# Create a "tallies.xml" file for the mesh tallies
tallies_file = openmc.Tallies([fission_rates, capture_rates])
cell_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
mat_mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()