Hi. I am trying to simulate a lead cooled fast reactor based on the core design of ALFRED (European LFR Demonstrator). The issue is that each time I run the code, I get the following error.
I have checked out the troubleshooting section on the OpenMC documentation but still cannot identify where the problem is in the geometry that is causing the particles to be lost during the simulation.
I’m using Jupyter Notebook to run the code.
Welcome to the community, thanks for posting
This error message normally mean that the particles are getting to a region of undefined or incorrect geometry.
The error message gives you a hint that surface 2 is where the particles are getting to before finding undefined geometry.
Feel free to post the geometry here and we can take a look.
As a new user, I am unable to upload attachments.
Here is my Google Drive link to the folder containing the .ipynb and .xml files.
Please have a look.
https://drive.google.com/drive/folders/1y1q-gex_QHRwm4rHxs2sX_bZhP-BgdBO?usp=share_link
Here is the example code extracted from the notebook and pasted here.
import openmc
# Gas gap composition
helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 1,'wo')
helium.temperature = 700
# Clad composition
t91 = openmc.Material(name='T91 Steel')
t91.set_density('g/cm3', 7.7)
t91.add_element('Fe', 0.884, 'wo')
t91.add_element('Cr', 0.09, 'wo')
t91.add_element('Mo', 0.01 , 'wo')
t91.add_element('Mn', 0.006, 'wo')
t91.add_element('Si', 0.005, 'wo')
t91.add_element('V', 0.002, 'wo')
t91.add_element('Ni', 0.002, 'wo')
t91.add_element('Nb', 0.001, 'wo')
t91.temperature = 700
# Coolant composition
lead = openmc.Material(name='Lead Coolant')
lead.add_element('Pb', 1)
lead.set_density('g/cm3', 10.58)
lead.temperature = 673
# Absorber material
b4c = openmc.Material(name='B4C')
b4c.add_nuclide('B10', 0.9*4)
b4c.add_nuclide('B11', 0.1*4)
b4c.add_nuclide('C12', 1)
b4c.set_density('g/cm3', 2.2)
b4c.temperature = 700
# Fuel composition for inner and outer zones
u235 = openmc.Material(name='U235')
u235.add_nuclide('U235', 1.0)
u235.set_density('g/cm3', 10.0)
u238 = openmc.Material(name='U238')
u238.add_nuclide('U238', 1.0)
u238.set_density('g/cm3', 10.0)
pu238 = openmc.Material(name='Pu238')
pu238.add_nuclide('Pu238', 1.0)
pu238.set_density('g/cm3', 10.0)
pu239 = openmc.Material(name='U235')
pu239.add_nuclide('Pu239', 1.0)
pu239.set_density('g/cm3', 10.0)
pu240 = openmc.Material(name='Pu240')
pu240.add_nuclide('Pu240', 1.0)
pu240.set_density('g/cm3', 10.0)
pu241 = openmc.Material(name='Pu241')
pu241.add_nuclide('Pu241', 1.0)
pu241.set_density('g/cm3', 10.0)
pu242 = openmc.Material(name='Pu242')
pu242.add_nuclide('Pu242', 1.0)
pu242.set_density('g/cm3', 10.0)
am241 = openmc.Material(name='Am241')
am241.add_nuclide('Am241', 1.0)
am241.set_density('g/cm3', 10.0)
o16 = openmc.Material(name='O16')
o16.add_nuclide('O16', 1.0)
o16.set_density('g/cm3', 10.0)
inner = openmc.Material.mix_materials(
[u235, u238, pu238, pu239, pu240, pu241, pu242, am241, o16],
[0.0019, 0.7509, 0.0046, 0.0612, 0.0383, 0.0106, 0.0134, 0.001, 0.1181],
'wo')
inner.temperature = 700 # This is equivalent to 427°C.
outer = openmc.Material.mix_materials(
[u235, u238, pu238, pu239, pu240, pu241, pu242, am241, o16],
[0.0018, 0.73, 0.0053, 0.0711, 0.0445, 0.0124, 0.0156, 0.0017, 0.1176],
'wo')
inner.temperature = 700
# Instantiate a materials collection and export it to XML
materials_file = openmc.Materials([b4c, helium, t91, lead, inner, outer])
materials_file.export_to_xml()
# Instantiate cylindrical surfaces fuel
# All lengths are in cm
# The ALFRED core active height is 90 cm
fuel_ir = openmc.ZCylinder(r=0.10, name='Fuel IR')
fuel_or = openmc.ZCylinder(r=0.45, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.465, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.525, name='Clad OR')
top = openmc.ZPlane(z0=45, name='top', boundary_type = 'reflective')
bottom = openmc.ZPlane(z0=-45, name='bottom', boundary_type = 'reflective')
# Instantiate cells and fill with materials
hole = openmc.Cell(name='Hole', region = -fuel_ir, fill = helium)
fuel_iz = openmc.Cell(name='Inner Zone Fuel', region = +fuel_ir & -fuel_or, fill = inner)
fuel_oz = openmc.Cell(name='Outer Zone Fuel', region = +fuel_ir & -fuel_or , fill = outer)
gap = openmc.Cell(name='Gap', region = +fuel_or & -clad_ir, fill = helium)
clad = openmc.Cell(name='Cladding', region = +clad_ir & -clad_or, fill = t91)
coolant = openmc.Cell(name='Coolant', region = +clad_or, fill = lead)
# Instantiate inner zone fuel universe
f_iz = openmc.Universe(cells=[hole, fuel_iz, gap, clad, coolant])
f_iz.plot(width=(2.0,2.0), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255), inner:(227,207,87),
outer:(205,51,51)})
# Instantiate outer zone fuel universe
f_oz = openmc.Universe(cells=[hole, fuel_oz, gap, clad, coolant])
f_oz.plot(width=(2.0,2.0), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255), inner:(227,207,87),
outer:(255,51,51)})
# Instantiate coolant universe
lead_cell = openmc.Cell(fill=lead)
lead_univ = openmc.Universe(cells=[lead_cell])
lead_univ.plot(width=(2.0,2.0), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255), inner:(227,207,87),
outer:(255,51,51)})
# Fuel assembly - inner zone
fa_iz = openmc.HexLattice()
fa_iz.center = (0., 0.)
fa_iz.pitch = (1.384,)
fa_iz.outer = lead_univ
fa_iz.orientation = 'x'
outer_ring = [f_iz]*36
ring_1 = [f_iz]*30
ring_2 = [f_iz]*24
ring_3 = [f_iz]*18
ring_4 = [f_iz]*12
ring_5 = [f_iz]*6
inner_ring = [f_iz]
fa_iz.universes = [outer_ring,
ring_1,
ring_2,
ring_3,
ring_4,
ring_5,
inner_ring]
print(fa_iz)
# Create the prism that will contain the lattice
outer_in_surface = openmc.model.hexagonal_prism(edge_length=7*1.384, orientation='x')
# Fill a cell with the lattice. This cell is filled with the lattice and contained within the prism.
main_in_assembly = openmc.Cell(fill=fa_iz, region=outer_in_surface & -top & +bottom)
# Fill a cell with a material that will surround the lattice
out_in_assembly = openmc.Cell(fill=lead, region=~outer_in_surface & -top & +bottom)
# Create a universe that contains both
main_in_u = openmc.Universe(cells=[main_in_assembly, out_in_assembly])
# Plot the universe
main_in_u.plot(origin = (0,0,0), pixels=(500, 500), width = (30.,30.), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255), inner:(227,207,87),
outer:(255,51,51)})
# Fuel assembly - outer zone
fa_oz = openmc.HexLattice()
fa_oz.center = (0., 0.)
fa_oz.pitch = (1.384,)
fa_oz.outer = lead_univ
fa_oz.orientation = 'x'
outer_ring = [f_oz]*36
ring_1 = [f_oz]*30
ring_2 = [f_oz]*24
ring_3 = [f_oz]*18
ring_4 = [f_oz]*12
ring_5 = [f_oz]*6
inner_ring = [f_oz]
fa_oz.universes = [outer_ring,
ring_1,
ring_2,
ring_3,
ring_4,
ring_5,
inner_ring]
print(fa_oz)
# Create the prism that will contain the lattice
outer_out_surface = openmc.model.hexagonal_prism(edge_length=7*1.384, orientation='x')
# Fill a cell with the lattice. This cell is filled with the lattice and contained within the prism.
main_out_assembly = openmc.Cell(fill=fa_oz, region=outer_out_surface)
# Fill a cell with a material that will surround the lattice
out_out_assembly = openmc.Cell(fill=lead, region=~outer_out_surface)
# Create a universe that contains both
main_out_u = openmc.Universe(cells=[main_out_assembly, out_out_assembly])
# Plot the universe
main_out_u.plot(origin = (0,0,0), pixels=(500, 500), width = (30.,30.), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255), inner:(227,207,87),
outer:(255,51,51)})
# Create a hexagonal lead cell
lead_assembly = openmc.model.hexagonal_prism(edge_length=7*1.384, orientation='x')
la_cell = openmc.Cell(fill=lead, region=lead_assembly)
out_la_cell = openmc.Cell(fill=lead, region=~lead_assembly)
lead_u = openmc.Universe(cells=[la_cell, out_la_cell])
lead_u.plot(origin = (0,0,0), pixels=(500, 500), width = (30.,30.), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255),
inner:(227,207,87),outer:(255,51,51)})
print(main_in_u, main_out_u, lead_u)
# We have 3 types of assemblies created.
# We can now make the entire reactor core by creating a lattice that is filled with the assemblies.
# Define the core lattice
core_lat = openmc.HexLattice(name='core')
core_lat.center = (0., 0.)
core_lat.pitch = (12*1.384,)
core_lat.outer = lead_univ
# Create rings of fuel universes that will fill the lattice
outer_ring = [lead_u, lead_u, main_out_u, main_out_u, main_out_u, main_out_u, lead_u]*6
ring_1 = [main_out_u, main_out_u, main_out_u, lead_u, main_out_u, main_out_u]*6
ring_2 = [lead_u, main_out_u, main_out_u, main_out_u, main_out_u]*6
ring_3 = [main_in_u]*24
ring_4 = [main_in_u]*18
ring_5 = [lead_u, main_in_u, main_in_u]*4
ring_6 = [main_in_u]*6
inner_ring = [lead_u]
core_lat.universes = [outer_ring,
ring_1,
ring_2,
ring_3,
ring_4,
ring_5,
ring_6,
inner_ring]
print(core_lat)
# Create the prism that will contain the lattice
outer_core_surface = openmc.model.hexagonal_prism(edge_length=180, boundary_type='reflective')
# Fill a cell with the lattice. This cell is filled with the lattice and contained within the prism.
core = openmc.Cell(fill=core_lat, region=outer_core_surface & -top & +bottom)
# Create a universe that contains the lattice
main_u = openmc.Universe(cells=[core])
main_u.plot(origin = (0,0,0), pixels=(1000, 1000), width = (400.,400.), color_by = 'material',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255),
inner:(227,207,87),outer:(255,51,51)})
main_u.plot(origin = (0,0,0), pixels=(1000, 1000), width = (400.,100.), color_by = 'material', basis='yz',
colors = {b4c:(50,150,0), t91:(25,25,112), helium:(255,255,255), lead:(0,100,255),
inner:(227,207,87),outer:(255,51,51)})
geometry = openmc.Geometry(main_u)
geometry.export_to_xml()
# OpenMC simulation parameters
batches = 5
inactive = 1
particles = 10000
##############################################################################
# Exporting to OpenMC settings.xml file #
##############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.temperature = {'method':'interpolation'}
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.energy_mode = 'continuous-energy'
settings_file.generations_per_batch = 1
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-10, -10, -10, 10, 10, 10]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.Source(space=uniform_dist)
settings_file.export_to_xml()
openmc.run()
error message
WARNING: After particle 7502 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 8752 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 2502 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 6252 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 5002 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 3752 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 1252 crossed surface 2 it could not be located in any
cell and it did not leak.
WARNING: After particle 2 crossed surface 2 it could not be located in any cell
and it did not leak.
WARNING: After particle 7503 crossed surface 2 it could not be located in any
cell and it did not leak.
...
RuntimeError: Maximum number of lost particles has been reached. ERROR: Maximum number of lost particles has been reached. ERROR: Maximum number of lost particles has been reached.
I think the option is to identify which surface is leaking, I tend to do this by setting the surface id values so that they are not automatically assigned. However I don’t know how to do this with latices.
It might also help to reduce the simulation script to a minimal reproducible example as the current script is quite large.
Additionally we could think about improving the error message in openmc so that leaking particles output more information. Perhaps the coordinates, material and cell that the particle is traveling from could be included.
Hi there,
It looks like the issue with your simulation might be caused by adding the same cells to multiple universes. Your code indicates that cells hole, gap, clad, and coolant are being added to both f_iz and f_oz universes. However, cells can only belong to a single universe, meaning that by adding these cells to the second universe f_oz, they are no longer associated with the first universe f_iz. This may explain the lost particles you are encountering in your simulation.
I hope this helps clarify the issue and resolve your problem. Let me know if you have any further questions!
Thanks a lot. The problem has been solved. It appears you’re right about the fact that the same cell, when added to more than one universes, results in OpenMC not being able to identify/locate particles crossing the cell’s surfaces. I re-instantiated the cells before adding them to f_oz and the k-effective eigenvalue simulation ran wonderfully without any issues.
Thanks again. God bless!
I met the same problem one day ago and now is solved.And i 'm using OpenMC to simulate a sodium cooled fast reactor(EBR-II)，we have quite a lot similar things,maybe we can help mutual in a certain.
Sure, feel free to inbox me.