Fix virtual surface crossing

include/openmc/cell.h

@@ -112,11 +112,20 @@ class Region {
 
   //! Determine if a particle is inside the cell for a complex cell.
   //!
-  //! Uses the comobination of half-spaces and binary operators to determine
+  //! Uses the combination of half-spaces and binary operators to determine
   //! if short circuiting can be used. Short cicuiting uses the relative and
   //! absolute depth of parentheses in the expression.
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
 
+  //! Find the oncoming boundary of this cell for a simple cell (only
+  //! intersection operators)
+  std::pair<double, int32_t> distance_simple(
+    Position r, Direction u, int32_t on_surface) const;
+
+  //! Find the oncoming boundary of this cell for a complex cell.
+  std::pair<double, int32_t> distance_complex(
+    Position r, Direction u, int32_t on_surface) const;
+
   //! BoundingBox if the paritcle is in a simple cell.
   BoundingBox bounding_box_simple() const;
 

src/cell.cpp

@@ -919,6 +919,18 @@ std::string Region::str() const
 
 std::pair<double, int32_t> Region::distance(
   Position r, Direction u, int32_t on_surface) const
+{
+  if (simple_) {
+    return distance_simple(r, u, on_surface);
+  } else {
+    return distance_complex(r, u, on_surface);
+  }
+}
+
+//==============================================================================
+
+std::pair<double, int32_t> Region::distance_simple(
+  Position r, Direction u, int32_t on_surface) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
@@ -947,6 +959,37 @@ std::pair<double, int32_t> Region::distance(
 
 //==============================================================================
 
+std::pair<double, int32_t> Region::distance_complex(
+  Position r, Direction u, int32_t on_surface) const
+{
+  double min_dist {INFTY};
+  int32_t i_surf {std::numeric_limits<int32_t>::max()};
+
+  for (int32_t token : expression_) {
+    // Ignore this token if it corresponds to an operator rather than a region.
+    if (token >= OP_UNION)
+      continue;
+
+    // Calculate the distance to this surface.
+    // Note the off-by-one indexing
+    bool coincident {std::abs(token) == std::abs(on_surface)};
+    double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+
+    // Check if this distance is the new minimum.
+    if (d < min_dist) {
+      if (min_dist - d >= FP_PRECISION * min_dist) {
+        if (!contains_complex(r + (d + TINY_BIT) * u, u, on_surface)) {
+          min_dist = d;
+          i_surf = -token;
+        }
+      }
+    }
+  }
+  return {min_dist, i_surf};
+}
+
+//==============================================================================
+
 bool Region::contains(Position r, Direction u, int32_t on_surface) const
 {
   if (simple_) {

tests/regression_tests/collision_track/case_1_Reactions/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/collision_track/case_8_2threads/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-5.642735E-02 1.494035E-02
+6.085380E-02 2.048137E-03

tests/regression_tests/complex_cell/results_true.dat

@@ -1,11 +1,11 @@
 k-combined:
-2.603220E-01 1.429366E-03
+2.589444E-01 9.409859E-04
 tally 1:
-2.624819E+00
-1.378200E+00
-2.730035E+00
-1.492361E+00
-1.013707E+00
-2.055807E-01
-1.123257E-01
-2.530233E-03
+2.578878E+00
+1.331090E+00
+2.702603E+00
+1.462070E+00
+9.847537E-01
+1.942139E-01
+1.084095E-01
+2.358175E-03

tests/regression_tests/filter_cellfrom/results_true.dat

@@ -1,53 +1,53 @@
 k-combined:
-9.035025E-02 2.654309E-03
+9.946255E-02 2.184346E-03
 tally 1:
-5.994069E+00
-3.594398E+00
+6.089246E+00
+3.712059E+00
 tally 2:
-4.559707E-04
-2.080739E-08
+4.646490E-04
+2.160545E-08
 tally 3:
-5.994525E+00
-3.594945E+00
+6.089710E+00
+3.712624E+00
 tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
-1.473892E+00
-2.187509E-01
+1.556619E+00
+2.452581E-01
 tally 6:
-5.223875E-05
-2.992861E-10
+0.000000E+00
+0.000000E+00
 tally 7:
-1.473945E+00
-2.187663E-01
+1.556619E+00
+2.452581E-01
 tally 8:
-1.885798E+01
-3.558423E+01
+1.001104E+02
+1.003264E+03
 tally 9:
-7.467961E+00
-5.580255E+00
+7.645865E+00
+5.850005E+00
 tally 10:
-5.082094E-04
-2.584432E-08
+4.646490E-04
+2.160545E-08
 tally 11:
-7.468470E+00
-5.581014E+00
+7.646330E+00
+5.850716E+00
 tally 12:
-1.885798E+01
-3.558423E+01
+1.001104E+02
+1.003264E+03
 tally 13:
 0.000000E+00
 0.000000E+00
 tally 14:
-2.739543E-04
-7.600983E-09
+3.197805E-04
+1.030696E-08
 tally 15:
-2.739543E-04
-7.600983E-09
+3.197805E-04
+1.030696E-08
 tally 16:
-7.881296E+01
-6.221087E+02
+0.000000E+00
+0.000000E+00
 tally 17:
-1.051397E+02
-1.106292E+03
+1.077571E+02
+1.162197E+03

tests/regression_tests/filter_mesh/results_true.dat

@@ -1 +1 @@
-e07ed2bc8893c69721abf61b123336f1f6128a3bee6ec63b84d1f549f31707a74a6ce885091ccc0eac6b7f16f7cab39ede4784584c08825829e108de878ea5fb
+62452fff4ffd89deba11463d67bc60a72aa2f6d65f773b113e9ed34b3c4b4ecc7c189401c9752d0f7a60c57bd62ba1a3447890f39e7217e2de5957b2ea16dd65