Cutting Distance Checks 574x Lost 104 Contacts — One Grain Sets the DEM Sector Size

A 574x Faster Search Lost 104 Contacts#

When aluminum particle agglomeration is solved with DEM, neighbor search is the first thing that slows down. I packed 2400 grains into one cluster and looked for contact pairs. Brute force takes 2,878,800 distance checks and finds 1539 contacts. Cutting sectors to a quarter of the grain diameter, then scanning only 27 cells, dropped that to 5015 checks. That is 574x fewer. The search found 1435 contacts.

104 pairs disappeared quietly. No error, no warning. One agglomerate simply split into two. Whether those pairs sat outside the stencil, and where the sector size has to stop, shows up in the cell count.

Drop s / d_max from 1.00 to 0.25 and the blue box shrinks with the distance checks, while red missed contacts appear.

Raising shell finds a setting where the red lines vanish. That setting is the later inequality $k\,s \ge d_{\max}$. If the slider moves and the red lines stay, that setting is not counting contacts.

The Cost of DEM Lives in a Double Loop, Not in Communication#

DPM (Discrete Parcel Method) couples particles to the flow field in both directions and nothing else. Grains never see each other. Cost scales with grain count.

DEM (Discrete Element Method) adds collisions, repulsion and attraction between grains. The structure is nearly identical to molecular dynamics. Because grain $i$ must know who it touches, a pairwise loop appears.

$N$ is the grain count. 2400 grains give 2.88 million pairs; 9310 grains give 43.3 million. The wall is this quadratic term, not MPI traffic. Even with good parallelization, each rank runs the same double loop locally. Adding cores leaves that loop in place. Cutting a global distance down to nearby cells is the same family of move as Fast Marching for the Eikonal wall distance.

What a Sector Grid Removes Is Distance Evaluations#

Cut the domain into cubic sectors of edge $s$. Register each grain in its own sector list and restrict the search to that sector and its neighbors. Scanning $2k+1$ cells along each axis sweeps the volume

where $k$ is the shell radius, the number of layers you look through. The expected number of grains inside is $N V_k / L^3$, with $L$ the domain edge. Distance evaluations drop to a term linear in $N$.

There Is Exactly One Safety Condition#

Grains $i$ and $j$ touch when their center distance falls below the sum of radii.

Here $r_{\max}$ is the largest radius in the domain and $d_{\max}$ the matching diameter. A $2k+1$ stencil, meanwhile, guarantees a search reach of $k\,s$. A grain can sit flush against its own sector boundary, so the reach you can rely on is $k\,s$, no more. The condition for missing nothing is therefore a single line.

With $s = d_{\max}$ you need $k=1$, so 27 cells suffice. Cutting down to $s = d_{\max}/4$ demands $k=4$, which is 729 cells. The experiment above kept 27 cells while shrinking $s$ by four, so the condition broke. Those 104 pairs were touching outside the stencil.

Six Faces, Twelve Edges, Eight Corners#

A $k=1$ stencil is one home sector plus 26 neighbors. Those 26 split into 6 that share a face, 12 that share an edge, and 8 that share only a corner. The slide deck I reviewed listed 12 faces and 6 edges. The two numbers are swapped.

The distinction matters because the 8 corner neighbors are the easiest to drop. Faces and edges alone give 19 cells. Running the same cluster with 19 cells missed 1 pair out of 1539. One pair looks harmless, but if it was the link holding an agglomerate together, one cluster gets reported as two. A corner neighbor only fires when both grains sit near the shared corner, so test cases almost never hit it. It hits halfway through a production run instead.

One Oversized Grain Sets the Grid for the Whole Domain#

What enters the inequality is $d_{\max}$, not the mean diameter. For monodisperse grains that makes no difference. Agglomeration problems are rarely monodisperse.

Idea 1 in the reviewed deck states as assumption 1 that the center grain has the largest diameter in the domain. The center grain attracts everything else, and the rest only feel repulsion among themselves. In that setup a single grain fixes $d_{\max}$.

In the experiment the center diameter was 0.060 and the largest of the rest was 0.030. 2399 grains sit at or below 0.030, yet the sector edge has to be 0.060 to stay safe. Every sector then holds eight times more volume than it needs. Distance checks scale with it: 269,450 as measured.

Cut Finer and Empty-Cell Lookups Become the New Cost#

So why not shrink $s$ and raise $k$? Substituting $s = d_{\max}/k$ turns the swept volume into

which is $27\,d_{\max}^3$ at $k=1$ and converges to $8\,d_{\max}^3$ as $k\to\infty$. However fine you cut, distance checks only ever fall by a factor of 27/8, that is 3.4x. The number of sectors to visit, meanwhile, grows as $(2k+1)^3 \approx 8k^3$. At $k=4$ that is 729 lookups per grain, so 1.75 million hash lookups for 2400 grains. Most of them hit empty sectors.

The measurement says exactly that. $k=4$ with $s=d_{\max}/4$ cut distance checks to 125,548, but sector lookups rose to 1,749,600 and the wall clock came out at 0.291 s. That is 3.5x slower than 0.082 s for $k=1$ with $s=d_{\max}$. Distance checks improved 2.1x while the total got worse.

Five Settings on the Same Cloud of 2400 Grains#

Assumptions and measurements live in one file. Standard library only.

import math, random
 
def pack_particles(n, seed=7):
    """One largest grain at the center; the rest injected as a Gaussian cloud"""
    rng = random.Random(seed)
    p = []
    for _ in range(n):
        d = 0.012 + 0.018 * rng.random() ** 2          # diameter 0.012 ~ 0.030
        c = [min(0.97, max(0.03, 0.5 + rng.gauss(0, 0.105))) for _ in range(3)]
        p.append((c[0], c[1], c[2], 0.5 * d))
    p[0] = (0.5, 0.5, 0.5, 0.030)                      # center: diameter 0.060
    return p
 
def brute_contacts(p):
    pairs, ev = set(), 0
    for i in range(len(p)):
        xi, yi, zi, ri = p[i]
        for j in range(i + 1, len(p)):
            xj, yj, zj, rj = p[j]
            ev += 1
            if (xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2 < (ri + rj)**2:
                pairs.add((i, j))
    return pairs, ev
 
def sector_bins(p, s, start=0):
    b = {}
    for i in range(start, len(p)):
        x, y, z, _ = p[i]
        b.setdefault((int(x/s), int(y/s), int(z/s)), []).append(i)
    return b
 
def neighbor_offsets(shell):
    r = range(-shell, shell + 1)
    return [(a, b, c) for a in r for b in r for c in r]
 
def binned_contacts(p, s, shell, start=0):
    bins, off = sector_bins(p, s, start), neighbor_offsets(shell)
    pairs, ev, look = set(), 0, 0
    for i in range(start, len(p)):
        xi, yi, zi, ri = p[i]
        cx, cy, cz = int(xi/s), int(yi/s), int(zi/s)
        for a, b, c in off:
            look += 1
            for j in bins.get((cx+a, cy+b, cz+c), ()):
                if j <= i:
                    continue
                xj, yj, zj, rj = p[j]
                ev += 1
                if (xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2 < (ri + rj)**2:
                    pairs.add((i, j))
    return pairs, ev, look
 
def split_big_contacts(p, s, shell):
    """Grain 0 (the largest) is handled by brute force; sectors sized by the rest"""
    pairs, ev, look = binned_contacts(p, s, shell, start=1)
    xb, yb, zb, rb = p[0]
    for j in range(1, len(p)):
        xj, yj, zj, rj = p[j]
        ev += 1
        if (xb-xj)**2 + (yb-yj)**2 + (zb-zj)**2 < (rb + rj)**2:
            pairs.add((0, j))
    return pairs, ev, look
 
N = 2400
P = pack_particles(N)
DMAX = 2 * max(q[3] for q in P)                 # largest diameter overall
DSUB = 2 * max(q[3] for q in P[1:])             # largest diameter excluding the center grain
truth, ev0 = brute_contacts(P)
print("N=%d  L=1.0  d_max=%.3f  d_max(sub)=%.3f" % (N, DMAX, DSUB))
print("all-pairs : contacts=%d  dist=%d" % (len(truth), ev0))
print()
print("%-20s %-7s %-6s %-9s %-9s %-5s %s" %
      ("case", "s/dmax", "cells", "dist", "lookups", "miss", "safe?"))
for f, shell in ((1.0, 1), (0.5, 1), (0.25, 1), (0.25, 4)):
    s = DMAX * f
    got, ev, look = binned_contacts(P, s, shell)
    print("%-20s %-7.2f %-6d %-9d %-9d %-5d %s" %
          ("sector shell=%d" % shell, f, (2*shell+1)**3, ev, look,
           len(truth - got), "yes" if shell*s >= DMAX - 1e-12 else "NO"))
got, ev, look = split_big_contacts(P, DSUB, 1)
print("%-20s %-7.2f %-6d %-9d %-9d %-5d %s" %
      ("big-split shell=1", DSUB/DMAX, 27, ev, look, len(truth - got), "yes"))
print()
print("required shell = ceil(d_max/s) :", [(f, math.ceil(1/f)) for f in (1.0, 0.5, 0.25)])

The output:

N=2400  L=1.0  d_max=0.060  d_max(sub)=0.030
all-pairs : contacts=1539  dist=2878800
 
case                 s/dmax  cells  dist      lookups   miss  safe?
sector shell=1       1.00    27     269450    64800     0     yes
sector shell=1       0.50    27     38914     64800     0     NO
sector shell=1       0.25    27     5015      64800     104   NO
sector shell=4       0.25    729    125548    1749600   0     yes
big-split shell=1    0.50    27     41229     64773     0     yes
 
required shell = ceil(d_max/s) : [(1.0, 1), (0.5, 2), (0.25, 4)]

Push s/d left and the blue curve (distance checks) falls, then flattens, while the orange curve (total) turns back up. Raising lookup / check cost slides the minimum to the right. Turn on pin grain at cell corner and the grain sits on a sector corner, which is why the red cells exist.

The dist column reads 269450, 38914, 5015, 125548, 41229. 5015 is 574 times smaller than the brute-force 2878800, and the contact count drops from 1539 to 1435, a loss of 104. Same 27 cells and 64800 lookups: s/dmax=1.00 is miss 0 and safe yes, 0.50 is miss 0 and safe NO, 0.25 is miss 104 and safe NO. shell=4 opens 729 cells, dist 125548, lookups 1749600, miss 0. Wall time goes from 0.082 s to 0.291 s, a factor of 3.5. big-split is dist 41229, lookups 64773, miss 0. 41229 against 269450 is the 6.5x in the text, and 0.023 s is 1/29 of the brute-force 0.668 s.

Nineteen cells, faces and edges only, drop 1 pair out of 1539. Smaller than 104, and still enough to split an agglomerate if that pair was the link. Diameters 0.060 and 0.030 are a factor of two in length, so sector volume jumps by 8. 2399 grains sit at or below 0.030, yet the cells follow 0.060. In $V_k$, $k=1$ is $27\,d_{\max}^3$ and a large $k$ only reaches $8\,d_{\max}^3$, so distance checks fall by 27/8, that is 3.4x. Buying that 3.4x with 729 cells came back as 0.291 s.

Look at the third row. s/dmax=0.50 with shell=1 violates the safety condition yet misses nothing. That is luck. Only one grain has diameter 0.060, and none of its contact partners happened to sit in a corner direction. Reseed the cloud or advance one step and the number changes. "I ran it and it was fine" is not evidence for a condition of this kind. A miss count of 0 with the inequality still NO is not a pass.

Splitting Off the Big Grain Buys Safety and Speed Together#

The last row is the fix. Handle the one center grain by brute force separately, and size sectors by the largest of the remaining 2399 grains, 0.030. The extra cost is 2399 distance checks, one $O(N)$ pass.

Distance checks fall to 41,229, which is 6.5x fewer than the 269,450 of the $s=d_{\max}$ setting. Sector lookups stay flat at 64,773. Nothing is missed, and this time the inequality guarantees it rather than luck. Among the remaining grains, $r_i + r_j \le 0.030 = 1 \cdot s$ holds by construction. Wall clock is 0.023 s, one twenty-ninth of the 0.668 s brute force.

The gap widens with scale. An agglomerate is a few large grains with thousands of small ones attached. Treating the handful in the tail of the size distribution separately needs less code than building a two-level grid. The point is to break the pattern where a few outliers set the cell size for everyone else.

What Those 104 Missing Contacts Change#

Miss one contact and the repulsion for that pair becomes zero. The two grains pass through each other, or stay overlapped and drop out of the surface-area update. In aluminum particle combustion the overlapped area of an agglomerate maps directly onto oxide-layer area. Burn rate comes out too high. Mass conservation never breaks, so no residual flags it. Same quiet failure as WENO smearing a small vortex while the average still looks fine.

So when you write a sector size into search code, three things are worth checking. Whether $k\,s \ge d_{\max}$ holds, whether that $d_{\max}$ is remeasured every step (merging and breakup change it), and how many grains are single-handedly pushing $d_{\max}$ up. If the answer to the third is a single digit, those grains are cheaper counted separately than binned.