ParticleUtils.H: pass geometry as a parameter to findParticlesInEachCell instead on relying on WarpX::GetInstance() (#5771)

This PR removes the call to `WarpX::GetInstance().Geom(lev);` inside
`findParticlesInEachCell`. Instead, the geometry is passed as a
parameter to the function.

Author: lucafedeli88

Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp

@@ -174,8 +174,8 @@ void DifferentialLuminosity::ComputeDiags (int step)
             ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
             ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
 
-            ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_1 );
-            ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_2 );
+            ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( warpx.Geom(lev), mfi, ptile_1 );
+            ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( warpx.Geom(lev), mfi, ptile_2 );
 
             // Species
             const auto soa_1 = ptile_1.getParticleTileData();

Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp

@@ -194,8 +194,8 @@ void DifferentialLuminosity2D::ComputeDiags (int step)
             ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
             ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
 
-            ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_1 );
-            ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_2 );
+            ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( warpx.Geom(lev), mfi, ptile_1 );
+            ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( warpx.Geom(lev), mfi, ptile_2 );
 
             // species 1
             const auto soa_1 = ptile_1.getParticleTileData();

Source/Evolve/WarpXEvolve.cpp

@@ -233,7 +233,7 @@ WarpX::Evolve (int numsteps)
         // Resample particles
         // +1 is necessary here because value of step seen by user (first step is 1) is different than
         // value of step in code (first step is 0)
-        mypc->doResampling(istep[0]+1, verbose);
+        mypc->doResampling(Geom(), istep[0]+1, verbose);
 
         if (evolve_scheme == EvolveScheme::Explicit) {
             applyMirrors(cur_time);

Source/Particles/Collision/BinaryCollision/BinaryCollision.H

@@ -305,14 +305,14 @@ public:
             global_debye_length_data = global_debye_length_fab.dataPtr();
         }
 
-        amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
-        auto const dV = AMREX_D_TERM(geom.CellSize(0), *geom.CellSize(1), *geom.CellSize(2));
+        amrex::Geometry const& geom_lev = WarpX::GetInstance().Geom(lev);
+        auto const dV = AMREX_D_TERM(geom_lev.CellSize(0), *geom_lev.CellSize(1), *geom_lev.CellSize(2));
 #if defined WARPX_DIM_RZ
         amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
         auto const lo = lbound(cbx);
         auto const hi = ubound(cbx);
         int const nz = hi.y - lo.y + 1;
-        auto const dr = geom.CellSize(0);
+        auto const dr = geom_lev.CellSize(0);
 #endif
 
         auto volume_factor = [=] AMREX_GPU_DEVICE(int i_cell) noexcept {
@@ -334,7 +334,7 @@ public:
 
             // Find the particles that are in each cell of this tile
             WARPX_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findParticlesInEachCell", prof_findParticlesInEachCell);
-            ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
+            ParticleBins bins_1 = findParticlesInEachCell( geom_lev, mfi, ptile_1 );
             WARPX_PROFILE_VAR_STOP(prof_findParticlesInEachCell);
 
             // Loop over cells, and collide the particles in each cell
@@ -551,8 +551,8 @@ public:
 
             // Find the particles that are in each cell of this tile
             WARPX_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findParticlesInEachCell", prof_findParticlesInEachCell);
-            ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
-            ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 );
+            ParticleBins bins_1 = findParticlesInEachCell( geom_lev, mfi, ptile_1 );
+            ParticleBins bins_2 = findParticlesInEachCell( geom_lev, mfi, ptile_2 );
             WARPX_PROFILE_VAR_STOP(prof_findParticlesInEachCell);
 
             // Loop over cells, and collide the particles in each cell

Source/Particles/MultiParticleContainer.H

@@ -31,6 +31,7 @@
 #include <AMReX_BLassert.H>
 #include <AMReX_Box.H>
 #include <AMReX_Config.H>
+#include <AMReX_Geometry.H>
 #include <AMReX_GpuControl.H>
 #include <AMReX_INT.H>
 #include <AMReX_MFIter.H>
@@ -184,9 +185,11 @@ public:
     /**
     * \brief This function loops over all species and performs resampling if appropriate.
     *
+    * @param[in] geom the geometry of all refinement levels
     * @param[in] timestep the current timestep.
+    * @param[in] verbose flag to control if information on the resampling process is printed out.
     */
-    void doResampling (int timestep, bool verbose);
+    void doResampling (const amrex::Vector<amrex::Geometry>& geom, int timestep, bool verbose);
 
 #ifdef WARPX_QED
     /** If Schwinger process is activated, this function is called at every

Source/Particles/MultiParticleContainer.cpp

@@ -1035,14 +1035,15 @@ MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt )
     collisionhandler->doCollisions(cur_time, dt, this);
 }
 
-void MultiParticleContainer::doResampling (const int timestep, const bool verbose)
+void MultiParticleContainer::doResampling (
+    const amrex::Vector<amrex::Geometry>& geom, const int timestep, const bool verbose)
 {
     for (auto& pc : allcontainers)
     {
         // do_resampling can only be true for PhysicalParticleContainers
         if (!pc->do_resampling){ continue; }
 
-        pc->resample(timestep, verbose);
+        pc->resample(geom, timestep, verbose);
     }
 }
 

Source/Particles/PhysicalParticleContainer.H

@@ -288,9 +288,12 @@ public:
     * \brief This function determines if resampling should be done for the current species, and
     * if so, performs the resampling.
     *
+    * @param[in] geom the geometry of all refinement levels
     * @param[in] timestep the current timestep.
+    * @param[in] verbose flag to control if information on the resampling process is printed out.
     */
-    void resample (int timestep, bool verbose=true) final;
+    void resample (
+        const amrex::Vector<amrex::Geometry>& geom, int timestep, bool verbose=true) final;
 
 #ifdef WARPX_QED
     //Functions decleared in WarpXParticleContainer.H

Source/Particles/PhysicalParticleContainer.cpp

@@ -3193,7 +3193,7 @@ PlasmaInjector* PhysicalParticleContainer::GetPlasmaInjector (int i)
     }
 }
 
-void PhysicalParticleContainer::resample (const int timestep, const bool verbose)
+void PhysicalParticleContainer::resample (const amrex::Vector<amrex::Geometry>& geom, const int timestep, const bool verbose)
 {
     // In heavily load imbalanced simulations, MPI processes with few particles will spend most of
     // the time at the MPI synchronization in TotalNumberOfParticles(). Having two profiler entries
@@ -3216,7 +3216,7 @@ void PhysicalParticleContainer::resample (const int timestep, const bool verbose
         {
             for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
             {
-                m_resampler(pti, lev, this);
+                m_resampler(geom[lev], pti, lev, this);
             }
         }
         deleteInvalidParticles();

Source/Particles/Resampling/LevelingThinning.H

@@ -11,6 +11,7 @@
 
 #include "Particles/WarpXParticleContainer_fwd.H"
 
+#include <AMReX_Geometry.H>
 #include <AMReX_REAL.H>
 
 #include <string>
@@ -47,11 +48,14 @@ public:
     /**
      * \brief A method that performs leveling thinning for the considered species.
      *
+     * @param[in] geom_lev the geometry of the current refinement level
      * @param[in] pti WarpX particle iterator of the particles to resample.
      * @param[in] lev the index of the refinement level.
      * @param[in] pc a pointer to the particle container.
      */
-    void operator() (WarpXParIter& pti, int lev, WarpXParticleContainer* pc) const final;
+    void operator() (
+        const amrex::Geometry& geom_lev, WarpXParIter& pti,
+        int lev, WarpXParticleContainer* pc) const final;
 
 private:
     amrex::Real m_target_ratio = amrex::Real(1.5);

Source/Particles/Resampling/LevelingThinning.cpp

@@ -64,8 +64,9 @@ void LevelingThinning::BackwardCompatibility (const std::string& species_name )
     );
 }
 
-void LevelingThinning::operator() (WarpXParIter& pti, const int lev,
-                                   WarpXParticleContainer * const pc) const
+void LevelingThinning::operator() (
+    const amrex::Geometry& geom_lev, WarpXParIter& pti,
+    const int lev, WarpXParticleContainer * const pc) const
 {
     using namespace amrex::literals;
 
@@ -79,7 +80,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev,
     // efficient to directly loop over the particles. Nevertheless, this structure with a loop over
     // the cells is more general and can be readily used to implement almost any other resampling
     // algorithm.
-    auto bins = ParticleUtils::findParticlesInEachCell(lev, pti, ptile);
+    auto bins = ParticleUtils::findParticlesInEachCell(geom_lev, pti, ptile);
 
     const auto n_cells = static_cast<int>(bins.numBins());
     auto *const indices = bins.permutationPtr();

Source/Particles/Resampling/Resampling.H

@@ -12,6 +12,7 @@
 #include "Particles/WarpXParticleContainer_fwd.H"
 
 #include <AMReX_REAL.H>
+#include <AMReX_Geometry.H>
 
 #include <memory>
 #include <string>
@@ -24,7 +25,9 @@ struct ResamplingAlgorithm
     /**
      * \brief Virtual operator() of the abstract ResamplingAlgorithm class
      */
-    virtual void operator() (WarpXParIter& /*pti*/, int /*lev*/, WarpXParticleContainer* /*pc*/) const = 0;
+    virtual void operator() (
+        const amrex::Geometry& /*geom_lev*/, WarpXParIter& /*pti*/,
+        int /*lev*/, WarpXParticleContainer* /*pc*/) const = 0;
 
     /**
      * \brief Virtual destructor of the abstract ResamplingAlgorithm class
@@ -76,11 +79,14 @@ public:
     /**
      * \brief A method that uses the ResamplingAlgorithm object to perform resampling.
      *
+     * @param[in] geom_lev the geometry of the current refinement level.
      * @param[in] pti WarpX particle iterator of the particles to resample.
      * @param[in] lev the index of the refinement level.
      * @param[in] pc a pointer to the particle container.
      */
-    void operator() (WarpXParIter& pti, int lev, WarpXParticleContainer* pc) const;
+    void operator() (
+        const amrex::Geometry& geom_lev, WarpXParIter& pti,
+        int lev, WarpXParticleContainer* pc) const;
 
 private:
     ResamplingTrigger m_resampling_trigger;

Source/Particles/Resampling/Resampling.cpp

@@ -38,7 +38,9 @@ bool Resampling::triggered (const int timestep, const amrex::Real global_numpart
     return m_resampling_trigger.triggered(timestep, global_numparts);
 }
 
-void Resampling::operator() (WarpXParIter& pti, const int lev, WarpXParticleContainer * const pc) const
+void Resampling::operator() (
+    const amrex::Geometry& geom_lev, WarpXParIter& pti,
+    const int lev, WarpXParticleContainer * const pc) const
 {
-    (*m_resampling_algorithm)(pti, lev, pc);
+    (*m_resampling_algorithm)(geom_lev, pti, lev, pc);
 }

Source/Particles/Resampling/VelocityCoincidenceThinning.H

@@ -15,6 +15,8 @@
 #include "Utils/ParticleUtils.H"
 
 #include <AMReX_Algorithm.H>
+#include <AMReX_Geometry.H>
+
 
 /**
  * \brief This class implements a particle merging scheme wherein particles
@@ -45,11 +47,14 @@ public:
     /**
      * \brief A method that performs merging for the considered species.
      *
+     * @param[in] geom_lev the geometry of the current refinement level
      * @param[in] pti WarpX particle iterator of the particles to resample.
      * @param[in] lev the index of the refinement level.
      * @param[in] pc a pointer to the particle container.
      */
-    void operator() (WarpXParIter& pti, int lev, WarpXParticleContainer* pc) const final;
+    void operator() (
+        const amrex::Geometry& geom_lev, WarpXParIter& pti,
+        int lev, WarpXParticleContainer* pc) const final;
 
     /**
      * \brief This merging routine requires functionality to sort a GPU vector

Source/Particles/Resampling/VelocityCoincidenceThinning.cpp

@@ -61,8 +61,9 @@ VelocityCoincidenceThinning::VelocityCoincidenceThinning (const std::string& spe
     }
 }
 
-void VelocityCoincidenceThinning::operator() (WarpXParIter& pti, const int lev,
-                                   WarpXParticleContainer * const pc) const
+void VelocityCoincidenceThinning::operator() (
+    const amrex::Geometry& geom_lev, WarpXParIter& pti,
+    const int lev, WarpXParticleContainer * const pc) const
 {
     using namespace amrex::literals;
 
@@ -85,7 +86,7 @@ void VelocityCoincidenceThinning::operator() (WarpXParIter& pti, const int lev,
     auto * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data();
 
     // Using this function means that we must loop over the cells in the ParallelFor.
-    auto bins = ParticleUtils::findParticlesInEachCell(lev, pti, ptile);
+    auto bins = ParticleUtils::findParticlesInEachCell(geom_lev, pti, ptile);
 
     const auto n_cells = static_cast<int>(bins.numBins());
     auto *const indices = bins.permutationPtr();

Source/Particles/WarpXParticleContainer.H

@@ -28,6 +28,7 @@
 
 #include <AMReX_Array.H>
 #include <AMReX_FArrayBox.H>
+#include <AMReX_Geometry.H>
 #include <AMReX_GpuAllocators.H>
 #include <AMReX_GpuContainers.H>
 #include <AMReX_INT.H>
@@ -485,7 +486,7 @@ public:
      * override the method for every derived class. Note that in practice this function is never
      * called because resample() is only called for PhysicalParticleContainers.
      */
-    virtual void resample (const int /*timestep*/, bool /*verbose*/) {}
+    virtual void resample (const amrex::Vector<amrex::Geometry>& /*geom*/, const int /*timestep*/, bool /*verbose*/) {}
 
     /**
      * When using runtime components, AMReX requires to touch all tiles

Source/Utils/ParticleUtils.H

@@ -11,6 +11,7 @@
 #include "Utils/WarpXConst.H"
 
 #include <AMReX_DenseBins.H>
+#include <AMReX_Geometry.H>
 #include <AMReX_Particles.H>
 
 #include <AMReX_BaseFwd.H>
@@ -24,12 +25,12 @@ namespace ParticleUtils {
      * particles of a given species present in each cell.
      * Note that this does *not* rearrange particle arrays.
      *
-     * @param[in] lev the index of the refinement level.
+     * @param[in] geom_lev the geometry of the current refinement level.
      * @param[in] mfi the MultiFAB iterator.
      * @param[in] ptile the particle tile.
      */
     amrex::DenseBins<typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType>
-    findParticlesInEachCell (int lev,
+    findParticlesInEachCell (amrex::Geometry const& geom_lev,
                              amrex::MFIter const & mfi,
                              WarpXParticleContainer::ParticleTileType & ptile);
 

Source/Utils/ParticleUtils.cpp

@@ -6,8 +6,6 @@
  */
 #include "ParticleUtils.H"
 
-#include "WarpX.H"
-
 #include <AMReX_Algorithm.H>
 #include <AMReX_Array.H>
 #include <AMReX_Box.H>
@@ -32,12 +30,11 @@ namespace ParticleUtils
     using ParticleTileType = typename WarpXParticleContainer::ParticleTileType;
     using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType;
     using ParticleBins = DenseBins<ParticleTileDataType>;
-    using index_type = typename ParticleBins::index_type;
 
     /* Find the particles and count the particles that are in each cell.
        Note that this does *not* rearrange particle arrays */
     ParticleBins
-    findParticlesInEachCell (int lev,
+    findParticlesInEachCell (Geometry const& geom_lev,
                              MFIter const & mfi,
                              ParticleTileType & ptile) {
 
@@ -46,11 +43,10 @@ namespace ParticleUtils
         auto ptd = ptile.getParticleTileData();
 
         // Extract box properties
-        Geometry const& geom = WarpX::GetInstance().Geom(lev);
         Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
         const auto lo = lbound(cbx);
-        const auto dxi = geom.InvCellSizeArray();
-        const auto plo = geom.ProbLoArray();
+        const auto dxi = geom_lev.InvCellSizeArray();
+        const auto plo = geom_lev.ProbLoArray();
 
         // Find particles that are in each cell;
         // results are stored in the object `bins`.