Maintain the high end of the 'roundoff domain' in both float and double precision (AMReX-Codes#2839)

atmyers · web-flow · commit 2d931f63cb4d · 2022-06-22T22:03:50.000+02:00
* 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
diff --git a/Src/Base/AMReX_Geometry.H b/Src/Base/AMReX_Geometry.H
@@ -67,6 +67,32 @@ public:
     int coord;
 };
 
+    namespace detail {
+        template <typename T>
+        T bisect_prob_hi (amrex::Real plo, amrex::Real phi, amrex::Real idx, int ilo, int ihi, amrex::Real tol) {
+            T hi = static_cast<T>(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<T>(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<T>(inside) - T(0.5);
+                                   }, static_cast<T>(tol));
+                return mid - static_cast<T>(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<Real,AMREX_SPACEDIM> RoundoffLoArray () const noexcept {
-        return {{AMREX_D_DECL(roundoff_domain.lo(0),roundoff_domain.lo(1),roundoff_domain.lo(2))}};
-    }
-
-    GpuArray<Real,AMREX_SPACEDIM> RoundoffHiArray () const noexcept {
-        return {{AMREX_D_DECL(roundoff_domain.hi(0),roundoff_domain.hi(1),roundoff_domain.hi(2))}};
+    GpuArray<ParticleReal,AMREX_SPACEDIM> 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,15 +430,15 @@ 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.
     *        All particles with positions inside the roundoff domain
     *        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<double, AMREX_SPACEDIM> roundoff_hi_d;
+    GpuArray<float , AMREX_SPACEDIM> roundoff_hi_f;
 
     //
     Box     domain;
diff --git 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<Real>(phi) - Real(0.5)*static_cast<Real>(deltax);
-        Real hi = static_cast<Real>(phi) + Real(0.5)*static_cast<Real>(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<Real>(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<float> (plo, phi, idx, ilo, ihi, ftol);
+        roundoff_hi_d[idim] = detail::bisect_prob_hi<double>(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
@@ -239,7 +239,7 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
     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<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>::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<ParticleReal, AMREX_SPACEDIM> 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<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
     Vector<std::map<int, int> > 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
@@ -517,7 +517,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 bool enforcePeriodic (P& p,
                       amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
                       amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& phi,
-                      amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& rhi,
+                      amrex::GpuArray<amrex::ParticleReal,AMREX_SPACEDIM> const& rhi,
                       amrex::GpuArray<int,AMREX_SPACEDIM> const& is_per) noexcept
 {
     bool shifted = false;
@@ -537,7 +537,7 @@ bool enforcePeriodic (P& p,
                 p.pos(idim) += static_cast<ParticleReal>(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<ParticleReal>(rhi[idim]);
             }
             shifted = true;
@@ -555,7 +555,7 @@ int
 partitionParticlesByDest (PTile& ptile, const PLocator& ploc, const ParticleBufferMap& pmap,
                           const GpuArray<Real,AMREX_SPACEDIM>& plo,
                           const GpuArray<Real,AMREX_SPACEDIM>& phi,
-                          const GpuArray<Real,AMREX_SPACEDIM>& rhi,
+                          const GpuArray<ParticleReal,AMREX_SPACEDIM>& rhi,
                           const GpuArray<int ,AMREX_SPACEDIM>& is_per,
                           int lev, int gid, int /*tid*/,
                           int lev_min, int lev_max, int nGrow, bool remove_negative)

Original file line number	Diff line number	Diff line change
`@@ -239,7 +239,7 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>`
`239`	`239`	`const auto& geom = Geom(0);`
`240`	`240`	`const auto plo = geom.ProbLoArray();`
`241`	`241`	`const auto phi = geom.ProbHiArray();`
`242`		`- const auto rhi = geom.RoundoffHiArray();`
	`242`	`+ const auto rhi = geom.ProbHiArrayInParticleReal();`
`243`	`243`	`const auto is_per = geom.isPeriodicArray();`
`244`	`244`
`245`	`245`	`return enforcePeriodic(p, plo, phi, rhi, is_per);`
`@@ -314,20 +314,21 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>::lo`
`314`	`314`
`315`	`315`	`if (! outside)`
`316`	`316`	`{`
`317`		`- if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(Real(p.pos(0)), Real(p.pos(1)), Real(p.pos(2)))))`
	`317`	`+ if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))))`
`318`	`318`	`{`
`319`		`- RealBox roundoff_domain = Geom(0).RoundoffDomain();`
	`319`	`+ RealBox prob_domain = Geom(0).ProbDomain();`
	`320`	`+ GpuArray<ParticleReal, AMREX_SPACEDIM> phi = Geom(0).ProbHiArrayInParticleReal();`
`320`	`321`	`for (int idim=0; idim < AMREX_SPACEDIM; ++idim)`
`321`	`322`	`{`
`322`		`- if (p.pos(idim) <= roundoff_domain.lo(idim)) {`
`323`		`- p.pos(idim) = std::nextafter((ParticleReal) roundoff_domain.lo(idim), (ParticleReal) roundoff_domain.hi(idim));`
	`323`	`+ if (p.pos(idim) <= prob_domain.lo(idim)) {`
	`324`	`+ p.pos(idim) = std::nextafter((ParticleReal) prob_domain.lo(idim), phi[idim]);`
`324`	`325`	`}`
`325`		`- if (p.pos(idim) >= roundoff_domain.hi(idim)) {`
`326`		`- p.pos(idim) = std::nextafter((ParticleReal) roundoff_domain.hi(idim), (ParticleReal) roundoff_domain.lo(idim));`
	`326`	`+ if (p.pos(idim) >= phi[idim]) {`
	`327`	`+ p.pos(idim) = std::nextafter(phi[idim], (ParticleReal) prob_domain.lo(idim));`
`327`	`328`	`}`
`328`	`329`	`}`
`329`	`330`
`330`		`- AMREX_ASSERT(! Geom(0).outsideRoundoffDomain(AMREX_D_DECL(Real(p.pos(0)), Real(p.pos(1)), Real(p.pos(2)))));`
	`331`	`+ AMREX_ASSERT(! Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))));`
`331`	`332`	`}`
`332`	`333`	`}`
`333`	`334`
`@@ -1233,7 +1234,7 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>`
`1233`	`1234`	`Vector<std::map<int, int> > new_sizes(num_levels);`
`1234`	`1235`	`const auto plo = Geom(0).ProbLoArray();`
`1235`	`1236`	`const auto phi = Geom(0).ProbHiArray();`
`1236`		`- const auto rhi = Geom(0).RoundoffHiArray();`
	`1237`	`+ const auto rhi = Geom(0).ProbHiArrayInParticleReal();`
`1237`	`1238`	`const auto is_per = Geom(0).isPeriodicArray();`
`1238`	`1239`	`for (int lev = lev_min; lev <= finest_lev_particles; ++lev)`
`1239`	`1240`	`{`