marquezj
commented
6 months ago I looked into this bug and it seems it is specific to RegularMesh(). I simplified the working example, and it looks that for very shallow angles of the source, combined with specific values of the grid in the perpendicular direction, it enters a loop:
import openmc
import numpy as np
outside_Sphere = openmc.Sphere(r=100, boundary_type='vacuum')
outside_cell = openmc.Cell(fill = None, region=-outside_Sphere)
uni = openmc.Universe(cells=[outside_cell])
geometry = openmc.Geometry(root=uni)
mesh_surface = openmc.RegularMesh() # <-- Important
mesh_surface.lower_left = [-30, -30, 30] # <-- Min z is important
mesh_surface.upper_right = [30, 30, 60] # <-- Max z is important
mesh_surface.dimension = [1, 1, 1]
reg_filter = openmc.MeshSurfaceFilter(mesh_surface)
mesh_surface_tally = openmc.Tally()
mesh_surface_tally.filters = [reg_filter]
mesh_surface_tally.scores = ['current']
tallies = openmc.Tallies([mesh_surface_tally])
settings = openmc.Settings()
source = openmc.IndependentSource()
polar = openmc.stats.Discrete([1.75e-07], [1]) # <--- Important, polar angle < aprox 1.75e-7
azimuthal = openmc.stats.Uniform(0.0, 2.0*np.pi)
source.angle = openmc.stats.PolarAzimuthal(polar, azimuthal)
settings.source = source
settings.run_mode = 'fixed source'
settings.particles = 10
settings.batches = 1
model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies)
sp_filename = model.run(threads=1)To be clear: this seems to be a bug of RegularMesh(). RectilinearMesh() works perfectly fine with the same grid.
Bug Description
With the following setup I run into an endless loop. A simple mono-energetic Source of quadratic-plane shape and a little farther away a rectangular RegularMesh. On this mesh the current is tallied with the openmc.MeshSurfaceFilter. fill = None, so no materials involved.
I'm happy to help finding the bug.
Steps to Reproduce
The following code reproduces the error. Note: with
settings.seed= 21it happens in batch 3. If seed is altered the occurrence time gets delayed (resp shifted)`import openmc import numpy as np
Materials
materials = openmc.Materials([]) materials.export_to_xml() geometry = openmc.Geometry([]) Big_sphere_radius = 80 # cm - for outside of the simulation
Definition of surfaces
outside_Sphere = openmc.Sphere(r=Big_sphere_radius, boundary_type='vacuum')
outside_cell = openmc.Cell(cell_id=999, fill = None, region=-outside_Sphere
& ~( rotated_detector_Plate)
uni = openmc.Universe(cells=[outside_cell]) geometry = openmc.Geometry(root=uni) geometry.export_to_xml()
Place Source
source_position = 0 offset_x = 0 offset_y = 0
x_source = openmc.stats.Uniform(-2+offset_x,2+offset_x) # source with no thickness y_source = openmc.stats.Uniform(-2+offset_y,2+offset_y) z_source = openmc.stats.Uniform(source_position, source_position)
source = openmc.IndependentSource() source.energy = openmc.stats.Discrete([3e6], [1.0]) # 3MeV source.space=openmc.stats.CartesianIndependent(x=x_source, y=y_source, z=z_source)
MESH SURFACE TALLY
mesh_surface = openmc.RegularMesh(mesh_id = 3) mesh_surface.lower_left = [-20, -20, 56] mesh_surface.upper_right = [20, 20, 60] mesh_surface.dimension = [10, 10, 1] reg_filter = openmc.MeshSurfaceFilter(mesh_surface, filter_id=4) mesh_surface_tally = openmc.Tally(tally_id=5, name='mesh_surface') mesh_surface_tally.filters = [reg_filter] mesh_surface_tally.scores = ['current']
tallies = openmc.Tallies([mesh_surface_tally])
Settings and Model
settings = openmc.Settings() settings.source = source settings.run_mode = 'fixed source' settings.particles = 100000 settings.batches = 50 settings.verbosity= 6 settings.seed= 21
model = openmc.Model() model.geometry = geometry model.settings = settings model.tallies = tallies
Run Simulation
sp_filename = model.run(output=True)`
Environment
Ubuntu 22.04 with OpenMC (0.14.1) (but happened also at the previous release)