diff --git a/src/Hipace.H b/src/Hipace.H index a9ba448252..bb5e686478 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -190,6 +190,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, @@ -206,6 +210,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 diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 5898568ca5..bbe425c1aa 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -92,6 +92,70 @@ Hipace::~Hipace () #endif } +void +Hipace::DefineSliceGDB (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm) +{ + std::map > 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 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 { @@ -111,8 +175,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 @@ -150,9 +218,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 diff --git a/src/fields/Fields.H b/src/fields/Fields.H index f533cffa1b..d75fbaa27f 100644 --- a/src/fields/Fields.H +++ b/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 m_poisson_solver; diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 139157ef63..996e29200c 100644 --- a/src/fields/Fields.cpp +++ b/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 @@ -52,47 +53,6 @@ Fields::AllocData (int lev, const amrex::BoxArray& ba, amrex::MFInfo().SetArena(amrex::The_Arena())); } - std::map > 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 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())); diff --git a/src/particles/PlasmaParticleContainer.H b/src/particles/PlasmaParticleContainer.H index 1fbdb24ef8..b16474e68d 100644 --- a/src/particles/PlasmaParticleContainer.H +++ b/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