From 2d931f63cb4d611d0d23d694726889647f8a482d Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 22 Jun 2022 15:03:50 -0500 Subject: [PATCH] Maintain the high end of the 'roundoff domain' in both float and double precision (#2839) * Maintain the high end of the 'roundoff domain' in both float and double precision * fix shadowing * fix warning * fix float conversion warning * fix logic * Update Src/Base/AMReX_Geometry.H * Update Src/Base/AMReX_Geometry.H --- Src/Base/AMReX_Geometry.H | 53 ++++++++++++++++++------- Src/Base/AMReX_Geometry.cpp | 48 ++++++++++------------ Src/Particle/AMReX_ParticleContainerI.H | 19 ++++----- Src/Particle/AMReX_ParticleUtil.H | 6 +-- 4 files changed, 74 insertions(+), 52 deletions(-) diff --git a/Src/Base/AMReX_Geometry.H b/Src/Base/AMReX_Geometry.H index 54a8b8630d3..4238793861d 100644 --- a/Src/Base/AMReX_Geometry.H +++ b/Src/Base/AMReX_Geometry.H @@ -67,6 +67,32 @@ public: int coord; }; + namespace detail { + template + T bisect_prob_hi (amrex::Real plo, amrex::Real phi, amrex::Real idx, int ilo, int ihi, amrex::Real tol) { + T hi = static_cast(phi - tol); + bool safe; + { + int i = int(Math::floor((hi - plo)*idx)) + 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/idx); + T mid = bisect(lo, hi, + [=] AMREX_GPU_HOST_DEVICE (T x) -> T + { + int i = int(Math::floor((x - plo)*idx)) + 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 @@ -168,8 +194,6 @@ public: //! Returns the problem domain. const RealBox& ProbDomain () const noexcept { return prob_domain; } - //! Returns the roundoff domain. - const RealBox& RoundoffDomain () const noexcept { return roundoff_domain; } //! Sets the problem domain. void ProbDomain (const RealBox& rb) noexcept { @@ -193,12 +217,12 @@ public: return {{AMREX_D_DECL(prob_domain.hi(0),prob_domain.hi(1),prob_domain.hi(2))}}; } - GpuArray RoundoffLoArray () const noexcept { - return {{AMREX_D_DECL(roundoff_domain.lo(0),roundoff_domain.lo(1),roundoff_domain.lo(2))}}; - } - - GpuArray RoundoffHiArray () const noexcept { - return {{AMREX_D_DECL(roundoff_domain.hi(0),roundoff_domain.hi(1),roundoff_domain.hi(2))}}; + GpuArray ProbHiArrayInParticleReal () const noexcept { +#ifdef AMREX_SINGLE_PRECISION_PARTICLES + return roundoff_hi_f; +#else + return roundoff_hi_d; +#endif } //! Returns the overall size of the domain by multiplying the ProbLength's together @@ -406,7 +430,7 @@ public: * are sure to be mapped to cells inside the Domain() box. Note that * the same need not be true for all points inside ProbDomain(). */ - bool outsideRoundoffDomain (AMREX_D_DECL(Real x, Real y, Real z)) const; + bool outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const; /** * \brief Returns true if a point is inside the roundoff domain. @@ -414,7 +438,7 @@ public: * are sure to be mapped to cells inside the Domain() box. Note that * the same need not be true for all points inside ProbDomain(). */ - bool insideRoundoffDomain (AMREX_D_DECL(Real x, Real y, Real z)) const; + bool insideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const; /** * \brief Compute the roundoff domain. Public because it contains an @@ -430,10 +454,11 @@ 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_domain" stores a phi - // that is very close to that in prob_domain, and for which all floating point numbers - // inside it according to a naive inequality check will map to a cell inside domain. - RealBox roundoff_domain; + // will map to a cell that is inside "domain". "roundoff_hi_d" and "roundoff_hi_f" each store + // a phi 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_hi_d; + GpuArray roundoff_hi_f; // Box domain; diff --git a/Src/Base/AMReX_Geometry.cpp b/Src/Base/AMReX_Geometry.cpp index 395f17e352b..1457db6b8d1 100644 --- a/Src/Base/AMReX_Geometry.cpp +++ b/Src/Base/AMReX_Geometry.cpp @@ -506,7 +506,6 @@ Geometry::computeRoundoffDomain () inv_dx[k] = 1.0_rt/dx[k]; } - roundoff_domain = prob_domain; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { int ilo = Domain().smallEnd(idim); @@ -516,40 +515,37 @@ Geometry::computeRoundoffDomain () Real idx = InvCellSize(idim); Real deltax = CellSize(idim); -#ifdef AMREX_SINGLE_PRECISION_PARTICLES - Real tolerance = std::max(1.e-4_rt*deltax, 2.e-7_rt*phi); -#else - Real tolerance = std::max(1.e-8_rt*deltax, 1.e-14_rt*phi); -#endif - // bisect the point at which the cell no longer maps to inside the domain - Real lo = static_cast(phi) - Real(0.5)*static_cast(deltax); - Real hi = static_cast(phi) + Real(0.5)*static_cast(deltax); - - Real mid = bisect(lo, hi, - [=] AMREX_GPU_HOST_DEVICE (Real x) -> Real - { - int i = int(Math::floor((x - plo)*idx)) + ilo; - bool inside = i >= ilo && i <= ihi; - return static_cast(inside) - Real(0.5); - }, tolerance); - roundoff_domain.setHi(idim, mid - tolerance); + 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); + + roundoff_hi_f[idim] = detail::bisect_prob_hi (plo, phi, idx, ilo, ihi, ftol); + roundoff_hi_d[idim] = detail::bisect_prob_hi(plo, phi, idx, ilo, ihi, dtol); } } bool -Geometry::outsideRoundoffDomain (AMREX_D_DECL(Real x, Real y, Real z)) const +Geometry::outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const { - bool outside = AMREX_D_TERM(x < roundoff_domain.lo(0) - || x >= roundoff_domain.hi(0), - || y < roundoff_domain.lo(1) - || y >= roundoff_domain.hi(1), - || z < roundoff_domain.lo(2) - || z >= roundoff_domain.hi(2)); +#ifdef AMREX_SINGLE_PRECISION_PARTICLES + bool outside = AMREX_D_TERM(x < prob_domain.lo(0) + || x >= roundoff_hi_f[0], + || y < prob_domain.lo(1) + || y >= roundoff_hi_f[1], + || z < prob_domain.lo(2) + || z >= roundoff_hi_f[2]); +#else + bool outside = AMREX_D_TERM(x < prob_domain.lo(0) + || x >= roundoff_hi_d[0], + || y < prob_domain.lo(1) + || y >= roundoff_hi_d[1], + || z < prob_domain.lo(2) + || z >= roundoff_hi_d[2]); +#endif return outside; } bool -Geometry::insideRoundoffDomain (AMREX_D_DECL(Real x, Real y, Real z)) const +Geometry::insideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const { return !outsideRoundoffDomain(AMREX_D_DECL(x, y, z)); } diff --git a/Src/Particle/AMReX_ParticleContainerI.H b/Src/Particle/AMReX_ParticleContainerI.H index f4ababb3a82..f6f51c572cf 100644 --- a/Src/Particle/AMReX_ParticleContainerI.H +++ b/Src/Particle/AMReX_ParticleContainerI.H @@ -239,7 +239,7 @@ ParticleContainer const auto& geom = Geom(0); const auto plo = geom.ProbLoArray(); const auto phi = geom.ProbHiArray(); - const auto rhi = geom.RoundoffHiArray(); + const auto rhi = geom.ProbHiArrayInParticleReal(); const auto is_per = geom.isPeriodicArray(); return enforcePeriodic(p, plo, phi, rhi, is_per); @@ -314,20 +314,21 @@ ParticleContainer::lo if (! outside) { - if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(Real(p.pos(0)), Real(p.pos(1)), Real(p.pos(2))))) + if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2)))) { - RealBox roundoff_domain = Geom(0).RoundoffDomain(); + RealBox prob_domain = Geom(0).ProbDomain(); + GpuArray phi = Geom(0).ProbHiArrayInParticleReal(); for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { - if (p.pos(idim) <= roundoff_domain.lo(idim)) { - p.pos(idim) = std::nextafter((ParticleReal) roundoff_domain.lo(idim), (ParticleReal) roundoff_domain.hi(idim)); + if (p.pos(idim) <= prob_domain.lo(idim)) { + p.pos(idim) = std::nextafter((ParticleReal) prob_domain.lo(idim), phi[idim]); } - if (p.pos(idim) >= roundoff_domain.hi(idim)) { - p.pos(idim) = std::nextafter((ParticleReal) roundoff_domain.hi(idim), (ParticleReal) roundoff_domain.lo(idim)); + if (p.pos(idim) >= phi[idim]) { + p.pos(idim) = std::nextafter(phi[idim], (ParticleReal) prob_domain.lo(idim)); } } - AMREX_ASSERT(! Geom(0).outsideRoundoffDomain(AMREX_D_DECL(Real(p.pos(0)), Real(p.pos(1)), Real(p.pos(2))))); + AMREX_ASSERT(! Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2)))); } } @@ -1233,7 +1234,7 @@ ParticleContainer Vector > new_sizes(num_levels); const auto plo = Geom(0).ProbLoArray(); const auto phi = Geom(0).ProbHiArray(); - const auto rhi = Geom(0).RoundoffHiArray(); + const auto rhi = Geom(0).ProbHiArrayInParticleReal(); const auto is_per = Geom(0).isPeriodicArray(); for (int lev = lev_min; lev <= finest_lev_particles; ++lev) { diff --git a/Src/Particle/AMReX_ParticleUtil.H b/Src/Particle/AMReX_ParticleUtil.H index 6732a271810..6623f353749 100644 --- a/Src/Particle/AMReX_ParticleUtil.H +++ b/Src/Particle/AMReX_ParticleUtil.H @@ -517,7 +517,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool enforcePeriodic (P& p, amrex::GpuArray const& plo, amrex::GpuArray const& phi, - amrex::GpuArray const& rhi, + amrex::GpuArray const& rhi, amrex::GpuArray const& is_per) noexcept { bool shifted = false; @@ -537,7 +537,7 @@ bool enforcePeriodic (P& p, p.pos(idim) += static_cast(phi[idim] - plo[idim]); } // clamp to avoid precision issues; - if (p.pos(idim) >= rhi[idim]) { + if (p.pos(idim) > rhi[idim]) { p.pos(idim) = static_cast(rhi[idim]); } shifted = true; @@ -555,7 +555,7 @@ int partitionParticlesByDest (PTile& ptile, const PLocator& ploc, const ParticleBufferMap& pmap, const GpuArray& plo, const GpuArray& phi, - const GpuArray& rhi, + const GpuArray& rhi, const GpuArray& is_per, int lev, int gid, int /*tid*/, int lev_min, int lev_max, int nGrow, bool remove_negative)