diff --git a/Src/Particle/AMReX_NeighborParticles.H b/Src/Particle/AMReX_NeighborParticles.H index ad5fedc25a8..8f94cffa967 100644 --- a/Src/Particle/AMReX_NeighborParticles.H +++ b/Src/Particle/AMReX_NeighborParticles.H @@ -233,6 +233,13 @@ public: /// void fillNeighbors (); + /// + /// Fill neighbor buffers, but only communicate particles within + /// physical distance search_radius of tile boundaries. + /// Reduces communication volume when search_radius << cell size. + /// + void fillNeighbors (amrex::Real search_radius); + /// /// This does an "inverse" fillNeighbors operation, meaning that it adds /// data from the ghost particles to the corresponding real ones. @@ -256,6 +263,14 @@ public: template void buildNeighborList (CheckPair const& check_pair, bool sort=false); + /// + /// Build a Neighbor List for each tile, using custom bin_size + /// instead of the mesh cell size. Useful when the interaction + /// radius is much smaller than dx. + /// + template + void buildNeighborList (CheckPair const& check_pair, amrex::Real bin_size, bool sort=false); + /// /// Build a Neighbor List for each tile /// @@ -407,6 +422,7 @@ protected: static constexpr int num_mask_comps = 3; //!< grid, tile, level size_t cdata_size; int m_num_neighbor_cells; + amrex::Real m_search_radius = amrex::Real(0); //!< physical search radius; 0 = use cell-based m_num_neighbor_cells Vector local_neighbors; Vector > mask_ptr; diff --git a/Src/Particle/AMReX_NeighborParticlesI.H b/Src/Particle/AMReX_NeighborParticlesI.H index 3a9f1a670ae..1617ab84488 100644 --- a/Src/Particle/AMReX_NeighborParticlesI.H +++ b/Src/Particle/AMReX_NeighborParticlesI.H @@ -280,7 +280,17 @@ NeighborParticleContainer { box.refine(ref_fac); } - box.grow(computeRefFac(0, lev)*m_num_neighbor_cells); + if (m_search_radius > amrex::Real(0)) { + const auto* dx_lev = this->Geom(lev).CellSize(); + IntVect grow_cells; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + grow_cells[d] = std::max(1, static_cast( + std::ceil(m_search_radius / dx_lev[d]))); + } + box.grow(grow_cells); + } else { + box.grow(computeRefFac(0, lev)*m_num_neighbor_cells); + } const Periodicity& periodicity = this->Geom(lev).periodicity(); const std::vector& pshifts = periodicity.shiftIntVect(); const BoxArray& ba = this->ParticleBoxArray(lev); @@ -550,6 +560,26 @@ getNeighborTags (Vector& tags, const ParticleType& p, Box shrink_box = pti.tilebox(); shrink_box.grow(-nGrow); + // Position-based early exit: when m_search_radius is set, skip particles + // that are farther than m_search_radius from all tile faces. + if (m_search_radius > amrex::Real(0)) { + const Box& tile_box = pti.tilebox(); + const int src_lev = src_tag.level; + const auto* dx = this->Geom(src_lev).CellSize(); + const auto* plo = this->Geom(src_lev).ProbLo(); + bool is_interior = true; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + amrex::Real face_lo = plo[d] + tile_box.smallEnd(d) * dx[d]; + amrex::Real face_hi = plo[d] + (tile_box.bigEnd(d) + 1) * dx[d]; + if ((p.pos(d) - face_lo) < m_search_radius || + (face_hi - p.pos(d)) < m_search_radius) { + is_interior = false; + break; + } + } + if (is_interior) { return; } + } + if (use_mask) { const BaseFab& mask = (*mask_ptr[src_tag.level])[src_tag.grid]; AMREX_ASSERT(this->finestLevel() == 0); @@ -604,14 +634,25 @@ getNeighborTags (Vector& tags, const ParticleType& p, Box tbx; for (int lev = 0; lev < this->numLevels(); ++lev) { - IntVect ref_fac = computeRefFac(0, lev); const Periodicity& periodicity = this->Geom(lev).periodicity(); const std::vector& pshifts = periodicity.shiftIntVect(); const BoxArray& ba = this->ParticleBoxArray(lev); const IntVect& iv = this->Index(p, lev); for (auto const& pshift : pshifts) { - Box pbox = amrex::grow(Box(iv, iv), ref_fac*nGrow) + pshift; + Box pbox; + if (m_search_radius > amrex::Real(0)) { + const auto* dx_lev = this->Geom(lev).CellSize(); + IntVect grow_cells; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + grow_cells[d] = std::max(1, static_cast( + std::ceil(m_search_radius / dx_lev[d]))); + } + pbox = amrex::grow(Box(iv, iv), grow_cells) + pshift; + } else { + IntVect ref_fac = computeRefFac(0, lev); + pbox = amrex::grow(Box(iv, iv), ref_fac*nGrow) + pshift; + } bool first_only = false; ba.intersections(pbox, isects, first_only, 0); for (const auto& isec : isects) @@ -638,6 +679,31 @@ template void NeighborParticleContainer ::fillNeighbors () { + m_search_radius = amrex::Real(0); +#ifdef AMREX_USE_GPU + fillNeighborsGPU(); +#else + fillNeighborsCPU(); +#endif + m_has_neighbors = true; +} + +template +void +NeighborParticleContainer +::fillNeighbors (amrex::Real search_radius) { + // Store search_radius; used by getNeighborTags and GetCommTagsBox + // for per-level nGrow = ceil(search_radius / dx[lev]) + m_search_radius = search_radius; + // m_num_neighbor_cells is still needed for the mask (single-level path) + { + const auto* dx = this->Geom(0).CellSize(); + amrex::Real dx_min = dx[0]; + for (int d = 1; d < AMREX_SPACEDIM; ++d) { + dx_min = std::min(dx_min, dx[d]); + } + m_num_neighbor_cells = std::max(1, static_cast(std::ceil(search_radius / dx_min))); + } #ifdef AMREX_USE_GPU fillNeighborsGPU(); #else @@ -718,7 +784,7 @@ buildNeighborList (CheckPair const& check_pair, bool /*sort*/) } #endif - auto& plev = this->GetParticles(lev); + auto& plev = this->GetParticles(lev); const auto& geom = this->Geom(lev); #ifdef AMREX_USE_OMP @@ -774,6 +840,115 @@ buildNeighborList (CheckPair const& check_pair, bool /*sort*/) } } +template +template +void +NeighborParticleContainer:: +buildNeighborList (CheckPair const& check_pair, amrex::Real bin_size, bool /*sort*/) +{ + AMREX_ASSERT(numParticlesOutOfRange(*this, m_num_neighbor_cells) == 0); + + BL_PROFILE("NeighborParticleContainer::buildNeighborList(bin_size)"); + + resizeContainers(this->numLevels()); + + for (int lev = 0; lev < this->numLevels(); ++lev) + { + m_neighbor_list[lev].clear(); + + for (MyParIter pti(*this, lev); pti.isValid(); ++pti) { + PairIndex index(pti.index(), pti.LocalTileIndex()); + m_neighbor_list[lev][index]; + } + +#ifndef AMREX_USE_GPU + neighbor_list[lev].clear(); + for (MyParIter pti(*this, lev); pti.isValid(); ++pti) { + PairIndex index(pti.index(), pti.LocalTileIndex()); + neighbor_list[lev][index]; + } +#endif + + auto& plev = this->GetParticles(lev); + const auto& geom = this->Geom(lev); + const auto* dx_arr = geom.CellSize(); + const auto* plo_arr = geom.ProbLo(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MyParIter pti(*this, lev); pti.isValid(); ++pti) + { + int gid = pti.index(); + int tid = pti.LocalTileIndex(); + auto index = std::make_pair(gid, tid); + + auto& ptile = plev[index]; + + if (ptile.numParticles() == 0) { continue; } + + Box tile_box = pti.tilebox(); + int ng = computeRefFac(0, lev).max() * m_num_neighbor_cells; + + // Build fine bins with bin_size instead of dx + int nbins[AMREX_SPACEDIM]; + amrex::Real phys_lo[AMREX_SPACEDIM], phys_hi[AMREX_SPACEDIM]; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + phys_lo[d] = plo_arr[d] + (tile_box.smallEnd(d) - ng) * dx_arr[d]; + phys_hi[d] = plo_arr[d] + (tile_box.bigEnd(d) + 1 + ng) * dx_arr[d]; + nbins[d] = std::max(1, static_cast( + std::ceil(static_cast(phys_hi[d] - phys_lo[d]) / bin_size))); + } + + GpuArray fine_dxi; + GpuArray fine_plo; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + fine_dxi[d] = static_cast(nbins[d]) / (phys_hi[d] - phys_lo[d]); + fine_plo[d] = phys_lo[d]; + } + + Dim3 fine_lo{0, 0, 0}; + Dim3 fine_hi{nbins[0]-1, AMREX_D_PICK(0, nbins[1]-1, nbins[1]-1), + AMREX_D_PICK(0, 0, nbins[2]-1)}; + + Gpu::DeviceVector off_bins_v; + Gpu::DeviceVector lo_v; + Gpu::DeviceVector hi_v; + Gpu::DeviceVector> dxi_v; + Gpu::DeviceVector> plo_v; + + off_bins_v.push_back(0); + off_bins_v.push_back(AMREX_D_TERM(nbins[0], * nbins[1], * nbins[2])); + lo_v.push_back(fine_lo); + hi_v.push_back(fine_hi); + dxi_v.push_back(fine_dxi); + plo_v.push_back(fine_plo); + + // num_cells=1: search 1 bin in each direction + // since bin_size ~ interaction radius, 1-cell search covers range + m_neighbor_list[lev][index].build(ptile, + check_pair, + off_bins_v, dxi_v, plo_v, lo_v, hi_v, 1); + +#ifndef AMREX_USE_GPU + const auto& counts = m_neighbor_list[lev][index].GetCounts(); + const auto& list = m_neighbor_list[lev][index].GetList(); + + int li = 0; + for (int i = 0; i < ptile.numParticles(); ++i) + { + auto cnt = counts[i]; + neighbor_list[lev][index].push_back(cnt); + for (size_t j = 0; j < cnt; ++j) + { + neighbor_list[lev][index].push_back(list[li++]+1); + } + } +#endif + } + } +} + template template void