reset particlecontainer dm and ba

src/Hipace.H

@@ -188,6 +188,10 @@ private:
     /** Pointer to current (and only) instance of class Hipace */
     static Hipace* m_instance;
 
+    amrex::Geometry m_slice_geom;
+    amrex::DistributionMapping m_slice_dm;
+    amrex::BoxArray m_slice_ba;
+
     /** \brief Predictor-corrector loop to calculate Bx and By.
      * 1. an initial Bx and By value is guessed.
      * 2. Using this Bx and By values, the plasma particles are advanced to the next slice,
@@ -204,6 +208,13 @@ private:
      */
     void PredictorCorrectorLoopToSolveBxBy (const int islice, const int lev);
 
+    /** \brief define Geometry, DistributionMapping and BoxArray for the slice.
+     * The slice MultiFAB and the plasma particles are defined on this GDB.
+     *
+     * \param[in] ba BoxArray of the whole domain
+     * \param[in] dm DistributionMappint of the whole domain
+     */
+    void DefineSliceGDB (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
 };
 
 #endif

src/Hipace.cpp

@@ -90,6 +90,70 @@ Hipace::~Hipace ()
 #endif
 }
 
+void
+Hipace::DefineSliceGDB (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm)
+{
+    std::map<int,amrex::Vector<amrex::Box> > boxes;
+    for (int i = 0; i < ba.size(); ++i) {
+        int rank = dm[i];
+        if (InSameTransverseCommunicator(rank)) {
+            boxes[rank].push_back(ba[i]);
+        }
+    }
+
+    // We assume each process may have multiple Boxes longitude direction, but only one Box in the
+    // transverse direction.  The union of all Boxes on a process is rectangular.  The slice
+    // BoxArray therefore has one Box per process.  The Boxes in the slice BoxArray have one cell in
+    // the longitude direction.  We will use the lowest longitude index in each process to construct
+    // the Boxes.  These Boxes do not have any overlaps. Transversely, there are no gaps between
+    // them.
+
+    amrex::BoxList bl;
+    amrex::Vector<int> procmap;
+    for (auto const& kv : boxes) {
+        int const iproc = kv.first;
+        auto const& boxes_i = kv.second;
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(boxes_i.size() > 0,
+                                         "We assume each process has at least one Box");
+        amrex::Box bx = boxes_i[0];
+        for (int j = 1; j < boxes_i.size(); ++j) {
+            amrex::Box const& bxj = boxes_i[j];
+            for (int idim = 0; idim < Direction::z; ++idim) {
+                AMREX_ALWAYS_ASSERT(bxj.smallEnd(idim) == bx.smallEnd(idim));
+                AMREX_ALWAYS_ASSERT(bxj.bigEnd(idim) == bx.bigEnd(idim));
+                if (bxj.smallEnd(Direction::z) < bx.smallEnd(Direction::z)) {
+                    bx = bxj;
+                }
+            }
+        }
+        bx.setBig(Direction::z, bx.smallEnd(Direction::z));
+        bl.push_back(bx);
+        procmap.push_back(iproc);
+    }
+
+    // Slice BoxArray
+    m_slice_ba = amrex::BoxArray(std::move(bl));
+
+    // Slice DistributionMapping
+    m_slice_dm = amrex::DistributionMapping(std::move(procmap));
+
+    // Slice Geometry
+    constexpr int lev = 0;
+    const int dir = AMREX_SPACEDIM-1;
+    // Set the lo and hi of domain and probdomain in the z direction
+    amrex::RealBox tmp_probdom = Geom(lev).ProbDomain();
+    amrex::Box tmp_dom = Geom(lev).Domain();
+    const amrex::Real dx = Geom(lev).CellSize(dir);
+    const amrex::Real hi = Geom(lev).ProbHi(dir);
+    const amrex::Real lo = hi - dx;
+    tmp_probdom.setLo(dir, lo);
+    tmp_probdom.setHi(dir, hi);
+    tmp_dom.setSmall(dir, 0);
+    tmp_dom.setBig(dir, 0);
+    m_slice_geom = amrex::Geometry(
+        tmp_dom, tmp_probdom, Geom(lev).Coord(), Geom(lev).isPeriodic());
+}
+
 bool
 Hipace::InSameTransverseCommunicator (int rank) const
 {
@@ -109,8 +173,12 @@ Hipace::InitData ()
     SetMaxGridSize(new_max_grid_size);
 
     AmrCore::InitFromScratch(0.0); // function argument is time
+    constexpr int lev = 0;
     m_beam_container.InitData(geom[0]);
-    m_plasma_container.InitData(geom[0]);
+    m_plasma_container.SetParticleBoxArray(lev, m_slice_ba);
+    m_plasma_container.SetParticleDistributionMap(lev, m_slice_dm);
+    m_plasma_container.SetParticleGeometry(lev, m_slice_geom);
+    m_plasma_container.InitData();
 }
 
 void
@@ -148,9 +216,8 @@ Hipace::MakeNewLevelFromScratch (
         dm.define(std::move(procmap));
     }
     SetDistributionMap(lev, dm); // Let AmrCore know
-
-    m_fields.AllocData(lev, ba, dm, Geom(lev));
-
+    DefineSliceGDB(ba, dm);
+    m_fields.AllocData(lev, ba, dm, Geom(lev), m_slice_ba, m_slice_dm);
 }
 
 void

src/fields/Fields.H

@@ -52,9 +52,13 @@ public:
      * \param[in] ba BoxArray for the simulation
      * \param[in] dm DistributionMapping for the whole simulation
      * \param[in] geom Geometry
+     * \param[in] slice_ba BoxArray for the slice
+     * \param[in] slice_dm DistributionMapping for the slice
      */
-    void AllocData (int lev, const amrex::BoxArray& ba,
-                    const amrex::DistributionMapping& dm, amrex::Geometry const& geom);
+    void AllocData (
+        int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
+        amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
+        const amrex::DistributionMapping& slice_dm);
 
     /** Class to handle transverse FFT Poisson solver on 1 slice */
     std::unique_ptr<FFTPoissonSolver> m_poisson_solver;

src/fields/Fields.cpp

@@ -15,8 +15,9 @@ Fields::Fields (Hipace const* a_hipace)
 }
 
 void
-Fields::AllocData (int lev, const amrex::BoxArray& ba,
-                   const amrex::DistributionMapping& dm, amrex::Geometry const& geom)
+Fields::AllocData (
+    int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
+    amrex::Geometry const& geom, const amrex::BoxArray& slice_ba, const amrex::DistributionMapping& slice_dm)
 {
     HIPACE_PROFILE("Fields::AllocData()");
     // Need at least 1 guard cell transversally for transverse derivative
@@ -32,47 +33,6 @@ Fields::AllocData (int lev, const amrex::BoxArray& ba,
                         amrex::MFInfo().SetArena(amrex::The_Arena()));
     }
 
-    std::map<int,amrex::Vector<amrex::Box> > boxes;
-    for (int i = 0; i < ba.size(); ++i) {
-        int rank = dm[i];
-        if (m_hipace->InSameTransverseCommunicator(rank)) {
-            boxes[rank].push_back(ba[i]);
-        }
-    }
-
-    // We assume each process may have multiple Boxes longitude direction, but only one Box in the
-    // transverse direction.  The union of all Boxes on a process is rectangular.  The slice
-    // BoxArray therefore has one Box per process.  The Boxes in the slice BoxArray have one cell in
-    // the longitude direction.  We will use the lowest longitude index in each process to construct
-    // the Boxes.  These Boxes do not have any overlaps. Transversely, there are no gaps between
-    // them.
-
-    amrex::BoxList bl;
-    amrex::Vector<int> procmap;
-    for (auto const& kv : boxes) {
-        int const iproc = kv.first;
-        auto const& boxes_i = kv.second;
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(boxes_i.size() > 0,
-                                         "We assume each process has at least one Box");
-        amrex::Box bx = boxes_i[0];
-        for (int j = 1; j < boxes_i.size(); ++j) {
-            amrex::Box const& bxj = boxes_i[j];
-            for (int idim = 0; idim < Direction::z; ++idim) {
-                AMREX_ALWAYS_ASSERT(bxj.smallEnd(idim) == bx.smallEnd(idim));
-                AMREX_ALWAYS_ASSERT(bxj.bigEnd(idim) == bx.bigEnd(idim));
-                if (bxj.smallEnd(Direction::z) < bx.smallEnd(Direction::z)) {
-                    bx = bxj;
-                }
-            }
-        }
-        bx.setBig(Direction::z, bx.smallEnd(Direction::z));
-        bl.push_back(bx);
-        procmap.push_back(iproc);
-    }
-
-    amrex::BoxArray slice_ba(std::move(bl));
-    amrex::DistributionMapping slice_dm(std::move(procmap));
-
     for (int islice=0; islice<(int) WhichSlice::N; islice++) {
         m_slices[lev][islice].define(slice_ba, slice_dm, FieldComps::nfields, m_slices_nguards,
                                      amrex::MFInfo().SetArena(amrex::The_Arena()));

src/particles/PlasmaParticleContainer.H

@@ -36,17 +36,22 @@ public:
 
     /** Allocate data for the beam particles and initialize particles with requested beam profile
      */
-    void InitData (const amrex::Geometry& geom);
+    void InitData ();
 
-    /** Initialize 1 xy slice of particles, with fixed number of particles per cell */
+    /** Initialize 1 xy slice of particles, with fixed number of particles per cell
+     *
+     * \param[in] a_num_particles_per_cell number of particles per cell in x and y
+     * \param[in] a_thermal_momentum_std standard deviation of the momentum distribution (3d)
+     * \param[in] a_thermal_momentum_mean average momentum (3d)
+     * \param[in] a_density plasma density (SI)
+     * \param[in] a_radius plasma radius. Only particles with x**2+y**2<a_radius**2 are injected
+     */
     void InitParticles (
         const amrex::IntVect& a_num_particles_per_cell,
         const amrex::RealVect& a_thermal_momentum_std,
         const amrex::RealVect& a_thermal_momentum_mean,
         const amrex::Real a_density,
-        const amrex::Real a_radius,
-        const amrex::Geometry& a_geom,
-        const amrex::RealBox& a_bounds);
+        const amrex::Real a_radius);
 
     amrex::Real m_density {0}; /**< Density of the plasma */
     /** maximum weighting factor gamma/(Psi +1) before particle is regarded as violating

src/particles/PlasmaParticleContainer.cpp

@@ -27,19 +27,10 @@ PlasmaParticleContainer::PlasmaParticleContainer (amrex::AmrCore* amr_core)
 }
 
 void
-PlasmaParticleContainer::InitData (const amrex::Geometry& geom)
+PlasmaParticleContainer::InitData ()
 {
     reserveData();
     resizeData();
 
-    const int dir = AMREX_SPACEDIM-1;
-    const amrex::Real dx = geom.CellSize(dir);
-    const amrex::Real hi = geom.ProbHi(dir);
-    const amrex::Real lo = hi - dx;
-
-    amrex::RealBox particleBox = geom.ProbDomain();
-    particleBox.setHi(dir, hi);
-    particleBox.setLo(dir, lo);
-
-    InitParticles(m_ppc,m_u_std, m_u_mean, m_density, m_radius, geom, particleBox);
+    InitParticles(m_ppc,m_u_std, m_u_mean, m_density, m_radius);
 }

src/particles/PlasmaParticleContainerInit.cpp

@@ -12,15 +12,14 @@ InitParticles (const IntVect& a_num_particles_per_cell,
                const amrex::RealVect& a_u_std,
                const amrex::RealVect& a_u_mean,
                const amrex::Real a_density,
-               const amrex::Real a_radius,
-               const Geometry& a_geom,
-               const RealBox& a_bounds)
+               const amrex::Real a_radius)
 {
     HIPACE_PROFILE("PlasmaParticleContainer::InitParticles");
 
     const int lev = 0;
-    const auto dx = a_geom.CellSizeArray();
-    const auto plo = a_geom.ProbLoArray();
+    const auto dx = ParticleGeom(lev).CellSizeArray();
+    const auto plo = ParticleGeom(lev).ProbLoArray();
+    const amrex::RealBox a_bounds = ParticleGeom(lev).ProbDomain();
 
     const int num_ppc = AMREX_D_TERM( a_num_particles_per_cell[0],
                                       *a_num_particles_per_cell[1],