Make NewtonInverse convergence tolerance precision-aware (#136)

src/SELF_Constants.f90

@@ -64,8 +64,22 @@ module SELF_Constants
 !                                                             !
 ! ------------------------------------------------------------!
   integer,parameter    :: maxInverseIters = 1000
-  real(prec),parameter :: newtonTolerance = 10.0**(-8)
+  ! Convergence tolerance on the reference-coordinate Newton step. This must be
+  ! precision-aware: the smallest achievable step is ~ (2/h)*|x|*epsilon, which
+  ! grows as the element size h shrinks. A fixed 1e-8 is unreachable in real32
+  ! on fine meshes, causing LocatePoints to reject interior points. Using
+  ! sqrt(epsilon) gives ~3.4e-4 (real32) / ~1.5e-8 (real64); the resulting
+  ! physical placement error ~ sqrt(epsilon)*h/2 is negligible in either case.
+  real(prec),parameter :: newtonTolerance = sqrt(epsilon(1.0_prec))
   integer,parameter    :: newtonMax = 500
+  ! Tolerance for accepting a located reference coordinate as lying inside the
+  ! bi-unit cube [-1,1]^d. A point on an element edge/corner inverts to xi = +/-1,
+  ! but its physical coordinate carries rounding ~|x|*epsilon which the Newton
+  ! map amplifies by ~2/h, so the recovered xi can sit ~ (2/h)*|x|*epsilon beyond
+  ! +/-1 on fine meshes. The acceptance window must exceed that, hence it is
+  ! precision-aware. max(1e-6, sqrt(epsilon)) leaves real64 behavior unchanged
+  ! (1e-6 wins) and widens the window to ~3.4e-4 in real32.
+  real(prec),parameter :: locateTolerance = max(1.0e-6_prec,sqrt(epsilon(1.0_prec)))
 
 !*************************************************************!
 ! ----------------- TIME STEPPING CONSTANTS ------------------!

src/SELF_Points_t.f90

@@ -208,7 +208,14 @@ subroutine LocatePoints_2D_Points_t(this,geometry)
     gMax(1) = maxval(bbMax(1,:))
     gMax(2) = maxval(bbMax(2,:))
     bbDiag = sqrt((gMax(1)-gMin(1))**2+(gMax(2)-gMin(2))**2)
-    slack = max(bbDiag*1.0e-8_prec,tiny(1.0_prec)*1.0e6_prec)
+    ! Geometric slack for the spatial-hash binning and the broad-phase AABB
+    ! reject. The relative factor must stay above the floating-point rounding of
+    ! the coordinates (~epsilon) so that a point sitting on an element edge or
+    ! corner is not spuriously rejected; a fixed 1e-8 is below epsilon in
+    ! single precision. max(1e-8, 100*epsilon) leaves double precision unchanged
+    ! and widens the slack just enough in single precision.
+    slack = max(bbDiag*max(1.0e-8_prec,100.0_prec*epsilon(1.0_prec)), &
+                tiny(1.0_prec)*1.0e6_prec)
     gMin = gMin-slack
     gMax = gMax+slack
 
@@ -299,8 +306,8 @@ subroutine LocatePoints_2D_Points_t(this,geometry)
         if(.not. converged) cycle
 
         ! Accept if reference coords lie in the (slightly inflated) bi-unit cube.
-        if(abs(xi(1)) <= 1.0_prec+1.0e-6_prec .and. &
-           abs(xi(2)) <= 1.0_prec+1.0e-6_prec) then
+        if(abs(xi(1)) <= 1.0_prec+locateTolerance .and. &
+           abs(xi(2)) <= 1.0_prec+locateTolerance) then
           this%elements(p) = iEl
           this%coordinates(p,1) = min(1.0_prec,max(-1.0_prec,xi(1)))
           this%coordinates(p,2) = min(1.0_prec,max(-1.0_prec,xi(2)))
@@ -379,7 +386,14 @@ subroutine LocatePoints_3D_Points_t(this,geometry)
     gMax(2) = maxval(bbMax(2,:))
     gMax(3) = maxval(bbMax(3,:))
     bbDiag = sqrt((gMax(1)-gMin(1))**2+(gMax(2)-gMin(2))**2+(gMax(3)-gMin(3))**2)
-    slack = max(bbDiag*1.0e-8_prec,tiny(1.0_prec)*1.0e6_prec)
+    ! Geometric slack for the spatial-hash binning and the broad-phase AABB
+    ! reject. The relative factor must stay above the floating-point rounding of
+    ! the coordinates (~epsilon) so that a point sitting on an element edge or
+    ! corner is not spuriously rejected; a fixed 1e-8 is below epsilon in
+    ! single precision. max(1e-8, 100*epsilon) leaves double precision unchanged
+    ! and widens the slack just enough in single precision.
+    slack = max(bbDiag*max(1.0e-8_prec,100.0_prec*epsilon(1.0_prec)), &
+                tiny(1.0_prec)*1.0e6_prec)
     gMin = gMin-slack
     gMax = gMax+slack
 
@@ -480,9 +494,9 @@ subroutine LocatePoints_3D_Points_t(this,geometry)
         call NewtonInverse_3D(geometry,iEl,N,xTarget,xi,converged)
         if(.not. converged) cycle
 
-        if(abs(xi(1)) <= 1.0_prec+1.0e-6_prec .and. &
-           abs(xi(2)) <= 1.0_prec+1.0e-6_prec .and. &
-           abs(xi(3)) <= 1.0_prec+1.0e-6_prec) then
+        if(abs(xi(1)) <= 1.0_prec+locateTolerance .and. &
+           abs(xi(2)) <= 1.0_prec+locateTolerance .and. &
+           abs(xi(3)) <= 1.0_prec+locateTolerance) then
           this%elements(p) = iEl
           this%coordinates(p,1) = min(1.0_prec,max(-1.0_prec,xi(1)))
           this%coordinates(p,2) = min(1.0_prec,max(-1.0_prec,xi(2)))
@@ -845,24 +859,32 @@ pure function ClampCell(rIdx,nCells) result(c)
   endfunction ClampCell
 
   subroutine BuildElementBBoxes_2D(geometry,N,nElem,bbMin,bbMax)
-    !! Axis-aligned bounding box of geometry%x%interior nodes for each element.
+    !! Axis-aligned bounding box of each element. The box is taken over the
+    !! element's boundary nodes (geometry%x%boundary), which are evaluated at
+    !! the element edges (reference coordinate = +/-1) and therefore reach the
+    !! true element boundary and corners. Using the interior nodes alone is not
+    !! sufficient: for Gauss quadrature the interior nodes are strictly inside
+    !! [-1,1], so a box built from them under-covers the element by a fixed
+    !! fraction of its size and rejects points lying on element edges/corners.
+    !! The boundary nodes give the correct extent for both Gauss and
+    !! Gauss-Lobatto node sets.
     implicit none
     type(SEMQuad),intent(in) :: geometry
     integer,intent(in) :: N,nElem
     real(prec),intent(out) :: bbMin(1:2,1:nElem),bbMax(1:2,1:nElem)
     ! Local
-    integer :: iEl,i,j
+    integer :: iEl,i,iSide
     real(prec) :: xv,yv
 
     do iEl = 1,nElem
       bbMin(1,iEl) = huge(1.0_prec)
       bbMin(2,iEl) = huge(1.0_prec)
       bbMax(1,iEl) = -huge(1.0_prec)
       bbMax(2,iEl) = -huge(1.0_prec)
-      do j = 1,N+1
+      do iSide = 1,4
         do i = 1,N+1
-          xv = geometry%x%interior(i,j,iEl,1,1)
-          yv = geometry%x%interior(i,j,iEl,1,2)
+          xv = geometry%x%boundary(i,iSide,iEl,1,1)
+          yv = geometry%x%boundary(i,iSide,iEl,1,2)
           if(xv < bbMin(1,iEl)) bbMin(1,iEl) = xv
           if(yv < bbMin(2,iEl)) bbMin(2,iEl) = yv
           if(xv > bbMax(1,iEl)) bbMax(1,iEl) = xv
@@ -874,12 +896,17 @@ subroutine BuildElementBBoxes_2D(geometry,N,nElem,bbMin,bbMax)
   endsubroutine BuildElementBBoxes_2D
 
   subroutine BuildElementBBoxes_3D(geometry,N,nElem,bbMin,bbMax)
+    !! Axis-aligned bounding box of each element, taken over the boundary nodes
+    !! (geometry%x%boundary) which are evaluated on the element faces (reference
+    !! coordinate = +/-1) and so reach the true element boundary, edges and
+    !! corners. See BuildElementBBoxes_2D for why the interior nodes alone are
+    !! insufficient (Gauss interior nodes under-cover the element).
     implicit none
     type(SEMHex),intent(in) :: geometry
     integer,intent(in) :: N,nElem
     real(prec),intent(out) :: bbMin(1:3,1:nElem),bbMax(1:3,1:nElem)
     ! Local
-    integer :: iEl,i,j,k
+    integer :: iEl,i,j,iSide
     real(prec) :: xv,yv,zv
 
     do iEl = 1,nElem
@@ -889,12 +916,12 @@ subroutine BuildElementBBoxes_3D(geometry,N,nElem,bbMin,bbMax)
       bbMax(1,iEl) = -huge(1.0_prec)
       bbMax(2,iEl) = -huge(1.0_prec)
       bbMax(3,iEl) = -huge(1.0_prec)
-      do k = 1,N+1
+      do iSide = 1,6
         do j = 1,N+1
           do i = 1,N+1
-            xv = geometry%x%interior(i,j,k,iEl,1,1)
-            yv = geometry%x%interior(i,j,k,iEl,1,2)
-            zv = geometry%x%interior(i,j,k,iEl,1,3)
+            xv = geometry%x%boundary(i,j,iSide,iEl,1,1)
+            yv = geometry%x%boundary(i,j,iSide,iEl,1,2)
+            zv = geometry%x%boundary(i,j,iSide,iEl,1,3)
             if(xv < bbMin(1,iEl)) bbMin(1,iEl) = xv
             if(yv < bbMin(2,iEl)) bbMin(2,iEl) = yv
             if(zv < bbMin(3,iEl)) bbMin(3,iEl) = zv

test/CMakeLists.txt

@@ -156,6 +156,7 @@ add_fortran_tests (
     "esatmo3d_motionless.f90"
     "esatmo3d_hydrostaticbalance.f90"
     "points_locate_2d.f90"
+    "points_locate_2d_fine.f90"
     "points_locate_2d_outside.f90"
     "points_locate_3d.f90"
     "points_eval_cache_2d.f90"

test/points_locate_2d_fine.f90

@@ -0,0 +1,184 @@
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+!
+! Maintainers : support@fluidnumerics.com
+! Official Repository : https://github.com/FluidNumerics/self/
+!
+! Copyright © 2024 Fluid Numerics LLC
+!
+! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+!
+! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+!
+! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
+!    the documentation and/or other materials provided with the distribution.
+!
+! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
+!    this software without specific prior written permission.
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+!
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+
+! Regression test for https://github.com/FluidNumerics/SELF/issues/136
+!
+! LocatePoints / NewtonInverse_2D must resolve interior points on a refined mesh
+! even in single precision. The Newton convergence test is on the reference-space
+! step delta = J^{-1} * res, with J ~ h/2 for an element of physical size h, so the
+! smallest achievable step is delta_floor ~ (2/h) * |x| * epsilon. This floor grows
+! as h shrinks; with a fixed absolute tolerance of 1e-8 the iteration could never
+! converge in real32 on fine meshes, and LocatePoints rejected interior points.
+!
+! Here we build a fine structured mesh (50 x 50 elements over [0,0.1]^2, so
+! h = 2e-3) and require that every strictly-interior point is located and that the
+! inverse map reconstructs the physical coordinate. In real32, delta_floor at
+! |x| ~ 0.05 is ~6e-6 -- far above the old 1e-8 tolerance (which this test would
+! have failed) but well below the precision-aware tolerance sqrt(epsilon) ~ 3.4e-4.
+program test
+
+  implicit none
+  integer :: exit_code
+
+  exit_code = points_locate_2d_fine()
+  if(exit_code /= 0) then
+    stop exit_code
+  endif
+
+contains
+  integer function points_locate_2d_fine() result(r)
+
+    use SELF_Constants
+    use SELF_Lagrange
+    use SELF_Mesh_2D
+    use SELF_Geometry_2D
+    use SELF_MappedScalar_2D
+    use SELF_Points
+
+    implicit none
+
+    integer,parameter :: controlDegree = 7
+    integer,parameter :: targetDegree = 16
+    integer,parameter :: nvar = 1
+    integer,parameter :: nPoints = 25
+    ! 50 x 50 elements (two 25-element tiles per direction) of width 2e-3 over
+    ! the domain [0, 0.1]^2, giving a small element size h = 2e-3.
+    integer,parameter :: nElemPerTile = 25
+    integer,parameter :: nTile = 2
+    real(prec),parameter :: dxElem = 2.0e-3_prec
+#ifdef DOUBLE_PRECISION
+    real(prec),parameter :: locateTol = 1.0e-6_prec
+    real(prec),parameter :: evalTol = 1.0e-5_prec
+#else
+    real(prec),parameter :: locateTol = 1.0e-4_prec
+    real(prec),parameter :: evalTol = 1.0e-3_prec
+#endif
+    real(prec),parameter :: lDomain = real(nElemPerTile*nTile,prec)*dxElem ! = 0.1
+    type(Lagrange),target :: interp
+    type(Mesh2D),target :: mesh
+    type(SEMQuad),target :: geometry
+    type(MappedScalar2D) :: f
+    type(Points) :: pts
+    real(prec) :: xIn(nPoints,2),fExact,xLoc,yLoc,err
+    real(prec) :: fSampled(nPoints,nvar)
+    real(prec) :: lS(0:controlDegree),lT(0:controlDegree)
+    real(prec) :: xReconstructed,yReconstructed
+    integer :: bcids(1:4)
+    integer :: p,iEl,i,j,found
+
+    call interp%Init(N=controlDegree, &
+                     controlNodeType=GAUSS, &
+                     M=targetDegree, &
+                     targetNodeType=UNIFORM)
+
+    ! Fine structured quad mesh: 50 x 50 elements over [0, 0.1]^2.
+    bcids(1:4) = [SELF_BC_PRESCRIBED, & ! South
+                  SELF_BC_PRESCRIBED, & ! East
+                  SELF_BC_PRESCRIBED, & ! North
+                  SELF_BC_PRESCRIBED] ! West
+    call mesh%StructuredMesh(nElemPerTile,nElemPerTile,nTile,nTile, &
+                             dxElem,dxElem,bcids)
+
+    call geometry%Init(interp,mesh%nElem)
+    call geometry%GenerateFromMesh(mesh)
+
+    call f%Init(interp,nvar,mesh%nelem)
+    call f%AssociateGeometry(geometry)
+    call f%SetEquation(1,'f = sin(3.14159265358979*x)*cos(3.14159265358979*y)')
+    call f%SetInteriorFromEquation(geometry,0.0_prec)
+
+    ! Pseudo-random sample points strictly inside the domain, kept away from the
+    ! origin so |x| is order 0.05 and the single-precision convergence floor is
+    ! firmly above the old 1e-8 tolerance.
+    do p = 1,nPoints
+      xIn(p,1) = lDomain*(0.1_prec+0.8_prec*real(modulo(7*p+1,nPoints),prec)/real(nPoints,prec))
+      xIn(p,2) = lDomain*(0.1_prec+0.8_prec*real(modulo(11*p+3,nPoints),prec)/real(nPoints,prec))
+    enddo
+
+    call pts%Init(nPoints,2)
+    call pts%SetPoints(xIn)
+    call pts%LocatePoints(geometry)
+
+    ! Every point must have been located (all are inside the mesh). Pre-fix, this
+    ! count dropped well below nPoints in single precision.
+    found = 0
+    do p = 1,nPoints
+      if(pts%elements(p) > 0) found = found+1
+    enddo
+    print*,"located ",found," of ",nPoints," points"
+    if(found /= nPoints) then
+      r = 1
+      return
+    endif
+
+    ! Verify inverse map: reconstruct x from (iEl, s, t) and compare to xIn.
+    err = 0.0_prec
+    do p = 1,nPoints
+      iEl = pts%elements(p)
+      lS = interp%CalculateLagrangePolynomials(pts%coordinates(p,1))
+      lT = interp%CalculateLagrangePolynomials(pts%coordinates(p,2))
+      xReconstructed = 0.0_prec
+      yReconstructed = 0.0_prec
+      do j = 1,controlDegree+1
+        do i = 1,controlDegree+1
+          xReconstructed = xReconstructed+lS(i-1)*lT(j-1)*geometry%x%interior(i,j,iEl,1,1)
+          yReconstructed = yReconstructed+lS(i-1)*lT(j-1)*geometry%x%interior(i,j,iEl,1,2)
+        enddo
+      enddo
+      err = max(err,abs(xReconstructed-xIn(p,1)),abs(yReconstructed-xIn(p,2)))
+    enddo
+    print*,"max inverse-map residual: ",err
+    if(err > locateTol) then
+      r = 1
+      return
+    endif
+
+    ! Verify field sampling against the analytic expression.
+    call pts%EvaluateScalar(f,fSampled)
+    err = 0.0_prec
+    do p = 1,nPoints
+      xLoc = xIn(p,1)
+      yLoc = xIn(p,2)
+      fExact = sin(3.14159265358979_prec*xLoc)*cos(3.14159265358979_prec*yLoc)
+      err = max(err,abs(fSampled(p,1)-fExact))
+    enddo
+    print*,"max scalar-sample error: ",err,"  tol=",evalTol
+    if(err > evalTol) then
+      r = 1
+      return
+    endif
+
+    call pts%Free()
+    call f%DissociateGeometry()
+    call f%Free()
+    call geometry%Free()
+    call mesh%Free()
+    call interp%Free()
+
+    r = 0
+
+  endfunction points_locate_2d_fine
+endprogram test