diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 6fb15c7821..419624b4e5 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -381,7 +381,7 @@ OpenPMDWriter::CopyBeams (MultiBeam& beams, const amrex::Vector< std::string > b if (np != 0) { // copy data from GPU to IO buffer - auto& soa = beam.getBeamSlice(WhichBeamSlice::This).GetStructOfArrays(); + auto& slice = beam.getBeamSlice(WhichBeamSlice::This); for (std::size_t idx=0; idx b ); } amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, - soa.GetIdCPUData().begin(), - soa.GetIdCPUData().begin() + np, + slice.GetIdCPUData().begin(), + slice.GetIdCPUData().begin() + np, m_uint64_beam_data[ibeam][idx].data() + m_offset[ibeam]); } AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - int(m_real_beam_data[ibeam].size()) == soa.NumRealComps(), + int(m_real_beam_data[ibeam].size()) == slice.NumRealComps(), "List of real names in openPMD Writer class does not match the beam"); for (std::size_t idx=0; idx b ); } amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, - soa.GetRealData(idx).begin(), - soa.GetRealData(idx).begin() + np, + slice.GetRealData(idx).begin(), + slice.GetRealData(idx).begin() + np, m_real_beam_data[ibeam][idx].data() + m_offset[ibeam]); } } diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H index 62b9575ed9..a4b3f90a88 100644 --- a/src/particles/beam/BeamParticleContainer.H +++ b/src/particles/beam/BeamParticleContainer.H @@ -28,7 +28,8 @@ struct BeamIdx w, // weight ux, uy, uz, // momentum real_nattribs_in_buffer, - real_nattribs=real_nattribs_in_buffer + real_nattribs=real_nattribs_in_buffer, + sx=real_nattribs, sy, sz }; enum { // no extra components stored in MultiBuffer, besides 64bit idcpu @@ -45,26 +46,7 @@ struct WhichBeamSlice { enum beam_slice : int { Next=0, This, N }; }; -using BeamTile = amrex::ParticleTile< - amrex::SoAParticle< - BeamIdx::real_nattribs, - BeamIdx::int_nattribs - >, - BeamIdx::real_nattribs, - BeamIdx::int_nattribs, - amrex::PolymorphicArenaAllocator - >; - -using BeamTileInit = amrex::ParticleTile< - amrex::SoAParticle< - BeamIdx::real_nattribs_in_buffer, - BeamIdx::int_nattribs_in_buffer - >, - BeamIdx::real_nattribs_in_buffer, - BeamIdx::int_nattribs_in_buffer, - // use PolymorphicArenaAllocator to either use Pinned or Device memory at runtime - amrex::PolymorphicArenaAllocator - >; +using BeamTile = amrex::ParticleTileRT; /** \brief Container for particles of 1 beam species. */ class BeamParticleContainer @@ -185,7 +167,7 @@ public: void resize (int which_slice, int num_particles, int num_slipped_particles); - BeamTileInit& getBeamInitSlice () { + BeamTile& getBeamInitSlice () { return m_init_slice; } @@ -215,7 +197,7 @@ private: std::array m_slices {}; std::array m_num_particles_without_slipped {}; std::array m_num_particles_with_slipped {}; - BeamTileInit m_init_slice {}; + BeamTile m_init_slice {}; BoxSorter m_init_sorter; uint64_t m_total_num_particles = 0; public: diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index 7f74fba361..36bb103afa 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -109,11 +109,22 @@ BeamParticleContainer::ReadParameters () queryWithParserAlt(pp, "spin_anom", m_spin_anom, pp_alt); } - getBeamInitSlice().define(m_do_spin_tracking ? 3 : 0, 0, nullptr, nullptr, - m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena()); + getBeamInitSlice().define( + BeamIdx::real_nattribs_in_buffer + (m_do_spin_tracking ? 3 : 0), + BeamIdx::int_nattribs_in_buffer, + nullptr, + nullptr, + m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena() + ); for (auto& beam_tile : m_slices) { - beam_tile.define(m_do_spin_tracking ? 3 : 0, 0, nullptr, nullptr, amrex::The_Arena()); + beam_tile.define( + BeamIdx::real_nattribs + (m_do_spin_tracking ? 3 : 0), + BeamIdx::int_nattribs, + nullptr, + nullptr, + amrex::The_Arena() + ); } } @@ -271,7 +282,7 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom) m_total_num_particles = getBeamInitSlice().size(); if (Hipace::HeadRank()) { m_init_sorter.sortParticlesByBox( - getBeamInitSlice().GetStructOfArrays().GetRealData(BeamIdx::z).dataPtr(), + getBeamInitSlice().GetRealData(BeamIdx::z).dataPtr(), getBeamInitSlice().size(), m_initialize_on_cpu, geom); } #else @@ -333,10 +344,10 @@ void BeamParticleContainer::TagByLevel (const int current_N_level, { HIPACE_PROFILE("BeamParticleContainer::TagByLevel()"); - auto& soa = getBeamSlice(which_slice).GetStructOfArrays(); - const amrex::Real * const pos_x = soa.GetRealData(BeamIdx::x).data(); - const amrex::Real * const pos_y = soa.GetRealData(BeamIdx::y).data(); - int * const p_mr_level = soa.GetIntData(BeamIdx::mr_level).data(); + auto& slice = getBeamSlice(which_slice); + const amrex::Real * const pos_x = slice.GetRealData(BeamIdx::x).data(); + const amrex::Real * const pos_y = slice.GetRealData(BeamIdx::y).data(); + int * const p_mr_level = slice.GetIntData(BeamIdx::mr_level).data(); const int lev1_idx = std::min(1, current_N_level-1); const int lev2_idx = std::min(2, current_N_level-1); @@ -395,9 +406,9 @@ BeamParticleContainer::initializeSlice (int slice, int which_slice) { ptd.rdata(BeamIdx::uy)[ip] = ptd_init.rdata(BeamIdx::uy)[idx_src]; ptd.rdata(BeamIdx::uz)[ip] = ptd_init.rdata(BeamIdx::uz)[idx_src]; if (do_spin_tracking) { - ptd.m_runtime_rdata[0][ip] = ptd_init.m_runtime_rdata[0][idx_src]; - ptd.m_runtime_rdata[1][ip] = ptd_init.m_runtime_rdata[1][idx_src]; - ptd.m_runtime_rdata[2][ip] = ptd_init.m_runtime_rdata[2][idx_src]; + ptd.rdata(BeamIdx::sx)[ip] = ptd_init.rdata(BeamIdx::sx)[idx_src]; + ptd.rdata(BeamIdx::sy)[ip] = ptd_init.rdata(BeamIdx::sy)[idx_src]; + ptd.rdata(BeamIdx::sz)[ip] = ptd_init.rdata(BeamIdx::sz)[idx_src]; } ptd.idcpu(ip) = ptd_init.idcpu(idx_src); ptd.idata(BeamIdx::nsubcycles)[ip] = 0; @@ -413,9 +424,9 @@ BeamParticleContainer::initializeSlice (int slice, int which_slice) { const amrex::RealVect initial_spin_norm = m_initial_spin / m_initial_spin.vectorLength(); amrex::ParallelFor(getNumParticles(which_slice), [=] AMREX_GPU_DEVICE (const int ip) { - ptd.m_runtime_rdata[0][ip] = initial_spin_norm[0]; - ptd.m_runtime_rdata[1][ip] = initial_spin_norm[1]; - ptd.m_runtime_rdata[2][ip] = initial_spin_norm[2]; + ptd.rdata(BeamIdx::sx)[ip] = initial_spin_norm[0]; + ptd.rdata(BeamIdx::sy)[ip] = initial_spin_norm[1]; + ptd.rdata(BeamIdx::sz)[ip] = initial_spin_norm[2]; } ); } @@ -442,60 +453,11 @@ BeamParticleContainer::ReorderParticles (int beam_slice, int step, amrex::Geomet HIPACE_PROFILE("BeamParticleContainer::ReorderParticles()"); const int np = getNumParticles(beam_slice); - const int np_total = getNumParticlesIncludingSlipped(beam_slice); auto& ptile = getBeamSlice(beam_slice); amrex::Gpu::DeviceVector perm; amrex::PermutationForDeposition(perm, np, ptile, slice_geom.Domain(), slice_geom, m_reorder_idx_type); - const unsigned int* permutations = perm.dataPtr(); - auto& soa = ptile.GetStructOfArrays(); - - { - typename BeamTile::SoA::IdCPU tmp_idcpu; - tmp_idcpu.setArena(amrex::The_Arena()); - tmp_idcpu.resize(np_total); - auto src = soa.GetIdCPUData().data(); - uint64_t* dst = tmp_idcpu.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = i < np ? src[permutations[i]] : src[i]; - }); - - amrex::Gpu::streamSynchronize(); - soa.GetIdCPUData().swap(tmp_idcpu); - } - - { // Create a scope for the temporary vector below - BeamTile::RealVector tmp_real; - tmp_real.setArena(amrex::The_Arena()); - tmp_real.resize(np_total); - for (int comp = 0; comp < soa.NumRealComps(); ++comp) { - auto src = soa.GetRealData(comp).data(); - amrex::ParticleReal* dst = tmp_real.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = i < np ? src[permutations[i]] : src[i]; - }); - - amrex::Gpu::streamSynchronize(); - soa.GetRealData(comp).swap(tmp_real); - } - } - - BeamTile::IntVector tmp_int; - tmp_int.setArena(amrex::The_Arena()); - tmp_int.resize(np_total); - for (int comp = 0; comp < soa.NumIntComps(); ++comp) { - auto src = soa.GetIntData(comp).data(); - int* dst = tmp_int.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = i < np ? src[permutations[i]] : src[i]; - }); - - amrex::Gpu::streamSynchronize(); - soa.GetIntData(comp).swap(tmp_int); - } + amrex::ReorderParticles(ptile, perm.dataPtr()); } } @@ -590,9 +552,9 @@ BeamParticleContainer::InSituComputeDiags (int islice) { const amrex::Real x = ptd.pos(0, ip); const amrex::Real y = ptd.pos(1, ip); - const amrex::Real sx = ptd.m_runtime_rdata[0][ip]; - const amrex::Real sy = ptd.m_runtime_rdata[1][ip]; - const amrex::Real sz = ptd.m_runtime_rdata[2][ip]; + const amrex::Real sx = ptd.rdata(BeamIdx::sx)[ip]; + const amrex::Real sy = ptd.rdata(BeamIdx::sy)[ip]; + const amrex::Real sz = ptd.rdata(BeamIdx::sz)[ip]; const amrex::Real w = ptd.rdata(BeamIdx::w)[ip]; if (!ptd.id(ip).is_valid() || x*x + y*y > insitu_radius_sq) { diff --git a/src/particles/beam/BeamParticleContainerInit.cpp b/src/particles/beam/BeamParticleContainerInit.cpp index 2b41163f18..ac5a2bd3e2 100644 --- a/src/particles/beam/BeamParticleContainerInit.cpp +++ b/src/particles/beam/BeamParticleContainerInit.cpp @@ -40,7 +40,7 @@ namespace */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddOneBeamParticle ( - const BeamTileInit::ParticleTileDataType& ptd, + const BeamTile::ParticleTileDataType& ptd, const amrex::Real& x, const amrex::Real& y, const amrex::Real& z, const amrex::Real& ux, const amrex::Real& uy, const amrex::Real& uz, const amrex::Real& sx, const amrex::Real& sy, const amrex::Real& sz, @@ -60,9 +60,9 @@ namespace ptd.rdata(BeamIdx::uy)[ip] = uyp; ptd.rdata(BeamIdx::uz)[ip] = uz; if (do_spin) { - ptd.m_runtime_rdata[0][ip] = sx; - ptd.m_runtime_rdata[1][ip] = sy; - ptd.m_runtime_rdata[2][ip] = sz; + ptd.rdata(BeamIdx::sx)[ip] = sx; + ptd.rdata(BeamIdx::sy)[ip] = sy; + ptd.rdata(BeamIdx::sz)[ip] = sz; } ptd.rdata(BeamIdx::w )[ip] = std::abs(weight); diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index b688f39285..5d1466f73a 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -88,12 +88,12 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { // Get particles SoA data - auto& soa1 = pti.GetStructOfArrays(); - amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data(); - amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data(); - amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data(); - const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data(); - const int* const ion_lev1 = soa1.GetIntData(PlasmaIdx::ion_lev).data(); + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + amrex::Real* const uy1 = ptile1.GetRealData(PlasmaIdx::uy_half_step).data(); + amrex::Real* const psi1 = ptile1.GetRealData(PlasmaIdx::psi_half_step).data(); + const amrex::Real* const w1 = ptile1.GetRealData(PlasmaIdx::w).data(); + const int* const ion_lev1 = ptile1.GetIntData(PlasmaIdx::ion_lev).data(); PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); amrex::Real q1 = species1.GetCharge(); @@ -154,12 +154,12 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { // Get particles SoA data for species 1 - auto& soa1 = pti.GetStructOfArrays(); - amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data(); - amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data(); - amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data(); - const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data(); - const int* const ion_lev1 = soa1.GetIntData(PlasmaIdx::ion_lev).data(); + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + amrex::Real* const uy1 = ptile1.GetRealData(PlasmaIdx::uy_half_step).data(); + amrex::Real* const psi1 = ptile1.GetRealData(PlasmaIdx::psi_half_step).data(); + const amrex::Real* const w1 = ptile1.GetRealData(PlasmaIdx::w).data(); + const int* const ion_lev1 = ptile1.GetIntData(PlasmaIdx::ion_lev).data(); PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); amrex::Real q1 = species1.GetCharge(); @@ -168,12 +168,11 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( // Get particles SoA data for species 2 auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); - auto& soa2 = ptile2.GetStructOfArrays(); - amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data(); - amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data(); - amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data(); - const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data(); - const int* const ion_lev2 = soa2.GetIntData(PlasmaIdx::ion_lev).data(); + amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); + amrex::Real* const uy2 = ptile2.GetRealData(PlasmaIdx::uy_half_step).data(); + amrex::Real* const psi2= ptile2.GetRealData(PlasmaIdx::psi_half_step).data(); + const amrex::Real* const w2 = ptile2.GetRealData(PlasmaIdx::w).data(); + const int* const ion_lev2 = ptile2.GetIntData(PlasmaIdx::ion_lev).data(); PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); amrex::Real q2 = species2.GetCharge(); @@ -265,11 +264,11 @@ CoulombCollision::doBeamPlasmaCoulombCollision ( for (PlasmaParticleIterator pti(species2); pti.isValid(); ++pti) { // // Get particles SoA data for species 1 - auto& soa1 = species1.getBeamSlice(WhichBeamSlice::This).GetStructOfArrays(); - amrex::Real* const ux1 = soa1.GetRealData(BeamIdx::ux).data(); - amrex::Real* const uy1 = soa1.GetRealData(BeamIdx::uy).data(); - amrex::Real* const psi1 = soa1.GetRealData(BeamIdx::uz).data(); - const amrex::Real* const w1 = soa1.GetRealData(BeamIdx::w).data(); + auto& ptile1 = species1.getBeamSlice(WhichBeamSlice::This); + amrex::Real* const ux1 = ptile1.GetRealData(BeamIdx::ux).data(); + amrex::Real* const uy1 = ptile1.GetRealData(BeamIdx::uy).data(); + amrex::Real* const psi1 = ptile1.GetRealData(BeamIdx::uz).data(); + const amrex::Real* const w1 = ptile1.GetRealData(BeamIdx::w).data(); BeamBins::index_type * const indices1 = bins1.permutationPtr(); BeamBins::index_type const * const offsets1 = bins1.offsetsPtr(); amrex::Real q1 = species1.GetCharge(); @@ -278,12 +277,12 @@ CoulombCollision::doBeamPlasmaCoulombCollision ( // Get particles SoA data for species 2 //auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); - auto& soa2 = pti.GetStructOfArrays(); - amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data(); - amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data(); - amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data(); - const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data(); - const int* const ion_lev2 = soa2.GetIntData(PlasmaIdx::ion_lev).data(); + auto& ptile2 = pti.GetParticleTile(); + amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); + amrex::Real* const uy2 = ptile2.GetRealData(PlasmaIdx::uy_half_step).data(); + amrex::Real* const psi2= ptile2.GetRealData(PlasmaIdx::psi_half_step).data(); + const amrex::Real* const w2 = ptile2.GetRealData(PlasmaIdx::w).data(); + const int* const ion_lev2 = ptile2.GetIntData(PlasmaIdx::ion_lev).data(); PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); amrex::Real q2 = species2.GetCharge(); diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index cfd61aef13..5cbdcb7c81 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -121,14 +121,12 @@ struct PlasmaIdx /** \brief Container for particles of 1 plasma species. */ class PlasmaParticleContainer - : public amrex::ParticleContainerPureSoA + : public amrex::ParticleContainerRTSoA { public: /** Constructor */ explicit PlasmaParticleContainer (std::string name) : - amrex::ParticleContainerPureSoA(), + amrex::ParticleContainerRTSoA(), m_name(name) {} @@ -273,6 +271,7 @@ public: amrex::Real m_charge = 0; /**< charge of each particle of this species, per Ion level */ int m_n_subcycles = 1; /**< number of subcycles in the plasma particle push */ bool m_do_push = true; /**< whether the plasma pusher is enabled */ + int m_ab5_permutation = 0; // ionization: @@ -323,17 +322,17 @@ private: amrex::Vector m_insitu_sum_idata; /** Prefix/path for the output files */ std::string m_insitu_file_prefix = ""; + /** If the runtime attributes of the AMReX particle container were already set up */ + bool m_components_allocated = false; }; /** \brief Iterator over boxes in a particle container */ -class PlasmaParticleIterator : public amrex::ParIterSoA< - PlasmaIdx::real_nattribs,PlasmaIdx::int_nattribs, amrex::PolymorphicArenaAllocator> +class PlasmaParticleIterator : public amrex::ParIterRTSoA { public: /** Constructor */ PlasmaParticleIterator (ContainerType& pc) - : amrex::ParIterSoA(pc, 0, DfltMfi) {} + : amrex::ParIterRTSoA(pc, 0, DfltMfi) {} }; #endif diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 23f085edb3..c7e2ab6b52 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -182,9 +182,22 @@ PlasmaParticleContainer::ReadParameters () void PlasmaParticleContainer::InitData (const amrex::Vector& geom3d) { - SetArena(amrex::The_Arena()); - reserveData(); - resizeData(); + if (!m_components_allocated) { + m_components_allocated = true; + + SetArena(amrex::The_Arena()); + + for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { + AddRealComp(); + } + + for (int j = 0; j < PlasmaIdx::int_nattribs; ++j) { + AddIntComp(); + } + + reserveData(); + resizeData(); + } if (!m_read_fine_patch) { m_read_fine_patch = true; @@ -330,12 +343,12 @@ PlasmaParticleContainer::TagByLevel (const int current_N_level, for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) { - auto& soa = pti.GetStructOfArrays(); + auto& ptile = pti.GetParticleTile(); const amrex::Real * const AMREX_RESTRICT pos_x = to_prev ? - soa.GetRealData(PlasmaIdx::x_prev).data() : soa.GetRealData(PlasmaIdx::x).data(); + ptile.GetRealData(PlasmaIdx::x_prev).data() : ptile.GetRealData(PlasmaIdx::x).data(); const amrex::Real * const AMREX_RESTRICT pos_y = to_prev ? - soa.GetRealData(PlasmaIdx::y_prev).data() : soa.GetRealData(PlasmaIdx::y).data(); - auto * AMREX_RESTRICT idcpup = soa.GetIdCPUData().data(); + ptile.GetRealData(PlasmaIdx::y_prev).data() : ptile.GetRealData(PlasmaIdx::y).data(); + auto * AMREX_RESTRICT idcpup = ptile.GetIdCPUData().data(); const int lev1_idx = std::min(1, current_N_level-1); const int lev2_idx = std::min(2, current_N_level-1); diff --git a/src/particles/pusher/BeamParticleAdvance.cpp b/src/particles/pusher/BeamParticleAdvance.cpp index f1a81ad8d7..d280385358 100644 --- a/src/particles/pusher/BeamParticleAdvance.cpp +++ b/src/particles/pusher/BeamParticleAdvance.cpp @@ -137,9 +137,9 @@ AdvanceBeamParticlesSlice ( amrex::RealVect spin {0._rt, 0._rt, 0._rt}; if (spin_tracking) { - spin[0] = ptd.m_runtime_rdata[0][ip]; - spin[1] = ptd.m_runtime_rdata[1][ip]; - spin[2] = ptd.m_runtime_rdata[2][ip]; + spin[0] = ptd.rdata(BeamIdx::sx)[ip]; + spin[1] = ptd.rdata(BeamIdx::sy)[ip]; + spin[2] = ptd.rdata(BeamIdx::sz)[ip]; } for (; i < n_subcycles; i++) { @@ -331,9 +331,9 @@ AdvanceBeamParticlesSlice ( ptd.rdata(BeamIdx::uz)[ip] = uz; if (spin_tracking) { - ptd.m_runtime_rdata[0][ip] = spin[0]; - ptd.m_runtime_rdata[1][ip] = spin[1]; - ptd.m_runtime_rdata[2][ip] = spin[2]; + ptd.rdata(BeamIdx::sx)[ip] = spin[0]; + ptd.rdata(BeamIdx::sy)[ip] = spin[1]; + ptd.rdata(BeamIdx::sz)[ip] = spin[2]; } }); } diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 76f9786b91..3c5256c2a5 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -64,6 +64,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, const bool can_ionize = plasma.m_can_ionize; const int n_subcycles = plasma.m_n_subcycles; + [[maybe_unused]] const int ab5_permutation = plasma.m_ab5_permutation; const auto enforceBC = EnforceBC(); const amrex::Real dz = gm[0].CellSize(2) / n_subcycles; @@ -225,11 +226,11 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp, Aabssqp, AabssqDxp, AabssqDyp, q_mass_clight_ratio); - ptd.rdata(PlasmaIdx::Fx1)[ip] = ux * psi_inv; - ptd.rdata(PlasmaIdx::Fy1)[ip] = uy * psi_inv; - ptd.rdata(PlasmaIdx::Fux1)[ip] = dz_ux; - ptd.rdata(PlasmaIdx::Fuy1)[ip] = dz_uy; - ptd.rdata(PlasmaIdx::Fpsi1)[ip] = dz_psi; + ptd.rdata(PlasmaIdx::Fx1 + ab5_permutation)[ip] = ux * psi_inv; + ptd.rdata(PlasmaIdx::Fy1 + ab5_permutation)[ip] = uy * psi_inv; + ptd.rdata(PlasmaIdx::Fux1 + ab5_permutation)[ip] = dz_ux; + ptd.rdata(PlasmaIdx::Fuy1 + ab5_permutation)[ip] = dz_uy; + ptd.rdata(PlasmaIdx::Fpsi1 + ab5_permutation)[ip] = dz_psi; const amrex::Real ab5_coeffs[5] = { ( 1901._rt / 720._rt ) * dz, // a1 times dz @@ -241,11 +242,15 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, HIPACE_LOOP_UNROLL for (int iab=0; iab<5; ++iab) { - xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + iab)[ip]; - yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + iab)[ip]; - ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + iab)[ip]; - uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + iab)[ip]; - psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + iab)[ip]; + int p = ab5_permutation + iab; + if (p >= 5) { + p -= 5; + } + xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + p)[ip]; + yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + p)[ip]; + ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + p)[ip]; + uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + p)[ip]; + psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + p)[ip]; } if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; @@ -268,36 +273,11 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, #endif } // loop over subcycles }); + } #ifdef HIPACE_USE_AB5_PUSH - if (!temp_slice && lev == current_N_level - 1) { - auto& rd = pti.GetStructOfArrays().GetRealData(); - - // shift force terms - rd[PlasmaIdx::Fx5].swap(rd[PlasmaIdx::Fx4]); - rd[PlasmaIdx::Fy5].swap(rd[PlasmaIdx::Fy4]); - rd[PlasmaIdx::Fux5].swap(rd[PlasmaIdx::Fux4]); - rd[PlasmaIdx::Fuy5].swap(rd[PlasmaIdx::Fuy4]); - rd[PlasmaIdx::Fpsi5].swap(rd[PlasmaIdx::Fpsi4]); - - rd[PlasmaIdx::Fx4].swap(rd[PlasmaIdx::Fx3]); - rd[PlasmaIdx::Fy4].swap(rd[PlasmaIdx::Fy3]); - rd[PlasmaIdx::Fux4].swap(rd[PlasmaIdx::Fux3]); - rd[PlasmaIdx::Fuy4].swap(rd[PlasmaIdx::Fuy3]); - rd[PlasmaIdx::Fpsi4].swap(rd[PlasmaIdx::Fpsi3]); - - rd[PlasmaIdx::Fx3].swap(rd[PlasmaIdx::Fx2]); - rd[PlasmaIdx::Fy3].swap(rd[PlasmaIdx::Fy2]); - rd[PlasmaIdx::Fux3].swap(rd[PlasmaIdx::Fux2]); - rd[PlasmaIdx::Fuy3].swap(rd[PlasmaIdx::Fuy2]); - rd[PlasmaIdx::Fpsi3].swap(rd[PlasmaIdx::Fpsi2]); - - rd[PlasmaIdx::Fx2].swap(rd[PlasmaIdx::Fx1]); - rd[PlasmaIdx::Fy2].swap(rd[PlasmaIdx::Fy1]); - rd[PlasmaIdx::Fux2].swap(rd[PlasmaIdx::Fux1]); - rd[PlasmaIdx::Fuy2].swap(rd[PlasmaIdx::Fuy1]); - rd[PlasmaIdx::Fpsi2].swap(rd[PlasmaIdx::Fpsi1]); - } -#endif + if (!temp_slice && lev == current_N_level - 1) { + plasma.m_ab5_permutation = (plasma.m_ab5_permutation + 4) % 5; } +#endif } diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index 05bb6f323d..abaefe8b75 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -402,9 +402,9 @@ SalameMultiplyBeamWeight (const amrex::Real W, Hipace* hipace) if (!beam.m_do_salame) continue; // For id and weights - auto& soa = beam.getBeamSlice(WhichBeamSlice::This).GetStructOfArrays(); - amrex::Real * const wp = soa.GetRealData(BeamIdx::w).data(); - auto * const idcpup = soa.GetIdCPUData().data(); + auto& slice = beam.getBeamSlice(WhichBeamSlice::This); + amrex::Real * const wp = slice.GetRealData(BeamIdx::w).data(); + auto * const idcpup = slice.GetIdCPUData().data(); amrex::ParallelFor( beam.getNumParticles(WhichBeamSlice::This), diff --git a/src/utils/AdaptiveTimeStep.cpp b/src/utils/AdaptiveTimeStep.cpp index 480cd001dd..ffa21ff91d 100644 --- a/src/utils/AdaptiveTimeStep.cpp +++ b/src/utils/AdaptiveTimeStep.cpp @@ -115,17 +115,17 @@ AdaptiveTimeStep::GatherMinUzSlice (MultiBeam& beams, const bool initial) // Extract particle properties // For momenta and weights if (initial) { - const auto& soa = beam.getBeamInitSlice().GetStructOfArrays(); + const auto& slice = beam.getBeamInitSlice(); num_particles = beam.getBeamInitSlice().size(); - uzp = soa.GetRealData(BeamIdx::uz).data(); - wp = soa.GetRealData(BeamIdx::w).data(); - idcpup = soa.GetIdCPUData().data(); + uzp = slice.GetRealData(BeamIdx::uz).data(); + wp = slice.GetRealData(BeamIdx::w).data(); + idcpup = slice.GetIdCPUData().data(); } else { - const auto& soa = beam.getBeamSlice(WhichBeamSlice::This).GetStructOfArrays(); + const auto& slice = beam.getBeamSlice(WhichBeamSlice::This); num_particles = beam.getNumParticles(WhichBeamSlice::This); - uzp = soa.GetRealData(BeamIdx::uz).data(); - wp = soa.GetRealData(BeamIdx::w).data(); - idcpup = soa.GetIdCPUData().data(); + uzp = slice.GetRealData(BeamIdx::uz).data(); + wp = slice.GetRealData(BeamIdx::w).data(); + idcpup = slice.GetIdCPUData().data(); } amrex::ReduceOps ReduceTuple diff --git a/src/utils/MultiBuffer.cpp b/src/utils/MultiBuffer.cpp index eaae5addcf..36678558c9 100644 --- a/src/utils/MultiBuffer.cpp +++ b/src/utils/MultiBuffer.cpp @@ -803,12 +803,12 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int for (int b = 0; b < m_nbeams; ++b) { auto& beam = beams.getBeam(b); const int num_particles = beam.getNumParticles(beam_slice); - auto& soa = beam.getBeamSlice(beam_slice).GetStructOfArrays(); + auto& ptile = beam.getBeamSlice(beam_slice); if (beam.communicateIdCpuComponent()) { // only pack idcpu component if it should be communicated memcpy_to_buffer(slice, bo.m_beam_idcpu[b].value(), // NOLINT(bugprone-unchecked-optional-access) - soa.GetIdCPUData().dataPtr(), + ptile.GetIdCPUData().dataPtr(), num_particles * sizeof(std::uint64_t)); } @@ -816,7 +816,7 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int // only pack real component if it should be communicated if (beam.communicateRealComponent(rcomp)) { memcpy_to_buffer(slice, bo.m_beam_real[b].at(rcomp), - soa.GetRealData(rcomp).dataPtr(), + ptile.GetRealData(rcomp).dataPtr(), num_particles * sizeof(amrex::Real)); } } @@ -825,7 +825,7 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int // only pack int component if it should be communicated if (beam.communicateIntComponent(icomp)) { memcpy_to_buffer(slice, bo.m_beam_int[b].at(icomp), - soa.GetIntData(icomp).dataPtr(), + ptile.GetIntData(icomp).dataPtr(), num_particles * sizeof(int)); } } @@ -857,16 +857,16 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i auto& beam = beams.getBeam(b); const int num_particles = get_metadata_location(slice)[b + 1]; beam.resize(beam_slice, num_particles, 0); - auto& soa = beam.getBeamSlice(beam_slice).GetStructOfArrays(); + auto& ptile = beam.getBeamSlice(beam_slice); if (beam.communicateIdCpuComponent()) { // only undpack idcpu component if it should be communicated memcpy_from_buffer(slice, bo.m_beam_idcpu[b].value(), // NOLINT(bugprone-unchecked-optional-access) - soa.GetIdCPUData().dataPtr(), + ptile.GetIdCPUData().dataPtr(), num_particles * sizeof(std::uint64_t)); } else { // if idcpu is not communicated, then we need to initialize it here - std::uint64_t* data_ptr = soa.GetIdCPUData().dataPtr(); + std::uint64_t* data_ptr = ptile.GetIdCPUData().dataPtr(); amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (int i) { amrex::ParticleIDWrapper{data_ptr[i]} = 1; amrex::ParticleCPUWrapper{data_ptr[i]} = 0; @@ -877,11 +877,11 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i if (beam.communicateRealComponent(rcomp)) { // only unpack real component if it should be communicated memcpy_from_buffer(slice, bo.m_beam_real[b].at(rcomp), - soa.GetRealData(rcomp).dataPtr(), + ptile.GetRealData(rcomp).dataPtr(), num_particles * sizeof(amrex::Real)); } else { // initialize per-slice-only real components to zero - amrex::Real* data_ptr = soa.GetRealData(rcomp).dataPtr(); + amrex::Real* data_ptr = ptile.GetRealData(rcomp).dataPtr(); amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (int i) { data_ptr[i] = amrex::Real(0.); }); @@ -892,11 +892,11 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i if (beam.communicateIntComponent(icomp)) { // only unpack int component if it should be communicated memcpy_from_buffer(slice, bo.m_beam_int[b].at(icomp), - soa.GetIntData(icomp).dataPtr(), + ptile.GetIntData(icomp).dataPtr(), num_particles * sizeof(int)); } else { // initialize per-slice-only int components to zero - int* data_ptr = soa.GetIntData(icomp).dataPtr(); + int* data_ptr = ptile.GetIntData(icomp).dataPtr(); amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (int i) { data_ptr[i] = 0; });