From 230a1a721ea630c54ca16e5967abd05cb09f5d01 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 22 Apr 2025 22:14:10 +0200 Subject: [PATCH 1/9] Test runtime only particle container --- cmake/dependencies/AMReX.cmake | 6 +- src/particles/beam/BeamParticleContainer.H | 27 +----- src/particles/beam/BeamParticleContainer.cpp | 97 ++++++++++--------- .../beam/BeamParticleContainerInit.cpp | 2 +- .../plasma/PlasmaParticleContainer.H | 8 +- .../plasma/PlasmaParticleContainer.cpp | 86 ++++++++-------- .../plasma/PlasmaParticleContainerInit.cpp | 6 +- src/particles/pusher/BeamParticleAdvance.cpp | 12 +-- 8 files changed, 119 insertions(+), 125 deletions(-) diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 321588741b..62d7e9fe5f 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -126,7 +126,7 @@ macro(find_amrex) mark_as_advanced(AMReX_TP_PROFILE) mark_as_advanced(USE_XSDK_DEFAULTS) - message(STATUS "AMReX: Using version '${AMREX_PKG_VERSION}' (${AMREX_GIT_VERSION})") + message(STATUS "AMReX: Using version '${AMREX_PKG_VERSION}' (${AMREX_GIT_VERSION})") else() message(STATUS "Searching for pre-installed AMReX ...") set(COMPONENT_PRECISION ${HiPACE_PRECISION} P${HiPACE_PRECISION}) @@ -146,10 +146,10 @@ set(HiPACE_amrex_src "" "Local path to AMReX source directory (preferred if set)") # Git fetcher -set(HiPACE_amrex_repo "https://github.com/AMReX-Codes/amrex.git" +set(HiPACE_amrex_repo "https://github.com/AlexanderSinn/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(HiPACE_amrex_internal)") -set(HiPACE_amrex_branch "development" +set(HiPACE_amrex_branch "Add_simpler_version_of_ParticleTile_using_2D_array" CACHE STRING "Repository branch for HiPACE_amrex_repo if(HiPACE_amrex_internal)") diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H index 1bf5630439..b59f58ff4d 100644 --- a/src/particles/beam/BeamParticleContainer.H +++ b/src/particles/beam/BeamParticleContainer.H @@ -27,7 +27,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 @@ -44,25 +45,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 - >; - -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::ParticleTile2; /** \brief Container for particles of 1 beam species. */ class BeamParticleContainer @@ -182,7 +165,7 @@ public: void resize (int which_slice, int num_particles, int num_slipped_particles); - BeamTileInit& getBeamInitSlice () { + BeamTile& getBeamInitSlice () { return m_init_slice; } @@ -212,7 +195,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 ddac67c776..39547f1408 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -92,25 +92,24 @@ BeamParticleContainer::ReadParameters () "Tilted beams and correlated energy spreads are only implemented for fixed weight beams"); } queryWithParserAlt(pp, "initialize_on_cpu", m_initialize_on_cpu, pp_alt); - auto& soa = getBeamInitSlice().GetStructOfArrays(); - soa.GetIdCPUData().setArena( - m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena()); - for (int rcomp = 0; rcomp < soa.NumRealComps(); ++rcomp) { - soa.GetRealData()[rcomp].setArena( - m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena()); - } - for (int icomp = 0; icomp < soa.NumIntComps(); ++icomp) { - soa.GetIntData()[icomp].setArena( - m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena()); - } queryWithParserAlt(pp, "do_spin_tracking", m_do_spin_tracking, pp_alt); if (m_do_spin_tracking) { getWithParserAlt(pp, "initial_spin", m_initial_spin, pp_alt); queryWithParserAlt(pp, "spin_anom", m_spin_anom, pp_alt); - for (auto& beam_tile : m_slices) { - // Use 3 real and 0 int runtime components - beam_tile.define(3, 0); - } + } + + getBeamInitSlice().define( + m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena(), + BeamIdx::real_nattribs_in_buffer, + BeamIdx::int_nattribs_in_buffer + ); + + for (auto& beam_tile : m_slices) { + beam_tile.define( + amrex::The_Arena(), + BeamIdx::real_nattribs + (m_do_spin_tracking ? 3 : 0), + BeamIdx::int_nattribs + ); } } @@ -396,9 +395,9 @@ BeamParticleContainer::initializeSlice (int slice, int which_slice) { 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]; } ); } @@ -431,48 +430,54 @@ BeamParticleContainer::ReorderParticles (int beam_slice, int step, amrex::Geomet 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(np_total); + amrex::Gpu::AsyncVector tmp_idcpu(np_total); - auto src = soa.GetIdCPUData().data(); + auto src = ptile.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]; + dst[i] = src[permutations[i]]; + }); + amrex::ParallelFor(np_total, + [=] AMREX_GPU_DEVICE (int i) { + src[i] = dst[i]; }); - - amrex::Gpu::streamSynchronize(); - soa.GetIdCPUData().swap(tmp_idcpu); } - { // Create a scope for the temporary vector below - BeamTile::RealVector tmp_real(np_total); - for (int comp = 0; comp < soa.NumRealComps(); ++comp) { - auto src = soa.GetRealData(comp).data(); + { + amrex::Gpu::AsyncVector tmp_real(np_total); + + for (int comp = 0; comp < ptile.NumRealComps(); ++comp) { + auto src = ptile.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]; + dst[i] = src[permutations[i]]; + }); + amrex::ParallelFor(np_total, + [=] AMREX_GPU_DEVICE (int i) { + src[i] = dst[i]; }); - - amrex::Gpu::streamSynchronize(); - soa.GetRealData(comp).swap(tmp_real); } } - BeamTile::IntVector tmp_int(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::AsyncVector tmp_int(np_total); - amrex::Gpu::streamSynchronize(); - soa.GetIntData(comp).swap(tmp_int); + for (int comp = 0; comp < ptile.NumIntComps(); ++comp) { + auto src = ptile.GetIntData(comp).data(); + int* dst = tmp_int.data(); + amrex::ParallelFor(np_total, + [=] AMREX_GPU_DEVICE (int i) { + dst[i] = src[permutations[i]]; + }); + amrex::ParallelFor(np_total, + [=] AMREX_GPU_DEVICE (int i) { + src[i] = dst[i]; + }); + } } } } @@ -570,9 +575,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 98d42cd6f2..e249fb705e 100644 --- a/src/particles/beam/BeamParticleContainerInit.cpp +++ b/src/particles/beam/BeamParticleContainerInit.cpp @@ -42,7 +42,7 @@ namespace */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddOneBeamParticle ( - const BeamTileInit::ParticleTileDataType& ptd, const amrex::Real& x, + 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& weight, const amrex::Long pid, const amrex::Long ip, const amrex::Real& speed_of_light, const EnforceBC& enforceBC) noexcept diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index b102f3ba22..95e6df8019 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -47,12 +47,12 @@ struct PlasmaIdx /** \brief Container for particles of 1 plasma species. */ class PlasmaParticleContainer - : public amrex::ParticleContainerPureSoA + : public amrex::ParticleContainerPureSoA2 { public: /** Constructor */ explicit PlasmaParticleContainer (std::string name) : - amrex::ParticleContainerPureSoA(), + amrex::ParticleContainerPureSoA2(), m_name(name) { ReadParameters(); @@ -250,12 +250,12 @@ private: }; /** \brief Iterator over boxes in a particle container */ -class PlasmaParticleIterator : public amrex::ParIterSoA +class PlasmaParticleIterator : public amrex::ParIterSoA2 { public: /** Constructor */ PlasmaParticleIterator (ContainerType& pc) - : amrex::ParIterSoA(pc, 0, DfltMfi) {} + : amrex::ParIterSoA2(pc, 0, DfltMfi) {} }; #endif diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 69a06b49b9..9ba78c3efa 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -179,6 +179,16 @@ PlasmaParticleContainer::ReadParameters () void PlasmaParticleContainer::InitData (const amrex::Geometry& geom) { + 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(); @@ -408,10 +418,8 @@ IonizationModule (const int lev, ptile_elec.resize(new_size); // Load electron soa and aos after resize - auto arrdata_ion = ptile_ion.GetStructOfArrays().realarray(); - auto arrdata_elec = ptile_elec.GetStructOfArrays().realarray(); - auto int_arrdata_elec = ptile_elec.GetStructOfArrays().intarray(); - auto idcpu_elec = ptile_elec.GetStructOfArrays().GetIdCPUData().data(); + auto ptd_ion = ptile_ion.getParticleTileData(); + auto ptd_elec = ptile_elec.getParticleTileData(); const int init_ion_lev = m_product_pc->m_init_ion_lev; @@ -428,30 +436,30 @@ IonizationModule (const int lev, // Copy ion data to new electron // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons - amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; - amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = lev; // current level - arrdata_elec[PlasmaIdx::x ][pidx] = arrdata_ion[PlasmaIdx::x ][ip]; - arrdata_elec[PlasmaIdx::y ][pidx] = arrdata_ion[PlasmaIdx::y ][ip]; - - arrdata_elec[PlasmaIdx::w ][pidx] = arrdata_ion[PlasmaIdx::w ][ip]; - arrdata_elec[PlasmaIdx::ux ][pidx] = 0._rt; - arrdata_elec[PlasmaIdx::uy ][pidx] = 0._rt; + ptd_elec.id(pidx) = 2; + ptd_elec.cpu(pidx) = lev; // current level + ptd_elec.rdata(PlasmaIdx::x )[pidx] = ptd_ion.rdata(PlasmaIdx::x)[ip]; + ptd_elec.rdata(PlasmaIdx::y )[pidx] = ptd_ion.rdata(PlasmaIdx::y)[ip]; + + ptd_elec.rdata(PlasmaIdx::w )[pidx] = ptd_ion.rdata(PlasmaIdx::w)[ip]; + ptd_elec.rdata(PlasmaIdx::ux )[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::uy )[pidx] = 0._rt; // Later we could consider adding a finite temperature to the ionized electrons - arrdata_elec[PlasmaIdx::psi ][pidx] = 1._rt; - arrdata_elec[PlasmaIdx::x_prev ][pidx] = arrdata_ion[PlasmaIdx::x_prev][ip]; - arrdata_elec[PlasmaIdx::y_prev ][pidx] = arrdata_ion[PlasmaIdx::y_prev][ip]; - arrdata_elec[PlasmaIdx::ux_half_step ][pidx] = 0._rt; - arrdata_elec[PlasmaIdx::uy_half_step ][pidx] = 0._rt; - arrdata_elec[PlasmaIdx::psi_half_step][pidx] = 1._rt; + ptd_elec.rdata(PlasmaIdx::psi )[pidx] = 1._rt; + ptd_elec.rdata(PlasmaIdx::x_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::ux_half_step )[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = 1._rt; #ifdef HIPACE_USE_AB5_PUSH #ifdef AMREX_USE_GPU #pragma unroll #endif for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - arrdata_elec[iforce][pidx] = 0._rt; + ptd_elec.rdata(iforce)[pidx] = 0._rt; } #endif - int_arrdata_elec[PlasmaIdx::ion_lev][pidx] = init_ion_lev; + ptd_elec.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } }); @@ -616,11 +624,8 @@ LaserIonization (const int islice, ptile_elec.resize(new_size); // Load electron soa and aos after resize - auto arrdata_ion = ptile_ion.GetStructOfArrays().realarray(); - auto arrdata_elec = ptile_elec.GetStructOfArrays().realarray(); - auto int_arrdata_elec = ptile_elec.GetStructOfArrays().intarray(); - auto idcpu_elec = ptile_elec.GetStructOfArrays().GetIdCPUData().data(); - auto idcpu_ion = ptile_ion.GetStructOfArrays().GetIdCPUData().data(); + auto ptd_ion = ptile_ion.getParticleTileData(); + auto ptd_elec = ptile_elec.getParticleTileData(); const int init_ion_lev = m_product_pc->m_init_ion_lev; @@ -699,29 +704,28 @@ LaserIonization (const int islice, const long pidx = pid + old_size; // Copy ion data to new electron // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons - amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; - amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = - amrex::ParticleCPUWrapper{idcpu_ion[pidx]}; // current level - arrdata_elec[PlasmaIdx::x ][pidx] = arrdata_ion[PlasmaIdx::x ][ip]; - arrdata_elec[PlasmaIdx::y ][pidx] = arrdata_ion[PlasmaIdx::y ][ip]; - arrdata_elec[PlasmaIdx::w ][pidx] = arrdata_ion[PlasmaIdx::w ][ip]; - arrdata_elec[PlasmaIdx::ux ][pidx] = ux * phys_const.c; - arrdata_elec[PlasmaIdx::uy ][pidx] = uy * phys_const.c; - arrdata_elec[PlasmaIdx::psi ][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz; //psi = gamma - uz - arrdata_elec[PlasmaIdx::x_prev ][pidx] = arrdata_ion[PlasmaIdx::x_prev][ip]; - arrdata_elec[PlasmaIdx::y_prev ][pidx] = arrdata_ion[PlasmaIdx::y_prev][ip]; - arrdata_elec[PlasmaIdx::ux_half_step ][pidx] = ux * phys_const.c; - arrdata_elec[PlasmaIdx::uy_half_step ][pidx] = uy * phys_const.c; - arrdata_elec[PlasmaIdx::psi_half_step][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz; + ptd_elec.id(pidx) = 2; // current level + ptd_elec.cpu(pidx) = ptd_ion.cpu(ip); + ptd_elec.rdata(PlasmaIdx::x )[pidx] = ptd_ion.rdata(PlasmaIdx::x )[ip]; + ptd_elec.rdata(PlasmaIdx::y )[pidx] = ptd_ion.rdata(PlasmaIdx::y )[ip]; + ptd_elec.rdata(PlasmaIdx::w )[pidx] = ptd_ion.rdata(PlasmaIdx::w )[ip]; + ptd_elec.rdata(PlasmaIdx::ux )[pidx] = ux * phys_const.c; + ptd_elec.rdata(PlasmaIdx::uy )[pidx] = uy * phys_const.c; + ptd_elec.rdata(PlasmaIdx::psi )[pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz; //psi = gamma - uz + ptd_elec.rdata(PlasmaIdx::x_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::ux_half_step )[pidx] = ux * phys_const.c; + ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = uy * phys_const.c; + ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz; #ifdef HIPACE_USE_AB5_PUSH #ifdef AMREX_USE_GPU #pragma unroll #endif for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - arrdata_elec[iforce][pidx] = 0._rt; + ptd_elec.rdata(iforce)[pidx] = 0._rt; } #endif - int_arrdata_elec[PlasmaIdx::ion_lev][pidx] = init_ion_lev; + ptd_elec.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } }); diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index 54b7140b0a..bac34f900f 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -176,8 +176,10 @@ InitParticles (const amrex::RealVect& a_u_std, scale_fac_fine /= 4.; } - auto& particles = GetParticles(lev); - auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + // auto& particles = GetParticles(lev); + // auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + + auto& particle_tile = DefineAndReturnParticleTile(0, mfi); auto old_size = particle_tile.size(); const auto new_size = old_size + total_num_particles; diff --git a/src/particles/pusher/BeamParticleAdvance.cpp b/src/particles/pusher/BeamParticleAdvance.cpp index 858f4a9099..023cf4662b 100644 --- a/src/particles/pusher/BeamParticleAdvance.cpp +++ b/src/particles/pusher/BeamParticleAdvance.cpp @@ -140,9 +140,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++) { @@ -328,9 +328,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]; } }); } From 00455fa404f24926effb065762b0e04951769665 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 22 Apr 2025 22:33:09 +0200 Subject: [PATCH 2/9] fix ab5 --- .../plasma/PlasmaParticleContainer.H | 1 + .../pusher/PlasmaParticleAdvance.cpp | 58 ++++++------------- 2 files changed, 20 insertions(+), 39 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 95e6df8019..6da0bff0df 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -196,6 +196,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: diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index a3bd479152..4cf2b1496e 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, clight_inv, q_mass_clight_ratio); - ptd.rdata(PlasmaIdx::Fx1)[ip] = clight_inv*(ux * psi_inv); - ptd.rdata(PlasmaIdx::Fy1)[ip] = clight_inv*(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] = clight_inv*(ux * psi_inv); + ptd.rdata(PlasmaIdx::Fy1 + ab5_permutation)[ip] = clight_inv*(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 @@ -243,11 +244,15 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, #pragma unroll #endif 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; @@ -270,36 +275,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 } From fe8c4ecc88c9ad741af9f6fbcb7421ec872d131e Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 22 Apr 2025 23:30:06 +0200 Subject: [PATCH 3/9] clean GetStructOfArrays --- src/diagnostics/OpenPMDWriter.cpp | 12 ++-- src/particles/beam/BeamParticleContainer.cpp | 10 ++-- src/particles/collisions/CoulombCollision.cpp | 57 +++++++++---------- .../plasma/PlasmaParticleContainer.cpp | 44 +++++++------- .../plasma/PlasmaParticleContainerInit.cpp | 3 - src/salame/Salame.cpp | 6 +- src/utils/AdaptiveTimeStep.cpp | 24 ++++---- src/utils/MultiBuffer.cpp | 22 +++---- 8 files changed, 85 insertions(+), 93 deletions(-) diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index fe778a28fe..161f4dd2b8 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -282,7 +282,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.cpp b/src/particles/beam/BeamParticleContainer.cpp index 39547f1408..96ad5eb167 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -267,7 +267,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 @@ -319,10 +319,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); diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index b67dc3557c..2ad257a6d8 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -89,12 +89,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(); @@ -155,12 +155,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(); @@ -169,12 +169,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(); @@ -268,11 +267,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(); @@ -281,12 +280,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.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 9ba78c3efa..e53d7c76aa 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -241,12 +241,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); @@ -315,8 +315,6 @@ IonizationModule (const int lev, mfi_ion.index(), mfi_ion.LocalTileIndex()); auto& ptile_ion = plevel_ion.at(index); - auto& soa_ion = ptile_ion.GetStructOfArrays(); // For momenta and weights - const amrex::Real clightsq = 1.0_rt / ( phys_const.c * phys_const.c ); // Calculation of E0 in SI units for denormalization const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * @@ -325,13 +323,13 @@ IonizationModule (const int lev, const amrex::Real E0 = Hipace::m_normalized_units ? wp * PhysConstSI::m_e * PhysConstSI::c / PhysConstSI::q_e : 1; - int * const ion_lev = soa_ion.GetIntData(PlasmaIdx::ion_lev).data(); - const amrex::Real * const x_prev = soa_ion.GetRealData(PlasmaIdx::x_prev).data(); - const amrex::Real * const y_prev = soa_ion.GetRealData(PlasmaIdx::y_prev).data(); - const amrex::Real * const uxp = soa_ion.GetRealData(PlasmaIdx::ux_half_step).data(); - const amrex::Real * const uyp = soa_ion.GetRealData(PlasmaIdx::uy_half_step).data(); - const amrex::Real * const psip =soa_ion.GetRealData(PlasmaIdx::psi_half_step).data(); - const auto * idcpup = soa_ion.GetIdCPUData().data(); + int * const ion_lev = ptile_ion.GetIntData(PlasmaIdx::ion_lev).data(); + const amrex::Real * const x_prev = ptile_ion.GetRealData(PlasmaIdx::x_prev).data(); + const amrex::Real * const y_prev = ptile_ion.GetRealData(PlasmaIdx::y_prev).data(); + const amrex::Real * const uxp = ptile_ion.GetRealData(PlasmaIdx::ux_half_step).data(); + const amrex::Real * const uyp = ptile_ion.GetRealData(PlasmaIdx::uy_half_step).data(); + const amrex::Real * const psip =ptile_ion.GetRealData(PlasmaIdx::psi_half_step).data(); + const auto * idcpup = ptile_ion.GetIdCPUData().data(); // Make Ion Mask and load ADK prefactors // Ion Mask is necessary to only resize electron particle tile once @@ -417,7 +415,7 @@ IonizationModule (const int lev, const auto new_size = old_size + num_new_electrons.dataValue(); ptile_elec.resize(new_size); - // Load electron soa and aos after resize + // Load electron after resize auto ptd_ion = ptile_ion.getParticleTileData(); auto ptd_elec = ptile_elec.getParticleTileData(); @@ -510,8 +508,6 @@ LaserIonization (const int islice, mfi_ion.index(), mfi_ion.LocalTileIndex()); auto& ptile_ion = plevel_ion.at(index); - auto& soa_ion = ptile_ion.GetStructOfArrays(); // for momenta and weights - const amrex::Real clightsq = 1.0_rt / ( phys_const.c * phys_const.c ); // Calcuation of E0 in SI units for denormalization const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * @@ -523,13 +519,13 @@ LaserIonization (const int islice, const amrex::Real omega0 = 2.0 * MathConst::pi * phys_const.c / lambda0; const bool linear_polarization = laser.LinearPolarization(); - int * const ion_lev = soa_ion.GetIntData(PlasmaIdx::ion_lev).data(); - const amrex::Real * const x_prev = soa_ion.GetRealData(PlasmaIdx::x_prev).data(); - const amrex::Real * const y_prev = soa_ion.GetRealData(PlasmaIdx::y_prev).data(); - const amrex::Real * const uxp = soa_ion.GetRealData(PlasmaIdx::ux_half_step).data(); - const amrex::Real * const uyp = soa_ion.GetRealData(PlasmaIdx::uy_half_step).data(); - const amrex::Real * const psip =soa_ion.GetRealData(PlasmaIdx::psi_half_step).data(); - const auto * idcpup = soa_ion.GetIdCPUData().data(); + int * const ion_lev = ptile_ion.GetIntData(PlasmaIdx::ion_lev).data(); + const amrex::Real * const x_prev = ptile_ion.GetRealData(PlasmaIdx::x_prev).data(); + const amrex::Real * const y_prev = ptile_ion.GetRealData(PlasmaIdx::y_prev).data(); + const amrex::Real * const uxp = ptile_ion.GetRealData(PlasmaIdx::ux_half_step).data(); + const amrex::Real * const uyp = ptile_ion.GetRealData(PlasmaIdx::uy_half_step).data(); + const amrex::Real * const psip =ptile_ion.GetRealData(PlasmaIdx::psi_half_step).data(); + const auto * idcpup = ptile_ion.GetIdCPUData().data(); // Make Ion Mask and load ADK prefactors // Ion Mask is necessary to only resize electron particle tile once @@ -623,7 +619,7 @@ LaserIonization (const int islice, const auto new_size = old_size + num_new_electrons.dataValue(); ptile_elec.resize(new_size); - // Load electron soa and aos after resize + // Load electron after resize auto ptd_ion = ptile_ion.getParticleTileData(); auto ptd_elec = ptile_elec.getParticleTileData(); diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index bac34f900f..121f20b36c 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -176,9 +176,6 @@ InitParticles (const amrex::RealVect& a_u_std, scale_fac_fine /= 4.; } - // auto& particles = GetParticles(lev); - // auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; - auto& particle_tile = DefineAndReturnParticleTile(0, mfi); auto old_size = particle_tile.size(); diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index 909b40981f..161d21c45f 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -414,9 +414,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 9e9ff57779..6d81e1ea72 100644 --- a/src/utils/AdaptiveTimeStep.cpp +++ b/src/utils/AdaptiveTimeStep.cpp @@ -117,17 +117,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 35e6d9b753..b95abe7d85 100644 --- a/src/utils/MultiBuffer.cpp +++ b/src/utils/MultiBuffer.cpp @@ -807,13 +807,13 @@ 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, get_buffer_offset(slice, offset_type::beam_idcpu, beams, laser, b, 0), - soa.GetIdCPUData().dataPtr(), + ptile.GetIdCPUData().dataPtr(), num_particles * sizeof(std::uint64_t)); } @@ -822,7 +822,7 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int if (beam.communicateRealComponent(rcomp)) { memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, beams, laser, b, rcomp), - soa.GetRealData(rcomp).dataPtr(), + ptile.GetRealData(rcomp).dataPtr(), num_particles * sizeof(amrex::Real)); } } @@ -832,7 +832,7 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int if (beam.communicateIntComponent(icomp)) { memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, beams, laser, b, icomp), - soa.GetIntData(icomp).dataPtr(), + ptile.GetIntData(icomp).dataPtr(), num_particles * sizeof(int)); } } @@ -861,17 +861,17 @@ 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, get_buffer_offset(slice, offset_type::beam_idcpu, beams, laser, b, 0), - 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; @@ -883,11 +883,11 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i // only unpack real component if it should be communicated memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, beams, laser, b, 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.); }); @@ -899,11 +899,11 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i // only unpack int component if it should be communicated memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, beams, laser, b, 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; }); From 25bd871aacb995c6c144081075f0cf90e712bbea Mon Sep 17 00:00:00 2001 From: Alexander Sinn <64009254+AlexanderSinn@users.noreply.github.com> Date: Sat, 17 May 2025 14:28:29 +0200 Subject: [PATCH 4/9] fix spin data --- src/particles/beam/BeamParticleContainer.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index ca4eb93411..f06f776525 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -383,9 +383,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; From 6b1410be17af425ad33d7911745ce3dad05628da Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sat, 26 Jul 2025 20:34:22 +0200 Subject: [PATCH 5/9] fix define --- src/particles/beam/BeamParticleContainer.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index 6590014e13..f61929f69b 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -101,16 +101,20 @@ BeamParticleContainer::ReadParameters () } getBeamInitSlice().define( - m_initialize_on_cpu ? amrex::The_Pinned_Arena() : amrex::The_Arena(), BeamIdx::real_nattribs_in_buffer + (m_do_spin_tracking ? 3 : 0), - BeamIdx::int_nattribs_in_buffer + 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( - amrex::The_Arena(), BeamIdx::real_nattribs + (m_do_spin_tracking ? 3 : 0), - BeamIdx::int_nattribs + BeamIdx::int_nattribs, + nullptr, + nullptr, + amrex::The_Arena() ); } } From 12b2b4f4cebcff6eb85721d18cb38a09c9d26b6b Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sat, 26 Jul 2025 21:05:21 +0200 Subject: [PATCH 6/9] change name --- src/particles/beam/BeamParticleContainer.H | 2 +- src/particles/plasma/PlasmaParticleContainer.H | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H index 8e44651192..ada855090d 100644 --- a/src/particles/beam/BeamParticleContainer.H +++ b/src/particles/beam/BeamParticleContainer.H @@ -45,7 +45,7 @@ struct WhichBeamSlice { enum beam_slice : int { Next=0, This, N }; }; -using BeamTile = amrex::ParticleTile2; +using BeamTile = amrex::ParticleTileRT; /** \brief Container for particles of 1 beam species. */ class BeamParticleContainer diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index f8772f5d87..774cc67c63 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -47,12 +47,12 @@ struct PlasmaIdx /** \brief Container for particles of 1 plasma species. */ class PlasmaParticleContainer - : public amrex::ParticleContainerPureSoA2 + : public amrex::ParticleContainerRTSoA { public: /** Constructor */ explicit PlasmaParticleContainer (std::string name) : - amrex::ParticleContainerPureSoA2(), + amrex::ParticleContainerRTSoA(), m_name(name) { ReadParameters(); @@ -251,12 +251,12 @@ private: }; /** \brief Iterator over boxes in a particle container */ -class PlasmaParticleIterator : public amrex::ParIterSoA2 +class PlasmaParticleIterator : public amrex::ParIterRTSoA { public: /** Constructor */ PlasmaParticleIterator (ContainerType& pc) - : amrex::ParIterSoA2(pc, 0, DfltMfi) {} + : amrex::ParIterRTSoA(pc, 0, DfltMfi) {} }; #endif From b16a6a17ab6d4ddb79cc5677dfd23cacd96cbf4b Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sat, 26 Jul 2025 21:29:44 +0200 Subject: [PATCH 7/9] use amrex ReorderParticles --- src/particles/beam/BeamParticleContainer.cpp | 52 +------------------- 1 file changed, 1 insertion(+), 51 deletions(-) diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index f61929f69b..7c394b6251 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -434,61 +434,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(); - - { - amrex::Gpu::AsyncVector tmp_idcpu(np_total); - - auto src = ptile.GetIdCPUData().data(); - uint64_t* dst = tmp_idcpu.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = src[permutations[i]]; - }); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - src[i] = dst[i]; - }); - } - - { - amrex::Gpu::AsyncVector tmp_real(np_total); - - for (int comp = 0; comp < ptile.NumRealComps(); ++comp) { - auto src = ptile.GetRealData(comp).data(); - amrex::ParticleReal* dst = tmp_real.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = src[permutations[i]]; - }); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - src[i] = dst[i]; - }); - } - } - - { - amrex::Gpu::AsyncVector tmp_int(np_total); - - for (int comp = 0; comp < ptile.NumIntComps(); ++comp) { - auto src = ptile.GetIntData(comp).data(); - int* dst = tmp_int.data(); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - dst[i] = src[permutations[i]]; - }); - amrex::ParallelFor(np_total, - [=] AMREX_GPU_DEVICE (int i) { - src[i] = dst[i]; - }); - } - } + amrex::ReorderParticles(ptile, perm.dataPtr()); } } From a0a00f251950f5135020682f3ca7f551b95a0ed0 Mon Sep 17 00:00:00 2001 From: Alexander Sinn <64009254+AlexanderSinn@users.noreply.github.com> Date: Fri, 3 Apr 2026 12:23:58 +0200 Subject: [PATCH 8/9] Use AMReX development --- cmake/dependencies/AMReX.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index b478cf2a21..5b47ed027b 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -146,10 +146,10 @@ set(HiPACE_amrex_src "" "Local path to AMReX source directory (preferred if set)") # Git fetcher -set(HiPACE_amrex_repo "https://github.com/AlexanderSinn/amrex.git" +set(HiPACE_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(HiPACE_amrex_internal)") -set(HiPACE_amrex_branch "Add_simpler_version_of_ParticleTile_using_2D_array" +set(HiPACE_amrex_branch "development" CACHE STRING "Repository branch for HiPACE_amrex_repo if(HiPACE_amrex_internal)") From ef95927e616a971d0efc7dd66c29e6618d27585f Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Mon, 20 Apr 2026 18:34:47 +0200 Subject: [PATCH 9/9] fix memory leak --- .../plasma/PlasmaParticleContainer.H | 2 ++ .../plasma/PlasmaParticleContainer.cpp | 22 +++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 971b42e121..5cbdcb7c81 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -322,6 +322,8 @@ 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 */ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 88179a63a9..c7e2ab6b52 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -182,18 +182,22 @@ PlasmaParticleContainer::ReadParameters () void PlasmaParticleContainer::InitData (const amrex::Vector& geom3d) { - SetArena(amrex::The_Arena()); + if (!m_components_allocated) { + m_components_allocated = true; - for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { - AddRealComp(); - } + SetArena(amrex::The_Arena()); - for (int j = 0; j < PlasmaIdx::int_nattribs; ++j) { - AddIntComp(); - } + for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { + AddRealComp(); + } - reserveData(); - resizeData(); + for (int j = 0; j < PlasmaIdx::int_nattribs; ++j) { + AddIntComp(); + } + + reserveData(); + resizeData(); + } if (!m_read_fine_patch) { m_read_fine_patch = true;