I’ve been attempting to generate a full cross section set using OpenMC and I’ve run into some strange issues with the physicality of cross sections. For some reason, the multigroup cross sections generated by the script below have a nonzero fission cross section for O-16 despite it being nonfissile. I checked the nuclear data to make sure it wasn’t incorrect, and surely enough there is no MT 18 reaction for O-16 in the O16.h5 file. In addition, the fission cross section for U-235 is zero.
It seems like there’s an issue with the way the fission cross section tallies are getting assigned to nuclides, because there are only 4 fissile isotopes (correctly corroborated by the chi cross sections), and there are 4 nonzero fission tallies, just for the wrong isotopes. I’m tempted to make this into an issue on GitHub, but wanted to verify that I’m not making any mistakes first.
from math import pi
import openmc
import openmc.deplete
import matplotlib.pyplot as plt
import numpy as np
import openmc.mgxs as mgxs
###############################################################################
# Define materials
###############################################################################
# Instantiate some Materials and register the appropriate Nuclides
uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)
uo2.add_element('Pu', 0.001)
helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)
zircaloy = openmc.Material(name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Sn', 0.014, 'wo')
zircaloy.add_element('Fe', 0.00165, 'wo')
zircaloy.add_element('Cr', 0.001, 'wo')
zircaloy.add_element('Zr', 0.98335, 'wo')
borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')
###############################################################################
# Create geometry
###############################################################################
# Define surfaces
pitch = 1.25984
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
# Define cells
fuel = openmc.Cell(fill=uo2, region=-fuel_or, name='fuel')
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir, name='gap')
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or, name='clad')
water = openmc.Cell(fill=borated_water, region=+clad_or & -box, name='water')
# Define overall geometry
geometry = openmc.Geometry([fuel, gap, clad, water])
###############################################################################
# Set volumes of depletable materials
###############################################################################
# Set material volume for depletion. For 2D simulations, this should be an area.
uo2.volume = pi * fuel_or.r**2
###############################################################################
# Transport calculation settings
###############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
# Create an initial uniform spatial source distribution over fissionable zones
# bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
# uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
# settings.source = openmc.IndependentSource(
# space=uniform_dist, constraints={'fissionable': True})
# entropy_mesh = openmc.RegularMesh()
# entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
# entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
# entropy_mesh.dimension = [10, 10, 1]
# settings.entropy_mesh = entropy_mesh
###############################################################################
# Initialize cross section tallies
###############################################################################
# Instantiate a 2-group EnergyGroups object
group_edges = np.array([0., 0.625, 20.0e6])
groups = mgxs.EnergyGroups(group_edges)
# Instantiate a few different sections
total = mgxs.TotalXS(domain=fuel, energy_groups=groups, name='total')
chi = mgxs.Chi(domain=fuel, energy_groups=groups, name='chi', by_nuclide=True)
fission = mgxs.FissionXS(domain=fuel, energy_groups=groups, name='fission', by_nuclide=True)
absorption = mgxs.AbsorptionXS(domain=fuel, energy_groups=groups, name='absorption')
scattering = mgxs.ScatterMatrixXS(domain=fuel, energy_groups=groups, name='scatter')
# Instantiate an empty Tallies object
tallies_file = openmc.Tallies()
tallies_file += chi.tallies.values()
tallies_file += fission.tallies.values()
# Add total tallies to the tallies file
tallies_file += total.tallies.values()
# Add absorption tallies to the tallies file
tallies_file += absorption.tallies.values()
# Add scattering tallies to the tallies file
tallies_file += scattering.tallies.values()
# Export to "tallies.xml"
tallies_file.export_to_xml()
###############################################################################
# Initialize and run depletion calculation
###############################################################################
model = openmc.Model(geometry=geometry, settings=settings)
# Create depletion "operator"
# openmc.config['chain_file'] = 'chain_simple.xml'
openmc.config['chain_file'] = '/projects/MCRE_studies/louime/data/depletion/chain_endfb80_sfr.xml'
op = openmc.deplete.CoupledOperator(model)
# Perform simulation using the predictor algorithm
time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days
power = 1740 # W/cm, for 2D simulations only (use W for 3D)
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
# integrator.integrate()
###############################################################################
# Read depletion calculation results
###############################################################################
# Open results file
results = openmc.deplete.Results("depletion_results.h5")
# Obtain K_eff as a function of time
time, keff = results.get_keff(time_units='d')
# Obtain U235 concentration as a function of time
_, n_U235 = results.get_atoms(uo2, 'U235')
# Obtain Xe135 capture reaction rate as a function of time
_, Xe_capture = results.get_reaction_rate(uo2, 'Xe135', '(n,gamma)')
###############################################################################
# Generate plots
###############################################################################
fig, ax = plt.subplots()
ax.errorbar(time, keff[:, 0], keff[:, 1], label="K-effective")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Keff")
plt.savefig('keff.png')
fig, ax = plt.subplots()
ax.plot(time, n_U235, label="U235")
ax.set_xlabel("Time [d]")
ax.set_ylabel("U235 atoms")
plt.savefig('U235.png')
fig, ax = plt.subplots()
ax.plot(time, Xe_capture, label="Xe135 capture")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Xe135 capture rate")
plt.savefig('XeCapture.png')
sp = openmc.StatePoint('openmc_simulation_n4.h5')
# Load the tallies from the statepoint into each MGXS object
total.load_from_statepoint(sp)
chi.load_from_statepoint(sp)
fission.load_from_statepoint(sp)
absorption.load_from_statepoint(sp)
scattering.load_from_statepoint(sp)
chi.print_xs()
fission.print_xs()
print(fission.get_xs(nuclides=['O16'], xs_type='micro'))