Fix raytrace infinite loop. (#3423)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: GuySten

src/mesh.cpp

@@ -989,14 +989,18 @@ 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)) {
           traveled_distance = distances[k].distance;
           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;
 }
 

tests/unit_tests/test_mesh.py

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