From 230a1a721ea630c54ca16e5967abd05cb09f5d01 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 22 Apr 2025 22:14:10 +0200 Subject: [PATCH 01/25] 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 02/25] 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 03/25] 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 04/25] 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 05/25] 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 06/25] 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 07/25] 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 08/25] 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 9135d5595d8255069b206e4a1ff8d7d1a7ba008a Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Tue, 14 Apr 2026 19:00:19 +0200 Subject: [PATCH 09/25] Dynamic Particle Components --- .github/workflows/linux.yml | 3 +- CMakeLists.txt | 11 - cmake/HiPACEFunctions.cmake | 7 - docs/source/building/building.rst | 1 - src/Hipace.cpp | 5 - .../plasma/PlasmaParticleContainer.H | 41 ++- .../plasma/PlasmaParticleContainer.cpp | 149 ++++++---- .../plasma/PlasmaParticleContainerInit.cpp | 71 +++-- .../pusher/PlasmaParticleAdvance.cpp | 268 ++++++++++-------- src/salame/Salame.cpp | 31 +- 10 files changed, 341 insertions(+), 246 deletions(-) diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 236f2d615e..e0ffe0a3f9 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -73,8 +73,7 @@ jobs: pip install -U -e ./tools cmake -S . -B build \ -DHiPACE_MPI=OFF \ - -DCMAKE_VERBOSE_MAKEFILE=ON \ - -DHiPACE_PUSHER=AB5 + -DCMAKE_VERBOSE_MAKEFILE=ON cmake --build build -j 2 - name: Run Tests run: ctest --test-dir build --output-on-failure diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d8e4c8277..314d8a4357 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,13 +77,6 @@ option(HiPACE_MPI "Multi-node support (message-passing)" ON) option(HiPACE_OPENPMD "openPMD I/O (HDF5, ADIOS)" ON) option(HiPACE_openpmd_mpi "parallel version of openPMD I/O" ${HiPACE_MPI}) -set(HiPACE_PUSHER_VALUES LEAPFROG AB5) -set(HiPACE_PUSHER LEAPFROG CACHE STRING "Plasma pusher (LEAPFROG/AB5)") -set_property(CACHE HiPACE_PUSHER PROPERTY STRINGS ${HiPACE_PUSHER_VALUES}) -if(NOT HiPACE_PUSHER IN_LIST HiPACE_PUSHER_VALUES) - message(FATAL_ERROR "HiPACE_PUSHER (${HiPACE_PUSHER}) must be one of ${HiPACE_PUSHER_VALUES}") -endif() - set(HiPACE_PRECISION_VALUES SINGLE DOUBLE) set(HiPACE_PRECISION DOUBLE CACHE STRING "Floating point precision (SINGLE/DOUBLE)") set_property(CACHE HiPACE_PRECISION PROPERTY STRINGS ${HiPACE_PRECISION_VALUES}) @@ -170,10 +163,6 @@ if(HiPACE_OPENPMD) target_link_libraries(HiPACE PUBLIC openPMD::openPMD) endif() -if(HiPACE_PUSHER STREQUAL "AB5") - target_compile_definitions(HiPACE PUBLIC HIPACE_USE_AB5_PUSH) -endif() - if(AMReX_LINEAR_SOLVERS) target_compile_definitions(HiPACE PUBLIC AMREX_USE_LINEAR_SOLVERS) endif() diff --git a/cmake/HiPACEFunctions.cmake b/cmake/HiPACEFunctions.cmake index 837f2867ba..fb264a2dcc 100644 --- a/cmake/HiPACEFunctions.cmake +++ b/cmake/HiPACEFunctions.cmake @@ -194,12 +194,6 @@ function(set_hipace_binary_name) set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".DEBUG") endif() - if(HiPACE_PUSHER STREQUAL "AB5") - set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".AB5") - else() - set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".LF") - endif() - # alias to the latest build, because using the full name is often confusing add_custom_command(TARGET HiPACE POST_BUILD COMMAND ${CMAKE_COMMAND} -E create_symlink @@ -274,6 +268,5 @@ function(hipace_print_summary) message(" MPI: ${HiPACE_MPI}") message(" OPENPMD: ${HiPACE_OPENPMD}") message(" PRECISION: ${HiPACE_PRECISION}") - message(" PUSHER: ${HiPACE_PUSHER}") message("") endfunction() diff --git a/docs/source/building/building.rst b/docs/source/building/building.rst index df523f6e1f..afd1f16966 100644 --- a/docs/source/building/building.rst +++ b/docs/source/building/building.rst @@ -190,7 +190,6 @@ or by providing arguments to the CMake call ``HiPACE_MPI`` **ON**/OFF Multi-node support (message-passing) ``HiPACE_PRECISION`` SINGLE/**DOUBLE** Floating point precision (single/double) ``HiPACE_OPENPMD`` **ON**/OFF openPMD I/O (HDF5, ADIOS2) - ``HiPACE_PUSHER`` **LEAPFROG**/AB5 Use leapfrog or fifth-order Adams-Bashforth plasma pusher ============================= ======================================== ========================================================= HiPACE++ can be configured in further detail with options from AMReX, which are documented in the `AMReX manual `__. diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 717e53d4fa..7f7c347c1b 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -285,11 +285,6 @@ Hipace::InitData () amrex::Print() << "using CUDA version " << __CUDACC_VER_MAJOR__ << "." << __CUDACC_VER_MINOR__ << "." << __CUDACC_VER_BUILD__ << "\n"; #endif -#ifdef HIPACE_USE_AB5_PUSH - amrex::Print() << "using the Adams-Bashforth plasma particle pusher\n"; -#else - amrex::Print() << "using the leapfrog plasma particle pusher\n"; -#endif m_multi_laser.InitData(); diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 971b42e121..c1ba0fe664 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -98,27 +98,46 @@ struct PlasmaIdx x=0, y, // position w, // weight, this will be returned by pos(2) ux, uy, // momentum - psi, // pseudo-potential at the particle position. ATTENTION what is stored is actually normalized psi+1 - x_prev, y_prev, // positions on the last non-temp slice - + psi, // pseudo-potential at the particle position. + // ATTENTION what is stored is actually normalized psi+1 ux_half_step, // momentum half a step behind the current slice for leapfrog pusher uy_half_step, // at the same step for AB5 pusher psi_half_step, // never effected by temp slice -#ifdef HIPACE_USE_AB5_PUSH - Fx1, Fx2, Fx3, Fx4, Fx5, // AB5 force terms + real_nattribs, + + // optional: + + aabssq, // Laser value at the particle position + real_nattribs_laser, + + x_prev=0, y_prev, // positions on the last non-temp slice + real_nattribs_temp_slice, + + Fx1=0, Fx2, Fx3, Fx4, Fx5, // AB5 force terms Fy1, Fy2, Fy3, Fy4, Fy5, // Fux1, Fux2, Fux3, Fux4, Fux5, // Fuy1, Fuy2, Fuy3, Fuy4, Fuy5, // Fpsi1, Fpsi2, Fpsi3, Fpsi4, Fpsi5, // -#endif - real_nattribs + real_nattribs_ab5_push }; enum { - ion_lev, // ionization level - int_nattribs + // optional: + + ion_lev=0, // ionization level + int_nattribs_ion_level }; }; +struct PlasmaComps { + bool use_laser = false; + bool use_temp_slice = false; + bool use_ab5_push = false; + bool use_ion_level = false; + + int offset_temp = -1; + int offset_ab5 = -1; +}; + /** \brief Container for particles of 1 plasma species. */ class PlasmaParticleContainer : public amrex::ParticleContainerRTSoA @@ -272,6 +291,10 @@ public: 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; + /** whether to use the AB5 or leapfrog pusher */ + bool m_use_ab5_push = false; + PlasmaComps m_comps; + bool m_is_on_temp_slice = false; // ionization: diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 88179a63a9..41c2a08530 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -49,14 +49,13 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_Da != 0, "Unknown Element"); } + queryWithParserAlt(pp, "use_ab5_push", m_use_ab5_push, pp_alt); queryWithParserAlt(pp, "n_subcycles", m_n_subcycles, pp_alt); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles >= 1, - "n_subcycles must be larger or equal to 1 sub-cycle (default is 1)"); -#ifdef HIPACE_USE_AB5_PUSH - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles == 1, - "Plasma subcycling only implemeted for leapfrog pusher!" - "Please set plasmas.n_subcycles = 1"); -#endif + "n_subcycles must be larger or equal to 1 sub-cycle (default is 1)"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_use_ab5_push || m_n_subcycles == 1, + "Plasma subcycling only implemeted for leapfrog pusher!" + "Please set plasmas.n_subcycles = 1"); queryWithParser(pp, "mass_Da", mass_Da); if(mass_Da != 0) { m_mass = phys_const.m_p * mass_Da / 1.007276466621; @@ -184,12 +183,46 @@ PlasmaParticleContainer::InitData (const amrex::Vector& geom3d) { SetArena(amrex::The_Arena()); + m_comps.use_laser = Hipace::m_use_laser; + m_comps.use_temp_slice = Hipace::GetInstance().m_multi_beam.AnySpeciesSalame() + || !Hipace::m_explicit; + m_comps.use_ab5_push = m_use_ab5_push; + m_comps.use_ion_level = m_can_ionize; + + int num_real_comps = 0; + for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { AddRealComp(); + ++num_real_comps; + } + + if (m_comps.use_laser) { + for (int j = 0; j < PlasmaIdx::real_nattribs_laser; ++j) { + AddRealComp(); + ++num_real_comps; + } } - for (int j = 0; j < PlasmaIdx::int_nattribs; ++j) { - AddIntComp(); + if (m_comps.use_temp_slice) { + m_comps.offset_temp = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_temp_slice; ++j) { + AddRealComp(); + ++num_real_comps; + } + } + + if (m_comps.use_ab5_push) { + m_comps.offset_ab5 = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_ab5_push; ++j) { + AddRealComp(); + ++num_real_comps; + } + } + + if (m_comps.use_ion_level) { + for (int j = 0; j < PlasmaIdx::int_nattribs_ion_level; ++j) { + AddIntComp(); + } } reserveData(); @@ -341,9 +374,11 @@ PlasmaParticleContainer::TagByLevel (const int current_N_level, { auto& ptile = pti.GetParticleTile(); const amrex::Real * const AMREX_RESTRICT pos_x = to_prev ? - ptile.GetRealData(PlasmaIdx::x_prev).data() : ptile.GetRealData(PlasmaIdx::x).data(); + ptile.GetRealData(PlasmaIdx::x_prev + m_comps.offset_temp).data() + : ptile.GetRealData(PlasmaIdx::x).data(); const amrex::Real * const AMREX_RESTRICT pos_y = to_prev ? - ptile.GetRealData(PlasmaIdx::y_prev).data() : ptile.GetRealData(PlasmaIdx::y).data(); + ptile.GetRealData(PlasmaIdx::y_prev + m_comps.offset_temp).data() + : ptile.GetRealData(PlasmaIdx::y).data(); auto * AMREX_RESTRICT idcpup = ptile.GetIdCPUData().data(); const int lev1_idx = std::min(1, current_N_level-1); @@ -455,9 +490,8 @@ IonizationModule (const int lev, if (!ptd_ion.id(ip).is_valid() || ptd_ion.cpu(ip) != lev) return; - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; // Define field at particle position reals amrex::ParticleReal ExmByp = 0., EypBxp = 0., Ezp = 0.; @@ -520,38 +554,46 @@ IonizationModule (const int lev, amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + auto comps = m_product_pc->m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice == m_comps.use_temp_slice); + AMREX_ALWAYS_ASSERT(comps.offset_temp == m_comps.offset_temp); + AMREX_ALWAYS_ASSERT(comps.use_ion_level == false); + // This kernel adds the new ionized electrons to the Plasma Particle Container amrex::ParallelFor(num_ions, [=] AMREX_GPU_DEVICE (long ip) { - if(p_ion_mask[ip] != 0) { + if (p_ion_mask[ip] != 0) { const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` 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 + // Later we could consider adding a finite temperature to the ionized electrons 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 - ptd_elec.rdata(PlasmaIdx::psi )[pidx] = 1._rt; // Assumes Aabssq == 0 - 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::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; + ptd_elec.rdata(PlasmaIdx::psi)[pidx] = 1._rt; // Assumes Aabssq == 0 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 - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd_elec.rdata(iforce)[pidx] = 0._rt; + + if (comps.use_temp_slice) { + ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd_elec.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } } -#endif - ptd_elec.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } }); @@ -646,9 +688,8 @@ LaserIonization (const int islice, [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, auto depos_order_xy) { - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; if (!ptd_ion.id(ip).is_valid() || !laser_bounds.contains(xp, yp)) return; @@ -719,6 +760,11 @@ LaserIonization (const int islice, amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + auto comps = m_product_pc->m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice == m_comps.use_temp_slice); + AMREX_ALWAYS_ASSERT(comps.offset_temp == m_comps.offset_temp); + AMREX_ALWAYS_ASSERT(comps.use_ion_level == false); + // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time. // It calculates the momentum of ionized electrons based on equations (B8) and (B9) // from the F. Massimo (2020) article and equation (12) from the P. Tomassini (2021) article. @@ -738,9 +784,8 @@ LaserIonization (const int islice, if(p_ion_mask[ip] != 0) { - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; if (!ptd_ion.id(ip).is_valid() || !laser_bounds.contains(xp, yp)) return; @@ -793,24 +838,28 @@ LaserIonization (const int islice, // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons ptd_elec.id(pidx) = 2; ptd_elec.cpu(pidx) = ptd_ion.cpu(ip); // 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] = ux; - ptd_elec.rdata(PlasmaIdx::uy )[pidx] = uy; - ptd_elec.rdata(PlasmaIdx::psi )[pidx] = psi; - 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::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; + ptd_elec.rdata(PlasmaIdx::uy)[pidx] = uy; + ptd_elec.rdata(PlasmaIdx::psi)[pidx] = psi; ptd_elec.rdata(PlasmaIdx::ux_half_step )[pidx] = ux; ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = uy; ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = psi; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd_elec.rdata(iforce)[pidx] = 0._rt; + + if (comps.use_temp_slice) { + ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd_elec.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } } -#endif - 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 02b6e2ca1d..a4686d0ef2 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -298,6 +298,9 @@ InitParticles (const amrex::RealVect& a_u_std, AMREX_ALWAYS_ASSERT(total_non_mirrored_particles == current_size); + const bool use_ab5_push = m_use_ab5_push; + auto comps = m_comps; + amrex::ParallelForRNG(current_size, [=] AMREX_GPU_DEVICE (unsigned int pidx, const amrex::RandomEngine& engine) { const amrex::Real x = ptd.rdata(PlasmaIdx::x)[pidx]; @@ -308,20 +311,27 @@ InitParticles (const amrex::RealVect& a_u_std, ptd.rdata(PlasmaIdx::ux)[pidx] = u[0]; ptd.rdata(PlasmaIdx::uy)[pidx] = u[1]; - ptd.rdata(PlasmaIdx::psi)[pidx] = plasma_psi(u[0], u[1], u[2], - /* Assumes Aabssq == 0 */ 0._rt); - ptd.rdata(PlasmaIdx::x_prev)[pidx] = x; - ptd.rdata(PlasmaIdx::y_prev)[pidx] = y; - ptd.rdata(PlasmaIdx::ux_half_step)[pidx] = u[0]; - ptd.rdata(PlasmaIdx::uy_half_step)[pidx] = u[1]; - ptd.rdata(PlasmaIdx::psi_half_step)[pidx] = ptd.rdata(PlasmaIdx::psi)[pidx]; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd.rdata(iforce)[pidx] = 0._rt; + ptd.rdata(PlasmaIdx::psi)[pidx] = + plasma_psi(u[0], u[1], u[2], /* Assumes Aabssq == 0 */ 0._rt); + + if (comps.use_laser) { + ptd.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + } + + if (comps.use_temp_slice) { + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = x; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = y; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } + } + + if (comps.use_ion_level) { + ptd.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } -#endif - ptd.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; }); if (m_do_symmetrize) { @@ -349,27 +359,38 @@ InitParticles (const amrex::RealVect& a_u_std, ptd.cpu(midx) = 0; // level 0 ptd.rdata(PlasmaIdx::x)[midx] = x_arr[imirror]; ptd.rdata(PlasmaIdx::y)[midx] = y_arr[imirror]; - ptd.rdata(PlasmaIdx::w)[midx] = ptd.rdata(PlasmaIdx::w)[pidx]; - ptd.rdata(PlasmaIdx::ux)[midx] = ptd.rdata(PlasmaIdx::ux)[pidx] * ux_arr[imirror]; - ptd.rdata(PlasmaIdx::uy)[midx] = ptd.rdata(PlasmaIdx::uy)[pidx] * uy_arr[imirror]; - ptd.rdata(PlasmaIdx::psi)[midx] = ptd.rdata(PlasmaIdx::psi)[pidx]; - ptd.rdata(PlasmaIdx::x_prev)[midx] = x_arr[imirror]; - ptd.rdata(PlasmaIdx::y_prev)[midx] = y_arr[imirror]; + ptd.rdata(PlasmaIdx::ux)[midx] = + ptd.rdata(PlasmaIdx::ux)[pidx] * ux_arr[imirror]; + ptd.rdata(PlasmaIdx::uy)[midx] = + ptd.rdata(PlasmaIdx::uy)[pidx] * uy_arr[imirror]; + ptd.rdata(PlasmaIdx::psi)[midx] = + ptd.rdata(PlasmaIdx::psi)[pidx]; ptd.rdata(PlasmaIdx::ux_half_step)[midx] = ptd.rdata(PlasmaIdx::ux_half_step)[pidx] * ux_arr[imirror]; ptd.rdata(PlasmaIdx::uy_half_step)[midx] = ptd.rdata(PlasmaIdx::uy_half_step)[pidx] * uy_arr[imirror]; ptd.rdata(PlasmaIdx::psi_half_step)[midx] = ptd.rdata(PlasmaIdx::psi_half_step)[pidx]; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd.rdata(iforce)[midx] = 0._rt; + + if (comps.use_laser) { + ptd.rdata(PlasmaIdx::aabssq)[midx] = 0._rt; } -#endif - ptd.idata(PlasmaIdx::ion_lev)[midx] = ptd.idata(PlasmaIdx::ion_lev)[pidx]; + if (comps.use_temp_slice) { + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[midx] = x_arr[imirror]; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[midx] = y_arr[imirror]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd.rdata(iforce + comps.offset_ab5)[midx] = 0._rt; + } + } + + if (comps.use_ion_level) { + ptd.idata(PlasmaIdx::ion_lev)[midx] = init_ion_lev; + } } }); } diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 3c5256c2a5..1a5c5d8c4b 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -79,17 +79,23 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, const amrex::Real clight_inv = 1._rt/phys_const.c; const amrex::Real charge_mass_clight_ratio = plasma.m_charge/(plasma.m_mass * phys_const.c); + const auto comps = plasma.m_comps; + const bool read_from_temp = plasma.m_is_on_temp_slice; + AMREX_ALWAYS_ASSERT(!(read_from_temp || temp_slice) || comps.use_temp_slice); + // Use OMP ParallelFor to use multiple threads when running on CPU omp::ParallelFor( amrex::TypeList< amrex::CompileTimeOptions<0, 1, 2, 3>, + amrex::CompileTimeOptions, amrex::CompileTimeOptions >{}, { Hipace::m_depos_order_xy, - Hipace::m_use_laser + Hipace::m_use_laser, + plasma.m_use_ab5_push }, int(pti.numParticles()), // int ParallelFor is 3-5% faster than amrex::Long version - [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser) { + [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser, auto use_ab5_push) { // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; @@ -106,10 +112,18 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; } - for (int i = 0; i < n_subcycles; i++) { + amrex::Real xp = 0._rt; + amrex::Real yp = 0._rt; + + if (!read_from_temp) { + xp = ptd.rdata(PlasmaIdx::x)[ip]; + yp = ptd.rdata(PlasmaIdx::y)[ip]; + } else { + xp = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + yp = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + } - amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev)[ip]; - amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev)[ip]; + for (int i = 0; i < n_subcycles; i++) { if (lev == 0 || lev_bounds.contains(xp, yp)) { ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt; @@ -134,150 +148,158 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, AabssqDyp *= 0.25_rt * laser_norm_ion; } -#ifndef HIPACE_USE_AB5_PUSH + if (!use_ab5_push.value) { - constexpr int nsub = 4; - const amrex::Real sdz = dz/nsub; + constexpr int nsub = 4; + const amrex::Real sdz = dz/nsub; - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; + amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; + amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; + amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; - // full push in momentum - // from t-1/2 to t+1/2 - // using the fields at t - for (int isub=0; isub= 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]; + } - } - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; -#else - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; - const amrex::Real psi_inv = 1._rt/psi; - - auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush( - ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp, - Aabssqp, AabssqDxp, AabssqDyp, q_mass_clight_ratio); - - 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 - ( -1387._rt / 360._rt ) * dz, // a2 times dz - ( 109._rt / 30._rt ) * dz, // a3 times dz - ( -637._rt / 360._rt ) * dz, // a4 times dz - ( 251._rt / 720._rt ) * dz // a5 times dz - }; - - HIPACE_LOOP_UNROLL - for (int iab=0; iab<5; ++iab) { - int p = ab5_permutation + iab; - if (p >= 5) { - p -= 5; + if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; + ptd.pos(0, ip) = xp; + ptd.pos(1, ip) = yp; + + if (!temp_slice) { + // update values of the last non temp slice + // the next push always starts from these + ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux; + ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy; + ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi; + ptd.rdata(PlasmaIdx::x_prev)[ip] = xp; + ptd.rdata(PlasmaIdx::y_prev)[ip] = yp; } - 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; - ptd.pos(0, ip) = xp; - ptd.pos(1, ip) = yp; - - if (!temp_slice) { - // update values of the last non temp slice - // the next push always starts from these - ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux; - ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy; - ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi; - ptd.rdata(PlasmaIdx::x_prev)[ip] = xp; - ptd.rdata(PlasmaIdx::y_prev)[ip] = yp; + ptd.rdata(PlasmaIdx::ux)[ip] = ux; + ptd.rdata(PlasmaIdx::uy)[ip] = uy; + ptd.rdata(PlasmaIdx::psi)[ip] = psi; } - - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; -#endif } // loop over subcycles + + ptd.pos(0, ip) = xp; + ptd.pos(1, ip) = yp; }); } -#ifdef HIPACE_USE_AB5_PUSH - if (!temp_slice && lev == current_N_level - 1) { + if (lev == current_N_level - 1) { + plasma.m_is_on_temp_slice = temp_slice; + } + + if (plasma.m_use_ab5_push && !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 abaefe8b75..7236c6e111 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -234,11 +234,20 @@ SalameGetJxJyFromBxBy (Hipace* hipace, const int lev) amrex::MultiFab& slicemf = hipace->m_fields.getSlices(lev); -#ifdef HIPACE_USE_AB5_PUSH - const amrex::Real dz = ( 1901._rt / 720._rt ) * hipace->m_3D_geom[lev].CellSize(Direction::z); -#else - const amrex::Real dz = 1.5_rt * hipace->m_3D_geom[lev].CellSize(Direction::z); -#endif + bool use_ab5_pusher = false; + for (int i=0; im_multi_plasma.GetNPlasmas(); ++i) { + if (i == 0) { + use_ab5_pusher = hipace->m_multi_plasma.m_all_plasmas[i].m_use_ab5_push; + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + hipace->m_multi_plasma.m_all_plasmas[i].m_use_ab5_push == use_ab5_pusher, + "All plasmas must use the same pusher when using SALAME" + ); + } + } + + const amrex::Real dz = (use_ab5_pusher ? 1901._rt / 720._rt : 1.5_rt) + * hipace->m_3D_geom[lev].CellSize(Direction::z); for ( amrex::MFIter mfi(slicemf, DfltMfiTlng); mfi.isValid(); ++mfi ){ @@ -280,7 +289,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) const amrex::Real dx_inv = gm.InvCellSize(0); const amrex::Real dy_inv = gm.InvCellSize(1); - const amrex::Real dz = gm.CellSize(2); + const amrex::Real dz = (plasma.m_use_ab5_push ? 1901._rt / 720._rt : 1.5_rt) + * gm.CellSize(2); // Offset for converting positions to indexes amrex::Real const x_pos_offset = GetPosOffset(0, gm, slice_fab.box()); @@ -314,13 +324,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) ptd.idata(PlasmaIdx::ion_lev)[ip] * charge_mass_c_ratio : charge_mass_c_ratio; -#ifdef HIPACE_USE_AB5_PUSH - ptd.rdata(PlasmaIdx::ux)[ip] = ( 1901._rt / 720._rt )*dz * q_m_c_ratio * Byp; - ptd.rdata(PlasmaIdx::uy)[ip] = -( 1901._rt / 720._rt )*dz * q_m_c_ratio * Bxp; -#else - ptd.rdata(PlasmaIdx::ux)[ip] = 1.5_rt*dz * q_m_c_ratio * Byp; - ptd.rdata(PlasmaIdx::uy)[ip] = -1.5_rt*dz * q_m_c_ratio * Bxp; -#endif + ptd.rdata(PlasmaIdx::ux)[ip] = dz * q_m_c_ratio * Byp; + ptd.rdata(PlasmaIdx::uy)[ip] = -dz * q_m_c_ratio * Bxp; }); } From 075ef7f28e25ab429fca8642a109cd08d5bbe818 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Wed, 15 Apr 2026 13:37:57 +0200 Subject: [PATCH 10/25] use laser component --- .../deposition/ExplicitDeposition.cpp | 9 +- .../deposition/PlasmaDepositCurrent.cpp | 58 ++++------- .../deposition/TemperatureDeposition.cpp | 56 +++-------- .../plasma/PlasmaParticleContainer.cpp | 41 ++------ .../plasma/PlasmaParticleContainerInit.cpp | 1 - .../pusher/PlasmaParticleAdvance.cpp | 98 ++++++++++--------- 6 files changed, 96 insertions(+), 167 deletions(-) diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index f281dbec29..5ab6f51515 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -52,7 +52,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real a_clight = pc.c; const amrex::Real clight_inv = 1._rt/pc.c; // The laser a0 is always normalized - const amrex::Real laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); + const amrex::Real a_laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); const amrex::Real charge_invvol_mu0 = plasma.m_charge * invvol * pc.mu0; const amrex::Real charge_mass_ratio = plasma.m_charge / plasma.m_mass; @@ -168,11 +168,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, amrex::Real Aabssqp = 0._rt; if (use_laser) { - // Its important that Aabssqp is first fully gathered and not used - // directly per cell like AabssqDxp and AabssqDyp - doLaserGatherShapeN(xp, yp, Aabssqp, arr, cache_idx[4], - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_fac * q_mass_ratio * q_mass_ratio; + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } // calculate gamma/psi for plasma particles @@ -203,6 +199,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, amrex::Real AabssqDyp = 0._rt; // Rename variables for NVCC lambda capture to work [[maybe_unused]] auto clight = a_clight; + [[maybe_unused]] auto laser_fac = a_laser_fac; if constexpr (use_laser) { // avoid going outside of domain if (shape_x * shape_y != 0._rt) { diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 418bad9e8f..e9037d1e17 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/src/particles/deposition/PlasmaDepositCurrent.cpp @@ -64,7 +64,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const int chi = deposit_chi ? Comps[which_slice]["chi"] : -1; const int rhomjz = deposit_rhomjz ? Comps[which_slice]["rhomjz"] : -1; const int n = deposit_n ? Comps[which_slice][n_str] : -1; - const int aabs = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; // Offset for converting positions to indexes const amrex::Real x_pos_offset = GetPosOffset(0, gm[lev], isl_fab.box()); @@ -84,8 +83,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const amrex::Real clight = pc.c; const amrex::Real charge_invvol = charge * invvol; const amrex::Real charge_mu0_mass_ratio = charge * pc.mu0 / mass; - const amrex::Real laser_norm = (charge/pc.q_e) * (pc.m_e/mass) - * (charge/pc.q_e) * (pc.m_e/mass); amrex::Gpu::DeviceScalar gpu_n_qsa_violation{}; int* AMREX_RESTRICT p_n_qsa_violation = nullptr; @@ -97,6 +94,9 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, &n_qsa_violation, &n_qsa_violation + 1, p_n_qsa_violation); } + const bool use_laser = Hipace::m_use_laser; + const bool can_ionize = plasma.m_can_ionize; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(isl_fab.box().ixType().cellCentered(), "jx, jy, jz, and rho must be cell centered in all directions."); @@ -104,13 +104,9 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, amrex::AnyCTO( // use compile-time options amrex::TypeList< - amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order - amrex::CompileTimeOptions, // can_ionize - amrex::CompileTimeOptions // use_laser + amrex::CompileTimeOptions<0, 1, 2, 3> // depos_order >{}, { - Hipace::m_depos_order_xy, - plasma.m_can_ionize, - Hipace::m_use_laser + Hipace::m_depos_order_xy }, // call deposition function // The three functions passed as arguments to this lambda @@ -118,28 +114,17 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, [&](auto is_valid, auto get_cell, auto deposit){ constexpr auto ctos = deposit.GetOptions(); constexpr int depos_order = ctos[0]; - constexpr int use_laser = ctos[2]; constexpr int stencil_size = depos_order + 1; - if constexpr (use_laser) { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{aabs}, - amrex::GpuArray{jx, jy, jz, rho, chi, rhomjz, n}); - } else { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{}, - amrex::GpuArray{jx, jy, jz, rho, chi, rhomjz, n}); - } + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{}, + amrex::GpuArray{jx, jy, jz, rho, chi, rhomjz, n}); }, // is_valid // return whether the particle is valid and should deposit [=] AMREX_GPU_DEVICE (int ip, auto ptd, - auto /*depos_order*/, - auto can_ionize, - auto /*use_laser*/) + auto /*depos_order*/) { // only deposit plasma currents on or below their according MR level return ptd.id(ip).is_valid() && @@ -150,9 +135,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, // get_cell // return the lowest cell index that the particle deposits into [=] AMREX_GPU_DEVICE (int ip, auto ptd, - auto depos_order, - auto /*can_ionize*/, - auto /*use_laser*/) -> amrex::IntVectND<2> + auto depos_order) -> amrex::IntVectND<2> { const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); @@ -161,7 +144,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; auto [shape_x, i] = shape_factor(xmid, 0); - auto [shape_y, j] = shape_factor(ymid, 0); return {i, j}; @@ -170,10 +152,8 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, // deposit the charge / current of one particle [=] AMREX_GPU_DEVICE (int ip, auto ptd, Array3 arr, - auto cache_idx, auto depos_idx, - auto depos_order, - auto can_ionize, - auto use_laser) noexcept + auto /*cache_idx*/, auto depos_idx, + auto depos_order) noexcept { const amrex::Real psi_inv = 1._rt/ptd.rdata(PlasmaIdx::psi)[ip]; const amrex::Real xp = ptd.pos(0, ip); @@ -184,22 +164,18 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, // calculate charge of the plasma particles amrex::Real q_invvol = charge_invvol * w; amrex::Real q_mu0_mass_ratio = charge_mu0_mass_ratio; - [[maybe_unused]] amrex::Real laser_norm_ion = laser_norm; - if constexpr (can_ionize) { + if (can_ionize) { const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); q_invvol *= p_ion_lev; q_mu0_mass_ratio *= p_ion_lev; - laser_norm_ion *= p_ion_lev * p_ion_lev; } const amrex::Real xmid = (xp - x_pos_offset) * dx_inv; const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; amrex::Real Aabssqp = 0._rt; - if constexpr (use_laser) { - doLaserGatherShapeN(xp, yp, Aabssqp, arr, cache_idx[0], - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + if (use_laser) { + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } // calculate gamma/psi for plasma particles diff --git a/src/particles/deposition/TemperatureDeposition.cpp b/src/particles/deposition/TemperatureDeposition.cpp index 603c21f82e..568e2c3bf4 100644 --- a/src/particles/deposition/TemperatureDeposition.cpp +++ b/src/particles/deposition/TemperatureDeposition.cpp @@ -51,23 +51,15 @@ DepositTemperature (PlasmaParticleContainer& plasma, // Extract box properties const amrex::Real dx_inv = gm[lev].InvCellSize(0); const amrex::Real dy_inv = gm[lev].InvCellSize(1); - // extract laser properties and boolean for the presence of the laser and for ionization - const PhysConst pc = get_phys_const(); - const int aabs = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; - const amrex::Real laser_norm = (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass) - * (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass); + const bool use_laser = Hipace::m_use_laser; // Loop over particles amrex::AnyCTO( // use compile-time options amrex::TypeList< - amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order - amrex::CompileTimeOptions, // can_ionize - amrex::CompileTimeOptions // use_laser + amrex::CompileTimeOptions<0, 1, 2, 3> // depos_order >{}, { - Hipace::m_temperature_depos_order, - plasma.m_can_ionize, - Hipace::m_use_laser + Hipace::m_temperature_depos_order }, // call deposition function // The three functions passed as arguments to this lambda @@ -75,29 +67,18 @@ DepositTemperature (PlasmaParticleContainer& plasma, [&](auto is_valid, auto get_start_cell, auto deposit){ constexpr auto ctos = deposit.GetOptions(); constexpr int depos_order = ctos[0]; - constexpr int use_laser = ctos[2]; constexpr int stencil_size = depos_order + 1; - if constexpr (use_laser) { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{aabs}, - amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); - } else { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{}, - amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); - } + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{}, + amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); }, // is_valid // return whether the particle is valid and should deposit [=] AMREX_GPU_DEVICE (int ip, auto ptd, - auto /*depos_order*/, - auto /*can_ionize*/, - auto /*use_laser*/) + auto /*depos_order*/) { // only deposit on or below their according MR level return ptd.id(ip).is_valid() && (lev == 0 || ptd.cpu(ip) >= lev); @@ -105,9 +86,7 @@ DepositTemperature (PlasmaParticleContainer& plasma, // get_start_cell // return the lowest cell index that the particle deposits into [=] AMREX_GPU_DEVICE (int ip, auto ptd, - auto depos_order, - auto /*can_ionize*/, - auto /*use_laser*/) -> amrex::IntVectND<2> + auto depos_order) -> amrex::IntVectND<2> { const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); @@ -124,24 +103,15 @@ DepositTemperature (PlasmaParticleContainer& plasma, // deposit of weight, momentum (ux, uy, uz) and their squares (uxsq, uysq, uzsq) [=] AMREX_GPU_DEVICE (int ip, auto ptd, Array3 arr, - auto cache_idx, auto depos_idx, - auto depos_order, - auto can_ionize, - auto use_laser) noexcept + auto /*cache_idx*/, auto depos_idx, + auto depos_order) noexcept { const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); amrex::Real Aabssqp = 0._rt; if (use_laser) { - amrex::Real laser_norm_ion = laser_norm; - if (can_ionize) { - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; - } - doLaserGatherShapeN<2>(xp, yp, Aabssqp, arr, cache_idx[0], - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } const amrex::Real uxp = ptd.rdata(PlasmaIdx::ux)[ip]; diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 41c2a08530..a7b1083786 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -469,7 +469,7 @@ IonizationModule (const int lev, const int max_ion_lev = m_max_ion_lev; long num_ions = ptile_ion.numParticles(); - + const bool use_laser = Hipace::m_use_laser; // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time // and calculates ionization probability. If ionization occurs, it increments @@ -511,8 +511,11 @@ IonizationModule (const int lev, const amrex::Real psi = ptd_ion.rdata(PlasmaIdx::psi_half_step)[ip]; // Compute probability of ionization p - const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, - /* Assumes Aabssq == 0 */ 0._rt); + amrex::Real Aabssqp = 0._rt; + if (use_laser) { + Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; + } + const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, Aabssqp); const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; if (ion_lev_loc >= max_ion_lev) { return; @@ -549,8 +552,6 @@ IonizationModule (const int lev, // Load electron after resize auto ptd_elec = ptile_elec.getParticleTileData(); - const int init_ion_lev = m_product_pc->m_init_ion_lev; - amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); @@ -712,9 +713,9 @@ LaserIonization (const int islice, const amrex::Real uy = ptd_ion.rdata(PlasmaIdx::uy_half_step)[ip]; const amrex::Real psi = ptd_ion.rdata(PlasmaIdx::psi_half_step)[ip]; + const amrex::Real Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; // Compute probability of ionization p - const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, - /* Assumes Aabssq == 0 */ 0._rt); + const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, Aabssqp); const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; if (ion_lev_loc >= max_ion_lev) { return; @@ -755,8 +756,6 @@ LaserIonization (const int islice, // Load electron after resize auto ptd_elec = ptile_elec.getParticleTileData(); - const int init_ion_lev = m_product_pc->m_init_ion_lev; - amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); @@ -833,6 +832,7 @@ LaserIonization (const int islice, const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` const long pidx = pid + old_size; + // No normalization required since these are always electrons const amrex::Real psi = plasma_psi(ux, uy, uz, amrex::abs(A*A)); // Copy ion data to new electron // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons @@ -885,22 +885,8 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) { // Loading the data const auto ptd = pti.GetParticleTile().getParticleTileData(); - amrex::Long const num_particles = pti.numParticles(); - - const PhysConst pc = get_phys_const(); const bool use_laser = Hipace::m_use_laser; - const amrex::Geometry& gm = Hipace::GetInstance().m_3D_geom[0]; - const int aabs_comp = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; - amrex::FArrayBox& isl_fab = Hipace::GetInstance().m_fields.getSlices(0)[pti]; - Array3 arr = isl_fab.array(); - const amrex::Real x_pos_offset = GetPosOffset(0, gm, isl_fab.box()); - const amrex::Real y_pos_offset = GetPosOffset(1, gm, isl_fab.box()); - const amrex::Real dx_inv = gm.InvCellSize(0); - const amrex::Real dy_inv = gm.InvCellSize(1); - const bool can_ionize = m_can_ionize; - const amrex::Real laser_norm = (m_charge/pc.q_e) * (pc.m_e/m_mass) - * (m_charge/pc.q_e) * (pc.m_e/m_mass); amrex::TypeMultiplier reduce_op; amrex::TypeMultiplier reduce_data(reduce_op); @@ -922,14 +908,7 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) amrex::Real Aabssqp = 0._rt; if (use_laser) { - amrex::Real laser_norm_ion = laser_norm; - if (can_ionize) { - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; - } - doLaserGatherShapeN<2>(x, y, Aabssqp, arr, aabs_comp, - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } // Particle's Lorentz factor diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index a4686d0ef2..f24768b442 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -298,7 +298,6 @@ InitParticles (const amrex::RealVect& a_u_std, AMREX_ALWAYS_ASSERT(total_non_mirrored_particles == current_size); - const bool use_ab5_push = m_use_ab5_push; auto comps = m_comps; amrex::ParallelForRNG(current_size, diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 1a5c5d8c4b..d6c0e802f5 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -107,9 +107,9 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::Real q_mass_clight_ratio = charge_mass_clight_ratio; amrex::Real laser_norm_ion = laser_norm; if (can_ionize) { - q_mass_clight_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip]; - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; + const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); + q_mass_clight_ratio *= p_ion_lev; + laser_norm_ion *= p_ion_lev * p_ion_lev; } amrex::Real xp = 0._rt; @@ -118,34 +118,49 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, if (!read_from_temp) { xp = ptd.rdata(PlasmaIdx::x)[ip]; yp = ptd.rdata(PlasmaIdx::y)[ip]; + if (temp_slice) { + // first temp slice + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip] = xp; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip] = yp; + } } else { xp = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; yp = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; } + amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; + amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; + amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; + + amrex::Real ux_half = ux; + amrex::Real uy_half = uy; + amrex::Real psi_half = psi; + for (int i = 0; i < n_subcycles; i++) { - if (lev == 0 || lev_bounds.contains(xp, yp)) { + if (i == 0 || lev == 0 || lev_bounds.contains(xp, yp)) { ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt; Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt; - Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; doGatherShapeN(xp, yp, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp, bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); if (use_laser.value) { + Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; + doLaserGatherShapeN(xp, yp, Aabssqp, AabssqDxp, AabssqDyp, slice_arr, aabs_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); + + Aabssqp *= 0.5_rt * laser_norm_ion; + AabssqDxp *= 0.25_rt * laser_norm_ion; + AabssqDyp *= 0.25_rt * laser_norm_ion; } ExmByp *= clight_inv; EypBxp *= clight_inv; Ezp *= clight_inv; - Aabssqp *= 0.5_rt * laser_norm_ion; - AabssqDxp *= 0.25_rt * laser_norm_ion; - AabssqDyp *= 0.25_rt * laser_norm_ion; } if (!use_ab5_push.value) { @@ -153,10 +168,6 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, constexpr int nsub = 4; const amrex::Real sdz = dz/nsub; - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; - // full push in momentum // from t-1/2 to t+1/2 // using the fields at t @@ -190,18 +201,9 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; - if (!temp_slice) { - // update values of the last non temp slice - // the next push always starts from these - ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux; - ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy; - ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi; - - if (comps.use_temp_slice) { - ptd.rdata(PlasmaIdx::x_prev)[ip] = xp; - ptd.rdata(PlasmaIdx::y_prev)[ip] = yp; - } - } + ux_half = ux; + uy_half = uy; + psi_half = psi; // half push in momentum // from t+1/2 to t+1 @@ -228,15 +230,9 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon; } - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; } else { - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; const amrex::Real psi_inv = 1._rt/psi; auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush( @@ -271,27 +267,39 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, } if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; - ptd.pos(0, ip) = xp; - ptd.pos(1, ip) = yp; - - if (!temp_slice) { - // update values of the last non temp slice - // the next push always starts from these - ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux; - ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy; - ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi; - ptd.rdata(PlasmaIdx::x_prev)[ip] = xp; - ptd.rdata(PlasmaIdx::y_prev)[ip] = yp; - } - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; + ux_half = ux; + uy_half = uy; + psi_half = psi; } } // loop over subcycles + if (use_laser.value) { + Aabssqp = 0._rt; + + doLaserGatherShapeN(xp, yp, + Aabssqp, slice_arr, aabs_comp, + dx_inv, dy_inv, x_pos_offset, y_pos_offset); + + Aabssqp *= laser_norm_ion; + + ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; + } + ptd.pos(0, ip) = xp; ptd.pos(1, ip) = yp; + + ptd.rdata(PlasmaIdx::ux)[ip] = ux; + ptd.rdata(PlasmaIdx::uy)[ip] = uy; + ptd.rdata(PlasmaIdx::psi)[ip] = psi; + + if (!temp_slice) { + // update values of the last non temp slice + // the next push always starts from these + ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux_half; + ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy_half; + ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi_half; + } }); } From 964d6708c28ac82dfc52e4f0287022add4f94910 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Wed, 15 Apr 2026 14:41:29 +0200 Subject: [PATCH 11/25] fix plasma init --- src/particles/plasma/PlasmaParticleContainer.H | 4 ++-- src/particles/plasma/PlasmaParticleContainerInit.cpp | 3 +++ src/particles/pusher/PlasmaParticleAdvance.cpp | 8 ++++---- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index c1ba0fe664..d463463ea3 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -107,8 +107,8 @@ struct PlasmaIdx // optional: - aabssq, // Laser value at the particle position - real_nattribs_laser, + aabssq=real_nattribs, // Laser value at the particle position + real_nattribs_laser=1, x_prev=0, y_prev, // positions on the last non-temp slice real_nattribs_temp_slice, diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index f24768b442..9d171287e3 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -312,6 +312,9 @@ InitParticles (const amrex::RealVect& a_u_std, ptd.rdata(PlasmaIdx::uy)[pidx] = u[1]; ptd.rdata(PlasmaIdx::psi)[pidx] = plasma_psi(u[0], u[1], u[2], /* Assumes Aabssq == 0 */ 0._rt); + ptd.rdata(PlasmaIdx::ux_half_step)[pidx] = u[0]; + ptd.rdata(PlasmaIdx::uy_half_step)[pidx] = u[1]; + ptd.rdata(PlasmaIdx::psi_half_step)[pidx] = ptd.rdata(PlasmaIdx::psi)[pidx]; if (comps.use_laser) { ptd.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index d6c0e802f5..87c29250b3 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -146,6 +146,10 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp, bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); + ExmByp *= clight_inv; + EypBxp *= clight_inv; + Ezp *= clight_inv; + if (use_laser.value) { Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; @@ -157,10 +161,6 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, AabssqDxp *= 0.25_rt * laser_norm_ion; AabssqDyp *= 0.25_rt * laser_norm_ion; } - - ExmByp *= clight_inv; - EypBxp *= clight_inv; - Ezp *= clight_inv; } if (!use_ab5_push.value) { From 6afa0de5279e0155fd8db88fd030292f8c3724b4 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Wed, 15 Apr 2026 16:22:35 +0200 Subject: [PATCH 12/25] Use CTO for plasma subcycling --- src/particles/pusher/PlasmaParticleAdvance.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 87c29250b3..ae9551c549 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -88,14 +88,17 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::TypeList< amrex::CompileTimeOptions<0, 1, 2, 3>, amrex::CompileTimeOptions, + amrex::CompileTimeOptions, amrex::CompileTimeOptions >{}, { Hipace::m_depos_order_xy, Hipace::m_use_laser, - plasma.m_use_ab5_push + plasma.m_use_ab5_push, + n_subcycles > 1 }, int(pti.numParticles()), // int ParallelFor is 3-5% faster than amrex::Long version - [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser, auto use_ab5_push) { + [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser, + auto use_ab5_push, auto use_subcycling) { // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; @@ -136,7 +139,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::Real uy_half = uy; amrex::Real psi_half = psi; - for (int i = 0; i < n_subcycles; i++) { + for (int i = 0; i < (use_subcycling ? n_subcycles : 1); i++) { if (i == 0 || lev == 0 || lev_bounds.contains(xp, yp)) { ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt; From 0efbe108d149aac16e304962334846a4e5dd3ac3 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 16 Apr 2026 13:07:53 +0200 Subject: [PATCH 13/25] gather laser from correct slice --- src/Hipace.cpp | 5 ++ .../deposition/TemperatureDeposition.cpp | 2 +- src/particles/plasma/MultiPlasma.H | 10 ++- src/particles/plasma/MultiPlasma.cpp | 8 +++ .../plasma/PlasmaParticleContainer.H | 11 +++ .../plasma/PlasmaParticleContainer.cpp | 69 +++++++++++++++++++ .../pusher/PlasmaParticleAdvance.cpp | 12 ---- 7 files changed, 103 insertions(+), 14 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 7f7c347c1b..9720e1d0c0 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -700,6 +700,11 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st // write laser aabs into fields MultiFab m_multi_laser.UpdateLaserAabs(islice, current_N_level, m_fields, m_3D_geom); + // interpolate laser aabs to plasma particles + for (int lev=0; lev(xmid, 0); auto [shape_y, j] = shape_factor(ymid, 0); - return {i-1, j-1}; + return {i, j}; }, // do_deposit // deposit of weight, momentum (ux, uy, uz) and their squares (uxsq, uysq, uzsq) diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index 452f7e8faa..cb9438e7c2 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -83,7 +83,7 @@ public: /** \brief Gather field values and push particles * - * \param[in,out] fields the general field class, modified by this function + * \param[in] fields the general field class * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] temp_slice if true, the temporary data (x_temp, ...) will be used * \param[in] lev MR level @@ -93,6 +93,14 @@ public: const Fields & fields, amrex::Vector const& gm, bool temp_slice, int lev, int const current_N_level); + /** \brief Gather Laser field to the particles + * + * \param[in] lev MR level + * \param[in] gm Geometry of the simulation, to get the cell size etc. + * \param[in] fields the general field class + */ + void GatherLaser (int lev, amrex::Geometry const& gm, const Fields & fields); + /** \brief Loop over plasma species and deposit their neutralizing background, if needed * * \param[in,out] fields the general field class, modified by this function diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index f685f10ed1..691c50a0ea 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -133,6 +133,14 @@ MultiPlasma::AdvanceParticles ( } } +void +MultiPlasma::GatherLaser (int lev, amrex::Geometry const& gm, const Fields & fields) +{ + for (int i=0; i @@ -868,6 +869,74 @@ LaserIonization (const int islice, } } +void +PlasmaParticleContainer:: +GatherLaser (const int lev, + const amrex::Geometry& geom, + const Fields& fields) +{ + if (!Hipace::m_use_laser) return; + HIPACE_PROFILE("PlasmaParticleContainer::GatherLaser()"); + + using namespace amrex::literals; + + const PhysConst phys_const = get_phys_const(); + + // Loop over particle boxes + for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) + { + // Extract field array from FabArray + const amrex::FArrayBox& slice_fab = fields.getSlices(lev)[pti]; + Array3 const slice_arr = slice_fab.const_array(); + const int aabs_comp = Comps[WhichSlice::This]["aabs"]; + + // Extract properties associated with physical size of the box + const amrex::Real dx_inv = geom.InvCellSize(0); + const amrex::Real dy_inv = geom.InvCellSize(1); + + // Offset for converting positions to indexes + amrex::Real const x_pos_offset = GetPosOffset(0, geom, slice_fab.box()); + const amrex::Real y_pos_offset = GetPosOffset(1, geom, slice_fab.box()); + + // loading the data + const auto ptd = pti.GetParticleTile().getParticleTileData(); + const bool can_ionize = m_can_ionize; + + const amrex::Real laser_norm = (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass) + * (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass); + + // Use OMP ParallelFor to use multiple threads when running on CPU + omp::ParallelFor( + amrex::TypeList< + amrex::CompileTimeOptions<0, 1, 2, 3> + >{}, { + Hipace::m_depos_order_xy + }, + pti.numParticles(), + [=] AMREX_GPU_DEVICE (int ip, auto depos_order) { + // only push plasma particles on their according MR level + if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; + + amrex::Real laser_norm_ion = laser_norm; + if (can_ionize) { + const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); + laser_norm_ion *= p_ion_lev * p_ion_lev; + } + + const amrex::Real xp = ptd.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd.rdata(PlasmaIdx::y)[ip]; + + amrex::Real Aabssqp = 0._rt; + doLaserGatherShapeN(xp, yp, + Aabssqp, slice_arr, aabs_comp, + dx_inv, dy_inv, x_pos_offset, y_pos_offset); + Aabssqp *= laser_norm_ion; + + ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; + }); + } +} + void PlasmaParticleContainer::InSituComputeDiags (int islice) { diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index ae9551c549..7e15dc6dc9 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -277,18 +277,6 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, } } // loop over subcycles - if (use_laser.value) { - Aabssqp = 0._rt; - - doLaserGatherShapeN(xp, yp, - Aabssqp, slice_arr, aabs_comp, - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - - Aabssqp *= laser_norm_ion; - - ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; - } - ptd.pos(0, ip) = xp; ptd.pos(1, ip) = yp; From f013fa87657c30d278bfa04a7b3010b1e41bdec8 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 17 Apr 2026 10:14:16 +0200 Subject: [PATCH 14/25] tests with aabssqdx --- .../deposition/ExplicitDeposition.cpp | 49 ++++++++++++------- .../plasma/PlasmaParticleContainer.H | 4 +- .../plasma/PlasmaParticleContainer.cpp | 18 ++++++- .../plasma/PlasmaParticleContainerInit.cpp | 4 ++ 4 files changed, 56 insertions(+), 19 deletions(-) diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index 5ab6f51515..9b4e451dcb 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -167,8 +167,23 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; amrex::Real Aabssqp = 0._rt; + amrex::Real AabssqDxp = 0._rt; + amrex::Real AabssqDyp = 0._rt; + // if (use_laser) { + // // Its important that Aabssqp is first fully gathered and not used + // // directly per cell like AabssqDxp and AabssqDyp + // doLaserGatherShapeN(xp, yp, Aabssqp, AabssqDxp, AabssqDyp, + // arr, cache_idx[4], + // dx_inv, dy_inv, x_pos_offset, y_pos_offset); + // Aabssqp *= a_laser_fac * q_mass_ratio * q_mass_ratio; + // AabssqDxp *= a_laser_fac * a_clight; + // AabssqDyp *= a_laser_fac * a_clight; + // } + if (use_laser) { Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; + AabssqDxp = ptd.rdata(PlasmaIdx::aabssqdx)[ip]; + AabssqDyp = ptd.rdata(PlasmaIdx::aabssqdy)[ip]; } // calculate gamma/psi for plasma particles @@ -195,23 +210,23 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real ExmBy_v = arr(i, j, cache_idx[2]); const amrex::Real EypBx_v = arr(i, j, cache_idx[3]); - amrex::Real AabssqDxp = 0._rt; - amrex::Real AabssqDyp = 0._rt; - // Rename variables for NVCC lambda capture to work - [[maybe_unused]] auto clight = a_clight; - [[maybe_unused]] auto laser_fac = a_laser_fac; - if constexpr (use_laser) { - // avoid going outside of domain - if (shape_x * shape_y != 0._rt) { - // need extra cells for gathering the laser - const amrex::Real xp1y00 = arr(i+1, j , cache_idx[4]); - const amrex::Real xm1y00 = arr(i-1, j , cache_idx[4]); - const amrex::Real x00yp1 = arr(i , j+1, cache_idx[4]); - const amrex::Real x00ym1 = arr(i , j-1, cache_idx[4]); - AabssqDxp = (xp1y00-xm1y00) * 0.5_rt * dx_inv * laser_fac * clight; - AabssqDyp = (x00yp1-x00ym1) * 0.5_rt * dy_inv * laser_fac * clight; - } - } + // amrex::Real AabssqDxp = 0._rt; + // amrex::Real AabssqDyp = 0._rt; + // // Rename variables for NVCC lambda capture to work + // [[maybe_unused]] auto clight = a_clight; + // [[maybe_unused]] auto laser_fac = a_laser_fac; + // if constexpr (use_laser) { + // // avoid going outside of domain + // if (shape_x * shape_y != 0._rt) { + // // need extra cells for gathering the laser + // const amrex::Real xp1y00 = arr(i+1, j , cache_idx[4]); + // const amrex::Real xm1y00 = arr(i-1, j , cache_idx[4]); + // const amrex::Real x00yp1 = arr(i , j+1, cache_idx[4]); + // const amrex::Real x00ym1 = arr(i , j-1, cache_idx[4]); + // AabssqDxp = (xp1y00-xm1y00) * 0.5_rt * dx_inv * laser_fac * clight; + // AabssqDyp = (x00yp1-x00ym1) * 0.5_rt * dy_inv * laser_fac * clight; + // } + // } amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), charge_density_mu0 * ( - shape_x * shape_y * ( diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 735492cc66..3db92780cd 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -108,7 +108,9 @@ struct PlasmaIdx // optional: aabssq=real_nattribs, // Laser value at the particle position - real_nattribs_laser=1, + aabssqdx, + aabssqdy, + real_nattribs_laser=3, x_prev=0, y_prev, // positions on the last non-temp slice real_nattribs_temp_slice, diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index b42f61a3af..9ae9e58f61 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -584,6 +584,12 @@ IonizationModule (const int lev, ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = 0._rt; ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = 1._rt; + if (comps.use_laser) { + ptd_elec.rdata(PlasmaIdx::aabssq)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdx)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdy)[ip] = 0._rt; + } + if (comps.use_temp_slice) { ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; @@ -849,6 +855,12 @@ LaserIonization (const int islice, ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = uy; ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = psi; + if (comps.use_laser) { + ptd_elec.rdata(PlasmaIdx::aabssq)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdx)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdy)[ip] = 0._rt; + } + if (comps.use_temp_slice) { ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; @@ -927,12 +939,16 @@ GatherLaser (const int lev, const amrex::Real yp = ptd.rdata(PlasmaIdx::y)[ip]; amrex::Real Aabssqp = 0._rt; + amrex::Real AabssqDxp = 0._rt; + amrex::Real AabssqDyp = 0._rt; doLaserGatherShapeN(xp, yp, - Aabssqp, slice_arr, aabs_comp, + Aabssqp, AabssqDxp, AabssqDyp, slice_arr, aabs_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); Aabssqp *= laser_norm_ion; ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; + ptd.rdata(PlasmaIdx::aabssqdx)[ip] = AabssqDxp; + ptd.rdata(PlasmaIdx::aabssqdy)[ip] = AabssqDyp; }); } } diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index 9d171287e3..da7c650ad4 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -318,6 +318,8 @@ InitParticles (const amrex::RealVect& a_u_std, if (comps.use_laser) { ptd.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + ptd.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; + ptd.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -377,6 +379,8 @@ InitParticles (const amrex::RealVect& a_u_std, if (comps.use_laser) { ptd.rdata(PlasmaIdx::aabssq)[midx] = 0._rt; + ptd.rdata(PlasmaIdx::aabssqdx)[midx] = 0._rt; + ptd.rdata(PlasmaIdx::aabssqdy)[midx] = 0._rt; } if (comps.use_temp_slice) { From 0054da8d92ba89978c2506df78e52bab9ccc8a6e Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 17 Apr 2026 14:43:04 +0200 Subject: [PATCH 15/25] fix prefactor --- .../plasma/PlasmaParticleContainer.cpp | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 9ae9e58f61..d81f5a28df 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -585,9 +585,9 @@ IonizationModule (const int lev, ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = 1._rt; if (comps.use_laser) { - ptd_elec.rdata(PlasmaIdx::aabssq)[ip] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdx)[ip] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdy)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -856,9 +856,9 @@ LaserIonization (const int islice, ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = psi; if (comps.use_laser) { - ptd_elec.rdata(PlasmaIdx::aabssq)[ip] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdx)[ip] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdy)[ip] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -914,8 +914,10 @@ GatherLaser (const int lev, const auto ptd = pti.GetParticleTile().getParticleTileData(); const bool can_ionize = m_can_ionize; - const amrex::Real laser_norm = (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass) + const amrex::Real laser_norm_qm = (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass) * (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass); + const amrex::Real laser_norm_c = phys_const.c * (phys_const.m_e/phys_const.q_e) + * (phys_const.m_e/phys_const.q_e); // Use OMP ParallelFor to use multiple threads when running on CPU omp::ParallelFor( @@ -929,10 +931,10 @@ GatherLaser (const int lev, // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; - amrex::Real laser_norm_ion = laser_norm; + amrex::Real laser_norm_qm_ion = laser_norm_qm; if (can_ionize) { const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); - laser_norm_ion *= p_ion_lev * p_ion_lev; + laser_norm_qm_ion *= p_ion_lev * p_ion_lev; } const amrex::Real xp = ptd.rdata(PlasmaIdx::x)[ip]; @@ -944,7 +946,9 @@ GatherLaser (const int lev, doLaserGatherShapeN(xp, yp, Aabssqp, AabssqDxp, AabssqDyp, slice_arr, aabs_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + Aabssqp *= laser_norm_qm_ion; + AabssqDxp *= laser_norm_c; + AabssqDyp *= laser_norm_c; ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; ptd.rdata(PlasmaIdx::aabssqdx)[ip] = AabssqDxp; From 4d4faec1bc53dc5de06728e21d2d3b5764dd3dd2 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 17 Apr 2026 15:15:58 +0200 Subject: [PATCH 16/25] prevent number of components accumulating each time step --- .../plasma/PlasmaParticleContainer.H | 1 + .../plasma/PlasmaParticleContainer.cpp | 68 ++++++++++--------- 2 files changed, 37 insertions(+), 32 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 3db92780cd..5263b85156 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -308,6 +308,7 @@ public: bool m_use_ab5_push = false; PlasmaComps m_comps; bool m_is_on_temp_slice = false; + bool m_components_allocated = false; // ionization: diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index d81f5a28df..983a61e577 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -182,52 +182,56 @@ PlasmaParticleContainer::ReadParameters () void PlasmaParticleContainer::InitData (const amrex::Vector& geom3d) { - SetArena(amrex::The_Arena()); + if (!m_components_allocated) { + m_components_allocated = true; - m_comps.use_laser = Hipace::m_use_laser; - m_comps.use_temp_slice = Hipace::GetInstance().m_multi_beam.AnySpeciesSalame() - || !Hipace::m_explicit; - m_comps.use_ab5_push = m_use_ab5_push; - m_comps.use_ion_level = m_can_ionize; + SetArena(amrex::The_Arena()); - int num_real_comps = 0; + m_comps.use_laser = Hipace::m_use_laser; + m_comps.use_temp_slice = Hipace::GetInstance().m_multi_beam.AnySpeciesSalame() + || !Hipace::m_explicit; + m_comps.use_ab5_push = m_use_ab5_push; + m_comps.use_ion_level = m_can_ionize; - for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { - AddRealComp(); - ++num_real_comps; - } + int num_real_comps = 0; - if (m_comps.use_laser) { - for (int j = 0; j < PlasmaIdx::real_nattribs_laser; ++j) { + for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { AddRealComp(); ++num_real_comps; } - } - if (m_comps.use_temp_slice) { - m_comps.offset_temp = num_real_comps; - for (int j = 0; j < PlasmaIdx::real_nattribs_temp_slice; ++j) { - AddRealComp(); - ++num_real_comps; + if (m_comps.use_laser) { + for (int j = 0; j < PlasmaIdx::real_nattribs_laser; ++j) { + AddRealComp(); + ++num_real_comps; + } } - } - if (m_comps.use_ab5_push) { - m_comps.offset_ab5 = num_real_comps; - for (int j = 0; j < PlasmaIdx::real_nattribs_ab5_push; ++j) { - AddRealComp(); - ++num_real_comps; + if (m_comps.use_temp_slice) { + m_comps.offset_temp = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_temp_slice; ++j) { + AddRealComp(); + ++num_real_comps; + } } - } - if (m_comps.use_ion_level) { - for (int j = 0; j < PlasmaIdx::int_nattribs_ion_level; ++j) { - AddIntComp(); + if (m_comps.use_ab5_push) { + m_comps.offset_ab5 = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_ab5_push; ++j) { + AddRealComp(); + ++num_real_comps; + } } - } - reserveData(); - resizeData(); + if (m_comps.use_ion_level) { + for (int j = 0; j < PlasmaIdx::int_nattribs_ion_level; ++j) { + AddIntComp(); + } + } + + reserveData(); + resizeData(); + } if (!m_read_fine_patch) { m_read_fine_patch = true; From 8d00c3f7bc07694a720122b6cfea4f0e6a6a0a19 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 17 Apr 2026 16:52:55 +0200 Subject: [PATCH 17/25] Add ResetPositions --- src/Hipace.cpp | 2 ++ .../deposition/ExplicitDeposition.cpp | 1 - src/particles/plasma/MultiPlasma.H | 3 ++ src/particles/plasma/MultiPlasma.cpp | 8 +++++ .../plasma/PlasmaParticleContainer.H | 2 ++ .../plasma/PlasmaParticleContainer.cpp | 29 +++++++++++++++++++ src/salame/Salame.cpp | 10 ++++--- 7 files changed, 50 insertions(+), 5 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 9720e1d0c0..217c58f253 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -1176,6 +1176,8 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice, const int current_N relative_Bfield_error_prev_iter = relative_Bfield_error; } // end of predictor corrector loop + m_multi_plasma.ResetPositions(); + if (relative_Bfield_error > 10. && m_predcorr_B_error_tolerance > 0.) { amrex::Print() << "WARNING: Predictor corrector loop may have diverged!\n" diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index 9b4e451dcb..65855b9c69 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -52,7 +52,6 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real a_clight = pc.c; const amrex::Real clight_inv = 1._rt/pc.c; // The laser a0 is always normalized - const amrex::Real a_laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); const amrex::Real charge_invvol_mu0 = plasma.m_charge * invvol * pc.mu0; const amrex::Real charge_mass_ratio = plasma.m_charge / plasma.m_mass; diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index cb9438e7c2..fde414b630 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -93,6 +93,9 @@ public: const Fields & fields, amrex::Vector const& gm, bool temp_slice, int lev, int const current_N_level); + /** \brief Reset positions after push to temp slice */ + void ResetPositions (); + /** \brief Gather Laser field to the particles * * \param[in] lev MR level diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index 691c50a0ea..f6816a21cd 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -133,6 +133,14 @@ MultiPlasma::AdvanceParticles ( } } +void +MultiPlasma::ResetPositions () +{ + for (int i=0; iExplicitMGSolveBxBy(lev, WhichSlice::This); } - } - if (hipace->m_N_level > 1) { - // tag to prev slice for push - hipace->m_multi_plasma.TagByLevel(current_N_level, hipace->m_3D_geom, true); + if (hipace->m_N_level > 1) { + // tag to prev slice for next push + hipace->m_multi_plasma.TagByLevel(current_N_level, hipace->m_3D_geom, true); + } } + + hipace->m_multi_plasma.ResetPositions(); } From 0a015ad78b1521b1e05bd7479a8e6bf987d58424 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Mon, 6 Jul 2026 11:08:48 +0200 Subject: [PATCH 18/25] fix --- src/particles/plasma/PlasmaParticleContainer.H | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 956fb49a2a..82d8ebfa30 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -310,6 +310,7 @@ public: bool m_use_ab5_push = false; PlasmaComps m_comps; bool m_is_on_temp_slice = false; + /** If the runtime attributes of the AMReX particle container were already set up */ bool m_components_allocated = false; // ionization: @@ -361,8 +362,6 @@ 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 */ From 5b832913ce30c419435a6e23220b0153cfcdc89e Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 11:06:33 +0200 Subject: [PATCH 19/25] dont store aabssqdx and aabssqdy per particle --- .../deposition/ExplicitDeposition.cpp | 52 +++++++------------ src/particles/particles_utils/FieldGather.H | 3 +- .../plasma/PlasmaParticleContainer.H | 4 +- .../plasma/PlasmaParticleContainer.cpp | 12 +---- .../plasma/PlasmaParticleContainerInit.cpp | 4 -- 5 files changed, 23 insertions(+), 52 deletions(-) diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index 65855b9c69..4a7244f9ec 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -52,6 +52,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real a_clight = pc.c; const amrex::Real clight_inv = 1._rt/pc.c; // The laser a0 is always normalized + const amrex::Real a_laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); const amrex::Real charge_invvol_mu0 = plasma.m_charge * invvol * pc.mu0; const amrex::Real charge_mass_ratio = plasma.m_charge / plasma.m_mass; @@ -166,23 +167,10 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; amrex::Real Aabssqp = 0._rt; - amrex::Real AabssqDxp = 0._rt; - amrex::Real AabssqDyp = 0._rt; - // if (use_laser) { - // // Its important that Aabssqp is first fully gathered and not used - // // directly per cell like AabssqDxp and AabssqDyp - // doLaserGatherShapeN(xp, yp, Aabssqp, AabssqDxp, AabssqDyp, - // arr, cache_idx[4], - // dx_inv, dy_inv, x_pos_offset, y_pos_offset); - // Aabssqp *= a_laser_fac * q_mass_ratio * q_mass_ratio; - // AabssqDxp *= a_laser_fac * a_clight; - // AabssqDyp *= a_laser_fac * a_clight; - // } - if (use_laser) { + // Its important that Aabssqp is first fully gathered and not used + // directly per cell like AabssqDxp and AabssqDyp Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; - AabssqDxp = ptd.rdata(PlasmaIdx::aabssqdx)[ip]; - AabssqDyp = ptd.rdata(PlasmaIdx::aabssqdy)[ip]; } // calculate gamma/psi for plasma particles @@ -209,23 +197,23 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real ExmBy_v = arr(i, j, cache_idx[2]); const amrex::Real EypBx_v = arr(i, j, cache_idx[3]); - // amrex::Real AabssqDxp = 0._rt; - // amrex::Real AabssqDyp = 0._rt; - // // Rename variables for NVCC lambda capture to work - // [[maybe_unused]] auto clight = a_clight; - // [[maybe_unused]] auto laser_fac = a_laser_fac; - // if constexpr (use_laser) { - // // avoid going outside of domain - // if (shape_x * shape_y != 0._rt) { - // // need extra cells for gathering the laser - // const amrex::Real xp1y00 = arr(i+1, j , cache_idx[4]); - // const amrex::Real xm1y00 = arr(i-1, j , cache_idx[4]); - // const amrex::Real x00yp1 = arr(i , j+1, cache_idx[4]); - // const amrex::Real x00ym1 = arr(i , j-1, cache_idx[4]); - // AabssqDxp = (xp1y00-xm1y00) * 0.5_rt * dx_inv * laser_fac * clight; - // AabssqDyp = (x00yp1-x00ym1) * 0.5_rt * dy_inv * laser_fac * clight; - // } - // } + amrex::Real AabssqDxp = 0._rt; + amrex::Real AabssqDyp = 0._rt; + // Rename variables for NVCC lambda capture to work + [[maybe_unused]] auto clight = a_clight; + [[maybe_unused]] auto laser_fac = a_laser_fac; + if constexpr (use_laser) { + // avoid going outside of domain + if (shape_x * shape_y != 0._rt) { + // need extra cells for gathering the laser + const amrex::Real xp1y00 = arr(i+1, j , cache_idx[4]); + const amrex::Real xm1y00 = arr(i-1, j , cache_idx[4]); + const amrex::Real x00yp1 = arr(i , j+1, cache_idx[4]); + const amrex::Real x00ym1 = arr(i , j-1, cache_idx[4]); + AabssqDxp = (xp1y00-xm1y00) * 0.5_rt * dx_inv * laser_fac * clight; + AabssqDyp = (x00yp1-x00ym1) * 0.5_rt * dy_inv * laser_fac * clight; + } + } amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), charge_density_mu0 * ( - shape_x * shape_y * ( diff --git a/src/particles/particles_utils/FieldGather.H b/src/particles/particles_utils/FieldGather.H index 40d9c39a9a..462aa96558 100644 --- a/src/particles/particles_utils/FieldGather.H +++ b/src/particles/particles_utils/FieldGather.H @@ -251,8 +251,7 @@ void doLaserGatherShapeN (const amrex::Real xp, auto [shape_y, j] = shape_factor(y, iy); auto [shape_x, i] = shape_factor(x, ix); - const amrex::Real x00y00 = slice_arr(i, j, aabs_comp); - Aabssqp += shape_x * shape_y * x00y00; + Aabssqp += shape_x * shape_y * slice_arr(i, j, aabs_comp); } } } diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 82d8ebfa30..377db39115 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -108,9 +108,7 @@ struct PlasmaIdx // optional: aabssq=real_nattribs, // Laser value at the particle position - aabssqdx, - aabssqdy, - real_nattribs_laser=3, + real_nattribs_laser=1, x_prev=0, y_prev, // positions on the last non-temp slice real_nattribs_temp_slice, diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index b0f156d5f1..b94de945f5 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -590,8 +590,6 @@ IonizationModule (const int lev, if (comps.use_laser) { ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -861,8 +859,6 @@ LaserIonization (const int islice, if (comps.use_laser) { ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; - ptd_elec.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -974,18 +970,12 @@ GatherLaser (const int lev, const amrex::Real yp = ptd.rdata(PlasmaIdx::y)[ip]; amrex::Real Aabssqp = 0._rt; - amrex::Real AabssqDxp = 0._rt; - amrex::Real AabssqDyp = 0._rt; doLaserGatherShapeN(xp, yp, - Aabssqp, AabssqDxp, AabssqDyp, slice_arr, aabs_comp, + Aabssqp, slice_arr, aabs_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); Aabssqp *= laser_norm_qm_ion; - AabssqDxp *= laser_norm_c; - AabssqDyp *= laser_norm_c; ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; - ptd.rdata(PlasmaIdx::aabssqdx)[ip] = AabssqDxp; - ptd.rdata(PlasmaIdx::aabssqdy)[ip] = AabssqDyp; }); } } diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index da7c650ad4..9d171287e3 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -318,8 +318,6 @@ InitParticles (const amrex::RealVect& a_u_std, if (comps.use_laser) { ptd.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; - ptd.rdata(PlasmaIdx::aabssqdx)[pidx] = 0._rt; - ptd.rdata(PlasmaIdx::aabssqdy)[pidx] = 0._rt; } if (comps.use_temp_slice) { @@ -379,8 +377,6 @@ InitParticles (const amrex::RealVect& a_u_std, if (comps.use_laser) { ptd.rdata(PlasmaIdx::aabssq)[midx] = 0._rt; - ptd.rdata(PlasmaIdx::aabssqdx)[midx] = 0._rt; - ptd.rdata(PlasmaIdx::aabssqdy)[midx] = 0._rt; } if (comps.use_temp_slice) { From c0a7e084b9eb9400afd49e5b0076bc16a1ab5c8d Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 11:29:43 +0200 Subject: [PATCH 20/25] fix ion level in Collisions --- src/particles/collisions/CoulombCollision.cpp | 14 +++++++++----- src/particles/plasma/PlasmaParticleContainer.cpp | 2 -- src/particles/pusher/PlasmaParticleAdvance.cpp | 6 +++--- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index 5d1466f73a..2e1c65cdfc 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -93,7 +93,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( 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(); + const int* const ion_lev1 = species1.m_can_ionize ? + ptile1.GetIntData(PlasmaIdx::ion_lev).data() : nullptr; PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); amrex::Real q1 = species1.GetCharge(); @@ -159,7 +160,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( 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(); + const int* const ion_lev1 = species1.m_can_ionize ? + ptile1.GetIntData(PlasmaIdx::ion_lev).data() : nullptr; PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); amrex::Real q1 = species1.GetCharge(); @@ -172,7 +174,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( 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(); + const int* const ion_lev2 = species2.m_can_ionize ? + ptile2.GetIntData(PlasmaIdx::ion_lev).data() : nullptr; PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); amrex::Real q2 = species2.GetCharge(); @@ -282,7 +285,8 @@ CoulombCollision::doBeamPlasmaCoulombCollision ( 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(); + const int* const ion_lev2 = species2.m_can_ionize ? + ptile2.GetIntData(PlasmaIdx::ion_lev).data() : nullptr; PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); amrex::Real q2 = species2.GetCharge(); @@ -332,7 +336,7 @@ CoulombCollision::doBeamPlasmaCoulombCollision ( ElasticCollisionPerez( cell_start1, cell_stop1, cell_start2, cell_stop2, indices1, indices2, - ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev2, ion_lev2, // passing ion_lev2 for beam particles, will never be used + ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, nullptr, ion_lev2, // passing nullptr for beam particles, will never be used q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2, dt, CoulombLog, inv_dV, clight, inv_c_SI, inv_c2_SI, normalized_units, background_density_SI, false, true, engine ); diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index b94de945f5..901d9421e1 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -945,8 +945,6 @@ GatherLaser (const int lev, const amrex::Real laser_norm_qm = (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass) * (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass); - const amrex::Real laser_norm_c = phys_const.c * (phys_const.m_e/phys_const.q_e) - * (phys_const.m_e/phys_const.q_e); // Use OMP ParallelFor to use multiple threads when running on CPU omp::ParallelFor( diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 7e15dc6dc9..02e6b4f74e 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -80,8 +80,8 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, const amrex::Real charge_mass_clight_ratio = plasma.m_charge/(plasma.m_mass * phys_const.c); const auto comps = plasma.m_comps; - const bool read_from_temp = plasma.m_is_on_temp_slice; - AMREX_ALWAYS_ASSERT(!(read_from_temp || temp_slice) || comps.use_temp_slice); + const bool read_from_prev = plasma.m_is_on_temp_slice; + AMREX_ALWAYS_ASSERT(!(read_from_prev || temp_slice) || comps.use_temp_slice); // Use OMP ParallelFor to use multiple threads when running on CPU omp::ParallelFor( @@ -118,7 +118,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::Real xp = 0._rt; amrex::Real yp = 0._rt; - if (!read_from_temp) { + if (!read_from_prev) { xp = ptd.rdata(PlasmaIdx::x)[ip]; yp = ptd.rdata(PlasmaIdx::y)[ip]; if (temp_slice) { From 99a056dfce910c2832a164452a84abb0df20c57c Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 12:57:38 +0200 Subject: [PATCH 21/25] test CI without Aabssqp in LaserIonization --- src/particles/plasma/PlasmaParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 901d9421e1..5b25bc361f 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -722,9 +722,9 @@ LaserIonization (const int islice, const amrex::Real uy = ptd_ion.rdata(PlasmaIdx::uy_half_step)[ip]; const amrex::Real psi = ptd_ion.rdata(PlasmaIdx::psi_half_step)[ip]; - const amrex::Real Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; + // const amrex::Real Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; // Compute probability of ionization p - const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, Aabssqp); + const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, 0._rt); // TODO Add Aabssqp const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; if (ion_lev_loc >= max_ion_lev) { return; From 648d8a04d58ca237f96efa671a6a594e65054db8 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 13:18:18 +0200 Subject: [PATCH 22/25] fix salame --- src/salame/Salame.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index 027c48d9ec..a0ee24d055 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -302,7 +302,7 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) const amrex::Real charge_mass_c_ratio = plasma.m_charge / (plasma.m_mass * get_phys_const().c); - const bool can_ionize = plasma.m_can_ionize; + auto comps = plasma.m_comps; omp::ParallelFor( amrex::TypeList>{}, @@ -312,8 +312,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; - const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev + comps.use_temp_slice)[ip]; + const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev + comps.use_temp_slice)[ip]; amrex::Real Bxp = 0._rt; amrex::Real Byp = 0._rt; @@ -322,7 +322,7 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) doBxByGatherShapeN(xp, yp, Bxp, Byp, slice_arr, bx_comp, by_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); - const amrex::Real q_m_c_ratio = can_ionize ? + const amrex::Real q_m_c_ratio = comps.use_ion_level ? ptd.idata(PlasmaIdx::ion_lev)[ip] * charge_mass_c_ratio : charge_mass_c_ratio; From 76b9d560a761fa1c506e3a92c0dd5cb127112704 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 13:54:35 +0200 Subject: [PATCH 23/25] try fix --- src/particles/plasma/PlasmaParticleContainer.cpp | 2 +- src/salame/Salame.cpp | 4 ++-- tests/beam_in_vacuum.SI.Serial.sh | 1 + tests/beam_in_vacuum.normalized.Serial.sh | 1 + tests/blowout_wake.Serial.sh | 1 + 5 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 5b25bc361f..0022d804ee 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -858,7 +858,7 @@ LaserIonization (const int islice, ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = psi; if (comps.use_laser) { - ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = amrex::abs(A*A); } if (comps.use_temp_slice) { diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index a0ee24d055..3648f54edf 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -312,8 +312,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; - const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev + comps.use_temp_slice)[ip]; - const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev + comps.use_temp_slice)[ip]; + const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; amrex::Real Bxp = 0._rt; amrex::Real Byp = 0._rt; diff --git a/tests/beam_in_vacuum.SI.Serial.sh b/tests/beam_in_vacuum.SI.Serial.sh index d667aa783a..c9bd631385 100755 --- a/tests/beam_in_vacuum.SI.Serial.sh +++ b/tests/beam_in_vacuum.SI.Serial.sh @@ -30,6 +30,7 @@ $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ hipace.tile_size = 8 \ hipace.depos_order_xy=0 \ diagnostic.field_data = all rho \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the result with theory diff --git a/tests/beam_in_vacuum.normalized.Serial.sh b/tests/beam_in_vacuum.normalized.Serial.sh index bc339714a5..8ae53d8154 100755 --- a/tests/beam_in_vacuum.normalized.Serial.sh +++ b/tests/beam_in_vacuum.normalized.Serial.sh @@ -30,6 +30,7 @@ $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \ hipace.tile_size = 8 \ hipace.depos_order_xy=0 \ diagnostic.field_data = all rho \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the result with theory diff --git a/tests/blowout_wake.Serial.sh b/tests/blowout_wake.Serial.sh index 0c9d46b408..94964bfe65 100755 --- a/tests/blowout_wake.Serial.sh +++ b/tests/blowout_wake.Serial.sh @@ -28,6 +28,7 @@ HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests # Run the simulation $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \ hipace.tile_size = 8 \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the results with checksum benchmark From 678124f2d59af66bc7d322265bc23f9cfb66d55e Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 15:18:12 +0200 Subject: [PATCH 24/25] Fix AB5 and reset CI because TemperatureDeposition now uses different shape factor for laser gather --- src/particles/plasma/PlasmaParticleContainer.cpp | 2 +- src/particles/pusher/PlasmaParticleAdvance.cpp | 12 +++++++----- src/salame/Salame.cpp | 1 + .../benchmarks_json/laser_ionization.1Rank.json | 10 +++++----- 4 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 0022d804ee..b81a2b7c04 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -790,7 +790,7 @@ LaserIonization (const int islice, [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, auto depos_order_xy) { - if(p_ion_mask[ip] != 0) { + if (p_ion_mask[ip] != 0) { const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 02e6b4f74e..cf246dc366 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -242,11 +242,12 @@ 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 + 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 int offset = ab5_permutation + comps.offset_ab5; + ptd.rdata(PlasmaIdx::Fx1 + offset)[ip] = ux * psi_inv; + ptd.rdata(PlasmaIdx::Fy1 + offset)[ip] = uy * psi_inv; + ptd.rdata(PlasmaIdx::Fux1 + offset)[ip] = dz_ux; + ptd.rdata(PlasmaIdx::Fuy1 + offset)[ip] = dz_uy; + ptd.rdata(PlasmaIdx::Fpsi1 + offset)[ip] = dz_psi; const amrex::Real ab5_coeffs[5] = { ( 1901._rt / 720._rt ) * dz, // a1 times dz @@ -262,6 +263,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, if (p >= 5) { p -= 5; } + p += comps.offset_ab5; 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]; diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index 3648f54edf..ad9c8b54ce 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -303,6 +303,7 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) const amrex::Real charge_mass_c_ratio = plasma.m_charge / (plasma.m_mass * get_phys_const().c); auto comps = plasma.m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice); omp::ParallelFor( amrex::TypeList>{}, diff --git a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json index e2b263b9d6..ce0d5e1cb6 100644 --- a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json +++ b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json @@ -31,11 +31,11 @@ "uy^2_ion": 0.0, "uy_elec": 0.0, "uy_ion": 0.0, - "uz^2_elec": 1.6568753834647e-06, - "uz^2_ion": 1.5025085926483e-20, - "uz_elec": 0.15044101093153, - "uz_ion": 4.8989912126984e-09, - "w_elec": 483046484939.32, + "uz^2_elec": 1.6613402997284e-06, + "uz^2_ion": 1.5330025597749e-20, + "uz_elec": 0.15061530702838, + "uz_ion": 4.9506114763176e-09, + "w_elec": 483046486891.44, "w_ion": 1400000000000.1, "|a^2|": 0.46877495720552 } From 9fcf70752f565f6251ac6671cd084ef3951cac5e Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Thu, 9 Jul 2026 16:26:20 +0200 Subject: [PATCH 25/25] Use old version of ExplicitDeposition --- src/particles/deposition/ExplicitDeposition.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index 4a7244f9ec..f281dbec29 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -52,7 +52,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const amrex::Real a_clight = pc.c; const amrex::Real clight_inv = 1._rt/pc.c; // The laser a0 is always normalized - const amrex::Real a_laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); + const amrex::Real laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e); const amrex::Real charge_invvol_mu0 = plasma.m_charge * invvol * pc.mu0; const amrex::Real charge_mass_ratio = plasma.m_charge / plasma.m_mass; @@ -170,7 +170,9 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, if (use_laser) { // Its important that Aabssqp is first fully gathered and not used // directly per cell like AabssqDxp and AabssqDyp - Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; + doLaserGatherShapeN(xp, yp, Aabssqp, arr, cache_idx[4], + dx_inv, dy_inv, x_pos_offset, y_pos_offset); + Aabssqp *= laser_fac * q_mass_ratio * q_mass_ratio; } // calculate gamma/psi for plasma particles @@ -201,7 +203,6 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, amrex::Real AabssqDyp = 0._rt; // Rename variables for NVCC lambda capture to work [[maybe_unused]] auto clight = a_clight; - [[maybe_unused]] auto laser_fac = a_laser_fac; if constexpr (use_laser) { // avoid going outside of domain if (shape_x * shape_y != 0._rt) {