Fix Roundoff Domain

There was a flaw in amrex::Geometry::computeRoundoffDomain. If probhi is
close to zero, std::nextafter towards negative infinity will be a very small
number (e.g., 1.e-300). So the function will be able to find the lowest
value of roundoff_hi that's outside the domain within a reasonable number of
iterations.

In the new implementation, we use bisection to find the lowest value of
roundoff_lo that's inside the domain and the lowest value fo roundoff_hi
that's outside the domain, up to a tolerance of epsilon times domain size.

Src/Base/AMReX_Geometry.cpp

@@ -522,9 +522,12 @@ Geometry::computeRoundoffDomain ()
     {
         int ilo = Domain().smallEnd(idim);
         int ihi = Domain().bigEnd(idim);
+        int ncells = ihi-ilo+1;
         Real plo = ProbLo(idim);
         Real phi = ProbHi(idim);
         Real dxinv = InvCellSize(idim);
+        auto epsilon = std::numeric_limits<ParticleReal>::epsilon()
+            * ParticleReal(phi-plo);
 
         // Check that the grid is well formed and that deltax > roundoff
         AMREX_ASSERT((plo + ihi*CellSize(idim)) < (plo + (ihi + 1)*CellSize(idim)));
@@ -533,29 +536,122 @@ Geometry::computeRoundoffDomain ()
         // roundoff_hi will be the lowest value that will be outside the domain
         roundoff_lo[idim] = static_cast<ParticleReal>(plo);
         roundoff_hi[idim] = static_cast<ParticleReal>(phi);
+        auto& rlo = roundoff_lo[idim];
+        auto& rhi = roundoff_hi[idim];
 
-        // Hopefully, the loops should never take more than 1 or 2 iterations.
-        // In obscure cases, more may be needed, for example if (phi - plo) is near round off
-        // compared to (phi + plo).
-        int iters = 0;
-        while (roundoff_lo[idim] < plo && iters < 20) {
-            roundoff_lo[idim] = std::nextafter(roundoff_lo[idim], roundoff_hi[idim]);
-            iters++;
+        auto is_outside = [=] (auto x) -> bool
+        {
+            auto idx = int(std::floor((x - plo)*dxinv));
+            return (idx < 0) || (idx >= ncells);
+        };
+
+        auto is_inside = [=] (auto x) -> bool
+        {
+            return !is_outside(x);
+        };
+
+        constexpr int maxiters = 100;
+        static int debug_nmax = 0;
+        int debug_nmax_this = 0;
+
+        ParticleReal rlo_out;
+        if (is_inside(rlo))
+        {
+            auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            int iters = 0;
+            while (is_inside(rlo-eps) && iters < maxiters) {
+                eps *= ParticleReal(2.);
+                ++iters;
+            }
+            rlo_out = rlo - eps;
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+        }
+        else
+        {
+            auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            int iters = 0;
+            while (is_outside(rlo+eps) && iters < maxiters) {
+                eps *= ParticleReal(2.);
+                ++iters;
+            }
+            rlo_out = rlo;
+            rlo += eps;
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+        }
+
+        {
+            // rlo_out is outside
+            // rlo     is inside
+            int iters = 0;
+            while (is_inside(rlo-epsilon) && iters < maxiters) {
+                auto rmid = ParticleReal(0.5)*(rlo_out+rlo);
+                if (rmid == rlo_out || rmid == rlo) {
+                    break;
+                }
+                if (is_inside(rmid)) {
+                    rlo = rmid;
+                } else {
+                    rlo_out = rmid;
+                }
+                ++iters;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+        }
+
+        ParticleReal rhi_in;
+        if (is_outside(rhi))
+        {
+            auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            int iters = 0;
+            while (is_outside(rhi-eps) && iters < maxiters) {
+                eps *= ParticleReal(2.);
+                ++iters;
+            }
+            rhi_in = rhi - eps;
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+        }
+        else
+        {
+            auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            int iters = 0;
+            while (is_inside(rhi+eps) && iters < maxiters) {
+                eps *= ParticleReal(2.);
+                ++iters;
+            }
+            rhi_in = rhi;
+            rhi += eps;
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
         }
-        ParticleReal rhi = roundoff_hi[idim];
-        while (iters < 20) {
-            rhi = std::nextafter(rhi, roundoff_lo[idim]);
-            if (int(std::floor((rhi - plo)*dxinv)) >= ihi + 1 - ilo) {
-                roundoff_hi[idim] = rhi;
-            } else {
-                break;
+
+        {
+            // rhi_in is inside
+            // rhi    is outside
+            int iters = 0;
+            while (is_outside(rhi-epsilon) && iters < maxiters) {
+                auto rmid = ParticleReal(0.5)*(rhi_in+rhi);
+                if (rmid == rhi_in || rmid == rhi) {
+                    break;
+                }
+                if (is_outside(rmid)) {
+                    rhi = rmid;
+                } else {
+                    rhi_in = rmid;
+                }
+                ++iters;
             }
-            iters++;
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
         }
 
-        // If iters is at the limit, then the grid is ill formed and reasonable values could not
-        // be found for the round off domain extent
-        AMREX_ASSERT(iters < 20);
+        if (debug_nmax_this > debug_nmax) {
+            debug_nmax = debug_nmax_this;
+            amrex::Print() << "xxxxx max iters = " << debug_nmax << "\n";
+        }
     }
 }