diff --git a/Src/Base/AMReX_Geometry.H b/Src/Base/AMReX_Geometry.H index 06948b075d7..75b1fb1fb24 100644 --- a/Src/Base/AMReX_Geometry.H +++ b/Src/Base/AMReX_Geometry.H @@ -67,56 +67,6 @@ public: int coord; }; - namespace detail { - template - T bisect_prob_lo (amrex::Real plo, amrex::Real /*phi*/, amrex::Real dxinv, int ilo, int ihi, amrex::Real tol) { - T lo = static_cast(plo + tol); - bool safe; - { - int i = int(std::floor((lo - plo)*dxinv)) + ilo; - safe = i >= ilo && i <= ihi; - } - if (safe) { - return lo; - } else { - // bisect the point at which the cell no longer maps to inside the domain - T hi = static_cast(plo + 0.5_rt/dxinv); - T mid = bisect(lo, hi, - [=] AMREX_GPU_HOST_DEVICE (T x) -> T - { - int i = int(std::floor((x - plo)*dxinv)) + ilo; - bool inside = i >= ilo && i <= ihi; - return static_cast(inside) - T(0.5); - }, static_cast(tol)); - return mid - static_cast(tol); - } - } - - template - T bisect_prob_hi (amrex::Real plo, amrex::Real phi, amrex::Real dxinv, int ilo, int ihi, amrex::Real tol) { - T hi = static_cast(phi - tol); - bool safe; - { - int i = int(std::floor((hi - plo)*dxinv)) + ilo; - safe = i >= ilo && i <= ihi; - } - if (safe) { - return hi; - } else { - // bisect the point at which the cell no longer maps to inside the domain - T lo = static_cast(phi - 0.5_rt/dxinv); - T mid = bisect(lo, hi, - [=] AMREX_GPU_HOST_DEVICE (T x) -> T - { - int i = int(std::floor((x - plo)*dxinv)) + ilo; - bool inside = i >= ilo && i <= ihi; - return static_cast(inside) - T(0.5); - }, static_cast(tol)); - return mid - static_cast(tol); - } - } - } - class Geometry : public CoordSys @@ -242,18 +192,10 @@ public: } [[nodiscard]] GpuArray ProbLoArrayInParticleReal () const noexcept { -#ifdef AMREX_SINGLE_PRECISION_PARTICLES - return roundoff_lo_f; -#else - return roundoff_lo_d; -#endif + return roundoff_lo; } [[nodiscard]] GpuArray ProbHiArrayInParticleReal () const noexcept { -#ifdef AMREX_SINGLE_PRECISION_PARTICLES - return roundoff_hi_f; -#else - return roundoff_hi_d; -#endif + return roundoff_hi; } //! Returns the overall size of the domain by multiplying the ProbLength's together @@ -489,11 +431,10 @@ private: RealBox prob_domain; // Due to round-off errors, not all floating point numbers for which plo >= x < phi - // will map to a cell that is inside "domain". "roundoff_{lo,hi}_{f,d}" each store - // a position that is very close to that in prob_domain, and for which all doubles and floats less than - // it will map to a cell inside domain. - GpuArray roundoff_lo_d, roundoff_hi_d; - GpuArray roundoff_lo_f, roundoff_hi_f; + // will map to a cell that is inside "domain". "roundoff_{lo,hi}" store + // positions that are very close to that in prob_domain, and for which all doubles and floats + // between them will map to a cell inside domain. + GpuArray roundoff_lo, roundoff_hi; // Box domain; diff --git a/Src/Base/AMReX_Geometry.cpp b/Src/Base/AMReX_Geometry.cpp index 235c7bb7674..5be69baa478 100644 --- a/Src/Base/AMReX_Geometry.cpp +++ b/Src/Base/AMReX_Geometry.cpp @@ -525,36 +525,49 @@ Geometry::computeRoundoffDomain () Real plo = ProbLo(idim); Real phi = ProbHi(idim); Real dxinv = InvCellSize(idim); - Real deltax = CellSize(idim); - Real ftol = std::max(1.e-4_rt*deltax, 2.e-7_rt*phi); - Real dtol = std::max(1.e-8_rt*deltax, 1.e-14_rt*phi); + // Check that the grid is well formed and that deltax > roundoff + AMREX_ASSERT((plo + ihi*CellSize(idim)) < (plo + (ihi + 1)*CellSize(idim))); + + // roundoff_lo will be the lowest value that will be inside the domain + // roundoff_hi will be the lowest value that will be outside the domain + roundoff_lo[idim] = static_cast(plo); + roundoff_hi[idim] = static_cast(phi); + + // 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++; + } + 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; + } + iters++; + } - roundoff_lo_f[idim] = detail::bisect_prob_lo (plo, phi, dxinv, ilo, ihi, ftol); - roundoff_lo_d[idim] = detail::bisect_prob_lo(plo, phi, dxinv, ilo, ihi, dtol); - roundoff_hi_f[idim] = detail::bisect_prob_hi (plo, phi, dxinv, ilo, ihi, ftol); - roundoff_hi_d[idim] = detail::bisect_prob_hi(plo, phi, dxinv, ilo, ihi, dtol); + // 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); } } bool Geometry::outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const { -#ifdef AMREX_SINGLE_PRECISION_PARTICLES - bool outside = AMREX_D_TERM(x < roundoff_lo_f[0] - || x >= roundoff_hi_f[0], - || y < roundoff_lo_f[1] - || y >= roundoff_hi_f[1], - || z < roundoff_lo_f[2] - || z >= roundoff_hi_f[2]); -#else - bool outside = AMREX_D_TERM(x < roundoff_lo_d[0] - || x >= roundoff_hi_d[0], - || y < roundoff_lo_d[1] - || y >= roundoff_hi_d[1], - || z < roundoff_lo_d[2] - || z >= roundoff_hi_d[2]); -#endif + bool outside = AMREX_D_TERM(x < roundoff_lo[0] + || x >= roundoff_hi[0], + || y < roundoff_lo[1] + || y >= roundoff_hi[1], + || z < roundoff_lo[2] + || z >= roundoff_hi[2]); return outside; }