Error when running depletion calculation

Hi,

I’m running into this error every time I try to run a depletion calculation with ‘diff_burnable_mats=True’:

ValueError: Number of material instances have not been determined. Call the Geometry.determine_paths() method.

How can I fix this?

I can’t attach the code, so I will paste it here, hope this is not a problem!

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 20 10:23:41 2022

@author: Marti
"""
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 18 12:58:18 2022

@author: Marti
"""
from math import pi, sin, cos, sqrt
import numpy as np
import matplotlib.pyplot as plt
import openmc
import openmc.deplete

 

# Materials definitions
openmc.Materials.cross_sections = "/Users/Marti/Downloads/endfb71_hdf5/cross_sections.xml"
inner_fuel = openmc.Material(name='Inner fuel')
inner_fuel.add_nuclide('U235', 0.19, 'wo')
inner_fuel.add_nuclide('U238', 75.09, 'wo')
inner_fuel.add_nuclide('Pu238', 0.46, 'wo')
inner_fuel.add_nuclide('Pu239', 6.12, 'wo')
inner_fuel.add_nuclide('Pu240', 3.83, 'wo')
inner_fuel.add_nuclide('Pu241', 1.06, 'wo')
inner_fuel.add_nuclide('Am241', 0.10, 'wo')
inner_fuel.add_nuclide('O16', 11.81, 'wo')
inner_fuel.set_density('kg/m3', 10661.4)

outer_fuel = openmc.Material(name='Outer fuel')
outer_fuel.add_nuclide('U235', 0.18, 'wo')
outer_fuel.add_nuclide('U238', 73.0, 'wo')
outer_fuel.add_nuclide('Pu238', 0.54, 'wo')
outer_fuel.add_nuclide('Pu239', 7.11, 'wo')
outer_fuel.add_nuclide('Pu240', 4.45, 'wo')
outer_fuel.add_nuclide('Pu241', 1.24, 'wo')
outer_fuel.add_nuclide('Am241', 0.17, 'wo')
outer_fuel.add_nuclide('O16', 11.75, 'wo')
outer_fuel.set_density('kg/m3', 10661.4)

clad = openmc.Material(name= 'Clad')
clad.add_nuclide('C0', 0.1, 'wo')
clad.add_element('Mn', 0.3, 'wo')
clad.add_element('P', 0.02, 'wo')
clad.add_element('Ni', 0.5, 'wo')
clad.add_element('Cr', 21.5, 'wo')
clad.add_element('Ti', 0.6, 'wo')
clad.add_element('Co', 0.3, 'wo')
clad.add_element('Cu', 0.15, 'wo')
clad.add_element('Al', 5.75, 'wo')
clad.add_nuclide('O16', 0.15, 'wo')
clad.add_element('Y', 0.28, 'wo')
clad.add_element('Fe', 70.08, 'wo')
clad.set_density('kg/m3', 7250.0)

sodium = openmc.Material(name = 'Sodium')
sodium.add_nuclide('Na23', 1.0)
sodium.set_density('kg/m3', 830.0)

boron_carbide = openmc.Material(name='Boron Carbide')
boron_carbide.add_nuclide('B10', 0.198, 'wo')
boron_carbide.add_nuclide('C0', 0.802, 'wo')
boron_carbide.set_density('kg/m3', 2450.0)

helium = openmc.Material(name = 'Helium')
helium.add_nuclide('He4', 1.0)
helium.set_density('kg/m3', 0.4282)

reflector_sodium = openmc.Material(name = 'Reflector Sodium')
reflector_sodium.add_nuclide('Na23', 1.0)
reflector_sodium.set_density('kg/m3', 833.0)

em10_steel = openmc.Material(name = 'Wrapper')
em10_steel.add_nuclide('C0', 0.1, 'wo')
em10_steel.add_element('Cr', 8.5, 'wo')
em10_steel.add_element('Mo', 1.0, 'wo')
em10_steel.add_element('Ni', 0.5, 'wo')
em10_steel.add_element('Mn', 0.5, 'wo')
em10_steel.add_element('Nb', 0.02, 'wo')
em10_steel.add_element('Ti', 0.02, 'wo')
em10_steel.add_element('W', 0.02, 'wo')
em10_steel.add_element('Si', 0.2, 'wo')
em10_steel.add_element('Fe', 89.0, 'wo')
em10_steel.set_density('kg/m3', 7753.0)

reflector_material = openmc.Material.mix_materials([sodium, em10_steel],[0.08, 0.92],'wo')

materials_file = openmc.Materials([inner_fuel, outer_fuel, sodium, em10_steel,
                                   reflector_material, reflector_sodium, helium,
                                   boron_carbide, clad])
materials_file.export_to_xml()

#Geometry definitions

sleeve_min_z = openmc.ZPlane(z0=-50, boundary_type = 'transmission')
sleeve_max_z = openmc.ZPlane(z0=+50, boundary_type = 'transmission')
axial_reflector_max_z = openmc.ZPlane(z0= +131, boundary_type = 'vacuum')
axial_reflector_min_z = openmc.ZPlane(z0= -80, boundary_type = 'transmission') 
gas_plenum_max_z = openmc.ZPlane(z0 = +61, boundary_type = 'transmission') 
gas_plenum_min_z = openmc.ZPlane (z0 = -171, boundary_type = 'vacuum')

fuel_radius_inner = openmc.ZCylinder(r=0.943/2)
clad_inner_radius_inner = openmc.ZCylinder(r=0.973/2)
clad_outer_radius_inner = openmc.ZCylinder(r=1.073/2)

fuel_region_inner = -fuel_radius_inner & -sleeve_max_z & +sleeve_min_z
gap_region_inner  = +fuel_radius_inner & -clad_inner_radius_inner  & -sleeve_max_z & +sleeve_min_z
clad_region_inner = +clad_inner_radius_inner & -clad_outer_radius_inner & -axial_reflector_max_z & +gas_plenum_min_z
axial_reflector_region_top_inner = -clad_inner_radius_inner & -axial_reflector_max_z & +gas_plenum_max_z
axial_reflector_region_bottom_inner = -clad_inner_radius_inner & +axial_reflector_min_z & -sleeve_min_z
gas_plenum_region_top_inner = -clad_inner_radius_inner & +sleeve_max_z & -gas_plenum_max_z
gas_plenum_region_bottom_inner = -clad_inner_radius_inner & +gas_plenum_min_z & -axial_reflector_min_z
surrounding_region_inner = +clad_outer_radius_inner & -axial_reflector_max_z & +gas_plenum_min_z

inner_fuel_cell = openmc.Cell(fill=inner_fuel, region=fuel_region_inner)
gap_cell_inner = openmc.Cell(fill=helium, region=gap_region_inner)
clad_cell_inner = openmc.Cell(fill=clad, region=clad_region_inner)
axial_reflector_top_cell_inner = openmc.Cell(fill=reflector_material,
                                             region=axial_reflector_region_top_inner)
axial_reflector_bottom_cell_inner = openmc.Cell(fill=reflector_material,
                                                region=axial_reflector_region_bottom_inner)
gas_plenum_cell_top_inner = openmc.Cell(fill = helium, region=gas_plenum_region_top_inner)
gas_plenum_cell_bottom_inner = openmc.Cell(fill = helium, 
                                           region=gas_plenum_region_bottom_inner) #
surrounding_cell_inner = openmc.Cell(fill=sodium, region=surrounding_region_inner)

inner_universe = openmc.Universe(cells=(inner_fuel_cell, gap_cell_inner, clad_cell_inner,
                                        axial_reflector_top_cell_inner, 
                                        axial_reflector_bottom_cell_inner, 
                                        gas_plenum_cell_top_inner,
                                        gas_plenum_cell_bottom_inner,
                                        surrounding_cell_inner))

fuel_radius_outer = openmc.ZCylinder(r=0.943/2)
clad_inner_radius_outer = openmc.ZCylinder(r=0.973/2)
clad_outer_radius_outer = openmc.ZCylinder(r=1.073/2)

fuel_region_outer = -fuel_radius_outer & -sleeve_max_z & +sleeve_min_z
gap_region_outer  = +fuel_radius_outer & -clad_inner_radius_outer  & -sleeve_max_z & +sleeve_min_z
clad_region_outer = +clad_inner_radius_outer & -clad_outer_radius_outer & -axial_reflector_max_z & +gas_plenum_min_z
axial_reflector_region_top_outer = -clad_inner_radius_outer & -axial_reflector_max_z & +gas_plenum_max_z
axial_reflector_region_bottom_outer = -clad_inner_radius_outer & +axial_reflector_min_z & -sleeve_min_z
gas_plenum_region_top_outer = -clad_inner_radius_outer & +sleeve_max_z & -gas_plenum_max_z
gas_plenum_region_bottom_outer = -clad_inner_radius_outer & +gas_plenum_min_z & -axial_reflector_min_z
surrounding_region_outer = +clad_outer_radius_outer & -axial_reflector_max_z & +gas_plenum_min_z

outer_fuel_cell = openmc.Cell(fill=outer_fuel, region=fuel_region_outer)
gap_cell_outer = openmc.Cell(fill=helium, region=gap_region_outer)
clad_cell_outer = openmc.Cell(fill=clad, region=clad_region_outer)
axial_reflector_top_cell_outer = openmc.Cell(fill=reflector_material,
                                             region=axial_reflector_region_top_outer)
axial_reflector_bottom_cell_outer = openmc.Cell(fill=reflector_material,
                                                region=axial_reflector_region_bottom_outer)
gas_plenum_cell_top_outer = openmc.Cell(fill = helium, region=gas_plenum_region_top_outer)
gas_plenum_cell_bottom_outer = openmc.Cell(fill = helium, 
                                           region=gas_plenum_region_bottom_outer) #
surrounding_cell_outer = openmc.Cell(fill=sodium, region=surrounding_region_outer)

outer_universe = openmc.Universe(cells=(outer_fuel_cell, gap_cell_outer, clad_cell_outer,
                                        axial_reflector_top_cell_outer, 
                                        axial_reflector_bottom_cell_outer, 
                                        gas_plenum_cell_top_outer,
                                        gas_plenum_cell_bottom_outer,
                                        surrounding_cell_outer))

sodium_fill_cell = openmc.Cell(fill=sodium)
sodium_universe = openmc.Universe(cells=(sodium_fill_cell,))

sodium_prism = openmc.hexagonal_prism(edge_length=11.91073605,
                                      orientation='x')
sodium_prism_region = sodium_prism & +gas_plenum_min_z & -axial_reflector_max_z
sodium_prism_out = ~sodium_prism & +gas_plenum_min_z & -axial_reflector_max_z
sodium_prism_cell = openmc.Cell(fill=sodium, region=sodium_prism_region)
sodium_prism_cell_out = openmc.Cell(fill=sodium, region=sodium_prism_out)
sodium_prism_universe = openmc.Universe(cells=(sodium_prism_cell, sodium_prism_cell_out))


reflector_radius = openmc.ZCylinder(r=1.9)
reflector_cladding_radius = openmc.ZCylinder(r=2.0062)

reflector_region = -reflector_radius & +gas_plenum_min_z & -axial_reflector_max_z
reflector_clad_region = (+reflector_radius & -reflector_cladding_radius & 
                         +gas_plenum_min_z & -axial_reflector_max_z)

reflector_cell = openmc.Cell(fill=reflector_material, region=reflector_region)
reflector_clad_cell = openmc.Cell(fill=em10_steel, region=reflector_clad_region)
reflector_out_cell = openmc.Cell(fill=sodium, region = +reflector_cladding_radius
                                 & +gas_plenum_min_z & -axial_reflector_max_z)
reflector_universe = openmc.Universe(cells = (reflector_cell, reflector_clad_cell,
                                              reflector_out_cell))

control_rod_radius = openmc.ZCylinder(r=0.915)
control_clad_inner_radius = openmc.ZCylinder(r=1.0415)
control_clad_outer_radius = openmc.ZCylinder(r=1.1412)

control_rod_region = (-control_rod_radius & +sleeve_max_z & -axial_reflector_max_z)
control_gap_region = (+control_rod_radius & -control_clad_inner_radius &
                      +sleeve_max_z & -axial_reflector_max_z)
control_clad_region = (+control_clad_inner_radius & -control_clad_outer_radius &
                       +gas_plenum_min_z & -axial_reflector_max_z)
follower_region = (-control_clad_outer_radius & +gas_plenum_min_z &
                   -sleeve_max_z)
control_outside_region = (+control_clad_outer_radius &
                          +gas_plenum_min_z & -axial_reflector_max_z)

control_rod_cell = openmc.Cell(fill=boron_carbide, region=control_rod_region)
control_gap_cell = openmc.Cell(fill=helium, region=control_gap_region)
control_clad_cell = openmc.Cell(fill=em10_steel, region=control_clad_region)
follower_cell = openmc.Cell(fill=reflector_material, region=follower_region)
control_surroundings_cell = openmc.Cell(fill=sodium, region=control_outside_region)

control_universe = openmc.Universe(cells=(control_rod_cell, control_gap_cell,
                                          control_clad_cell, follower_cell,
                                          control_surroundings_cell))

inner_lattice = openmc.HexLattice()
inner_lattice.center = (0., 0.)
inner_lattice.pitch = (21.08/17,)
inner_lattice.orientation = 'x'
inner_lattice.outer = sodium_universe

inner_ring_one = [inner_universe] * 48
inner_ring_two = [inner_universe] * 42
inner_ring_three = [inner_universe] * 36
inner_ring_four = [inner_universe] * 30
inner_ring_five = [inner_universe] * 24
inner_ring_six = [inner_universe] * 18
inner_ring_seven = [inner_universe] * 12
inner_ring_eight = [inner_universe] * 6
inner_ring_nine = [inner_universe]

inner_lattice.universes = [inner_ring_one, inner_ring_two, inner_ring_three,
                           inner_ring_four, inner_ring_five, inner_ring_six,
                           inner_ring_seven, inner_ring_eight, inner_ring_nine]

inner_prism = openmc.model.hexagonal_prism(edge_length=11.65092843, orientation='x')
wrapper_tube_inner = openmc.model.hexagonal_prism(edge_length=11.91073605,
                                                        orientation='x')
wrapper_tube_region_inner = (~inner_prism & wrapper_tube_inner & 
                             +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_inner = (~wrapper_tube_inner & -axial_reflector_max_z & 
                              +gas_plenum_min_z)

main_in_assembly_cell = openmc.Cell(fill=inner_lattice,
                                    region=inner_prism & -axial_reflector_max_z & 
                                    +gas_plenum_min_z)

wrapper_tube_inner_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_inner)
wrapper_tube_inner_outside_cell = openmc.Cell(fill=sodium,
                                              region=outside_wrapper_tube_inner)
inner_assembly_universe = openmc.Universe(cells=[main_in_assembly_cell, 
                                                 wrapper_tube_inner_cell,
                                                 wrapper_tube_inner_outside_cell])

outer_lattice = openmc.HexLattice()
outer_lattice.center = (0., 0.)
outer_lattice.pitch = (21.08/17,)
outer_lattice.orientation = 'x'
outer_lattice.outer = sodium_universe

out_ring_zero = [outer_universe] * 54
outer_ring_one = [outer_universe] * 48
outer_ring_two = [outer_universe] * 42
outer_ring_three = [outer_universe] * 36
outer_ring_four = [outer_universe] * 30
outer_ring_five = [outer_universe] * 24
outer_ring_six = [outer_universe] * 18
outer_ring_seven = [outer_universe] * 12
outer_ring_eight = [outer_universe] * 6
outer_ring_nine = [outer_universe]

outer_lattice.universes = [outer_ring_one, outer_ring_two, outer_ring_three,
                           outer_ring_four, outer_ring_five, outer_ring_six,
                           outer_ring_seven, outer_ring_eight, outer_ring_nine]

outer_prism = openmc.model.hexagonal_prism(edge_length=11.65092843, orientation='x')
wrapper_tube_outer = openmc.model.hexagonal_prism(edge_length=11.91073605,
                                                        orientation='x')
wrapper_tube_region_outer = (~outer_prism & wrapper_tube_outer & 
                             +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_outer = (~wrapper_tube_outer & -axial_reflector_max_z & 
                              +gas_plenum_min_z)

main_out_assembly_cell = openmc.Cell(fill=outer_lattice,
                                    region=outer_prism & -axial_reflector_max_z & 
                                    +gas_plenum_min_z)
wrapper_tube_outer_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_outer)
wrapper_tube_outer_outside_cell = openmc.Cell(fill=sodium,
                                              region=outside_wrapper_tube_outer)
outer_assembly_universe = openmc.Universe(cells=[main_out_assembly_cell,
                                                 wrapper_tube_outer_cell,
                                                 wrapper_tube_outer_outside_cell])

reflector_lattice = openmc.HexLattice()
reflector_lattice.center = (0., 0.)
reflector_lattice.pitch = (4.315,)
reflector_lattice.orientation = 'x'
reflector_lattice.outer = sodium_universe

reflector_inner_ring_one = [reflector_universe] * 12
reflector_inner_ring_two = [reflector_universe] * 6
reflector_inner_ring_three = [reflector_universe]

reflector_lattice.universes = [reflector_inner_ring_one, reflector_inner_ring_two, reflector_inner_ring_three]
reflector_prism = openmc.model.hexagonal_prism(edge_length=11.65092843, orientation='x')
reflector_wrapper =openmc.model.hexagonal_prism(edge_length=11.91073605, orientation='x')

wrapper_tube_region_ref = (~reflector_prism & reflector_wrapper & 
                             +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_ref = (~reflector_wrapper & -axial_reflector_max_z &
                              +gas_plenum_min_z)

wrapper_tube_ref_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_ref)
wrapper_tube_outside_ref_cell = openmc.Cell(fill=sodium,
                                              region=outside_wrapper_tube_ref)
main_reflector_assembly_cell = openmc.Cell(fill=reflector_lattice,
                                    region=reflector_prism & -axial_reflector_max_z &
                                    +gas_plenum_min_z)
reflector_assembly_universe = openmc.Universe(cells=[main_reflector_assembly_cell,
                                                     wrapper_tube_ref_cell,
                                                     wrapper_tube_outside_ref_cell])

control_lattice = openmc.HexLattice()
control_lattice.center = (0., 0.)
control_lattice.pitch = (2.43,)
control_lattice.orientation = 'x'
control_lattice.outer = sodium_universe
control_ring_one = [control_universe] * 18
control_ring_two = [control_universe] * 12
control_ring_three = [control_universe] * 6
control_ring_four = [control_universe]

control_lattice.universes = [control_ring_one, control_ring_two, control_ring_three, control_ring_four]
control_wrapper_1 = openmc.hexagonal_prism(edge_length=8.775724092, orientation='x')
control_wrapper_2 = openmc.hexagonal_prism(edge_length=9.006664199, orientation='x')
control_wrapper_3 = openmc.hexagonal_prism(edge_length=11.65092843, orientation='x')
control_wrapper_4 = openmc.hexagonal_prism(edge_length=11.91073605, orientation='x')

control_assembly_cell = openmc.Cell(fill=control_lattice, region=control_wrapper_1
                                    & +gas_plenum_min_z & -axial_reflector_max_z)
control_wrapper_1_cell = openmc.Cell(fill=em10_steel, region = ~control_wrapper_1
                                     & control_wrapper_2 & +gas_plenum_min_z
                                     & -axial_reflector_max_z)
control_wrapper_gap_cell = openmc.Cell(fill=sodium, region=~control_wrapper_2
                                       & control_wrapper_3 & +gas_plenum_min_z
                                       & -axial_reflector_max_z)
control_wrapper_2_cell = openmc.Cell(fill=em10_steel, region = ~control_wrapper_3
                                     & control_wrapper_4 & +gas_plenum_min_z
                                     & -axial_reflector_max_z)
control_outside_cell = openmc.Cell(fill=sodium, region=~control_wrapper_4 & +gas_plenum_min_z
                                   & -axial_reflector_max_z)

control_assembly_universe = openmc.Universe(cells=[control_assembly_cell,
                                                   control_wrapper_1_cell,
                                                   control_wrapper_gap_cell,
                                                   control_wrapper_2_cell,
                                                   control_outside_cell])

main_lattice = openmc.HexLattice()
main_lattice.center = (0., 0.)
main_lattice.pitch = (21.08,)
main_lattice.orientation = 'y'
main_lattice.outer = sodium_universe

reflector_ring_one = [reflector_assembly_universe] * 96
reflector_ring_two = [reflector_assembly_universe] * 90
reflector_ring_three = [reflector_assembly_universe] * 84
reflector_ring_four = ([reflector_assembly_universe] * 5 + [outer_assembly_universe] * 4 + 
                       [reflector_assembly_universe] * 4) * 6
reflector_ring_five = ([reflector_assembly_universe] + [outer_assembly_universe] * 11) * 6
out_ring_one = [outer_assembly_universe] * 66
out_ring_two = ([outer_assembly_universe] * 5 + [control_assembly_universe] + 
                ([outer_assembly_universe] * 9 + [control_assembly_universe]) * 5 + [outer_assembly_universe] * 4)
out_ring_three = ([outer_assembly_universe] * 2 + ([control_assembly_universe] + [inner_assembly_universe] * 4 +
                  [control_assembly_universe] + [outer_assembly_universe] * 3) * 5 +
                  [control_assembly_universe] + [inner_assembly_universe] * 4  
                  + [control_assembly_universe] + [outer_assembly_universe])
in_ring_one = [inner_assembly_universe] * 48
in_ring_two = [inner_assembly_universe] * 42
in_ring_three = ([inner_assembly_universe]*3 + [sodium_prism_universe]) * 9
in_ring_four = [inner_assembly_universe] * 30
in_ring_five = [inner_assembly_universe] * 24
in_ring_six = ([control_assembly_universe] + [inner_assembly_universe] * 2) * 6
in_ring_seven = [inner_assembly_universe] * 12
in_ring_eight = [inner_assembly_universe] * 6
in_ring_nine = [reflector_assembly_universe] * 1

main_lattice.universes = [reflector_ring_one, reflector_ring_two, reflector_ring_three,
                          reflector_ring_four, reflector_ring_five,
                          out_ring_one,out_ring_two,out_ring_three,in_ring_one,in_ring_two,
                          in_ring_three,in_ring_four,in_ring_five,in_ring_six,
                          in_ring_seven,in_ring_eight,in_ring_nine]

core_surface = openmc.model.hexagonal_prism(edge_length=347.82, boundary_type='vacuum')
core = openmc.Cell(fill=main_lattice,
                   region = core_surface & -axial_reflector_max_z & 
                   +gas_plenum_min_z)
outside_core = openmc.Cell(fill=sodium,
                           region = ~core_surface & -axial_reflector_max_z & 
                           +gas_plenum_min_z)
core_universe = openmc.Universe(cells=[core, outside_core])
geom = openmc.Geometry(inner_assembly_universe)
geom.export_to_xml()





point = openmc.stats.Point((0, 0, 0))
source = openmc.Source(space=point)
settings = openmc.Settings()
settings.source = source
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
settings.export_to_xml()
#openmc.run()

#depletion calculation

inner_fuel.volume = np.pi * ((0.943/2) ** 2) * 217 * 225
outer_fuel.volume = np.pi * ((0.943/2) ** 2) * 217 * 228

chain = openmc.deplete.Chain.from_xml("/home/mpf/support_files/chain_casl_sfr.xml")
model = openmc.Model(geometry=geom, settings=settings, materials=materials_file)
operator = openmc.deplete.CoupledOperator(model, "/home/mpf/support_files/chain_casl_sfr.xml", diff_burnable_mats=True)

power = 16000000

time_steps = [30] * 2
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps,
                                                power, timestep_units = 'd')
integrator.integrate()

Hi @marti. All you should need to do is to call:

geom.determine_paths()

after your geom variable has been created.

Hi Paul, thanks for the quick answer. I added geom.determine_paths() after I define geom and it still returns the same error.

@marti It appears that the problem is as follows:

  • When it tries to differentiate burnable materials, it’s finding that outer_fuel doesn’t have a number of instances determined (which should usually happen with geom.determine_paths())
  • Your Geometry object has a root universe of inner_assembly_universe
  • However, nowhere in inner_assembly_universe is the outer_fuel material, so the call to geom.determine_paths() doesn’t actually count the number of instances for that material.
  • You’ve explicitly given a list of materials to the Model class that includes outer_fuel, and this ends up confusing the diff_burnable_mats option.

To fix this problem, simply omit the list of materials when you’re creating a Model:

model = openmc.Model(geometry=geom, settings=settings)

Also, with that in place, you don’t need to explicitly call geom.determine_paths() as that’s already done under the hood – sorry for misleading you before.

Hi Paul, when I do that I get the error:
DataError: Cross sections were not specified in Model.materials and the OPENMC_CROSS_SECTIONS environment variable is not set.

This is probably because I haven’t actually set an OPENMC_CROSS_SECTIONS. I am using a MacBook and I’m not sure how to do that, so I downloaded the ENDFB71 cross sections file.

Hi again @marti. You can set environment variables either directly from a terminal (export OPENMC_CROSS_SECTIONS=...) or it may be more convenient to set them permanently in your terminal profile. See here for further instructions on that.