From 7ef97e25b96dd3ec54c8815306601e06fdf2b2c2 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 6 Apr 2023 14:28:57 -0700 Subject: [PATCH] Rework handling of roundoff domain extent (#3199) ## Summary This reworks previous changes that created the round off domain which was done to avoid round off errors when locating particles in the domain, specifically when handling periodic boundaries. This PR changes how this is done so that the roundoff domain makes as small a change as possible. ## Additional background Previous PRs #2839 and #2950 created a round off domain that has a slightly smaller extent that the actual domain. Any particles that ended up outside the round off domain but inside the actual domain were shifted to be on the boundary of the round off domain. The main purpose was to handle problems with periodic boundary with mixed precision (single precision particles and double precision grid). In that case, it is possible for a particle to have `z < lo` and `z+(hi-lo) > hi`. It is also possible for a particle to have `z < hi` but `(z-lo)*dzinv > ihi` putting it in a grid cell beyond the domain. Unfortunately, the roundoff domains was being set further inside that was needed and would affect valid particles that would not have those issues. While the problem is small, the effects can be noticeable when examining simulation properties, for example charge conservation in WarpX. This PR changes the roundoff domain so that only problematic particles are affected. The code now uses `std::nextafter` to find the exact place where the problems occurs. At the lower end, it ensures that the casted position will be greater than ProbLo. At the upper end, it ensures that `(z-lo)*dzinv <= ihi`. One comment is that this assumes that particle grid coordinates will always be calculated at the higher precision, with expressions like `iz = (int)(z - lo)`, without `lo` being cast to `ParticleReal`. Co-authored-by: atmyers --- Src/Base/AMReX_Geometry.H | 71 ++++--------------------------------- Src/Base/AMReX_Geometry.cpp | 57 +++++++++++++++++------------ 2 files changed, 41 insertions(+), 87 deletions(-) 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; }