Fix Increment, OK, and EnforcePeriodicWhere for terrain-fitted partic…

…les. (AMReX-Codes#3556) Follow-on to AMReX-Codes#3499. The proposed changes: - [x] fix a bug or incorrect behavior in AMReX - [ ] add new capabilities to AMReX - [ ] changes answers in the test suite to more than roundoff level - [ ] are likely to significantly affect the results of downstream AMReX users - [ ] include documentation in the code and/or rst files, if appropriate Co-authored-by: Ann Almgren <[email protected]>
ajnonaka · Sep 18, 2023 · cbf064a · cbf064a
1 parent e9f4435
commit cbf064a
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 20 deletions.
diff --git a/Src/Particle/AMReX_ParIter.H b/Src/Particle/AMReX_ParIter.H
@@ -111,13 +111,14 @@ protected:
 };
 
 template <typename T_ParticleType, int NArrayReal=0, int NArrayInt=0,
-          template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
+          template<class> class Allocator=DefaultAllocator, class T_CellAssignor=DefaultAssignor>
 class ParIter_impl
-    : public ParIterBase_impl<false, T_ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
+    : public ParIterBase_impl<false, T_ParticleType, NArrayReal, NArrayInt, Allocator, T_CellAssignor>
 {
 public:
 
     using ParticleType=T_ParticleType;
+    using CellAssignor=T_CellAssignor;
     static constexpr int NStructReal = ParticleType::NReal;
     static constexpr int NStructInt = ParticleType::NInt;
 
@@ -139,13 +140,14 @@ public:
 };
 
 template <typename T_ParticleType, int NArrayReal=0, int NArrayInt=0,
-          template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
+          template<class> class Allocator=DefaultAllocator, class T_CellAssignor=DefaultAssignor>
 class ParConstIter_impl
-    : public ParIterBase_impl<true,T_ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
+    : public ParIterBase_impl<true,T_ParticleType, NArrayReal, NArrayInt, Allocator, T_CellAssignor>
 {
 public:
 
     using ParticleType     = T_ParticleType;
+    using CellAssignor     = T_CellAssignor;
     using ContainerType    = ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
     using ParticleTileType = typename ContainerType::ParticleTileType;
     using AoS              = typename ContainerType::AoS;

diff --git a/Src/Particle/AMReX_ParticleContainer.H b/Src/Particle/AMReX_ParticleContainer.H
@@ -1188,7 +1188,7 @@ public:
     void AssignCellDensitySingleLevel (int rho_index, MultiFab& mf, int level,
                                        int ncomp=1, int particle_lvl_offset = 0) const;
 
-    template <typename P>
+    template <typename P, typename Assignor=CellAssignor>
     IntVect Index (const P& p, int lev) const;
 
     /**

diff --git a/Src/Particle/AMReX_ParticleContainerI.H b/Src/Particle/AMReX_ParticleContainerI.H
@@ -66,7 +66,7 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
 
 template <typename ParticleType, int NArrayReal, int NArrayInt,
           template<class> class Allocator, class CellAssignor>
-template <typename P>
+template <typename P, typename Assignor>
 IntVect
 ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>::Index (const P& p, int lev) const
 {
@@ -75,7 +75,7 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
     const auto plo = geom.ProbLoArray();
     const auto dxi = geom.InvCellSizeArray();
 
-    return CellAssignor{}(p, plo, dxi, domain);
+    return Assignor{}(p, plo, dxi, domain);
 }
 
 template <typename ParticleType, int NArrayReal, int NArrayInt,
@@ -196,11 +196,11 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
             AMREX_ASSERT(ba.ixType().cellCentered());
 
             if (local_grid < 0) {
-                iv = Index(p_prime, lev);
+                iv = Index<amrex::Particle<0, 0>, DefaultAssignor>(p_prime, lev);
                 ba.intersections(Box(iv, iv), isects, true, 0);
                 grid = isects.empty() ? -1 : isects[0].first;
             } else {
-                iv = Index(p_prime, lev);
+                iv = Index<amrex::Particle<0, 0>, DefaultAssignor>(p_prime, lev);
                 if (ba[local_grid].contains(iv))
                 {
                     grid = local_grid;
@@ -586,7 +586,8 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
                               amrex::Array4<amrex::Real> const& count)
         {
             const auto& p = make_particle<ConstParticleType>{}(ptd, ip);
-            auto iv = getParticleCell(p, plo, dxi, domain);
+            CellAssignor assignor;
+            IntVect iv = assignor(p, plo, dxi, domain);
             amrex::Gpu::Atomic::AddNoRet(&count(iv), 1.0_rt);
         }, false);
 }

diff --git a/Src/Particle/AMReX_ParticleUtil.H b/Src/Particle/AMReX_ParticleUtil.H
@@ -75,11 +75,9 @@ numParticlesOutOfRange (Iterator const& pti, IntVect nGrow)
     {
         auto p = make_particle<ParticleType>{}(ptd,i);
         if ((p.id() < 0)) { return false; }
-        IntVect iv(
-            AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
-                         int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
-                         int(amrex::Math::floor((p.pos(2)-plo[2])*dxi[2]))));
-        iv += domain.smallEnd();
+        using AssignorType = typename Iterator::CellAssignor;
+        AssignorType assignor;
+        IntVect iv = assignor(p, plo, dxi, domain);
         return !box.contains(iv);
     });
     int hv = amrex::get<0>(reduce_data.value(reduce_op));
@@ -113,11 +111,9 @@ numParticlesOutOfRange (Iterator const& pti, IntVect nGrow)
     {
         ParticleType p(tile_data,i);
         if ((p.id() < 0)) { return false; }
-        auto iv = IntVect(
-            AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
-                         int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
-                         int(amrex::Math::floor((p.pos(2)-plo[2])*dxi[2]))));
-        iv += domain.smallEnd();
+        using AssignorType = typename Iterator::CellAssignor;
+        AssignorType assignor;
+        IntVect iv = assignor(p, plo, dxi, domain);
         return !box.contains(iv);
     });
     int hv = amrex::get<0>(reduce_data.value(reduce_op));