Implementation — GPU Optimization & Reference Codebases

This document is the practical engineering layer: how to make MPM fast on a GPU and the reference implementations worth studying or building on. The product questions — is MPM shipping in commercial games, how novel would an MPM-native game be, how to talk about that — live in their own doc, Games.

1. The central performance problem: the P2G scatter

Recall the four-phase step (Theory §4): reset → particle-to-grid (P2G) → grid update → grid-to-particle (G2P). Profiling almost always lands on P2G as the bottleneck, for one structural reason:

Many particles write to the same grid node at the same time. A node’s mass and momentum are sums over every particle in its B-spline stencil. On a GPU, thousands of particle-threads racing to accumulate into shared node memory means atomic write contention — the dominant cost.

G2P, by contrast, is a gather (each particle reads from its stencil and writes only to itself) — embarrassingly parallel, no contention. So optimization effort concentrates on P2G.

2. GPU optimization techniques

The authoritative reference is Gao et al., GPU Optimization of Material Point Methods, SIGGRAPH Asia 2018 (PDF, code: kuiwuchn/GPUMPM). The techniques below come from that paper and the surrounding GPU-MPM literature.

2.1 Tame the atomic contention

• Spatial sorting + tiling. Sort particles by the grid cell they occupy so that particles writing to the same nodes are processed together. This both reduces scattered random writes and enables shared-memory staging.
• Shared-memory (block-local) accumulation. Stage a tile of the grid in fast on-chip shared memory; threads in a block accumulate there (cheap atomics or warp reductions), then flush the tile to global memory once. Turns many expensive global atomics into a few.
• Warp-level reduction. Particles mapped to the same node within a warp combine their contributions via warp shuffle intrinsics before a single atomic, cutting contention by the warp width.
• Fixed-point atomics. GPU floating-point atomic adds are slower (and non-deterministic in ordering, which perturbs results). Encoding the accumulator as fixed-point integers allows fast integer atomicAdd and makes the scatter order-independent / deterministic — valuable for reproducibility. This is a known, widely-used trick — the canonical citation is from GPU molecular dynamics (Le Grand, Götz & Walker, SPFP, 2013; also OpenMM), and EA’s PB-MPM WebGPU code uses fixed-point atomicAdd for P2G. Two caveats: (a) determinism from fixed-point scatter is same-platform run-to-run reproducibility — cross-platform bitwise determinism is a much harder property (see the engine-determinism links in References Determinism & reproducibility); (b) the scatter is only one nondeterminism source — any other race in the pipeline (e.g. the order in which per-node slots or lists are claimed) must be neutralized separately, typically by sorting the racy structure by a stable key before any floating-point math consumes it.

The worked shape of “sort + tile + stage in shared memory.” The three bullets above are usually given abstractly. EA’s PB-MPM main branch is the most legible complete instance of all of them at once, so it is worth reading as the reference implementation of this subsection rather than as a PB-MPM curiosity (§3.3 covers the algorithm; this is only the memory strategy).

It bins particles into 6×6-cell tiles (“bukkits”) each substep via three small kernels (bukkitCount → bukkitAllocate → bukkitInsert) that build an index table — particles themselves never move or reorder, which is what keeps §2.1’s determinism and index-stability intact. The entire iteration body is then ONE kernel (g2p2g.wgsl), one workgroup per 64-particle slice of a bukkit:

1. Threads cooperatively load the tile’s grid nodes + a 1-node halo into workgroup shared memory, applying the grid update during the load — mass-normalize, collide, boundary — so the read is the grid pass.
2. Barrier. Each thread gathers its particle from shared memory, does its per-particle work, and scatters back into shared-memory atomics.
3. Barrier. The tile flushes to the global grid — one atomicAdd per node per channel per tile (halo overlap is why it must still be an add).

Global atomic traffic per particle drops from 27 nodes × channels to a per-tile flush — one-to-two orders of magnitude less contention, and global atomics are exactly what is slow on Apple GPUs. Three properties are worth extracting:

• Determinism survives the whole scheme. Shared tiles accumulate the same fixed-point integers; the tile partition is position-determined, not race-determined; integer adds commute within and across tiles. The bin insertion order is racy, but it only permutes commutative sums and per-particle writes. Tiling is therefore compatible with a bit-reproducibility invariant — rare among optimizations, and worth knowing before profiling ever asks for it.
• The shared-memory budget closes in 3D. A 4³-cell tile + 1-cell halo = 6³ = 216 nodes × 4 channels × 4 B × 2 tiles ≈ 7 KB, comfortable against a 32 KB workgroup limit, with 64 threads each staging ~4 nodes.
• The grid clear can ride along. Rotate three grid buffers per iteration — read, atomic-accumulate, and next-write, which the fused kernel zeroes with plain non-atomic stores while it is already touching that memory. No clear dispatch exists in the loop at all. The separate-dispatch form of the same economy is to have the grid-update kernel zero each node’s accumulator in place right after decoding it; its one precondition is that each accumulator has exactly one reader, which a pass that reads neighboring nodes’ accumulators during that same dispatch violates.

2.2 Data layout & memory coalescing

• AoSoA (Array-of-Structs-of-Arrays). Group particle attributes in small batches (e.g. 32 = one warp) of contiguous per-attribute arrays. This balances the cache-locality of AoS against the coalesced access of SoA, so a warp reads consecutive addresses in one transaction. Gao et al. use this layout.
• Sort particles for memory coherence. Cell-sorted particles access nearby grid nodes, maximizing cache hits and coalescing on both P2G and G2P.

2.3 Sparse grids

Most of the background grid is empty (material occupies a fraction of the domain). Sparse grid structures — notably SPGrid (Sparse Paged Grid, Setaluri et al.) — store only active blocks, slashing memory and letting high-resolution simulations fit on-GPU. Active-block bookkeeping is updated as particles move.

2.4 Reported results

Gao et al. report ~20–50× speedups over an optimized multicore CPU implementation, scaling to millions of particles — bringing offline-quality MPM toward interactive rates and making large-scale visual-effects (VFX) and engineering work practical.

Real-time vs. throughput. The Gao-style optimizations maximize throughput (particles/second). For real-time games, the harder wall is the timestep: explicit MPM’s CFL (Courant–Friedrichs–Lewy) stability limit forces many small steps per frame. Two complementary attacks: MLS-MPM (Moving Least Squares MPM — do each step ~2× cheaper) and PB-MPM (position-based MPM — take larger, stable steps). See §3.

3. Reference implementations

A curated, build-on-able set — ordered roughly from most readable to most production-optimized. Always open the actual LICENSE in each repo before relying on it; terms can change.

3.1 nialltl’s incremental_mpm — the best starting point

A readable MLS-MPM implementation in Unity using data-oriented High-Performance C#, with a detailed companion guide (nialltl.neocities.org/articles/mpm_guide). The headline demo is a real-time ~40,000-particle fluid sim. The guide builds up in three parts:

1. Part 1 — basic MLS-MPM movement using APIC (affine particle-in-cell) transfers, no deformation model. Establishes the four-phase loop, the 3×3 quadratic-B-spline stencil, and the affine momentum matrix C.

2. Part 2 — elastic solids via a Neo-Hookean constitutive model, introducing the deformation gradient F and its update. The coded stress (from the guide), with elastic_mu as μ and elastic_lambda as λ, is:

   $$\begin{aligned} J &= \det \mathbf{F} \\ \mathbf{F}^{-T} &= \left(\mathbf{F}^{T}\right)^{-1} \\ \mathbf{P} &= \mu\left(\mathbf{F} - \mathbf{F}^{-T}\right) + \lambda \log(J)\,\mathbf{F}^{-T} \\ \boldsymbol{\sigma} &= \frac{1}{J}\,\mathbf{P}\,\mathbf{F}^{T} \end{aligned}$$

3. Part 3 — real-time fluids (credited partly to Grant Kot), via three tweaks that buy a much larger timestep:

   • Recompute particle volume every frame from a density estimate in the style of Smoothed Particle Hydrodynamics (SPH) (scatter mass to grid, gather back) instead of integrating it — more memory reads, but far larger stable Δt.
   • A Newtonian-fluid constitutive model with a Tait equation of state for weakly-compressible fluids (as in weakly-compressible SPH, WCSPH).
   • Softened, predictive boundary conditions acting on particles directly, to curb high-speed tunneling through the domain walls.

   The credit is explicit in the guide — “credits to Grant Kot for the altered volume calculations and boundary conditions tricks” — and the first of those traces back to Kot’s own 2010 applet note: “Instead of integrating the density over time (which is what most of the MPM papers do), I do a density summation every frame. Because this is not dependent on previously calculated values of density, there is no accumulated error.” Origin and the rest of his work: References §3.

3.2 taichi_mpm — the 88-line MLS-MPM

The canonical compact MLS-MPM reference from Hu et al. (the MLS-MPM authors), written in the Taichi language for high performance with CPU/CUDA backends. The famous 88-line version is the clearest end-to-end MLS-MPM in existence (a more-readable variant by David Medina also circulates). MIT-licensed and explicitly commercial-OK. Cross- reference it against the SIGGRAPH 2016 course notes for the theory.

3.3 EA SEED’s PB-MPM — Position-Based MPM (real-time)

Presented at SIGGRAPH 2024 by Chris Lewin (EA SEED). PB-MPM is a semi-implicit compliant-constraint (Position-Based-Dynamics-style) reformulation of MPM that is stable at essentially any timestep while staying “as easy to implement as an explicit integrator.” This directly attacks MPM’s single biggest real-time barrier — the small-timestep stability requirement — and is explicitly framed for real-time applications where the hardest requirement is stability under autonomous, potentially violent player input.

Algorithmically, PB-MPM differs from explicit MPM by only a small number of changes: it wraps the ParticleToGrid → GridUpdate → GridToParticle cycle in an iteration loop with an added SolveConstraints step (the lines marked in red in the paper’s Algorithm 1).

Repository: github.com/electronicarts/pbmpm

• WebGPU implementation; ~64% JavaScript, ~33% WGSL (WebGPU Shading Language — the GPU compute shaders).
• License: BSD-3-Clause (commercial-friendly).
• No formal releases — a reference/demo codebase, not a drop-in engine.
• Runs in Chrome-family browsers; Firefox unsupported, mobile not properly supported.
• Demo scenes: blockCrusher, coiling, colliders, splashing.
• Its README points to nialltl’s incremental_mpm as a learning resource.

The rest of §3.3 is a close reading of both branches.

3.3.1 Conventions decoder

Their code and standard velocity-form MLS-MPM describe the same discretization in different clothes. Translation table, for reading their shaders — the right-hand column is the convention of a typical explicit velocity-form 3D MLS-MPM (normalized unit-cube domain, fixed-point scatter):

Their fixed-point discipline is the same as the right-hand column’s, for the same reason (WebGPU, like Metal, has no float atomics; integer sums commute → deterministic scatter — §2.1). Two small deltas: their multiplier is runtime-tunable (a precision dial worth remembering), and their encode truncates where the right-hand column rounds — rounding is the better half of that trade, since truncation carries a toward-zero bias.

3.3.2 The algorithm as implemented

The loop. The particle carries displacement and deformationDisplacement (D) instead of velocity and C. Per substep (sim.js, both branches; defaults dt = 1/240, iterationCount = 5):

emit / recycle particles
for iteration in 0..iterationCount:            # default 5 (UI: 2..100)
    particleUpdate     # constraint projection: nudge D toward the material's target
    gridZero
    particleToGrid     # scatter mass + momentum(d + D·offset) — NO stress term
    gridUpdate         # mass-normalize, collide, guardian — all in displacement form
    gridToParticle     # gather d and D (= B·4)
particleIntegrate      # x += d;  F ← (I+D)F;  plasticity; gravity → next d;  once

The entire transfer cycle runs inside the constraint loop. That is the whole trick: particleUpdate projects each particle’s D toward a constitutive target shape, and the P2G→grid→G2P pass in between is what makes neighboring particles negotiate — Jacobi-style position-based dynamics where the MPM grid is the constraint coupling. Nothing anywhere applies a force scaled by stiffness·dt², so nothing can explode; stiffness is expressed as a convergence rate (elasticRelaxation, liquidRelaxation, default 1.5 — over-relaxed). Stability is unconditional; what varies with iteration count is how stiff materials manage to look.

Warm start, for free: d and D persist across substeps as the next loop’s initial guess, and the integrate step adds gravity into the next substep’s initial displacement. Preserve this in any port — zeroing D between substeps silently costs iterations.

Materials as constraint projections (particleUpdatePBMPM.wgsl). Each material answers “what shape does F want to be, and how fast do we relax toward it”:

Their 2×2 SVD is closed-form (Gimeno’s rotation-angle construction, matrix.inc.wgsl) so they can afford one per particle per iteration. That luxury does not survive 3D: the elastic projection needs only $\mathbf{U}\mathbf{V}^{T}$, which is the polar rotation R, so a Newton polar iteration replaces it in the loop and the full SVD is paid once per substep at integrate time. Sand’s in-loop unilateral $\Sigma \ge 1$ clamp genuinely needs singular values and has no such escape.

Integration and plasticity (particleIntegrate.wgsl) runs once per substep. For liquids, F is never stored as a matrix at all — only a scalar liquidDensity, integrated by the first-order identity $\det((\mathbf{I}+\mathbf{D})\mathbf{F}) \approx (1 + \mathrm{tr}\,\mathbf{D})\det \mathbf{F}$. Solids do $\mathbf{F} \leftarrow (\mathbf{I}+\mathbf{D})\mathbf{F}$, then SVD, then:

• Safety clamp on singular values (≈ [0.1–0.2, 10⁴]) — their explosion backstop. Note they sanitize silently; a clamp that hides its own firing rate is how an invariant rots, so count first.
• Sand: Drucker-Prager return mapping in log-strain space — the Klár 2016 scheme, with logJp carrying the volume-debt role ($\beta = 0.5$). Expansion forgets all deformation ($\Sigma \leftarrow 1$).
• Visco: clamp Σ into $[1/y,\; y]$ where $y = \exp(1 - \mathrm{plasticity})$, then rescale Σ to restore the original determinant — plastic shear that conserves volume. Three lines.

Displacement-form boundaries. gridUpdate works on predicted position: it asks “would this node’s displaced position penetrate?” — colliders subtract only the radial component of displacement past contact (with a gap term so approach up to the surface is free), and the guardian projects the displaced node position back into bounds. Position-predictive BCs are naturally tunnel-resistant; the velocity-form equivalent is a sign test, which only knows “moving inward,” not “will cross this substep.” Worth remembering when next fighting a fast-material leak. (The main branch adds borderFriction mixing at the guardian; the sig branch zeroes.)

Grid-recovered volume. gridToParticle.wgsl carries a comment worth quoting in paraphrase: standard MPM volume integration slowly loses liquid volume under shear, so they recover an objective volume measure from the grid — scatter particle rest volume (their 4th grid channel; emitted as $1/\mathrm{ppc}$ per particle so a rest-packed node sums to ~1), gather it back, reciprocate, and blend into the tracked density only where the grid reads compressed (at a free surface the estimate reads spuriously expanded — fewer neighbors, not more volume — so tension is never trusted). It is an SPH-flavored density estimate used as a drift corrector, and it transplants cleanly into an explicit solver.

3.3.3 GPU-resident particle lifecycle

particleCount lives in a GPU buffer; a one-thread setIndirectArgs kernel refreshes indirect dispatch args; drains push freed indices onto an atomic-cursor free list; emitters pop it. Emission is grid-dispatched — one thread per cell decides whether its cell spawns this frame, so rate is per-area and collision-free by construction. Tombstone schemes (dead particles kept in place) cost dispatch threads forever and never reuse the slot; when a scene needs hoses, waterfalls, or blast-debris cleanup, this is the pattern, and it composes with the tiling in §2.1 because dead particles simply never get binned.

3.3.4 Read and rejected

Recorded so nobody re-litigates them from scratch:

• Per-particle positional push-out of colliders (particleIntegrate pushes penetrating particles out along the normal). This is an energy-injecting positional projection; even their comment concedes the grid should have done the work. A grid BC plus a bounded-penetration probe is the better pairing.
• Per-frame createBuffer in the JS driver (grid buffers allocated every update, GC’d). Demo-grade resource management; nothing to learn.
• Truncating fixed-point encode (i32(x·mult)) — toward-zero bias; round() is strictly better (3.3.1).
• Silent sanitization — their Σ clamps and density floors fix states without counting them.
• Mouse-grab sets displacement directly — cute for a demo; a kinematic grab in a game wants to go through the collider/impulse path, not teleport displacement.

3.3.5 Not in their scope (don’t look for it there)

3D. Multi-body contact (single velocity field, so two elastic bodies weld on contact — they simply do not have separable bodies). Cutting/CPIC. Two-way rigid coupling and articulation. Snow. Determinism as a verified property — their scheme is deterministic by the fixed-point argument of §2.1, but nothing tests it. Any renderer beyond a point splat.

3.3.6 File map (siggraph2024 branch)

3.3.7 License & origin

The repo is BSD-3-Clause (Electronic Arts, 2024 — EA’s stock variant, no added terms). Redistribution of their code, including transliterated shaders, requires retaining the copyright notice + conditions + disclaimer, a notices entry in binary distributions, and no use of EA’s name for promotion. Algorithms and ideas are not expression: a technique re-derived and re-written against a different architecture carries no notice obligation — but if anyone ever ports their WGSL wholesale, add the LICENSE.md text to a third-party notices file that ships with builds. Requested citation either way:

Chris Lewin. A Position Based Material Point Method. ACM SIGGRAPH 2024 Talks. ea.com/seed · paper PDF

See References for the full citation block.

What PB-MPM buys, and what it does not. Unconditional stability, not accuracy. Stiffness is a function of iteration count and dt — materials are “relaxation rates,” not Young’s moduli, and cannot be calibrated to physical units; the relaxation is inherently dissipative. It is the right solver for game-feel materials at display-rate substeps on mobile; an explicit MLS-MPM remains the right one for calibrated or violent dynamics.

4. Suggested path forward (for building with MPM)

1. Learn the algorithm from nialltl’s incremental examples (most readable, MIT, Unity-native).
2. Cross-reference Yuanming Hu’s 88-line MLS-MPM and the SIGGRAPH 2016 course (mpm.graphics) for deeper theory.
3. Scale to GPU with the Gao et al. techniques (spatial sort, sparse grid, fixed-point atomics, AoSoA) when particle counts grow.
4. If you hit the timestep/stability wall in a shipping context, study EA’s PB-MPM (siggraph2024 branch first), purpose-built to solve it.

Building this for a game specifically? The product-side questions — where MPM ships today, whether any commercial game runs it as core physics, and what makes an MPM-native game newly viable — are covered in Games.

5. Engineering rules that hold for any GPU MPM

Every GPU MPM implementation ends up rediscovering these. They follow from the structure of the method and of GPU execution, not from any one codebase, and each is stated with its reason, because the reason is what tells you when it applies.

1. Fixed-point integer atomics are what make the scatter deterministic. Floating-point addition is not associative and GPU thread order is not reproducible, so a float-atomic P2G produces a different grid on every run. Accumulate into scaled integers with InterlockedAdd, which is associative, and convert once in the grid update. Metal has no float atomics at all, so this is also the portable choice.
2. Particle arrays must stay index-stable. Sorting particles for memory locality, or compacting out dead ones, reorders the accumulation and destroys the determinism just bought. Kill a particle by relabeling and parking it — a dead group id, a position outside the domain — never by compaction. If spatial sorting is needed for performance, sort an index buffer, never the particle buffer itself.
3. Determinism is same-platform only. Same binary, same device, same seed → byte-identical replay. A different GPU architecture, shader compiler, or FMA-contraction policy → statistically equivalent at best. Golden-file regression tests work; cross-platform golden files do not. Corollary: any edit inside a hot kernel recompiles it, and the compiler may contract the same arithmetic differently, so a replay digest moves for a change that altered no logic. A digest is a same-code invariant, not a “nothing happened” check.
4. Consume-and-zero in place is free — only where each accumulator is read exactly once. The grid clear can be folded into the grid update when every node accumulator has one reader. A pass that reads neighboring nodes’ accumulators in the same dispatch — a multi-field contact pass, for instance — must keep its explicit clear.
5. dt is not a free performance knob. The timestep is entangled with oscillation, plastic yield rates, and boundary friction; raising it to buy frames changes the result, not just the cost. Treat (dt, gravity, substeps per frame) as one tuned set, and re-tune the others when one moves.
6. There are two timestep limits, and the second depends on the running state. The acoustic limit, $\Delta t \le \mathrm{CFL} \cdot \Delta x / c_p$, is a property of the material and can be read off a preset before it runs. The advective limit — material must not cross a cell in a step, $\Delta t \le \Delta x / \max_p\!\left(\lVert \mathbf{v}_p \rVert + 1.5\,\Delta x\,\lVert \mathbf{C}_p \rVert_F\right)$ for an APIC/MLS transfer, where the affine term counts because it advects too — is a property of the current state and has to be monitored, not precomputed. The smaller one sets the step.
7. Gravity and damping are environmental, not material. A material is the same on Earth and on the Moon. Keep gravity and damping on the solver, set per scene, and keep the material record to the constitutive law alone (E, ν, hardening, yield, cohesion, friction angle). Give each quantity one owner, or the two drift apart and the bug reads as “the material changed” when nothing did. A fixed per-substep velocity blend is a damper whose rate is $(1-\beta)/\Delta t$; expose a rate in 1/s and derive the per-step factor from it.
8. Keep synchronous readback off the per-frame path. A synchronous buffer read costs however long the GPU happens to be backed up, not the size of the copy, because it blocks until the queued work drains — “a few floats, once a second” can freeze a frame for seconds. Read back asynchronously, and keep any per-tick state readback O(groups), never O(particles).
9. Fixed-size GPU buffers truncate silently. An append buffer, collider array, or query result that overflows its capacity caps with no error and produces physics that is subtly, deterministically wrong. Size for the densest scene you can construct, and carry an overflow counter you actually check.
10. A dispatch has a per-axis group limit, and exceeding it fails silently and partially. A 256³ field cleared as 262,144 one-dimensional groups exceeds the 65,535-groups-per-axis limit of Direct3D and of the engines built on it. The clear partly does not happen, and with a freshly allocated buffer the first frame still looks correct, because zeroed memory and a successful clear are indistinguishable. Dispatch 3D fields with 3D groups.
11. A fixed-point readback has a quantum, and its error scales with particle count, not time. Summing a fixed-point per-particle quantity measures through a quantizer. The rounding errors cancel while the motion is lively and stop canceling when it goes quiet: thousands of particles sharing a similar small velocity round the same way, and the error becomes a bias. A conservation check that drifts only at rest is usually measuring its own quantum.

Sources for this document

• Gao, Wang, Wu, Pradhana, Sifakis, Yuksel, Jiang. GPU Optimization of Material Point Methods, SIGGRAPH Asia 2018 — PDF · code.
• nialltl. MLS-MPM guide & incremental_mpm — guide · code (MIT).
• Hu et al. MLS-MPM / taichi_mpm — project · code (MIT).
• Lewin (EA SEED). A Position Based Material Point Method, SIGGRAPH 2024 — code (BSD-3) · SEED publications.
• For the games / commercial-applications angle, see Games.
• See References for full citations and more.