Hello all,
I am encountering an issue in collecting MGXS data. This happens after I run a depletion simulation to generate a depletion results file and then run separate k-eigenvalue simulations reading in material compositions at various burnup indeces from those depletion results. The reason I am running things this way is because I want summary and statepoint files for each depletion step (so if someone knows an easier way to do this, that would also be very helpful). The error I am encountering is a LookupError when I try to load data from the statepoint file. I only get this error when I use a non-zero burnup index, which makes me think that perhaps altering the materials.xml file is affecting the MGXS tallies (though the fission tally I included appears to work fine)? I’ve included a simplified script here based off the pincell depletion example that demonstrates the specific issue. Ideally, I would like to be able to run this script with do_depletion=True initially, then set it to False and run do_eigen=True for each burnup index with no errors. Any help would be much appreciated.
Thank you,
Luke Seifert
from math import pi
import openmc
import openmc.deplete
import matplotlib.pyplot as plt
import numpy as np
import openmc.mgxs as mgxs
# Modify these values
burnup_index = 0
chain_file = '../data/chain/chain_endfb80_sfr.xml'
xs_xml = "/home/seifluke/xs/endfb-viii.0-hdf5/cross_sections.xml"
do_depletion = True
do_keigen = True
###############################################################################
# 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.)
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)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
# 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
###############################################################################
# Add tallies
###############################################################################
tallies_file = openmc.Tallies()
energy_list=[0.0, 748.5, 5531.0, 24790.0, 497900.0, 2231000.0, 20000000.0]
energy_groups = mgxs.EnergyGroups(group_edges=energy_list)
mgxs_lib = mgxs.Library(geometry)
mgxs_lib.energy_groups = energy_groups
mgxs_lib.legendre_order = 1
mgxs_lib.num_delayed_groups = 6
mgxs_lib.mgxs_types = ['transport','total','capture','absorption','fission','nu-fission',
'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix',
'nu-scatter matrix', 'scatter probability matrix',
'consistent scatter matrix', 'consistent nu-scatter matrix',
'inverse-velocity', 'chi', 'chi-prompt', 'prompt-nu-fission',
'chi-delayed', 'delayed-nu-fission','beta', 'decay-rate',
'prompt-nu-fission matrix', 'delayed-nu-fission matrix',
'diffusion-coefficient', 'nu-diffusion-coefficient']
mgxs_lib.domain_type = 'material'
mgxs_lib.domains = geometry.get_all_materials().values()
mgxs_lib.by_nuclide = True
mgxs_lib.xs_type = 'micro'
library = openmc.data.DataLibrary.from_xml(xs_xml)
nuclide_data_list = [i['materials'][0] for i in library.libraries][:556]
mgxs_lib.nuclides = nuclide_data_list
mgxs_lib.build_library()
mgxs_lib.add_to_tallies_file(tallies_file, merge=False)
mesh = openmc.RegularMesh()
mesh.dimension = [1, 1, 1]
mesh.lower_left = np.array([-0.39218, -0.39218, -1.e50] )
mesh.upper_right = np.array([0.39218, 0.39218, 1.e50] )
mesh_filter = openmc.MeshFilter(mesh)
score = 'fission'
name = 'fission'
tally = openmc.Tally(name=name)
tally.filters = [mesh_filter]
tally.scores = [score]
tallies_file.append(tally)
tallies_file.export_to_xml()
###############################################################################
# Initialize and run depletion calculation
###############################################################################
model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies_file)
# Create depletion "operator"
op = openmc.deplete.CoupledOperator(model, chain_file)
# Perform simulation using the predictor algorithm
time_steps = [10, 10, 10] # days
power = 174 # W/cm, for 2D simulations only (use W for 3D)
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
if do_depletion:
integrator.integrate()
###############################################################################
# Run keigenvalue solve at a specified burnup
###############################################################################
# Open results file
results = openmc.deplete.Results("depletion_results.h5")
# Update materials with depleted material composition
materials = results.export_to_materials(burnup_index=burnup_index)
model = openmc.Model.from_xml()
model.materials = materials
model.materials.export_to_xml()
if do_keigen:
openmc.run()
# Show that non-MGXS tally is working
fiss_tally = sp.get_tally(scores=['fission'])
print(fiss_tally.get_values())
###############################################################################
# Evaluate MGXS data
###############################################################################
delayed_groups=list(range(1,7))
material_xml_path='./materials.xml'
geometry_xml_path='./geometry.xml'
statepoint_path = f'./statepoint.100.h5'
sp = openmc.StatePoint(statepoint_path)
geom = openmc.Geometry.from_xml(geometry_xml_path, material_xml_path)
domain = openmc.Materials().from_xml(material_xml_path)[0]
energy_groups = openmc.mgxs.EnergyGroups(group_edges=energy_list)
one_group = mgxs.EnergyGroups(group_edges=np.array([energy_groups.group_edges[0], energy_groups.group_edges[-1]]))
beta = mgxs.Beta(domain=domain, energy_groups=energy_groups, delayed_groups=delayed_groups, by_nuclide=True)
try:
beta.load_from_statepoint(sp)
print('Pass')
except LookupError:
print('Fail')