From b9ee9a3284f91a5cd16dc8920c39c03baa3fe23c Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 17 Aug 2023 01:23:10 -0400 Subject: [PATCH] Local Move of Lost Particles --- src/ImpactX.cpp | 15 +- src/particles/CollectLost.cpp | 157 +++++++++++++++++---- src/particles/ImpactXParticleContainer.cpp | 5 - 3 files changed, 141 insertions(+), 36 deletions(-) diff --git a/src/ImpactX.cpp b/src/ImpactX.cpp index 281218fa6..c67f1fbb7 100644 --- a/src/ImpactX.cpp +++ b/src/ImpactX.cpp @@ -33,7 +33,8 @@ namespace impactx { ImpactX::ImpactX () : AmrCore(initialization::init_amr_core()), - m_particle_container(std::make_unique(this)) + m_particle_container(std::make_unique(this)), + m_particles_lost(std::make_unique(this)) { // todo: if amr.n_cells is provided, overwrite/redefine AmrCore object @@ -80,6 +81,18 @@ namespace impactx // init blocks / grids & MultiFabs AmrCore::InitFromScratch(0.0); amrex::Print() << "boxArray(0) " << boxArray(0) << std::endl; + + // alloc particle containers + // the lost particles have an extra runtime attribute: s when it was lost + bool comm = true; + m_particles_lost->AddRealComp(comm); + + // have to resize here, not in the constructor because grids have not + // been built when constructor was called. + m_particle_container->reserveData(); + m_particle_container->resizeData(); + m_particles_lost->reserveData(); + m_particles_lost->resizeData(); } void ImpactX::evolve () diff --git a/src/particles/CollectLost.cpp b/src/particles/CollectLost.cpp index 6252ce570..d3d03fa40 100644 --- a/src/particles/CollectLost.cpp +++ b/src/particles/CollectLost.cpp @@ -11,48 +11,145 @@ #include #include +#include +#include +#include namespace impactx { void collect_lost_particles (ImpactXParticleContainer& source, ImpactXParticleContainer& dest) { + BL_PROFILE("impactX::collect_lost_particles"); + using SrcData = ImpactXParticleContainer::ParticleTileType::ConstParticleTileDataType; + using DstData = ImpactXParticleContainer::ParticleTileType::ParticleTileDataType; - // copy all particles marked with a negative ID from source to destination - bool const local = true; - dest.copyParticles( - source, - [=] AMREX_GPU_HOST_DEVICE(const SrcData &src, int ip) { - return src.id(ip) < 0; - }, - local - ); + RefPart const ref_part = source.GetRefParticle(); + auto const s_lost = ref_part.s; - // flip IDs back to positive in destination - int const nLevel = dest.finestLevel(); + // have to resize here, not in the constructor because grids have not + // been built when constructor was called. + dest.reserveData(); + dest.resizeData(); + + // copy all particles marked with a negative ID from source to destination + int const nLevel = source.finestLevel(); for (int lev = 0; lev <= nLevel; ++lev) { // loop over all particle boxes using ParIt = ImpactXParticleContainer::iterator; + auto& plevel = source.GetParticles(lev); + + // TODO: is it safe to add OpenMP parallelism here? + for (ParIt pti(source, lev); pti.isValid(); ++pti) { + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + if (plevel.find(index) == plevel.end()) continue; + + auto& ptile_source = plevel.at(index); + auto const np = ptile_source.numParticles(); + if (np == 0) continue; // no particles in source tile + + // we will copy particles that were marked as lost, with a negative id + auto predicate = [] AMREX_GPU_HOST_DEVICE (const SrcData& src, int ip, const amrex::RandomEngine& /*engine*/) noexcept + { + return src.id(ip) < 0; + }; + + // FIXME This seems to segfault + auto& ptile_dest = dest.DefineAndReturnParticleTile( + lev, pti.index(), pti.LocalTileIndex()); + + // count how many particles we will copy + amrex::ReduceOps reduce_op; + amrex::ReduceData reduce_data(reduce_op); + { + const auto src_data = ptile_source.getConstParticleTileData(); + + const amrex::RandomEngine rng{}; // unused + reduce_op.eval(np, reduce_data, [=] AMREX_GPU_HOST_DEVICE (int ip) { + return predicate(src_data, ip, rng) ? 1 : 0; + }); + } + int const np_to_move = amrex::get<0>(reduce_data.value()); + if (np_to_move == 0) continue; // no particles to move from source tile + + // allocate memory in destination + int const dst_index = ptile_dest.numParticles(); + ptile_dest.resize(dst_index + np_to_move); + + // copy particles + // skipped in loop below: integer compile-time or runtime attributes + AMREX_ALWAYS_ASSERT(SrcData::NAI == 0); + AMREX_ALWAYS_ASSERT(ptile_source.NumRuntimeIntComps() == 0); + + // first runtime attribute in destination is s position when particle got lost + int const s_index = dest.NumRuntimeRealComps() - 1; + auto copy_and_mark_negative = [&s_index, &s_lost](DstData& dst, const SrcData& src, int src_ip, int dst_ip) noexcept + { + dst.m_aos[dst_ip] = src.m_aos[src_ip]; + + for (int j = 0; j < SrcData::NAR; ++j) + dst.m_rdata[j][dst_ip] = src.m_rdata[j][src_ip]; + for (int j = 0; j < src.m_num_runtime_real; ++j) + dst.m_runtime_rdata[j][dst_ip] = src.m_runtime_rdata[j][src_ip]; + + // unused: integer compile-time or runtime attributes + //for (int j = 0; j < SrcData::NAI; ++j) + // dst.m_idata[j][dst_ip] = src.m_idata[j][src_ip]; + //for (int j = 0; j < src.m_num_runtime_int; ++j) + // dst.m_runtime_idata[j][dst_ip] = src.m_runtime_idata[j][src_ip]; + + // flip id to positive in destination + dst.id(dst_ip) = amrex::Math::abs(dst.id(dst_ip)); + + // remember the current s of the ref particle when lost + dst.m_runtime_rdata[s_index][dst_ip] = s_lost; + }; + + amrex::filterAndTransformParticles( + ptile_dest, + ptile_source, + predicate, + copy_and_mark_negative, + 0, + dst_index + ); + + // remove particles with negative ids in source + { + int n_removed = 0; + auto ptile_src_data = ptile_source.getParticleTileData(); + for (int ip = 0; ip < np; ++ip) + { + if (ptile_source.id(ip) < 0) + n_removed++; + else + { + if (n_removed > 0) + { + // move down + int const new_index = ip - n_removed; + + ptile_src_data.m_aos[new_index] = ptile_src_data.m_aos[ip]; + + for (int j = 0; j < SrcData::NAR; ++j) + ptile_src_data.m_rdata[j][new_index] = ptile_src_data.m_rdata[j][ip]; + for (int j = 0; j < ptile_src_data.m_num_runtime_real; ++j) + ptile_src_data.m_runtime_rdata[j][new_index] = ptile_src_data.m_runtime_rdata[j][ip]; + + // unused: integer compile-time or runtime attributes + //for (int j = 0; j < SrcData::NAI; ++j) + // dst.m_idata[j][new_index] = src.m_idata[j][ip]; + //for (int j = 0; j < ptile_src_data.m_num_runtime_int; ++j) + // dst.m_runtime_idata[j][new_index] = src.m_runtime_idata[j][ip]; + } + } + } + AMREX_ALWAYS_ASSERT(np_to_move == n_removed); + ptile_source.resize(np - n_removed); + } -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (ParIt pti(dest, lev); pti.isValid(); ++pti) { - const int np = pti.numParticles(); - - // preparing access to particle data: AoS - using PType = ImpactXParticleContainer::ParticleType; - auto &aos = pti.GetArrayOfStructs(); - PType *AMREX_RESTRICT aos_ptr = aos().dataPtr(); - - amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) { - PType &p = aos_ptr[i]; - p.id() = -p.id(); - }); - } - } - - // TODO remove particles with negative ids in source + } // particle tile loop + } // lev } } // namespace impactx diff --git a/src/particles/ImpactXParticleContainer.cpp b/src/particles/ImpactXParticleContainer.cpp index ebdb63b2c..9973ea029 100644 --- a/src/particles/ImpactXParticleContainer.cpp +++ b/src/particles/ImpactXParticleContainer.cpp @@ -100,11 +100,6 @@ namespace impactx // number of particles to add int const np = x.size(); - // have to resize here, not in the constructor because grids have not - // been built when constructor was called. - reserveData(); - resizeData(); - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); /* Create a temporary tile to obtain data from simulation. This data