Fix lost particles after virtual surface crossings in complex regions

src/cell.cpp

@@ -999,6 +999,7 @@ std::pair<double, int32_t> Region::distance_to_nearest_surface(Position r,
 std::pair<double, int32_t> Region::distance_complex(
   Position r, Direction u, int32_t on_surface) const
 {
+  const Position r_initial {r};
   const bool in_region = contains_complex(r, u, on_surface);
   double total_distance {0.0};
 
@@ -1009,12 +1010,23 @@ std::pair<double, int32_t> Region::distance_complex(
       return {INFTY, std::numeric_limits<int32_t>::max()};
     }
 
-    // Move to the candidate surface and determine which side of it the ray is
-    // entering. The surface normal is used instead of evaluating the surface
-    // equation because accumulated roundoff may place the point slightly to
-    // the wrong side of a curved surface.
-    r += distance * u;
+    // Move to the candidate surface. The position is recomputed from the
+    // starting point with the same arithmetic as GeometryState::move_distance
+    // rather than accumulated step by step. When the true boundary is found,
+    // region membership has therefore been evaluated at the position where the
+    // particle will actually be transported. Positions at virtual crossings
+    // are used only within this search.
+    const double previous_distance = total_distance;
     total_distance += distance;
+    if (total_distance == previous_distance) {
+      // The step was lost to roundoff, so the position would not move. Take
+      // the smallest representable step so that the search always advances.
+      total_distance = std::nextafter(total_distance, INFTY);
+    }
+    r = r_initial + total_distance * u;
+
+    // Determine which side of the surface the ray is entering. Use the surface
+    // normal, because roundoff may place the point on the wrong side.
     i_surf = std::abs(i_surf);
     const auto& surf {*model::surfaces[i_surf - 1]};
     if (u.dot(surf.normal(r)) <= 0.0) {

tests/cpp_unit_tests/test_region.cpp

@@ -100,6 +100,65 @@ class CoincidentSurfaceFixture {
   }
 };
 
+// Helper class for testing accumulated distance through virtual crossings
+class AccumulatedDistanceFixture {
+public:
+  AccumulatedDistanceFixture()
+  {
+    pugi::xml_document doc;
+    auto plane = doc.append_child("surface");
+    plane.append_attribute("id") = 1;
+    plane.append_attribute("type") = "x-plane";
+    plane.append_attribute("coeffs") = 0.0;
+
+    const double offsets[] {-57452.33336021505, 43403.32187479845,
+      62125.04976135607, 62125.04976135607};
+    for (int i = 0; i < 4; ++i) {
+      plane.attribute("id") = i + 1;
+      plane.attribute("coeffs") = offsets[i];
+      openmc::model::surfaces.push_back(
+        std::make_unique<openmc::SurfaceXPlane>(plane));
+      openmc::model::surface_map[i + 1] = i;
+    }
+  }
+
+  ~AccumulatedDistanceFixture()
+  {
+    openmc::model::surfaces.clear();
+    openmc::model::surface_map.clear();
+  }
+};
+
+// Helper class for testing coincident planes within a single region
+class CoincidentPlaneFixture {
+public:
+  CoincidentPlaneFixture()
+  {
+    pugi::xml_document doc;
+    auto plane = doc.append_child("surface");
+    plane.append_attribute("id") = 1;
+    plane.append_attribute("type") = "x-plane";
+    plane.append_attribute("coeffs") = 0.0;
+
+    // Surfaces 1 and 2 are the same plane with different IDs. Surface 3 is far
+    // away and is never reached.
+    const double offsets[] {2718.775540449465, 2718.775540449465, 1.0e9};
+    for (int i = 0; i < 3; ++i) {
+      plane.attribute("id") = i + 1;
+      plane.attribute("coeffs") = offsets[i];
+      openmc::model::surfaces.push_back(
+        std::make_unique<openmc::SurfaceXPlane>(plane));
+      openmc::model::surface_map[i + 1] = i;
+    }
+  }
+
+  ~CoincidentPlaneFixture()
+  {
+    openmc::model::surfaces.clear();
+    openmc::model::surface_map.clear();
+  }
+};
+
 } // anonymous namespace
 
 TEST_CASE("Test region simplification")
@@ -223,3 +282,60 @@ TEST_CASE("Ignore roundoff-scale virtual surface crossings")
   REQUIRE(distance == Catch::Approx(603.9161175466262));
   REQUIRE(surface == 1);
 }
+
+TEST_CASE("Maintain position consistency through virtual surface crossings")
+{
+  AccumulatedDistanceFixture fixture;
+  openmc::Region source("-1 | -2 | -3", 0);
+  openmc::Region destination("4", 0);
+
+  // Planes 1 and 2 are virtual crossings and plane 3 is the true boundary.
+  // Surface 4 is the same plane as seen from the destination cell.
+  //
+  // The coordinates are deliberate: a position accumulated step by step lands
+  // short of surface 3, and rounder values do not show the difference.
+  openmc::Position r {-85626.45049221347, 0.0, 0.0};
+  openmc::Direction u {0.9182609440159194, 0.39597580569397484, 0.0};
+
+  auto [distance, surface] = source.distance(r, u, 0);
+  r += distance * u;
+
+  REQUIRE(distance == Catch::Approx(160903.60938949839));
+  REQUIRE(surface == 3);
+  REQUIRE(destination.contains(r, u, surface));
+}
+
+TEST_CASE("Find boundary at coincident surfaces within one region")
+{
+  CoincidentPlaneFixture fixture;
+
+  // After about 166 m the first candidate crossing lands short of both planes,
+  // and the remaining step is lost to roundoff. The search must still return a
+  // distance that carries the particle across both planes.
+  const openmc::Position r_start {-9948.187586907628, 0.0, 0.0};
+  const openmc::Direction u {0.7623218202875262, 0.6471981476437588, 0.0};
+
+  SECTION("Leaving a union")
+  {
+    openmc::Region region("-1 | -2", 0);
+    auto [distance, surface] = region.distance(r_start, u, 0);
+    openmc::Position r = r_start;
+    r += distance * u;
+
+    REQUIRE(distance == Catch::Approx(16616.2935262426));
+    REQUIRE((surface == 1 || surface == 2));
+    REQUIRE_FALSE(region.contains(r, u, surface));
+  }
+
+  SECTION("Entering an intersection")
+  {
+    openmc::Region region("(1 2) | 3", 0);
+    auto [distance, surface] = region.distance(r_start, u, 0);
+    openmc::Position r = r_start;
+    r += distance * u;
+
+    REQUIRE(distance == Catch::Approx(16616.2935262426));
+    REQUIRE((surface == 1 || surface == 2));
+    REQUIRE(region.contains(r, u, surface));
+  }
+}

tests/regression_tests/filter_mesh/results_true.dat

@@ -1 +1 @@
-e54a31cf702974deda984a8a7d699ed0291b16a222f88c2b91d6295e4973f972c9fce7f90bfd2c18faea14f59247237f29ded153a1b0ed96f28afaf32174cb3c
+f3ea866c98528e5a06df31d664bfbbeb8ba2cceb0c813352d880c6038dcc4913f26af7b5c8ab4d12983964c30c8b0c1bf46d0e8a48fcf923940bdd337686e7a4