Add Tests/RoundoffDomain

Also fix an assertion issue in AMReX-Codes#3247.
WeiqunZhang · Apr 10, 2023 · 49f9c4e · 49f9c4e
1 parent bb47ec9
commit 49f9c4e
Show file tree

Hide file tree

Showing 6 changed files with 121 additions and 2 deletions.
diff --git a/Src/Particle/AMReX_ParticleUtil.H b/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;

diff --git a/Tests/CMakeLists.txt b/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)

diff --git a/Tests/RoundoffDomain/CMakeLists.txt b/Tests/RoundoffDomain/CMakeLists.txt
@@ -0,0 +1,6 @@
+set(_sources main.cpp)
+
+setup_test(_sources)
+
+unset(_sources)
+
diff --git a/Tests/RoundoffDomain/GNUmakefile b/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
diff --git a/Tests/RoundoffDomain/Make.package b/Tests/RoundoffDomain/Make.package
@@ -0,0 +1 @@
+CEXE_sources += main.cpp
diff --git a/Tests/RoundoffDomain/main.cpp b/Tests/RoundoffDomain/main.cpp
@@ -0,0 +1,72 @@
+#include <AMReX.H>
+#include <AMReX_Print.H>
+#include <AMReX_Random.H>
+#include <AMReX_Geometry.H>
+
+using namespace amrex;
+
+int main(int argc, char* argv[])
+{
+    amrex::Initialize(argc,argv);
+    for (int icell = 0; icell < 10000; ++icell)
+    {
+        int ncells = amrex::Random_int(102400) + 4;
+        Box domain(IntVect(0),IntVect(ncells-1));
+
+        for (int ieps = 0; ieps < 10000; ++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() * 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.0);
+                auto rlom = rlo[idim] - epsilon;
+                epsilon = std::numeric_limits<ParticleReal>::epsilon()
+                    * std::max(ParticleReal(geom.CellSize(idim)),std::abs(rhi[idim]))
+                    * ParticleReal(2.0);
+                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::Print() << "Geometry: " << geom << std::endl;
+                    amrex::Print().SetPrecision(17)
+                        << " rlo[" << idim << "] = " << rlo[idim]
+                        << " rhi[" << idim << "] = " << rhi[idim]
+                        << " rlo_minus = " << rlom
+                        << " rhi_plus  = " << rhip << "\n";
+                    amrex::Print() << " ilo = " << index(rlo[idim])
+                                   << " ihi = " << index(rhi[idim])
+                                   << " ilo-1 = " << index(rlom)
+                                   << " ihi+1 = " << index(rhip)
+                                   << "\n";
+                    amrex::Abort("Failed");
+                }
+            }
+        }
+    }
+    amrex::Finalize();
+}