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 iteration, 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.
AMReX-Codes · Apr 10, 2023 · eeaf026 · eeaf026
1 parent 8cae802
commit eeaf026
Showing 1 changed file with 138 additions and 18 deletions.
diff --git a/Src/Base/AMReX_Geometry.cpp b/Src/Base/AMReX_Geometry.cpp
@@ -522,6 +522,7 @@ 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);
@@ -533,29 +534,148 @@ 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++;
+        constexpr auto rlowest = std::numeric_limits<ParticleReal>::lowest();
+
+        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 = 200;
+        static int debug_nmax = 0;
+        int debug_nmax_this = 0;
+
+        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;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+
+            auto rlo_out = rlo - eps; // rlo_out is outside. rlo is inside.
+            iters = 0;
+            eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            while (is_inside(rlo-eps) && 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 = 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;
+        else
+        {
+            auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            int iters = 0;
+            while (is_outside(rlo+eps) && iters < maxiters) {
+                eps *= ParticleReal(2.);
+                ++iters;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+
+            auto rlo_out = rlo; // rlo_out is outside
+            rlo += eps;         // rlo     is inside
+            iters = 0;
+            eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
+            while (is_inside(rlo-eps) && 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;
             }
-            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 (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;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+
+            auto rhi_in = rhi - eps; // rhi_in is inside. rhi is outside.
+            iters = 0;
+            while (is_outside(std::nextafter(rhi,rlowest)) && 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;
+            }
+            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;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+
+            auto rhi_in = rhi; // rhi_in is inside
+            rhi += eps;        // rhi    is outside
+            iters = 0;
+            while (is_outside(std::nextafter(rhi,rlowest)) && 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;
+            }
+            AMREX_ALWAYS_ASSERT(iters < maxiters);
+            debug_nmax_this = std::max(debug_nmax_this, iters);
+        }
+
+        if (debug_nmax_this > debug_nmax) {
+            debug_nmax = debug_nmax_this;
+            amrex::Print() << "xxxxx max iters = " << debug_nmax << "\n";
+        }
     }
 }