Fix numerical cancellation in RectLattice::distance for large pitch values (#3853)

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

include/openmc/math_functions.h

@@ -233,5 +233,17 @@ 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"
@@ -260,25 +261,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
   // Determine the oncoming edge.
   double x0 {copysign(0.5 * pitch_[0], u.x)};
   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 the distance to each oncoming edge independently. Comparing these
+  // distances directly (rather than reconstructing the crossing position)
+  // avoids the floating-point cancellation that occurs for large pitches.
+  double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
+  double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
+  double dz = 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);
+    double z0 {copysign(0.5 * pitch_[2], u.z)};
+    dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY;
   }
 
-  // Determine which lattice boundaries are being crossed
+  // The distance to the nearest lattice boundary is the smallest axial
+  // distance.
+  double d = std::min({dx, dy, dz});
+
+  // Determine which lattice boundaries are being crossed. The axis attaining
+  // the minimum is exactly equal to d, so at least one translation is always
+  // set for a finite crossing; a near-equal second axis indicates a corner
+  // crossing.
   array<int, 3> lattice_trans = {0, 0, 0};
-  if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
+  if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION))
     lattice_trans[0] = copysign(1, u.x);
-  if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
+  if (isclose(d, dy, FP_COINCIDENT, 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);
-  }
+  if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION))
+    lattice_trans[2] = copysign(1, u.z);
 
   return {d, lattice_trans};
 }

src/math_functions.cpp

@@ -959,4 +959,12 @@ 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

tests/cpp_unit_tests/test_math.cpp

@@ -355,3 +355,24 @@ TEST_CASE("Test broaden_wmp_polynomials")
     REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val));
   }
 }
+
+TEST_CASE("Test isclose")
+{
+  using openmc::isclose;
+
+  // Identical values are always close, regardless of tolerances.
+  REQUIRE(isclose(1.0, 1.0, 0.0, 0.0));
+  REQUIRE(isclose(0.0, 0.0, 0.0, 0.0));
+
+  // Absolute tolerance governs comparisons near zero.
+  REQUIRE(isclose(0.0, 1e-15, 0.0, 1e-14));
+  REQUIRE_FALSE(isclose(0.0, 1e-13, 0.0, 1e-14));
+
+  // Relative tolerance scales with the magnitude of the operands.
+  REQUIRE(isclose(1.0e12, 1.0e12 + 1.0, 1e-12, 0.0));
+  REQUIRE_FALSE(isclose(1.0e12, 1.0e12 + 10.0, 1e-12, 0.0));
+
+  // The looser of the two tolerances wins.
+  REQUIRE(isclose(1.0, 1.0 + 1e-13, 0.0, 1e-12));
+  REQUIRE(isclose(1.0e6, 1.0e6 + 1e-4, 1e-9, 0.0));
+}

tests/regression_tests/lattice_large_pitch/inputs_true.dat

@@ -0,0 +1,51 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density value="0.001" units="g/cm3"/>
+      <nuclide name="N14" ao="1.0"/>
+    </material>
+    <material id="2">
+      <density value="7.0" units="g/cm3"/>
+      <nuclide name="Fe56" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="2" universe="1"/>
+    <cell id="2" material="1" universe="2"/>
+    <cell id="3" fill="3" region="1 -2 3 -4" universe="4"/>
+    <lattice id="3">
+      <pitch>100.0 100.0</pitch>
+      <outer>2</outer>
+      <dimension>3 3</dimension>
+      <lower_left>-150.0 -150.0</lower_left>
+      <universes>
+1 2 1 
+2 1 2 
+1 2 1 </universes>
+    </lattice>
+    <surface id="1" name="minimum x" type="x-plane" boundary="vacuum" coeffs="-150.0"/>
+    <surface id="2" name="maximum x" type="x-plane" boundary="vacuum" coeffs="150.0"/>
+    <surface id="3" name="minimum y" type="y-plane" boundary="vacuum" coeffs="-150.0"/>
+    <surface id="4" name="maximum y" type="y-plane" boundary="vacuum" coeffs="150.0"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>6 6</dimension>
+      <lower_left>-150.0 -150.0</lower_left>
+      <upper_right>150.0 150.0</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1</filters>
+      <scores>flux</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/lattice_large_pitch/results_true.dat

@@ -0,0 +1,73 @@
+tally 1:
+1.749265E+01
+3.275925E+01
+4.159430E+01
+1.799568E+02
+7.027454E+01
+4.967568E+02
+6.789716E+01
+4.652258E+02
+4.272935E+01
+1.861129E+02
+2.073342E+01
+4.374279E+01
+4.028297E+01
+1.641642E+02
+6.734263E+01
+4.627171E+02
+1.037657E+02
+1.078682E+03
+1.108394E+02
+1.240927E+03
+7.236770E+01
+5.302034E+02
+4.356361E+01
+1.921990E+02
+6.731536E+01
+4.615977E+02
+9.850684E+01
+9.810272E+02
+5.164863E+02
+2.670336E+04
+5.097770E+02
+2.604770E+04
+1.064847E+02
+1.137536E+03
+7.052986E+01
+5.003220E+02
+7.065046E+01
+5.036723E+02
+1.044890E+02
+1.103250E+03
+5.121544E+02
+2.632389E+04
+5.065331E+02
+2.570752E+04
+1.044794E+02
+1.101249E+03
+6.569142E+01
+4.361120E+02
+4.739812E+01
+2.313289E+02
+7.402284E+01
+5.591689E+02
+1.066905E+02
+1.149521E+03
+1.038665E+02
+1.086813E+03
+7.160008E+01
+5.217744E+02
+4.219705E+01
+1.816100E+02
+2.089838E+01
+4.464899E+01
+4.154271E+01
+1.745866E+02
+7.081253E+01
+5.049997E+02
+6.673273E+01
+4.509038E+02
+4.184624E+01
+1.771546E+02
+1.853898E+01
+3.538608E+01

tests/regression_tests/lattice_large_pitch/test.py

@@ -0,0 +1,70 @@
+"""Regression test for rectangular lattices with large pitch values.
+
+Large pitches used to trigger a segmentation fault in ``RectLattice::distance``
+because the boundary-crossing check compared a reconstructed crossing position
+against an absolute tolerance that did not scale with the geometry (see #3852).
+This test transports particles across a lattice with a large pitch to ensure the
+crossing logic remains robust.
+
+"""
+
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+    model = openmc.Model()
+
+    # Large pitch that previously caused floating-point cancellation
+    pitch = 100.0
+    n = 3
+
+    air = openmc.Material()
+    air.set_density('g/cm3', 0.001)
+    air.add_nuclide('N14', 1.0)
+    metal = openmc.Material()
+    metal.set_density('g/cm3', 7.0)
+    metal.add_nuclide('Fe56', 1.0)
+
+    metal_cell = openmc.Cell(fill=metal)
+    metal_uni = openmc.Universe(cells=[metal_cell])
+    air_cell = openmc.Cell(fill=air)
+    air_uni = openmc.Universe(cells=[air_cell])
+
+    lattice = openmc.RectLattice()
+    lattice.lower_left = (-pitch*n/2, -pitch*n/2)
+    lattice.pitch = (pitch, pitch)
+    lattice.outer = air_uni
+    lattice.universes = [
+        [metal_uni, air_uni, metal_uni],
+        [air_uni, metal_uni, air_uni],
+        [metal_uni, air_uni, metal_uni],
+    ]
+
+    box = openmc.model.RectangularPrism(pitch*n, pitch*n, boundary_type='vacuum')
+    root_cell = openmc.Cell(region=-box, fill=lattice)
+    model.geometry = openmc.Geometry([root_cell])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 10
+    model.settings.particles = 1000
+
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (6, 6)
+    mesh.lower_left = (-pitch*n/2, -pitch*n/2)
+    mesh.upper_right = (pitch*n/2, pitch*n/2)
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally = openmc.Tally(tally_id=1)
+    tally.filters = [mesh_filter]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+
+    return model
+
+
+def test_lattice_large_pitch(model):
+    harness = PyAPITestHarness('statepoint.10.h5', model)
+    harness.main()