I discovered that the only way I could use the CellBornFilter was by manually defining each pin as an individual cell. For example, when modeling a single assembly of a reactor core, I had to explicitly define each fuel pin as a separate cell instead of using commands that generate new instances of a base cell. This approach makes geometry construction more cumbersome, but it seems to work.
However, the problem arises again when I try to model a system where the fuel pins contain a dispersion of TRISO particles. In this case, even if I manually define the pin cells, the pins themselves are filled with a lattice of TRISO particles. As a result, I cannot extract tally information using the cellFilter or CellBornFilter at the pin level. I suspect that the elementary cells in this case are the individual TRISO particles, which makes data extraction particularly inconvenient.
Does anyone have any suggestions or guidance on how to handle this issue and, specifically, on how to extract the score of interest at the ‘pin level’ instead of the ‘TRISO-particle’ level? Probably, the more correct and general question is: Is there a way to use the CellBornFilter also for cells that are filled with universes or lattices? Below, I’ve included a simple case to better illustrate my problem.
def hex_assembly_builder(pitch=1.0, center=[0,0], orientation='y', nrings=1):
dic_info = {}
dic_info['central'] = {}
dic_info['central']["pin:1"] = (center[0],center[1])
for ring in range(nrings):
dic_info[f'ring:{ring+1}'] = {}
hex_side = pitch*(ring+1)
#print(hex_side)
starting_points = [
(center[0] + hex_side*np.cos(np.pi/2) , center[1] + hex_side*np.sin(np.pi/2)),
(center[0] + hex_side*np.cos(np.pi/6) , center[1] + hex_side*np.sin(np.pi/6)),
(center[0] + hex_side*np.cos(-np.pi/6) , center[1] + hex_side*np.sin(-np.pi/6)),
(center[0] + hex_side*np.cos(-np.pi/2) , center[1] + hex_side*np.sin(-np.pi/2)),
(center[0] + hex_side*np.cos(np.pi*7/6) , center[1] + hex_side*np.sin(np.pi*7/6)),
(center[0] + hex_side*np.cos(np.pi*5/6) , center[1] + hex_side*np.sin(np.pi*5/6)),
]
if orientation == 'x':
theta = np.radians(30)
R = np.array([[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
starting_points = np.dot(starting_points,R.T)
t = np.linspace(0, 1, ring+2)
xc=[]
yc=[]
for i in range(6):
x0, y0 = starting_points[i]
x1, y1 = starting_points[(i + 1)%6]
if i == 0:
x_values = x0 + (x1 - x0) * t
y_values = y0 + (y1 - y0) * t
elif i==5:
x_values = x0 + (x1 - x0) * t[1:-1]
y_values = y0 + (y1 - y0) * t[1:-1]
else:
x_values = x0 + (x1 - x0) * t[1:]
y_values = y0 + (y1 - y0) * t[1:]
xc.append(x_values)
yc.append(y_values)
xc = np.concatenate(xc)
yc= np.concatenate(yc)
#print(len(xc))
for ival in range(len(xc)):
dic_info[f"ring:{ring+1}"][f"pin:{ival+1}"] = (xc[ival],yc[ival])
return dic_info
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 10000
inactive = 4000
particles = 20000
TRISO_temp = 900 # fuel particle temperature in [K]
mod_temp = 600 # moderator and reactor bulk temperature in [K]
################################################################################
# Exporting to OpenMC materials.xml file
################################################################################
openmc.Materials.cross_sections = r"INSERT HERE YOUR cross_sections.xml FILE PATH"
fuel_25 = openmc.Material(material_id=125, name='UCO Fuel') # from AGR-2 Fuel Compact Report
fuel_25.add_element('U' , 1.0, enrichment=15.00)
fuel_25.add_nuclide('O16', 1.428)
fuel_25.add_element('C' , 0.392)
fuel_25.set_density('g/cm3', 10.97)
fuel_25.temperature = TRISO_temp
buff_mat = openmc.Material(material_id=2, name='buff_mat')
buff_mat.add_element('C', 1.0)
buff_mat.set_density('g/cm3', 0.95)
# buff_mat.add_s_alpha_beta('c_Graphite')
buff_mat.temperature = TRISO_temp
pyc1 = openmc.Material(material_id=3, name='Inner PyC')
pyc1.add_element('C', 1.0)
pyc1.set_density('g/cm3', 1.90)
# pyc1.add_s_alpha_beta('c_Graphite')
pyc1.temperature = TRISO_temp
sic = openmc.Material(material_id=4, name = 'sic')
sic.add_element('Si', 0.5)
sic.add_element('C', 0.5)
sic.set_density('g/cm3', 3.18)
sic.temperature = TRISO_temp
pyc2 = openmc.Material(material_id=5, name='Outer PyC')
pyc2.add_element('C', 1.0)
pyc2.set_density('g/cm3', 1.90)
# pyc2.add_s_alpha_beta('c_Graphite')
pyc2.temperature = TRISO_temp
graph = openmc.Material(material_id=6, name='Graphite')
graph.add_element('C', 1.0)
graph.set_density('g/cm3', 1.70)
# graph.add_s_alpha_beta('c_Graphite')
graph.temperature = mod_temp
graph_cnt = openmc.Material(material_id=60, name='Graphite Cnt')
graph_cnt.add_element('C', 1.0)
graph_cnt.set_density('g/cm3', 1.70)
# graph_cnt.add_s_alpha_beta('c_Graphite')
graph_cnt.temperature = mod_temp
graph_ne = openmc.Material(material_id=600, name='Graphite NE')
graph_ne.add_element('C', 1.0)
graph_ne.set_density('g/cm3', 1.70)
# graph_ne.add_s_alpha_beta('c_Graphite')
graph_ne.temperature = mod_temp
# Instantiate a Materials collection and export to xml
materials_list = [fuel_25, fuel_35, buff_mat, pyc1, sic, pyc2, graph, graph_cnt, graph_ne]
materials_file = openmc.Materials(materials_list)
materials_file.export_to_xml()
################################################################################
# Exporting to OpenMC geometry.xml file
################################################################################
# Geometry variables
r_pellet = 1.75/2 # fuel compact radius [cm]
# Define TRISO universes
kernelsph = openmc.Sphere(r=250e-4) # try 212.5um and 250um
buffsph = openmc.Sphere(r=350e-4)
pyc1sph = openmc.Sphere(r=390e-4)
sicsph = openmc.Sphere(r=425e-4)
pyc2sph = openmc.Sphere(r=470e-4)
triso_outer_radius = 465e-4
layers = [kernelsph, buffsph, pyc1sph, sicsph, pyc2sph]
triso_mats_25 = [fuel_25, buff_mat, pyc1, sic, pyc2]
triso_cells_25 = []
triso_cells_25.append(openmc.Cell(fill=triso_mats_25[0], region=-layers[0]))
triso_cells_25.append(openmc.Cell(fill=triso_mats_25[1], region=+layers[0] & -layers[1]))
triso_cells_25.append(openmc.Cell(fill=triso_mats_25[2], region=+layers[1] & -layers[2]))
triso_cells_25.append(openmc.Cell(fill=triso_mats_25[3], region=+layers[2] & -layers[3]))
triso_cells_25.append(openmc.Cell(fill=triso_mats_25[4], region=+layers[3] & -layers[4]))
triso_universe_25 = openmc.Universe(cells=triso_cells_25)
dic_pin_xy = hex_assembly_builder(pitch = pin_pitch, nrings=1)
pin_surfaces_list = []
pin_regions_list = []
pin_cells_list = []
pin_cells_outer_list = []
TRISO_particles_list = []
pin_counter = 0
for i_key_ring,key_ring in enumerate(list(dic_pin_xy.keys())):
for i_key_pin,key_pin in enumerate(list(dic_pin_xy[key_ring].keys())):
r_pin = openmc.ZCylinder(r=r_pellet, x0=dic_pin_xy[key_ring][key_pin][0], y0=dic_pin_xy[key_ring][key_pin][1])
maxz_it = openmc.ZPlane(z0=+core_height/2, boundary_type='reflective')
minz_it = openmc.ZPlane(z0=-core_height/2, boundary_type='reflective')
pin_inside_region = -r_pin & +minz_it & -maxz_it
lower_left, upp_right = pin_inside_region.bounding_box
shape = (4, 4, 4)
pitch = (upp_right - lower_left)/shape
ceneters_25_for_iteration = openmc.model.pack_spheres(radius=triso_outer_radius, region=pin_inside_region, pf=0.25)
#print("\n", "centers moved: ", ceneters_25_for_iteration)
triso_particles_25_for_iteration = [openmc.model.TRISO(triso_outer_radius, fill=triso_universe_25, center=c) for c in ceneters_25_for_iteration]
triso_latt_25_for_iteration = openmc.model.create_triso_lattice(triso_particles_25_for_iteration, lower_left, pitch, shape, graph)
triso_latt_25_for_iteration.fill = triso_latt_25_for_iteration
fuel_cell = openmc.Cell(cell_id=pin_counter,fill=triso_latt_25_for_iteration, region = pin_inside_region )
pin_counter+=1
TRISO_particles_list = TRISO_particles_list + triso_particles_25_for_iteration
pin_surfaces_list.append(r_pin)
pin_regions_list.append(pin_inside_region)
pin_universe = openmc.Universe(cells=(fuel_cell,))
pin_cells_list.append(fuel_cell)
outer_radial_surface = openmc.ZCylinder(r=15,boundary_type='vacuum')
graphite_matrix_region = -outer_radial_surface
for reg in pin_regions_list:
graphite_matrix_region = graphite_matrix_region & ~reg #+reg & -minz & +maxz
graphite_cell = openmc.Cell(fill=graph,region=graphite_matrix_region)
root_universe = openmc.Universe(cells = pin_cells_list + [graphite_cell] )
# Define root universe
geometry = openmc.Geometry(root = root_universe)
geometry.export_to_xml()
p=assembly_pitch
s=p/np.sqrt(3)
core = openmc.Plot(plot_id=5)
core.filename = 'core'
core.basis = 'xy'
core.level = 1006
core.origin = [0, 0, 0]
core.width = [20,20]
core.pixels = [500, 500]
core.color_by = 'material'
core.colors = {fuel_25:'green', fuel_35:'yellow', buff_mat:'lightgrey', pyc1:'lightgrey', sic:'lightgrey',
pyc2:'lightgrey', graph:'blue'}
plots = openmc.Plots([core])
plots.export_to_xml()
openmc.plot_geometry()
### DEFINE TALLIES ###
cell_filter = openmc.CellFilter(pin_cells_list)
cellborn_filter = openmc.CellBornFilter(pin_cells_list)
universe_system_filter = openmc.UniverseFilter(root_universe)
tally_pfm = openmc.Tally()
tally_pfm.filters = [cell_filter,cellborn_filter]
tally_pfm.scores = ['nu-fission']
tallies = openmc.Tallies([tally_pfm])
tallies.export_to_xml()
### SETTINGS ###
settings = openmc.Settings()
settings.batches = 25
settings.inactive = 5
settings.particles = 10000
settings.source = openmc.Source(space=openmc.stats.Point((0., 0., 0.)))
settings.export_to_xml()
openmc.run()`
`