IBM particle-laden: force exchange ships ~96% zeros; a ghost-particle exchange would be ~3000x smaller

[sbryngelson]

Static read of master while inspecting the multi-rank IBM for a large particle-laden weak-scaling study. The immersed-boundary force/torque exchange moves a great deal more data than the problem requires, and a standard ghost-particle exchange looks like it would replace it at a tiny fraction of the cost. Writing up the arithmetic so the assumptions can be checked. This is static analysis, not profiling — please push back where the model is wrong.

Reference configuration throughout: 7500 spherical IBs per rank, particle diameter 10 cells, ib_neighborhood_radius = 1 (the default), $300^3$ cells per rank, 3D, TVD-RK3, 8 ranks/node. wp is compile-time selectable in m_precision_select.f90; numbers below assume double, so halve them for a single-precision build.

1. Every pass sends the whole neighborhood array

s_communicate_ib_forces in src/simulation/m_ibm.fpp packs all tracked IBs on each pass:

do i = 1, num_ibs
    send_ids(i) = patch_ib(i)%gbl_patch_id
    send_ft(1:3,i) = forces(i,:)
    send_ft(4:6,i) = torques(i,:)
end do

num_ibs is the neighborhood count, not the owned count — m_start_up.fpp filters with f_neighborhood_ranks_own_location, so at R=1 it spans a 3x3x3 block of subdomains: 27 x 7500 = 202,500 entries, giving a 10.5 MB message. There are 4R passes per dimension (2R accumulate + 2R back-propagate), 12 per RK substep, and s_propagate_immersed_boundaries(s) runs per substep:

one message	10.5 MB
per rank per RK substep	126 MB
per rank per timestep (RK3)	379 MB
per node per timestep	3.0 GB
entries with a nonzero local contribution	~7,500 of 202,500 (3.7%)

A rank only has contributions for IBs overlapping its own subdomain. About 96% of the payload is zeros, re-sent 36 times per timestep. This is comparable to or larger than the fluid halo exchange for the same configuration, though that comparison depends on sys_size and buff_size so I am not putting a precise ratio on it.

2. A ghost-particle exchange would be dramatically smaller

The force on a particle only needs summing over the ranks whose subdomains it actually intersects — at most 2 per direction for a 10-cell particle in a 300-cell subdomain. It does not need to reach every rank in the neighborhood.

With centroids uniform, a particle is wholly interior with probability $(1 - 10/300)^3 = 0.903$. So per rank:

• ~6,775 of 7,500 particles need no communication at all
• ~725 straddle a boundary; expected neighbour-sends ~775 (crossing k dimensions implies $2^k - 1$ neighbours)

approach	volume per substep	rounds
current sweep	126 MB	12
straddlers to all 26 neighbours (conservative)	0.98 MB	1
straddlers to only the ranks they intersect	40 KB	1

That is ~130x for the lazy version and ~3100x for the targeted one, in a single MPI_Neighbor_alltoallv instead of 12 sequential rounds. Per timestep: 379 MB becomes ~121 KB.

This is the ghost/halo particle pattern standard in parallel DEM and MD. The pieces already exist in the tree — f_local_rank_owns_location is used exactly this way in m_collisions.fpp to make one rank responsible for each collision pair.

This also reframes the O(R) argument. The sweep is O(R) in message count versus O(R^d) for reaching every neighbourhood rank directly, which is its motivation. But reaching every neighbourhood rank is not required; only the ranks a given particle intersects are. Once the exchange is sparse, the quantity being minimised is no longer the binding cost.

3. Host round-trip on every pass

forces and torques are local automatic arrays (m_ibm.fpp:981), not in a persistent device region, and the unpack loop carries

$:GPU_PARALLEL_LOOP(..., copy='[forces, torques, recv_forces_snap, recv_torques_snap]')

so four arrays of num_ibs * 3 reals move host<->device on every pass: roughly 39 MB x 12 passes x 3 substeps ~ 1.4 GB per rank per timestep. Cray MPICH supports GPU-aware MPI (MPICH_GPU_SUPPORT_ENABLED=1); device pointers could go straight to MPI_SENDRECV. Several copy= clauses also look stronger than needed and want copyin/copyout.

4. Rank placement is not topology-aware

MPI_CART_CREATE in src/common/m_mpi_common.fpp passes reorder = .false., so the Cartesian grid maps row-major onto MPI_COMM_WORLD. Neighbour separation at 8 ranks/node:

ranks	z-neighbour	y-neighbour	x-neighbour
4096 (16^3)	same node	2 nodes	32 nodes
32768 (32^3)	same node	4 nodes	128 nodes
74088 (42^3)	same node	5.25 nodes	220 nodes

At scale the x-direction sweep crosses ~220 nodes of a dragonfly on each of its passes. This grows monotonically with rank count, which is the right shape for degradation that only shows up at large node counts.

5. Ownership handoff is latency-bound, not bandwidth-bound

Worth separating from the above, because it is often assumed to be the expensive part. s_handoff_ib_ownership broadcasts only newly-owned patches, so for particles that move slowly relative to a subdomain the payload is nearly empty. What it does cost is 26 neighbour rounds per substep, so ~78 latency exposures per timestep, which interacts badly with item 4.

It also does allocate(send_buf(buf_size), recv_bufs(buf_size, max_nbrs)) on every call, with buf_size sized for num_local_ibs_max and 26 receive buffers. Hoisting that allocation looks free.

6. Global reduction each step under cfl_adap_dt

m_time_steppers.fpp calls s_mpi_allreduce_min(dt_local, dt) in s_compute_dt, dispatched every step under cfl_adap_dt (only at t_step == 0 under cfl_const_dt). At tens of thousands of ranks that is a global min-reduction and a hard synchronisation point per step, which also exposes load imbalance.

Suggested order

1. Sparse ghost-particle exchange (item 2) — the large one, and it subsumes much of the rest.
2. GPU-aware MPI (item 3) — no algorithm change.
3. Topology-aware rank mapping (item 4) — reorder = .true. or MPICH_RANK_REORDER_METHOD.
4. Fuse the IB exchange into the existing fluid halo exchange — same neighbours, same cadence, removes a set of separate latency rounds.
5. Persistent neighbour collectives — the pattern is static between ownership changes, so MPI_Neighbor_alltoallv_init can amortise setup.
6. Hoist the handoff allocations (item 5), and relax the dt reduction (item 6) where CFL margin allows.

Each is independently measurable.

Related observations

• The back-propagation loop uses min(2*ib_neighborhood_radius, num_procs_x - 1), the same count as accumulation... each rank appears to end with a window offset by R. Withdrawn — this was wrong, see the correction below. The 2R count is correct and exactly minimal; I had missed that updates are gated per particle by s_get_neighborhood_idx, which prevents a particle's data from propagating outside its own neighbourhood. Nothing in items 1-6 depends on this.
• s_detect_ib_collisions (the ghost-point detector) is commented out in src/simulation/m_collisions.fpp in favour of s_detect_ib_collisions_n2. At neighbourhood-scale num_ibs the quadratic path is very expensive for any collisional run. (This one stands.)

Related: #1536 also concerns ib_neighborhood_radius. See also #1704 for a lubrication correction proposal.

---

Edited after filing: added the ghost-particle comparison in item 2, which changes the recommended order; separated the ownership handoff (item 5) as latency- rather than bandwidth-bound; added the precision caveat; removed a specific fluid-halo ratio I could not pin down; and struck the back-propagation observation, which was incorrect.

[sbryngelson]

cc @danieljvickers — this came out of looking at the multi-rank IBM, so you are much better placed than I am to say which parts of the model above are right.

Five things I could not settle from the source alone, and where a wrong guess would change the conclusions:

1. Which branch produced the large Frontier runs? Everything above is read off master. If the runs used a branch that differs in s_communicate_ib_forces or re-enables the ghost-point collision detector, several items shrink or disappear.
2. Was collision_model nonzero in the scaling case? It defaults to 0. If collisions were off, the quadratic-detector concern does not apply to those runs, though it still applies to collisional ones.
3. cfl_adap_dt or cfl_const_dt? That decides whether the global s_mpi_allreduce_min in item 6 fires every step or only at step 0.
4. Is num_ibs really ~27 x 7500 in that configuration? The whole volume estimate scales linearly with it. I inferred the 3x3x3 neighbourhood from f_neighborhood_ranks_own_location plus R=1, but if the packing loop is bounded by something tighter in practice the numbers come down proportionally.
5. The min(2*ib_neighborhood_radius, ...) count in the back-propagation loop. I traced it as ending one neighbourhood off-centre past 2R+1 ranks per direction, but that is a paper trace, not a run. A case on four or more ranks per direction with a particle straddling a boundary, comparing patch_ib%force on the two ranks that share it, would settle it either way in a few minutes.

Item 2 in the issue is the one I would most like your read on. If the sparse ghost-particle exchange is sound, it makes most of the rest moot, and it is close to what the collision path already does with f_local_rank_owns_location. If there is a reason the sweep has to reach every rank in the neighbourhood that I have missed, that would be very useful to know.

[sbryngelson]

Retracting the first of the two "related observations" — the back-propagation loop count. I was wrong, and @danieljvickers was right to push back.

I claimed the back-propagation phase runs 2R passes where R are needed, and that past 2R+1 ranks per direction each rank ends with its neighbourhood window shifted by R. That is incorrect. The loop count is right, and 2R is not merely sufficient but exactly minimal.

What I got wrong. I modelled the sweep as a scalar sliding-window sum over ranks. The code is not that, because every update is gated per particle:

call s_get_neighborhood_idx(recv_ids(i), j)
if (j > 0) then     ! only for particles this rank tracks

forces and torques are per-particle arrays, and a particle only exists in the arrays of ranks that track it. So its partial sums cannot propagate outside its own neighbourhood: ranks beyond it never hold an entry, never send one, and never overwrite anything. The window physically cannot slide off the particle, which is the failure mode my scalar model produced. The extra passes are idempotent — they re-deliver the same correct value.

And the count is minimal. After accumulation, the rank at the +R edge of a particle's neighbourhood holds the complete sum; carrying it to the -R edge is 2R hops. Re-running the trace with the tracking filter included, at R=1 with 42 ranks, R back-propagation passes leave 42 of 126 tracked entries stale, 2R is correct, and anything beyond 2R is correct but wasteful. Same pattern at R=2 and R=3. So the R I proposed would have introduced precisely the bug I thought I was reporting.

Consequently the "25% more communication than needed" in that bullet is also withdrawn, and the 8-in-2D / 12-in-3D counts in the associated write-up are correct.

What stands. Nothing else in the issue depends on this. Items 1 through 6 concern message volume, host round-trips, rank placement, the handoff, and the dt reduction, none of which touch the loop count. The second related observation, that s_detect_ib_collisions is commented out in favour of the quadratic path, is unaffected.

One presentational note that is worth something. The reason I misread it is that the pseudocode in the accompanying paper writes the force and torque as single values, with no per-particle indexing and no sign of the tracking filter. Read literally it is a scalar sliding window. Subscripting them per particle, or adding one line stating that updates apply only to particles in the receiving rank's neighbourhood, would prevent the same misreading — and would make the correctness argument visible, since once the filter is stated the 2R count is self-evidently right.

Apologies for the noise, and thanks for catching it quickly.