# On the Deviation of Results in Using OpenMC to Produce Multi group Cross Sections for Openmoc in the Presence of Vacuum Boundaries

I need to use OpenMC to create a cross-section library of a single three-dimensional fuel rod for use by OpenMOC. When all boundary conditions are reflection, I can use the cross-section library produced to calculate good results in OpenMOC. However, when the upper and lower surface boundary conditions are vacuum, the error reaches several thousand pcm. Can someone tell me the reason for this.Here is the code,

import numpy as np
import matplotlib.pyplot as plt
plt.style.use(‘seaborn-dark’)

import openmoc

import openmc
import openmc.mgxs as mgxs
import openmc.model as model
import openmc.data
from openmc.openmoc_compatible import get_openmoc_geometry

# create a model object to tie together geometry, materials, settings, and tallies

#model = openmc.model()
model = model.Model()

# 1.6% enriched fuel

fuel = openmc.Material(name=‘3.1% Fuel’)
fuel.set_density(‘g/cm3’, 10.257)

# borated water

water = openmc.Material(name=‘Water’)
water.set_density(‘g/cm3’, 0.743)

# zircaloy

gap = openmc.Material(name=‘gap’)
gap.set_density(‘g/cm3’, 0.0001786)

# zircaloy

zircaloy.set_density(‘g/cm3’, 6.56)

# Instantiate a Materials collection

model.materials = openmc.Materials([fuel, gap ,water, zircaloy])

# Create box to surround the geometry

#box = openmc.model.rectangular_prism(1.26, 1.26, boundary_type=‘reflective’)
min_x = openmc.XPlane(x0=-0.63, boundary_type=‘reflective’)
max_x = openmc.XPlane(x0=+0.63, boundary_type=‘reflective’)
min_y = openmc.YPlane(y0=-0.63, boundary_type=‘reflective’)
max_y = openmc.YPlane(y0=+0.63, boundary_type=‘reflective’)
min_z = openmc.ZPlane(z0=-182.88, boundary_type=‘vacuum’)
max_z = openmc.ZPlane(z0=+182.88, boundary_type=‘vacuum’)

# Create a Universe to encapsulate a fuel pin

pin_cell_universe = openmc.Universe(name=‘1.6% Fuel Pin’)

# Create fuel Cell

fuel_cell = openmc.Cell(name=‘1.6% Fuel’)
fuel_cell.fill = fuel

# Create a gap Cell

gap_cell.fill = gap

# Create a moderator Cell

moderator_cell = openmc.Cell(name=‘1.6% Moderator’)
moderator_cell.fill = water

# Create root Cell

root_cell = openmc.Cell(name=‘root cell’)
root_cell.fill = pin_cell_universe

root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z

# Create root Universe

root_universe = openmc.Universe(universe_id=0, name=‘root universe’)

# Create Geometry and set root Universe

model.geometry = openmc.Geometry(root_universe)

# OpenMC simulation parameters

batches = 500
inactive = 50
particles = 10000

# Instantiate a Settings object

settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.output = {‘tallies’: False}

# Create an initial uniform spatial source distribution over fissionable zones

bounds = [-1.26, -1.26, -1.26, 1.26, 1.26, 1.26]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)

# Activate tally precision triggers

settings.trigger_active = True
settings.trigger_max_batches = settings.batches * 4

model.settings = settings

# Instantiate a “fine” 8-group EnergyGroups object

fine_groups = mgxs.EnergyGroups([0., 0.058, 0.14, 0.28,
0.625, 4.0, 5.53e3, 821.0e3, 20.0e6])

# Extract all Cells filled by Materials

openmc_cells = model.geometry.get_all_material_cells().values()

xs_library = {}

# Instantiate 8-group cross sections for each cell

for cell in openmc_cells:
xs_library[cell.id] = {}
xs_library[cell.id][‘total’] = mgxs.TotalXS(groups=fine_groups)
xs_library[cell.id][‘transport’] = mgxs.TransportXS(groups=fine_groups)
xs_library[cell.id][‘fission’] = mgxs.FissionXS(groups=fine_groups)
xs_library[cell.id][‘nu-fission’] = mgxs.FissionXS(groups=fine_groups, nu=True)
xs_library[cell.id][‘nu-scatter’] = mgxs.ScatterMatrixXS(groups=fine_groups, nu=True)
xs_library[cell.id][‘chi’] = mgxs.Chi(groups=fine_groups)

# Create a tally trigger for +/- 0.01 on each tally used to compute the multi-group cross sections

tally_trigger = openmc.Trigger(‘std_dev’, 1e-2)

# Add the tally trigger to each of the multi-group cross section tallies

for cell in openmc_cells:
for mgxs_type in xs_library[cell.id]:
xs_library[cell.id][mgxs_type].tally_trigger = tally_trigger

# Instantiate an empty Tallies object

tallies = openmc.Tallies()

# Iterate over all cells and cross section types

for cell in openmc_cells:
for rxn_type in xs_library[cell.id]:

``````    # Set the cross sections domain to the cell
xs_library[cell.id][rxn_type].domain = cell

# Tally cross sections by nuclide
xs_library[cell.id][rxn_type].by_nuclide = True

# Add OpenMC tallies to the tallies file for XML generation
for tally in xs_library[cell.id][rxn_type].tallies.values():
tallies.append(tally, merge=True)
``````

model.tallies = tallies

# Run OpenMC

sp_file = model.run()

# Load the last statepoint file

sp = openmc.StatePoint(sp_file)
#sp = openmc.StatePoint(‘statepoint.2000.h5’)

# Iterate over all cells and cross section types

for cell in openmc_cells:
for rxn_type in xs_library[cell.id]:

# Create an OpenMOC Geometry from the OpenMC Geometry

openmoc_geometry = get_openmoc_geometry(sp.summary.geometry)

# Get all OpenMOC cells in the gometry

openmoc_cells = openmoc_geometry.getRootUniverse().getAllCells()

# Inject multi-group cross sections into OpenMOC Materials

for cell_id, cell in openmoc_cells.items():

``````# Ignore the root cell
if cell.getName() == 'root cell':
continue

# Get a reference to the Material filling this Cell
openmoc_material = cell.getFillMaterial()

# Set the number of energy groups for the Material
openmoc_material.setNumEnergyGroups(fine_groups.num_groups)

# Extract the appropriate cross section objects for this cell
transport = xs_library[cell_id]['transport']
nufission = xs_library[cell_id]['nu-fission']
nuscatter = xs_library[cell_id]['nu-scatter']
chi = xs_library[cell_id]['chi']

# Inject NumPy arrays of cross section data into the Material
# NOTE: Sum across nuclides to get macro cross sections needed by OpenMOC
openmoc_material.setSigmaT(transport.get_xs(nuclides='sum').flatten())
openmoc_material.setNuSigmaF(nufission.get_xs(nuclides='sum').flatten())
openmoc_material.setSigmaS0(nuscatter.get_xs(nuclides='sum').flatten())
openmoc_material.setChi(chi.get_xs(nuclides='sum').flatten())
``````

# Generate tracks for OpenMOC

track_generator = openmoc.TrackGenerator3D(openmoc_geometry, num_azim=16, num_polar=6,
azim_spacing=0.1, z_spacing=0.75)
track_generator.generateTracks()

# Run OpenMOC

solver = openmoc.CPUSolver(track_generator)
solver.computeEigenvalue()
solver.setConvergenceThreshold(1e-6)

# Print report of keff and bias with OpenMC

openmoc_keff = solver.getKeff()
#openmc_keff = sp.keff.n
openmc_keff = sp.k_combined.nominal_value
bias = (openmoc_keff - openmc_keff) * 1e5

print(‘openmc keff = {0:1.6f}’.format(openmc_keff))
print(‘openmoc keff = {0:1.6f}’.format(openmoc_keff))
print(‘bias [pcm]: {0:1.1f}’.format(bias))

The results of running the above code,
openmc_keff = 1.198678
openmoc_keff = 1.127340
bias [pcm]: -7133.8