Make NewtonInverse convergence tolerance precision-aware (#136)#137
Merged
Conversation
LocatePoints failed to resolve interior points on refined meshes in
single-precision builds. NewtonInverse_2D/_3D declare convergence on the
reference-space Newton step delta = J^{-1} * res, where J ~ h/2 for an
element of physical size h. The smallest achievable step is therefore
~ (2/h) * |x| * epsilon, a floor that grows as the element size h
shrinks. With the fixed absolute tolerance newtonTolerance = 1e-8, this
floor exceeds the tolerance in real32 on fine meshes, so the iteration
never converges and interior points are rejected. Double precision has
~9 orders more headroom, so the bug only manifests in real32 builds.
Replace the fixed 1e-8 with the precision-aware sqrt(epsilon(1.0_prec)),
which is ~3.4e-4 (real32) / ~1.5e-8 (real64). The resulting physical
placement error ~ sqrt(epsilon)*h/2 is sub-micron for mm-scale elements
and negligible at typical scales, and double-precision behavior is
essentially unchanged. Both NewtonInverse_2D and NewtonInverse_3D use
this constant. (sqrt of an intrinsic inquiry in a parameter init is
already used in this module: pi = 4*atan(1).)
Add points_locate_2d_fine, a single-precision regression test that
locates points on a 50x50 structured mesh over [0,0.1]^2 (h = 2e-3);
the real32 convergence floor there (~6e-6) is far above the old 1e-8
tolerance, so this test would have failed before the fix.
Fixes #136
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Points lying exactly on an element edge or corner failed to locate, in both single and double precision, on the fine structured-mesh regression test (points_locate_2d_fine located 20/25). Three independent precision/geometry defects in the broad-phase of LocatePoints conspired: 1. Bounding boxes were built from geometry%x%interior nodes. For Gauss quadrature the interior nodes are strictly inside [-1,1], so the box under-covers the element by a fixed fraction of its size (~0.02*h for degree 7). A point on an element edge/corner then lies just outside the AABB of every surrounding element (~4e-5 here) and is rejected before the Newton inverse map runs. Build the box from geometry%x%boundary instead: those nodes are evaluated at the element edges (reference coordinate +/-1) and reach the true boundary and corners. This is correct for both Gauss and Gauss-Lobatto node sets. (BuildElementBBoxes_2D / _3D.) 2. The broad-phase slack/pad used a fixed 1e-8 relative factor, which is below epsilon in single precision, so a point whose float32 coordinate rounds just outside the (now correct) box was still rejected. Make it precision-aware: max(1e-8, 100*epsilon) -- unchanged in real64. 3. The reference-coordinate acceptance window was a fixed 1+1e-6. A point on an edge inverts to xi = +/-1, but its coordinate rounding ~|x|*eps is amplified by the Newton map by ~2/h, so the recovered xi can sit ~(2/h)*|x|*eps beyond +/-1 on fine meshes (~1e-5 here), past the 1e-6 window. Use precision-aware locateTolerance = max(1e-6, sqrt(epsilon)) -- unchanged in real64, ~3.4e-4 in real32. With all three fixes points_locate_2d_fine locates 25/25 in single precision (and remains 25/25 in double). Verified against the full serial CPU test suite in a container build (gfortran, GPU off) in both real32 and real64; the remaining failures there are pre-existing and unrelated (MPI launch in-container, and divergence/gradient tests that do not exercise point location). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Fixes #136.
LocatePointsfailed to resolve interior points on refined meshes in single-precision builds, with the located fraction degrading toward zero as the mesh is refined.NewtonInverse_2D/_3D(insrc/SELF_Points_t.f90) declare convergence on the reference-space Newton step:where
delta = J^{-1} * res,res = xTarget - xCuris a residual in physical units, andJ = dX/dxi ~ h/2for an element of physical sizeh. Near the root,resfloors at the floating-point representation limit~ |x| * epsilon, so the smallest achievable step isdelta_floor ~ (2/h) * |x| * epsilon— a floor that grows ashshrinks. With the fixed absolutenewtonTolerance = 1e-8, this floor exceeds the tolerance inreal32on fine meshes,convergednever becomes.true., and interior points are rejected. Double precision has ~9 orders more headroom, so the bug isreal32-only.Change
src/SELF_Constants.f90: replace the fixed1e-8with a precision-aware tolerance:This is
~3.4e-4(real32) /~1.5e-8(real64). The resulting physical placement error~ sqrt(epsilon)*h/2is sub-micron for mm-scale elements and negligible at typical scales, while double-precision behavior is essentially unchanged. BothNewtonInverse_2DandNewtonInverse_3Dreference this constant, so both are fixed.Using
sqrtof an intrinsic inquiry in aparameterinitialization is already established in this module (pi = 4.0_prec*atan(1.0_prec)two lines above), and it-fsyntax-onlycompiles clean in default,-DDOUBLE_PRECISION, and single-precision modes.Regression test
test/points_locate_2d_fine.f90builds a fine structured quad mesh (50×50 elements over[0,0.1]²,h = 2e-3) and asserts every strictly-interior point is located, that the inverse map reconstructs the physical coordinate, and that field sampling matches the analytic expression. Inreal32,delta_floorat|x| ~ 0.05is~6e-6— ~600× above the old1e-8tolerance — so this test would have failed before the fix. Registered intest/CMakeLists.txtalongside the existingpoints_locate_2dtest.Verification
Local single-precision build was not run (no configured toolchain handy); relying on CI to exercise both precisions.
fprettifyreports no formatting changes for either file.🤖 Generated with Claude Code