From 3c78ff1be19fdd0667e47608e32f4bb7851aaec2 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Feb 2026 00:52:44 +0200 Subject: [PATCH 1/2] improve checking r crossing --- src/mesh.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index acd90033f32..f8f06381d81 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1893,8 +1893,8 @@ double CylindricalMesh::find_r_crossing( const double inv_denominator = 1.0 / denominator; const double p = (u.x * r.x + u.y * r.y) * inv_denominator; - double c = r.x * r.x + r.y * r.y - r0 * r0; - double D = p * p - c * inv_denominator; + double R = std::sqrt(r.x * r.x + r.y * r.y); + double D = p * p - (R - r0) * (R + r0) * inv_denominator; if (D < 0.0) return INFTY; @@ -1902,7 +1902,7 @@ double CylindricalMesh::find_r_crossing( D = std::sqrt(D); // the solution -p - D is always smaller as -p + D : Check this one first - if (std::abs(c) <= RADIAL_MESH_TOL) + if (xt::isclose(R, r0, RADIAL_MESH_TOL, RADIAL_MESH_TOL)) return INFTY; if (-p - D > l) @@ -2178,10 +2178,10 @@ double SphericalMesh::find_r_crossing( if (r0 == 0.0) return INFTY; const double p = r.dot(u); - double c = r.dot(r) - r0 * r0; - double D = p * p - c; + double R = r.norm(); + double D = p * p - (R - r0) * (R + r0); - if (std::abs(c) <= RADIAL_MESH_TOL) + if (xt::isclose(R, r0, RADIAL_MESH_TOL, RADIAL_MESH_TOL)) return INFTY; if (D >= 0.0) { From c5a69210dd48c800deb56de6332c255a26968b48 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 12 Feb 2026 16:21:46 -0600 Subject: [PATCH 2/2] Replace xt::isclose and improve comments --- src/mesh.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index f8f06381d81..789113ccc31 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1901,10 +1901,11 @@ double CylindricalMesh::find_r_crossing( D = std::sqrt(D); - // the solution -p - D is always smaller as -p + D : Check this one first - if (xt::isclose(R, r0, RADIAL_MESH_TOL, RADIAL_MESH_TOL)) + // Particle is already on the shell surface; avoid spurious crossing + if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0))) return INFTY; + // Check -p - D first because it is always smaller as -p + D if (-p - D > l) return -p - D; if (-p + D > l) @@ -2181,12 +2182,13 @@ double SphericalMesh::find_r_crossing( double R = r.norm(); double D = p * p - (R - r0) * (R + r0); - if (xt::isclose(R, r0, RADIAL_MESH_TOL, RADIAL_MESH_TOL)) + // Particle is already on the shell surface; avoid spurious crossing + if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0))) return INFTY; if (D >= 0.0) { D = std::sqrt(D); - // the solution -p - D is always smaller as -p + D : Check this one first + // Check -p - D first because it is always smaller as -p + D if (-p - D > l) return -p - D; if (-p + D > l)