Material Point Method: Step by Step

Every simulation method is a bet about what is hard. Rigid bodies bet that shape never changes. Cloth bets that connectivity never changes. The Material Point Method bets that connectivity is the thing you should never have had in the first place — and pays for that bet by renting a grid, one timestep at a time.

That single trade explains almost everything about how MPM code looks, why it is fast on a GPU, why it is slow on a clock, and why your snow sometimes behaves like jello. This page builds that intuition from arithmetic you can do on paper, with a worked transfer you can drag around, and ends by mapping the phases onto the kernels of a real solver.

1 · The bargain

There are two classical ways to simulate a deforming material, and each one fails in a characteristic, visible way.

A Lagrangian method glues its sample points to the material. Think of a mesh whose vertices ride along with the stuff. This is great bookkeeping: mass is exactly conserved, and every point remembers its own history. It is also brittle. Push the material hard enough and elements shear, flatten, and eventually turn inside out, at which point the math that computes forces from element shape starts returning nonsense.

An Eulerian method fixes its sample points in space and watches material flow past. Think of a fixed grid of cells. Nothing can tangle, because nothing moves. But every step you have to advect — copy each cell's contents into its neighbors — and that copy is a blur. Do it a few hundred times and a crisp block of snow has smeared into a fog with no edges and no memory of which grain came from where.

MPM's answer is to refuse to choose. Material lives on particles that carry mass, velocity and deformation history and move freely through space, with no notion of who their neighbors are. Once per step, every particle dumps its state onto a background grid, the grid does the one job a grid is genuinely good at — computing spatial derivatives and resolving contact — and the answer is handed straight back to the particles. Then the grid is thrown away.

What that buys, in the language of things you would otherwise have to build:

And the bill, which arrives later in this page: a small timestep, a lot of memory traffic, contact that is too sticky, and a resolution floor that means a thin object simply cannot be simulated no matter how you tune it.

2 · Two data structures

Before the algorithm, the two records it moves data between. Everything else in MPM is a consequence of what is in these two structs.

Note what is not on the particle: any reference to another particle. No neighbor list, no spring, no element index. Two particles that are about to collide have no idea the other exists. They will find out through the grid, and section 4 shows exactly how.

Note also what is not on the node: anything material. A node has no idea whether the mass sitting on it is snow or water. It is a place where momentum is added up and divided by mass.

3 · The four-phase loop

One MPM timestep is always the same four phases, in the same order. Every variant of MPM you will read about — APIC, MLS-MPM, PB-MPM — changes what is carried across these arrows, never the shape of the loop.

Two remarks that will pay off later. First, phase 2 is where every particle writes into nodes that other particles are also writing into, simultaneously. On a GPU that is a pile of atomic adds into contended addresses, and it is the single biggest performance problem in the method. Second, phase 3 is the only place where two separate pieces of material can influence each other. That is the whole contact model.

4 · Do one step by hand

This is the section that builds the intuition. We will take two particles moving toward each other, transfer them to a grid, solve, and transfer back — with real numbers at every stage. Nothing here is simplified away; these are the same formulas the compute shader runs, in 2D with $\mathrm{dx} = 1$ so the arithmetic is legible.

4.1 · Which nodes does a particle touch?

A particle at position x does not talk to "the cell it is in". It talks to a 3×3 block of nodes (3×3×3 in 3D) chosen so the particle sits near the middle. The solver finds that block with two lines:

int3  base = int3(floor(p.x * invdx - 0.5));   // lowest node of the 3×3×3 block
float3 fx  = p.x * invdx - float3(base);       // where in the block, in [0.5, 1.5)

fx is the particle's offset from that lowest node, measured in cells. The −0.5 is what centers the block: it guarantees fx always lands in [0.5, 1.5), so the particle is always between the first and third node with the middle node nearest. Take particle A at $x = 2.3$: $\mathrm{base} = \lfloor 1.8 \rfloor = 1$ and $\mathrm{fx} = 1.3$. It touches nodes 1, 2 and 3 — and it is closest to node 3, so node 3 will get the most of it.

4.2 · How much of the particle does each node get?

The share each node gets is its weight. MPM uses a quadratic B-spline, which in code is three lines with no branches:

w0 = 0.5 * (1.5 - fx) * (1.5 - fx);
w1 = 0.75 - (fx - 1.0) * (fx - 1.0);
w2 = 0.5 * (fx - 0.5) * (fx - 0.5);   // w0 + w1 + w2 == 1, always

For A's $\mathrm{fx} = 1.3$: $w = (0.02,\; 0.66,\; 0.32)$. They sum to exactly 1 — that property is called partition of unity, and it is what makes the transfer conserve mass. In 2D the weight for a node is the product of its x-weight and its y-weight; in 3D, all three. Nine numbers per particle in 2D, twenty-seven in 3D.

4.3 · Scatter, solve, gather

Now the whole cycle, with two particles heading at each other at ±1 unit/s. Step through it, drag them closer or further apart, and watch what happens at the nodes they share. Neither particle ever reads the other.

Walk the phases and read the numbers. The four things worth noticing, in order:

1 — the shared node is a mass-weighted average

At the default spacing, node (3, 3) receives weight 0.189 from A (moving right) and weight 0.437 from B (moving left). Its momentum is $0.189\cdot(+1) + 0.437\cdot(-1) = -0.248$ and its mass is 0.626, so its velocity is −0.396. A node cannot hold two velocities. Divide momentum by mass and the disagreement is resolved, permanently, in favor of whichever side brought more mass. That single division is MPM's entire collision response.

2 — the particles change velocity without ever meeting

In phase 4 each particle gathers $\sum w\,v_\mathrm{node}$ over its own nine nodes. A's stencil now contains nodes that B polluted, so A comes back slower than $+1$. B likewise. Neither did a distance check, a raycast, or a neighbor query. This is what "self-collision is automatic" actually means at the arithmetic level.

3 — momentum is conserved exactly; energy is not

The verdict bar shows total momentum before and after: identical to the last digit. That is a theorem, not a coincidence — it follows directly from partition of unity, and it holds for every MPM transfer. Total kinetic energy, meanwhile, drops hard.

Two different things are behind that drop, and it is worth separating them because only one is a defect. The first is real physics we have deliberately switched off: this lab scatters bare momentum, with no stress term, so two masses meeting head-on simply merge their motion — a perfectly inelastic collision. In the full solver that energy is not lost; it is stored in each particle's deformation gradient and pushed back out as stress on the following step. The second is transfer dissipation, and it is unrecoverable: any motion the grid's one-vector-per-node representation cannot express — rotation, and any structure finer than a cell — averages to nothing on the way in and does not come back. That second part is the defining flaw of the scheme, and section 6 is about how much of it can be clawed back.

4 — separate them far enough and nothing happens

Push the gap up and watch the shared-node count step down: six shared nodes, then three, then at a gap of 4.5 — about 2.2 cells of separation — zero, and the exchange stops dead. The particles pass without ever affecting each other. The interaction range of MPM is exactly the stencil width — three cells, from the shape function. It is not a tunable radius and not a material property. A grid cell is the length scale of everything.

4.4 · What the real kernel adds

The lab above scatters plain momentum, $w\,m\,v$. The production kernel scatters $w\,m\,(v + C\,\mathrm{dpos})$, where dpos is the offset from particle to node and C is a 3×3 matrix. Two separate things are folded into that one matrix:

• the particle's affine velocity field — its local spin and stretch, which plain averaging would throw away (section 6);
• the particle's internal stress, converted to a force. This is the MLS-MPM trick: instead of evaluating shape-function gradients, the stress is scaled by $-\tfrac{4}{dx^{2}}\,dt$ and simply added to C, so one scatter delivers momentum and force together.

// an MLS-MPM P2G kernel (HLSL), the lines that matter
float3x3 stress = (-4.0 * invdx * invdx * _Dt) * PF;   // force, as an affine term
float3x3 C      = stress + p.C;                        // + the particle's own affine field
...
float3 mom = particleMass * (p.v + mul(C, dpos));      // one scatter carries both

That fusion is why modern real-time MPM starts from MLS-MPM rather than classical MPM: it removes an entire gradient evaluation from the innermost loop, at no cost in behavior.

5 · Where the weights come from

The three-line weight function is doing more work than it looks. Here is the same B-spline drawn across a row of nodes, so you can see the property that matters: at every position, exactly three bumps are non-zero and they always sum to one.

5.1 · Why not the obvious linear weights?

You could split a particle between the two nodes it sits between, linearly. That works, it is cheaper, and it produces a notorious artifact called cell-crossing instability. The reason is not in the weight — it is in the weight's derivative, which is what force depends on.

6 · What the transfer forgets

In the lab you saw total kinetic energy drop on every round trip. Here is the mechanism, and it is worth seeing because it explains three decades of method names.

Consider four particles rotating around a single node. Their velocities are tangential, in four different directions. Scatter them onto the node and they cancel: the node's momentum is zero, so its velocity is zero. Gather back and all four particles come home dead. The rotation is gone, and there was never a moment where anything "lost" it — the node simply has no slot in which spin could be stored.

The family tree, which is just a sequence of answers to "what should we transfer back?":

7 · F — the memory in a particle

Everything so far would give you a fluid with no memory. What makes MPM able to simulate a rubber ball, packed snow, and a sand pile is one 3×3 matrix per particle: the deformation gradient, F.

F answers one question: if I drew a tiny arrow inside this material when it was at rest, what would that arrow look like now? It is the local linear map from rest shape to current shape. That is all.

7.1 · From F to force

A constitutive model is a pure function from F to stress. That is the entire definition of "material" in MPM — swap the function, get a different substance, change nothing else in the solver. The elastic model graphics MPM has favored for snow and rubber since the Disney paper is fixed corotated, and the intuition behind it is worth a sentence:

float3x3 R  = polar_rotation(F);           // the nearest pure rotation to F
PF  = (2.0 * mu) * mul(F - R, transpose(F));   // shear: penalize deviation from that rotation
PF += lambda * (J - 1.0) * J * I;              // volume: penalize J drifting from 1

R is the rotation part of F. Subtracting it means rotation is free — spin a rubber block and it feels nothing, which is exactly right and is the whole reason this model is preferred over the naive alternatives. mu and lambda are the two stiffness constants derived from Young's modulus E and Poisson's ratio nu; mu resists shape change, lambda resists volume change.

7.2 · Plasticity is a clamp

Elastic materials spring back. Snow packs, sand piles, and mud slumps — they keep the dent. The mechanism is startlingly simple: split F into an elastic part that makes stress and a plastic part that does not, and enforce the split by clamping.

Factor F with an SVD into $\mathbf{U}\Sigma\mathbf{V}^{T}$. The diagonal of $\Sigma$ holds the singular values — the stretch factors along the principal axes. Snow plasticity clamps each one into $[1-\theta_c,\; 1+\theta_s]$. Whatever you clamped away is now permanent deformation. Sand does the same thing against a cone rather than a box.

8 · The material zoo

With F, a stress function, and a yield surface in place, a material is just a row of constants. Here is a worked palette from one production solver — the names are its own, but the shape of the table is what every MPM material set looks like.

8.1 · Water is a different shape of code

Notice that Water has no Young's modulus. A weakly-compressible fluid does not resist shape change at all — only volume change — so the entire deviatoric branch, and with it the per-particle SVD, disappears:

// fluid: pressure only, no shear term, no SVD on this path
float pres = bulkK * (pow(max(J, 1e-4), -eosGamma) - 1.0);
float dd   = -pres * J;
PF = float3x3(dd,0,0, 0,dd,0, 0,0,dd);   // purely isotropic

That is why the fluid path is meaningfully cheaper per particle than the solid path, and why bulkModulus for a fluid is doing double duty: it sets both how crisp a splash looks and how strongly a floating object is pushed up.

9 · Tuning a scenario

This is where intuition earns its keep, because most MPM parameters do not do what their name suggests, and several of them are entangled. The following is the field-tested version.

9.1 · Which knob owns what

9.2 · The five traps

9.3 · Symptom → knob

9.4 · How to tune, as a process

Nobody has calibrated intuition for MPM parameters — not even people who have shipped one. Twiddling a dial and asking "is that better?" does not converge. Two rules make it converge:

1. Physics first on a deliberately ugly render. Draw particles as big soft overlapping blobs, lighting off, iso-surface off, and tune motion until the material reads as the right substance. Human vision is biased toward motion, and a pretty render hides motion behind silhouette and shading artifacts. Only once the motion is right do you switch on the real renderer and tune the look — against locked physics. The debug renderer is standing tooling, not scaffolding to delete.
2. Side-by-side, one axis per round. Run three or four deterministic sims simultaneously that differ in exactly one parameter and ask "which of these is more X" — a question vision answers far better than "is this good". Keep the previous champion on screen so you notice when round 2 gains something and quietly loses what round 1 had. Bracket wide first to find the regime boundaries on both sides, then bisect. Lock the scene seed across resets, or every variant is judged on a different scene and the comparison is noise that looks like signal.

9.5 · The performance model, in one figure

The last thing to internalize before touching a config: what actually costs time. It is not the particle count in memory. It is the number of times per second you stream all of that state through the two transfers, which the CFL condition puts a floor under.

10 · Reading a solver

Everything above maps onto a small number of places in any MLS-MPM codebase, whatever the language. Find these eight first and the rest of the file reads itself. An explicit MLS-MPM solver is the one to read first; PB-MPM is a different bargain and should be read second.

10.1 · Two implementation details worth knowing up front

A deterministic scatter uses integer atomics, on purpose. Floating-point addition is not associative, so the order in which particles land on a node would change the result — and GPU thread order is not reproducible. A scatter that accumulates into fixed-point uint values with InterlockedAdd, which is associative, buys bit-identical replays on the same machine. It is also the only portable option: Metal has no float atomics at all.

The grid clear can be folded into the grid update. If each node's accumulator is read exactly once per step, in the grid-update kernel, it can be zeroed right there and the separate clear dispatch disappears. This works only because of the "exactly once" property — a pass that reads neighboring nodes' accumulators during the same dispatch (a multi-field contact pass, for instance) has to keep its explicit clear.

11 · Other MPM tutorials

Tutorials and guides only — no papers. Each entry says what it is actually good for, because they are not interchangeable.

Written guides & courses

Read-the-code walkthroughs

These are implementations small enough to be read as tutorials, which for this method is often the faster path than prose.