Question about getNeighborTags search box scaling in multi-level NeighborParticleContainer

[desjardi]

I'm using NeighborParticleContainer with fillNeighbors() for DEM collisions in a multi-level setup (max_level=4, ref_ratio=2). I'm seeing severe performance issues and crashes in cacheNeighborInfo, and I believe it traces to the search box construction in getNeighborTags.

In the multi-level path (use_mask = false), in Src/Particle/AMReX_NeighborParticlesI.H around line 607:

IntVect ref_fac = computeRefFac(0, lev);

This scales the search box relative to level 0, so calling fillNeighbors() with nGrow=2 on a 5-level hierarchy produces a ±32-cell search box at the finest level (since computeRefFac(0, 4) = 16).

I think computeRefFac(src_tag.level, lev) might be the intended behavior? The physical reasoning: m_num_neighbor_cells defines a search radius in terms of cells on the particle's home level. When searching for neighbors on a different level, the number of cells should scale by the ratio between those two levels, not between level 0 and the destination. A particle living on level 4 and searching level 4 should need nGrow cells, not 16 * nGrow.

This would also be consistent with the single-level use_mask = true path, which uses nGrow directly with no scaling.

The same pattern appears in GetCommTagsBox around line 283 of the same file.

Is this the intended behavior, or is the 0 a bug? Am I misunderstanding how m_num_neighbor_cells is meant to be interpreted across levels?

[desjardi]

I realize that a key part of the question is what m_num_neighbor_cells is meant to represent. It's a single global integer, but it gets scaled by computeRefFac(0, lev) in the multi-level path. So is it meant to specify a number of cells on each level (in which case the scaling shouldn't be there), or a physical distance in level-0 cells (in which case it works as written, but arguably does not make sense with deep hierarchies)? A more standard DEM usage may call for a distance at this finest level instead?

[desjardi]

Quick follow-up: I tested replacing computeRefFac(0, lev) with just the unscaled value in both locations:

In getNeighborTags (line 614):

// was: Box pbox = amrex::grow(Box(iv, iv), ref_fac*nGrow) + pshift;
Box pbox = amrex::grow(Box(iv, iv), nGrow) + pshift;

In GetCommTagsBox (line 283):

// was: box.grow(computeRefFac(0, lev)*m_num_neighbor_cells);
box.grow(m_num_neighbor_cells);

This treats m_num_neighbor_cells as a constant number of cells on each level. With these changes, my multi-level DEM case (max_level=4, ref_ratio=2) runs correctly and at expected speed.

[desjardi]

Following up again on the broader issue behind my original question (sorry for spamming! I'm getting more familiar with the code so my thinking is maturing). In my application of interest here (DEM collisions with AMR), I'd like to be able to use coarser Eulerian meshes without being overly penalized on the collision/neighbor side. Right now, when dx >> rcrit, two things scale badly:

1. buildNeighborList bins particles at the mesh cell resolution. When dx is much larger than the interaction radius, the bins are too coarse and the pair search becomes unnecessarily expensive. This one seems straightforward to fix — allow passing a custom bin size (or rcrit) to NeighborList::build() so the binning is independent of the mesh.

2. fillNeighbors communicates all particles within no * dx of tile boundaries, which is a much larger shell than the physical interaction radius requires. A position-based fill — fillNeighbors(Real rcrit) — would only communicate particles actually within rcrit of boundaries, which ties directly into my earlier question about m_num_neighbor_cells semantics.

I have working implementations of both on my side and would be happy to contribute PRs. Would there be interest in this, and is there a preferred API direction?

[atmyers]

Hi @desjardi,

Thanks for your question and your interest in AMReX!

I think the original intention behind IntVect ref_fac = computeRefFac(0, lev); was that, if the particles have fixed size and the cell spacing on the base level has been chosen so that you have to search one cell to find all possible collision partners, you'd have to search ref_ratio cells on level 1, ref_ratio^2 on level 2, etc. This is kind of a hacky way of achieving a fixed search distance across multiple levels. However, as you point out, this doesn't allow enough flexibility, because you may well have set things up so that one cell is enough on the finest level, or you may want to use a different bin radius than dx altogether.

The improvements you mention sound good and I would be happy to review and accept PRs along the lines of your suggestions. Maybe for the API, you could add another constructor that takes the physical search distance instead of num_neighbor_cells, which would instead be computed?

[desjardi]

Thanks Andrew! That confirms my understanding of the intent.

Rather than a new constructor, I went with method overloads that take the physical search distance at call time:

1. fillNeighbors(Real search_radius) — stores search_radius in a new m_search_radius member, computes m_num_neighbor_cells from it for the mask path, and threads it through getNeighborTags / GetCommTagsBox. In getNeighborTags, when m_search_radius > 0, the search box uses ceil(search_radius / dx[lev]) per level instead of ref_fac * nGrow. It also does an early-out for particles that are farther than search_radius from all tile faces, which cuts communication volume significantly when search_radius << dx.

2. buildNeighborList(check_pair, Real bin_size) — builds fine bins of size bin_size (typically ≈ interaction radius) instead of dx, then searches ±1 bin. This decouples the pair search resolution from the mesh.

The overload approach seemed more flexible than a constructor since the search radius can change between calls (e.g., different radii for different physics, or adapting to particle size distributions). The original fillNeighbors() and buildNeighborList(check_pair) remain unchanged.

I have a working CPU implementation — happy to clean it up and submit a PR. Does the overload approach sound reasonable, or would you prefer the constructor route?

[atmyers]

Hi @desjardi,

Yes, your approach sounds reasonable, I'd be happy to take a look at your PRs.

Thanks,
Andrew

[desjardi]

PR submitted: #5262

The position-based filtering in getNeighborTags is in the CPU path; GPU correctness is maintained via m_num_neighbor_cells. Happy to iterate on any feedback, thanks for your time!