bagnold

Granular matter on the GPU, loadable through kernels: MLS-MPM with Drucker-Prager sand plasticity (Klar et al., SIGGRAPH 2016). The reference standards are the constitutive model's own return map, reproduced to machine precision, and the granular-flow laws of the experimental literature. Named for Ralph Bagnold (The Physics of Blown Sand and Desert Dunes, 1941).

Simulating sand by discrete elements stores per-contact friction history, order 200 bytes per grain, which caps particle counts long before a scene is interesting. MPM takes friction from the constitutive model's yield surface instead: state is per-particle, the grid is scratch, and nothing is stored per contact pair, so a hundred million grains fit on one card and run at over a billion particle-updates per second, with piles, silos, and collapses that track the published experimental laws.

A 390k-particle column released and slumping to its repose angle, simulated with this kernel (media/make_collapse.py), each particle path-traced as a sand grain in Cycles (media/render_blender.py).

Usage

from kernels import get_kernel

bag = get_kernel("phanerozoic/bagnold", version=1, trust_remote_code=True)

sim = bag.Sim(grid=(512, 512, 256), dx=1.25e-3,
              material=bag.sand(phi_deg=35.0))
sim.add(positions, ppc=8)          # particles per cell fixes their volume
sim.run(1000)                      # timestep from the elastic CFL
rho = sim.density()                # kg/m3, per particle

col = bag.Collider.from_sdf(phi, normals, dx, mu=0.3)   # phi < 0 inside solid
sim = bag.Sim(grid=g, dx=dx, material=bag.sand(), collider=col)

version selects the release branch; trust_remote_code is required by kernels for publishers without the trusted-publisher mark. Sim owns the grid scratch, picks a stable step, and keeps particles sorted for the shared-memory scatter; mpm_step exposes every array and parameter directly.

API

• mu=None is a stick boundary, under which a tall column hangs on its walls rather than loading the outlet; mu ~ 0.3 is a finite wall.
• E sets the elastic wave speed and therefore the timestep; plastic response is governed by the friction angle, so E can be reduced for a larger step while elastic strains stay small.
• sand(hardening=True) makes the friction angle track accumulated plastic strain per Klar et al. sec. 4; off by default.

Method

MLS-MPM transfers particle momentum to a background grid, applies forces and boundary conditions there, and transfers back, with a Drucker-Prager return map in principal log-strain space deciding how much deformation is elastic and how much is permanent. The p2g scatter is 85% of a step and its atomics 93% of that, so particles are sorted by grid block every 32 steps: one CUDA block owns a 4^3 node patch, accumulates into a 6^3 shared tile, and flushes once (Gao et al. 2018), cutting p2g from 22.5 to 7.1 ms/step. Drifted particles fall back to global atomics, so the sort interval is a performance knob, not a correctness one. g2p stays on the global path, where L2 serves 97% of its 27 reads per particle.

Measured

Throughput on one RTX 6000 Ada (48 GB), sand cube, dense grid:

The per-particle rate is flat, so cost scales with particle count and VRAM is the limit; past the point where the grid fills, sort and grid work start to dominate (221M particles run at 164 Mparticle/s) and indices are 64-bit throughout.

Head-to-head, same scene, same material, same GPU in each row:

taichi_elements loses most of its step rebuilding particle-in-block lists every substep (46 of 73 ms); bagnold amortizes its sort over 32 steps and wins with per-step sorting as well (~17 ms). Against Warp both pay the same dense cubic grid, and the margin is the sorted shared-memory scatter against global atomics.

Validation

The constitutive model is exact: bench/constitutive_probe.py replicates the kernel's return map and drives a single material element to critical state in triaxial compression:

The model returns the fed angle to machine precision; the fixed-corotated stress used in p2g against the Hencky-derived return map is bounded at 0.1 degree.

Angle of repose converges toward the fed angle under refinement:

Monotone across the 16x refinement; fitting theta(dx) brackets the limit to ~33-37 degrees, consistent with the fed 35 and not claimed as converged. Four surface estimators differ by 0.3-0.6 degree and share the trend, and runout is unchanged for floor mu >= tan(35 deg), above which the sand shears internally instead of sliding.

Collapse runout reproduces the two-regime empirical law (Lube 2004; Lajeunesse 2004), breaking at aspect ratio a ~ 1.7 with the correct morphology sequence:

Magnitudes run 60-90% of experiment (fitted 0.76a and 0.84 a^0.79 against 1.24a and 1.6 a^0.5), the known signature of rate-independent Drucker-Prager: no rate weakening, so the thin flowing front refreezes early. Resolution and basal friction are ruled out (bench/runout_conv.py). Treat absolute runout as conservative by tens of percent; the scaling structure, transition point, and morphology are quantitative.

Silo discharge is head-independent and Beverloo-scaled: on a 4.2M-particle silo the rate is flat to a few percent while the head falls:

The rates fit Beverloo's form at rms 2.5% with C = 0.514 against a literature 0.55-0.65, carrying a 4.44 mm aperture deficit that examples/orifice_profile.py localizes to a dead annulus of ~2.5 grid cells at the rim, shrinking under refinement while the core exit speed holds. Measured separately, exit velocity, effective aperture, and discharge density multiply to the observed W ~ D^3.31, and integrating the measured local density reproduces the silo's 103.1 g/s to 0.1.

Density is recoverable from particle state: bagnold.density(pF, pJp) = rho0 / (det(F) Jp) reads a dropped block at a mean 2155 kg/m3, with grid mass and summed particle volumes agreeing to 5.4% as independent checks.

Demonstrations

Each demo is one file in demos/, runs in seconds to minutes on one card, and checks itself. calibrate_friction.py recovers a hidden friction angle from a pile to 0.04 degree in 34 simulations by golden-section search. design_hopper.py sizes a hopper orifice for a target discharge (21.1 mm for 134 g/s in 8 simulations). learn_surrogate.py trains an MLP on simulated (friction angle -> pile height) pairs and inverts it (32.4 degrees recovered).

Requirements and limits

• float32, CUDA, compute capability 8.0+, Linux x86_64; on Windows load_local.py JIT-builds the same source.
• Dense background grid: memory is nx*ny*nz*4 floats regardless of how much holds material.
• Explicit integration, dt <= ~0.3 dx / sqrt(E/rho): a second of hourglass is tens of thousands of steps.
• Discharged material must be removed or it piles beneath the outlet and blocks it (examples/hourglass.py parks it outside the grid).
• Absolute collapse runout is conservative by tens of percent; scaling structure and morphology are quantitative.

References

Klar et al., "Drucker-Prager Elastoplasticity for Sand Animation" (SIGGRAPH 2016); Hu et al., MLS-MPM (2018); Gao et al., GPU MPM optimization (2018); Lube et al. (2004) and Lajeunesse et al. (2004) on column collapse; Beverloo et al. (1961) on silo discharge.

License

Apache-2.0.