Fix numerical cancellation in RectLattice::distance for large pitch values

include/openmc/math_functions.h

@@ -233,5 +233,18 @@ void get_energy_index(
 
 double standard_normal_cdf(double z);
 
+//==============================================================================
+//! Return true if two floating-point values are approximately equal within a
+//! combined relative and absolute tolerance.
+//!
+//! \param a first floating point value
+//! \param b second floating point value
+//! \param rel_tol relative tolerance
+//! \param abs_tol absolute tolerance
+//! \return true if a and b are approximately equal, false otherwise
+//==============================================================================
+
+bool isclose(double a, double b, double rel_tol, double abs_tol);
+
 } // namespace openmc
 #endif // OPENMC_MATH_FUNCTIONS_H

src/lattice.cpp

@@ -10,6 +10,7 @@
 #include "openmc/geometry.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/string_utils.h"
 #include "openmc/vector.h"
 #include "openmc/xml_interface.h"
@@ -262,27 +263,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
   double y0 {copysign(0.5 * pitch_[1], u.y)};
   double z0;
 
-  double d = std::min(
-    u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
+  // Evaluate axial distances.
+  double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
+  double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
+  double dz = is_3d_ ? (u.z != 0.0 ? (z0 - z) / u.z : INFTY) : INFTY;
+
   if (is_3d_) {
-    z0 = copysign(0.5 * pitch_[2], u.z);
-    d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
+    z0 = std::copysign(0.5 * pitch_[2], u.z);
+    dz = (u.z != 0.0) ? (z0 - z) / u.z : INFTY;
   }
 
-  // Determine which lattice boundaries are being crossed
+  // Minimum distance
+  double d = std::min({dx, dy, dz});
+
   array<int, 3> lattice_trans = {0, 0, 0};
-  if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
-    lattice_trans[0] = copysign(1, u.x);
-  if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
-    lattice_trans[1] = copysign(1, u.y);
-  if (is_3d_) {
-    if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
-      lattice_trans[2] = copysign(1, u.z);
-  }
+
+  // Determine which directions are crossed
+  if (isclose(d, dx, 1e-12, FP_PRECISION))
+    lattice_trans[0] = std::copysign(1, u.x);
+
+  if (isclose(d, dy, 1e-12, FP_PRECISION))
+    lattice_trans[1] = std::copysign(1, u.y);
+
+  if (is_3d_ && isclose(d, dz, 1e-12, FP_PRECISION))
+    lattice_trans[2] = std::copysign(1, u.z);
 
   return {d, lattice_trans};
 }
-
 //==============================================================================
 
 void RectLattice::get_indices(

src/math_functions.cpp

@@ -959,4 +959,11 @@ void get_energy_index(
   }
 }
 
+// Return true if two floating-point values are approximately equal within a
+// combined relative and absolute tolerance.
+bool isclose(double a, double b, double rel_tol, double abs_tol)
+{
+  return std::abs(a - b) <=
+         std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
+}
 } // namespace openmc

src/mesh.cpp

@@ -28,6 +28,7 @@
 #include "openmc/geometry.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/material.h"
+#include "openmc/math_functions.h"
 #include "openmc/memory.h"
 #include "openmc/message_passing.h"
 #include "openmc/openmp_interface.h"
@@ -1861,13 +1862,13 @@ double CylindricalMesh::find_r_crossing(
   // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0
 
   const double r0 = grid_[0][shell];
-  if (r0 == 0.0)
+  if (isclose(r0, 0.0, 0.0, FP_PRECISION))
     return INFTY;
 
   const double denominator = u.x * u.x + u.y * u.y;
 
   // Direction of flight is in z-direction. Will never intersect r.
-  if (std::abs(denominator) < FP_PRECISION)
+  if (isclose(denominator, 0.0, 0.0, FP_PRECISION))
     return INFTY;
 
   // inverse of dominator to help the compiler to speed things up
@@ -1883,7 +1884,7 @@ double CylindricalMesh::find_r_crossing(
   D = std::sqrt(D);
 
   // Particle is already on the shell surface; avoid spurious crossing
-  if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
+  if (isclose(R, r0, RADIAL_MESH_TOL, FP_PRECISION))
     return INFTY;
 
   // Check -p - D first because it is always smaller as -p + D
@@ -2157,14 +2158,14 @@ double SphericalMesh::find_r_crossing(
   // solve |r+s*u| = r0
   // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !)
   const double r0 = grid_[0][shell];
-  if (r0 == 0.0)
+  if (isclose(r0, 0.0, 0.0, FP_PRECISION))
     return INFTY;
   const double p = r.dot(u);
   double R = r.norm();
   double D = p * p - (R - r0) * (R + r0);
 
   // Particle is already on the shell surface; avoid spurious crossing
-  if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
+  if (isclose(R, r0, RADIAL_MESH_TOL, FP_PRECISION))
     return INFTY;
 
   if (D >= 0.0) {