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()