From 0cf2735e67011af83a186505a2462e0919442c5d Mon Sep 17 00:00:00 2001
From: MaxThevenet <maxence.thevenet@desy.de>
Date: Tue, 20 Jul 2021 14:06:59 +0200
Subject: [PATCH] enforceBC works for all levels (#566)

* enforceBC is level-agnostic

* doc
---
 src/particles/pusher/BeamParticleAdvance.cpp  |  6 ++--
 src/particles/pusher/GetAndSetPosition.H      | 14 ++++++---
 src/particles/pusher/GetDomainLev.H           | 31 +++++++++++++++++++
 .../pusher/PlasmaParticleAdvance.cpp          |  5 ++-
 4 files changed, 48 insertions(+), 8 deletions(-)
 create mode 100644 src/particles/pusher/GetDomainLev.H

diff --git a/src/particles/pusher/BeamParticleAdvance.cpp b/src/particles/pusher/BeamParticleAdvance.cpp
index 810b6f4462..acc0603d92 100644
--- a/src/particles/pusher/BeamParticleAdvance.cpp
+++ b/src/particles/pusher/BeamParticleAdvance.cpp
@@ -4,6 +4,7 @@
 #include "utils/Constants.H"
 #include "GetAndSetPosition.H"
 #include "utils/HipaceProfilerWrapper.H"
+#include "GetDomainLev.H"
 
 void
 AdvanceBeamParticlesSlice (BeamParticleContainer& beam, Fields& fields, amrex::Geometry const& gm,
@@ -66,8 +67,9 @@ AdvanceBeamParticlesSlice (BeamParticleContainer& beam, Fields& fields, amrex::G
 
     const auto getPosition = GetParticlePosition<BeamParticleContainer>(beam, offset);
     const auto setPosition = SetParticlePosition<BeamParticleContainer>(beam, offset);
-    const auto enforceBC = EnforceBC<BeamParticleContainer>(beam, lev, offset);
-
+    const auto enforceBC = EnforceBC<BeamParticleContainer>(
+        beam, GetDomainLev(gm, box, 1, lev), GetDomainLev(gm, box, 0, lev),
+        gm.isPeriodicArray(), offset);
     const amrex::Real zmin = xyzmin[2];
 
     // Declare a DenseBins to pass it to doDepositionShapeN, although it will not be used.
diff --git a/src/particles/pusher/GetAndSetPosition.H b/src/particles/pusher/GetAndSetPosition.H
index 418d97a005..ca8d7c308d 100644
--- a/src/particles/pusher/GetAndSetPosition.H
+++ b/src/particles/pusher/GetAndSetPosition.H
@@ -146,15 +146,19 @@ struct EnforceBC
 
     /** Constructor.
      * \param a_ptile tile containing the macroparticles
-     * \param lev level of MR
+     * \param plo lower end of the physical domain where particles should stay
+     * \param phi upper end of the physical domain where particles should stay
+     * \param is_per periodicity
      * \param a_offset offset to apply to the particle indices
      */
-    EnforceBC (T_ParTile& a_ptile, const int lev, int a_offset = 0) noexcept
+    EnforceBC (T_ParTile& a_ptile, const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>& plo,
+               const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>& phi,
+               const amrex::GpuArray<int,AMREX_SPACEDIM>& is_per, int a_offset = 0) noexcept
     {
 
-        m_plo    = Hipace::GetInstance().Geom(lev).ProbLoArray();
-        m_phi    = Hipace::GetInstance().Geom(lev).ProbHiArray();
-        m_is_per = Hipace::GetInstance().Geom(lev).isPeriodicArray();
+        m_plo    = plo;
+        m_phi    = phi;
+        m_is_per = is_per;
         AMREX_ALWAYS_ASSERT(m_is_per[0] == m_is_per[1]);
 
         m_periodicity = {true, true, false};
diff --git a/src/particles/pusher/GetDomainLev.H b/src/particles/pusher/GetDomainLev.H
new file mode 100644
index 0000000000..76b55fa778
--- /dev/null
+++ b/src/particles/pusher/GetDomainLev.H
@@ -0,0 +1,31 @@
+#ifndef HIPACE_GETDOMAINLEV_H_
+#define HIPACE_GETDOMAINLEV_H_
+
+#include <AMReX.H>
+
+namespace
+{
+    amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> GetDomainLev (
+        const amrex::Geometry& gm, const amrex::Box& box, bool is_lo, int lev)
+    {
+        if (lev == 0) {
+            return is_lo ? gm.ProbLoArray() : gm.ProbHiArray();
+        } else {
+            const auto dx = gm.CellSizeArray();
+            const amrex::IntVect& small = box.smallEnd();
+            const amrex::IntVect& big = box.bigEnd();
+            const auto plo = gm.ProbLoArray();
+            if (is_lo) {
+                return {{AMREX_D_DECL(plo[0]+(small[0]+0.5)*dx[0],
+                                      plo[1]+(small[1]+0.5)*dx[1],
+                                      plo[2]+(small[2]+0.5)*dx[2])}};
+            } else {
+                return {{AMREX_D_DECL(plo[0]+(big[0]+0.5)*dx[0],
+                                      plo[1]+(big[1]+0.5)*dx[1],
+                                      plo[2]+(big[2]+0.5)*dx[2])}};
+            }
+        }
+    }
+}
+
+#endif // HIPACE_GETANDSETPOSITION_H_
diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp
index 624c033ed6..0048d48f42 100644
--- a/src/particles/pusher/PlasmaParticleAdvance.cpp
+++ b/src/particles/pusher/PlasmaParticleAdvance.cpp
@@ -1,6 +1,7 @@
 #include "PlasmaParticleAdvance.H"
 
 #include "particles/PlasmaParticleContainer.H"
+#include "GetDomainLev.H"
 #include "FieldGather.H"
 #include "PushPlasmaParticles.H"
 #include "UpdateForceTerms.H"
@@ -109,7 +110,9 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, Fields & fields,
         using PTileType = PlasmaParticleContainer::ParticleTileType;
         const auto getPosition = GetParticlePosition<PTileType>(pti.GetParticleTile());
         const auto SetPosition = SetParticlePosition<PTileType>(pti.GetParticleTile());
-        const auto enforceBC = EnforceBC<PTileType>(pti.GetParticleTile(), lev);
+        const auto enforceBC = EnforceBC<PTileType>(
+            pti.GetParticleTile(), GetDomainLev(gm, pti.tilebox(), 1, lev),
+            GetDomainLev(gm, pti.tilebox(), 0, lev), gm.isPeriodicArray());
         const amrex::Real zmin = xyzmin[2];
         const amrex::Real dz = dx[2];