Skip to content

Add physical search radius to NeighborParticleContainer #5262

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
desjardi wants to merge 1 commit into
base: development
Choose a base branch
from
+195 −4
Open

Add physical search radius to NeighborParticleContainer #5262

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions Src/Particle/AMReX_NeighborParticles.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -256,6 +263,14 @@ public:
template <class CheckPair>
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 <class CheckPair>
void buildNeighborList (CheckPair const& check_pair, amrex::Real bin_size, bool sort=false);

///
/// Build a Neighbor List for each tile
///
Expand Down Expand Up @@ -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<NeighborCommTag> local_neighbors;
Vector<std::unique_ptr<iMultiFab> > mask_ptr;

Expand Down
183 changes: 179 additions & 4 deletions Src/Particle/AMReX_NeighborParticlesI.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,17 @@ NeighborParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>
{
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<int>(
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<IntVect>& pshifts = periodicity.shiftIntVect();
const BoxArray& ba = this->ParticleBoxArray(lev);
Expand Down Expand Up @@ -550,6 +560,26 @@ getNeighborTags (Vector<NeighborCopyTag>& 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<int>& mask = (*mask_ptr[src_tag.level])[src_tag.grid];
AMREX_ASSERT(this->finestLevel() == 0);
Expand Down Expand Up @@ -604,14 +634,25 @@ getNeighborTags (Vector<NeighborCopyTag>& 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<IntVect>& 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<int>(
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)
Expand All @@ -638,6 +679,31 @@ template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
void
NeighborParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>
::fillNeighbors () {
m_search_radius = amrex::Real(0);
#ifdef AMREX_USE_GPU
fillNeighborsGPU();
#else
fillNeighborsCPU();
#endif
m_has_neighbors = true;
}

template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
void
NeighborParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>
::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<int>(std::ceil(search_radius / dx_min)));
}
#ifdef AMREX_USE_GPU
fillNeighborsGPU();
#else
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -774,6 +840,115 @@ buildNeighborList (CheckPair const& check_pair, bool /*sort*/)
}
}

template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
template <class CheckPair>
void
NeighborParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>::
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<int>(
std::ceil(static_cast<amrex::Real>(phys_hi[d] - phys_lo[d]) / bin_size)));
}

GpuArray<Real, AMREX_SPACEDIM> fine_dxi;
GpuArray<Real, AMREX_SPACEDIM> fine_plo;
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
fine_dxi[d] = static_cast<Real>(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<int> off_bins_v;
Gpu::DeviceVector<Dim3> lo_v;
Gpu::DeviceVector<Dim3> hi_v;
Gpu::DeviceVector<GpuArray<Real,AMREX_SPACEDIM>> dxi_v;
Gpu::DeviceVector<GpuArray<Real,AMREX_SPACEDIM>> 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 <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
template <class CheckPair, class OtherPCType>
void
Expand Down
Loading