Add Tests/RoundoffDomain

Src/Base/AMReX_Geometry.cpp

@@ -556,23 +556,26 @@ Geometry::computeRoundoffDomain ()
         {
             auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
             int iters = 0;
-            while (is_inside(rlo-eps) && iters < maxiters) {
+            rlo_out = rlo - eps;
+            while (is_inside(rlo_out) && iters < maxiters) {
                 eps *= ParticleReal(2.);
+                rlo_out = rlo - eps;
                 ++iters;
             }
-            rlo_out = rlo - eps;
             numiters = std::max(numiters,iters);
         }
         else
         {
             auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
             int iters = 0;
-            while (is_outside(rlo+eps) && iters < maxiters) {
+            auto rtmp = rlo + eps;
+            while (is_outside(rtmp) && iters < maxiters) {
                 eps *= ParticleReal(2.);
+                rtmp = rlo + eps;
                 ++iters;
             }
             rlo_out = rlo;
-            rlo += eps;
+            rlo = rtmp;
             numiters = std::max(numiters,iters);
         }
 
@@ -606,23 +609,29 @@ Geometry::computeRoundoffDomain ()
         {
             auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
             int iters = 0;
-            while (is_inside(rhi+eps) && iters < maxiters) {
+            rhi_out = rhi + eps;
+            while (is_inside(rhi_out) && iters < maxiters) {
                 eps *= ParticleReal(2.);
+                rhi_out = rhi + eps;
                 ++iters;
             }
-            rhi_out = rhi + eps;
             numiters = std::max(numiters,iters);
         }
         else
         {
             auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
             int iters = 0;
-            while (is_outside(rhi-eps) && iters < maxiters) {
+            // Yes, we have to write it this way for Intel compiler.
+            // is_outside(rhi-eps) could be different from is_outside(rtmp),
+            // where rtmp = rhs-eps.
+            auto rtmp = rhi - eps;
+            while (is_outside(rtmp) && iters < maxiters) {
                 eps *= ParticleReal(2.);
+                rtmp = rhi - eps;
                 ++iters;
             }
             rhi_out = rhi;
-            rhi -= eps;
+            rhi = rtmp;
             numiters = std::max(numiters,iters);
         }
 

Src/Particle/AMReX_ParticleUtil.H

@@ -589,7 +589,7 @@ bool enforcePeriodic (P& p,
             }
             shifted = true;
         }
-        AMREX_ASSERT( (p.pos(idim) >= plo[idim] ) && ( p.pos(idim) <= phi[idim] ));
+        AMREX_ASSERT( (p.pos(idim) >= rlo[idim] ) && ( p.pos(idim) <= rhi[idim] ));
     }
 
     return shifted;

Tests/CMakeLists.txt

@@ -1,7 +1,7 @@
 #
 # List of subdirectories to search for CMakeLists.
 #
-set( AMREX_TESTS_SUBDIRS AsyncOut MultiBlock Reinit Amr CLZ Parser CTOParFor)
+set( AMREX_TESTS_SUBDIRS AsyncOut MultiBlock Reinit Amr CLZ Parser CTOParFor RoundoffDomain)
 
 if (AMReX_PARTICLES)
    list(APPEND AMREX_TESTS_SUBDIRS Particles)

Tests/RoundoffDomain/CMakeLists.txt

@@ -0,0 +1,8 @@
+set(_sources main.cpp)
+set(_input_files)
+
+setup_test(_sources _input_files)
+
+unset(_sources)
+unset(_input_files)
+

Tests/RoundoffDomain/GNUmakefile

@@ -0,0 +1,40 @@
+AMREX_HOME ?= ../../../amrex
+
+PREC ?= 1
+
+ifeq ($(PREC),1)
+  PRECISION = DOUBLE
+  USE_SINGLE_PRECISION_PARTICLES = FALSE
+else ifeq ($(PREC),2)
+  PRECISION = FLOAT
+  USE_SINGLE_PRECISION_PARTICLES = FALSE
+else ifeq ($(PREC),3)
+  PRECISION = DOUBLE
+  USE_SINGLE_PRECISION_PARTICLES = TRUE
+else ifeq ($(PREC),4)
+  PRECISION = FLOAT
+  USE_SINGLE_PRECISION_PARTICLES = TRUE
+endif
+
+DEBUG	= FALSE
+
+DIM	= 3
+
+COMP    = gcc
+
+USE_MPI   = FALSE
+USE_OMP   = FALSE
+USE_CUDA  = FALSE
+USE_HIP   = FALSE
+USE_SYCL  = FALSE
+
+BL_NO_FORT = TRUE
+
+TINY_PROFILE = FALSE
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ./Make.package
+include $(AMREX_HOME)/Src/Base/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules

Tests/RoundoffDomain/Make.package

@@ -0,0 +1 @@
+CEXE_sources += main.cpp

Tests/RoundoffDomain/main.cpp

@@ -0,0 +1,82 @@
+#include <AMReX.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_Print.H>
+#include <AMReX_Random.H>
+#include <random>
+
+using namespace amrex;
+
+int main(int argc, char* argv[])
+{
+    amrex::Initialize(argc,argv);
+    ULong seed;
+    {
+        std::random_device rd; // non-deterministic random numbers
+        std::uniform_int_distribution<ULong> dist(0,std::numeric_limits<ULong>::max());
+        seed = dist(rd);
+        amrex::ResetRandomSeed(seed);
+    }
+    for (int icell = 0; icell < 10000; ++icell)
+    {
+        int ncells = int(amrex::Random_int(102400)) + 4;
+        Box domain(IntVect(0),IntVect(ncells-1));
+
+        for (int ieps = 0; ieps < 1000; ++ieps)
+        {
+            std::array<Real,AMREX_SPACEDIM> rblo{AMREX_D_DECL(Real(0.),Real(-1.),Real(-0.3))};
+            std::array<Real,AMREX_SPACEDIM> rbhi{AMREX_D_DECL(Real(1.),Real( 0.),Real( 0.5))};
+            if (ieps % 100 != 0) {
+                auto eps = (amrex::Random() - Real(0.5)) * Real(1.e-4);
+                AMREX_D_TERM(rblo[0] += eps;,
+                             rblo[1] -= eps;,
+                             rblo[2] += eps);
+                AMREX_D_TERM(rbhi[0] -= eps;,
+                             rbhi[1] += eps;,
+                             rbhi[2] -= eps);
+            }
+
+            RealBox rb(rblo, rbhi);
+            Geometry geom(domain, rb, 0, {AMREX_D_DECL(0,0,0)});
+
+            auto rlo = geom.ProbLoArrayInParticleReal();
+            auto rhi = geom.ProbHiArrayInParticleReal();
+            auto plo = geom.ProbLoArray();
+            auto dxinv = geom.InvCellSizeArray();
+            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                auto index = [&] (ParticleReal x) -> int
+                {
+                    return int(std::floor((x - plo[idim])*dxinv[idim]));
+                };
+                auto epsilon = std::numeric_limits<ParticleReal>::epsilon()
+                    * std::max(ParticleReal(geom.CellSize(idim)),std::abs(rlo[idim]))
+                    * ParticleReal(2.5);
+                auto rlom = rlo[idim] - epsilon;
+                epsilon = std::numeric_limits<ParticleReal>::epsilon()
+                    * std::max(ParticleReal(geom.CellSize(idim)),std::abs(rhi[idim]))
+                    * ParticleReal(2.5);
+                auto rhip = rhi[idim] + epsilon;
+                bool pass = (index(rlom)      == -1)
+                    &&      (index(rlo[idim]) == 0 )
+                    &&      (index(rhi[idim]) == ncells-1)
+                    &&      (index(rhip)      == ncells);
+                if (!pass) {
+                    amrex::AllPrint().SetPrecision(17)
+                        << "Random seed = " << seed << "\n"
+                        << "RealBox: " << rb << "\n"
+                        << "Geometry: " << geom << "\n"
+                        << " rlo[" << idim << "] = " << rlo[idim]
+                        << " rhi[" << idim << "] = " << rhi[idim]
+                        << " rlo_minus = " << rlom
+                        << " rhi_plus  = " << rhip << "\n"
+                        << " ilo = " << index(rlo[idim])
+                        << " ihi = " << index(rhi[idim])
+                        << " ilo-1 = " << index(rlom)
+                        << " ihi+1 = " << index(rhip)
+                        << "\n";
+                    amrex::Abort("Failed");
+                }
+            }
+        }
+    }
+    amrex::Finalize();
+}