Resize grid based on beam extent

src/ImpactX.H

@@ -84,6 +84,9 @@ namespace impactx
         //! Delete level data
         void ClearLevel (int lev) override;
 
+        //! Resize the mesh, based on the extent of the bunch of particle
+        void ResizeMesh ();
+
         /** these are the physical/beam particles of the simulation */
         std::unique_ptr<ImpactXParticleContainer> m_particle_container;
 

src/ImpactX.cpp

@@ -97,6 +97,24 @@ namespace impactx
         amrex::ignore_unused(lev);
     }
 
+    void ImpactX::ResizeMesh () {
+        // Extract the min and max of the particle positions
+        auto const [x_min, x_max, y_min, y_max, z_min, z_max] = m_particle_container->MinAndMaxPositions();
+        // Resize the domain size
+        // The box is expanded slightly beyond the min and max of particles.
+        // This controlled by the variable `frac` below.
+        const amrex::Real frac=0.1;
+        amrex::RealBox rb(
+            {x_min-frac*(x_max-x_min), y_min-frac*(y_max-y_min), z_min-frac*(z_max-z_min)}, // Low bound
+            {x_max+frac*(x_max-x_min), y_max+frac*(y_max-y_min), z_max+frac*(z_max-z_min)}); // High bound
+        amrex::Geometry::ResetDefaultProbDomain(rb);
+        for (int lev = 0; lev <= max_level; ++lev) {
+            amrex::Geometry g = Geom(lev);
+            g.ProbDomain(rb);
+            amrex::AmrMesh::SetGeometry(lev, g);
+        }
+    }
+
     void ImpactX::evolve (int num_steps)
     {
         BL_PROFILE("ImpactX::evolve");
@@ -106,6 +124,13 @@ namespace impactx
             BL_PROFILE("ImpactX::evolve::step");
             amrex::Print() << " ++++ Starting step=" << step << "\n";
 
+            // Note: The following operation assume that
+            // the particles are in x, y, z coordinates.
+            // Resize the mesh, based on `m_particle_container` extent
+            ResizeMesh();
+            // Redistribute particles in the new mesh
+            m_particle_container->Redistribute();
+
             // push all particles
             Push(*m_particle_container, m_lattice);
 

src/particles/ImpactXParticleContainer.H

@@ -14,6 +14,8 @@
 
 #include <AMReX_Vector.H>
 
+#include <tuple>
+
 
 namespace impactx
 {
@@ -83,6 +85,16 @@ namespace impactx
                        amrex::Vector<amrex::ParticleReal> const & y,
                        amrex::Vector<amrex::ParticleReal> const & z);
 
+        /** Compute the min and max of the particle position in each dimension
+        *
+        * @returns x_min, x_max, y_min, y_max, z_min, z_max
+        */
+        std::tuple<
+            amrex::ParticleReal, amrex::ParticleReal,
+            amrex::ParticleReal, amrex::ParticleReal,
+            amrex::ParticleReal, amrex::ParticleReal>
+        MinAndMaxPositions ();
+
     }; // ImpactXParticleContainer
 
 } // namespace impactx

src/particles/ImpactXParticleContainer.cpp

@@ -82,4 +82,43 @@ namespace impactx
 //        Redistribute();
 
     }
+
+    std::tuple<
+            amrex::ParticleReal, amrex::ParticleReal,
+            amrex::ParticleReal, amrex::ParticleReal,
+            amrex::ParticleReal, amrex::ParticleReal>
+    ImpactXParticleContainer::MinAndMaxPositions ()
+    {
+        using PType = ImpactXParticleContainer::SuperParticleType;
+
+        // Get min and max for the local rank
+        amrex::ReduceOps<amrex::ReduceOpMin, amrex::ReduceOpMin, amrex::ReduceOpMin, amrex::ReduceOpMax, amrex::ReduceOpMax, amrex::ReduceOpMax> reduce_ops;
+        auto r = amrex::ParticleReduce<amrex::ReduceData<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>>(
+            *this,
+            [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
+            {
+                amrex::ParticleReal x = p.pos(0);
+                amrex::ParticleReal y = p.pos(1);
+                amrex::ParticleReal z = p.pos(2);
+                return {x, y, z, x, y, z};
+            },
+            reduce_ops);
+
+        // Get min and max across all ranks
+        std::vector<amrex::ParticleReal> xyz_min = {
+            amrex::get<0>(r),
+            amrex::get<1>(r),
+            amrex::get<2>(r)
+        };
+        amrex::ParallelDescriptor::ReduceRealMin(xyz_min.data(), xyz_min.size());
+        std::vector<amrex::ParticleReal> xyz_max = {
+            amrex::get<3>(r),
+            amrex::get<4>(r),
+            amrex::get<5>(r)
+        };
+        amrex::ParallelDescriptor::ReduceRealMax(xyz_max.data(), xyz_max.size());
+
+        return {xyz_min[0], xyz_min[1], xyz_min[2], xyz_max[0], xyz_max[1], xyz_max[2]};
+    }
+
 } // namespace impactx