Fix raytrace infinite loop. (#3423)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten 2025-06-12 01:54:04 +03:00 committed by GitHub
parent 2eeba89992
commit e678b1a057
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
2 changed files with 40 additions and 2 deletions

View file

@ -989,7 +989,7 @@ void StructuredMesh::raytrace_mesh(
// For all directions outside the mesh, find the distance that we need
// to travel to reach the next surface. Use the largest distance, as
// only this will cross all outer surfaces.
int k_max {0};
int k_max {-1};
for (int k = 0; k < n; ++k) {
if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
(distances[k].distance > traveled_distance)) {
@ -997,6 +997,10 @@ void StructuredMesh::raytrace_mesh(
k_max = k;
}
}
// Assure some distance is traveled
if (k_max == -1) {
traveled_distance += TINY_BIT;
}
// If r1 is not inside the mesh, exit here
if (traveled_distance >= total_distance)
@ -1011,7 +1015,7 @@ void StructuredMesh::raytrace_mesh(
}
// If inside the mesh, Tally inward current
if (in_mesh)
if (in_mesh && k_max >= 0)
tally.surface(ijk, k_max, !distances[k_max].max_surface, true);
}
}
@ -1207,6 +1211,7 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(
d.next_index--;
d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
}
return d;
}

View file

@ -618,3 +618,36 @@ def test_mesh_material_volumes_serialize():
assert new_volumes.by_element(1) == [(None, 1.0)]
assert new_volumes.by_element(2) == [(2, 0.5), (1, 0.5)]
assert new_volumes.by_element(3) == [(2, 1.0)]
def test_raytrace_mesh_infinite_loop():
# Create a model with one large spherical cell
sphere = openmc.Sphere(r=100, boundary_type='vacuum')
cell = openmc.Cell(region=-sphere)
model = openmc.Model()
model.geometry = openmc.Geometry([cell])
# Create a regular mesh and associated tally
mesh_surface = openmc.RegularMesh()
mesh_surface.lower_left = (-30, -30, 30)
mesh_surface.upper_right = (30, 30, 60)
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']
model.tallies = [mesh_surface_tally]
# Define a source such that the z position is on a mesh boundary with a very
# small directional cosine in the z direction
polar = openmc.stats.delta_function(1.75e-7)
azimuthal = openmc.stats.Uniform(0.0, 2.0*pi)
model.settings.source = openmc.IndependentSource(
angle=openmc.stats.PolarAzimuthal(polar, azimuthal)
)
model.settings.run_mode = 'fixed source'
model.settings.particles = 10
model.settings.batches = 1
# Run the model; this should not cause an infinite loop
model.run()