Solver Internals

Three lines of arithmetic

The whole transfer basis is one function, in shaders/common/mpm_utils.glsl:

Node	Weight
w0 (left / below)	0.5 * (1.5 - fx) * (1.5 - fx)
w1 (centre)	0.75 - (fx - 1.0) * (fx - 1.0)
w2 (right / above)	0.5 * (fx - 0.5) * (fx - 0.5)

That is the quadratic B-spline. It has three properties the solver depends on and none of them are optional:

• The three weights sum to 1 for every fx. This is what makes the transfer conserve mass. Any bug that breaks it leaks or invents material.
• It is C¹ continuous. The weights and their first derivatives vary smoothly as a particle crosses a cell boundary, so nothing pops when it does. Linear weights would be cheaper and would visibly stutter.
• It has compact support of exactly 3 cells. Nine nodes in 2D, always. Not a neighbour search, not a variable count — a fixed double loop.

The - 0.5 that centres it

The stencil origin is computed like this, identically in five shaders:

cell_pos = (pos - grid_origin) / dx
base     = floor(cell_pos - 0.5)
fx       = cell_pos - base

Without the - 0.5, base would be the cell the particle is standing in and the 3×3 block would be lopsided — one node behind, two ahead. With it, base + 1 is the nearest node, so the block is centred and fx lands in [0.5, 1.5] rather than [0, 1). That shifted range is exactly the domain over which the three weight expressions above are the correct B-spline pieces. The offset and the weight formulas are one design, not two.

One constant instead of a matrix inverse

APIC gives each particle an affine velocity matrix C so that momentum transfers carry local rotation and shear instead of losing them. The general form needs an inertia-like matrix D inverted per particle per step.

For a quadratic B-spline, D is always (1/4) · dx² · I. It does not depend on where the particle is inside its cell. So D⁻¹ is the compile-time constant 4 / dx², and the inverse disappears from the inner loop entirely. That is the MLS-MPM identity (Hu et al., 2018), and it is why you see this literal 4.0 in both directions:

// p2g
affine = (-dt * p_vol * 4.0 * inv_dx * inv_dx) * stress + mass * C;

// g2p
new_C += 4.0 * inv_dx * inv_dx * weight * outer(g_vel, dpos);

Change the interpolation order and that constant is wrong. It is not a tuning parameter.

Stress and APIC in one term

Look at the P2G line again. It fuses two things that are conceptually separate:

• -dt · p_vol · D⁻¹ · stress — the internal force the material exerts, integrated over the step.
• mass · C — the APIC affine momentum.

Both end up multiplying the same dpos in the scatter, so MLS-MPM adds them once and multiplies once:

momentum = weight * (mass * vel + affine * dpos) * z_occ

The practical consequence is that the constitutive model is the only thing a material author writes. Every material branch in mpm_p2g.comp — all sixty-odd of them — exists to produce one mat2 stress. The transfer machinery below it never changes. That is why a new material is a branch and not a subsystem.

p_vol is p_spacing², set once on the CPU at mpm_solver.cpp:1011.

Order matters, and it is deliberate

mpm_grid_op.comp runs once per grid cell, and its ordering encodes several decisions that are easy to get wrong:

Why the grid is thrown away

The grid holds no state between substeps. It is zeroed, scattered into, operated on, gathered from, and discarded. Nothing about a node survives to the next step.

This is what makes MPM behave like a mesh-free method despite having a mesh: the grid is a scratchpad for resolving contact and internal forces, not a representation of the material. Material state lives entirely on particles — F, C, temperature, and the per-material lanes.

It also means the grid resolution is a purely numerical knob. Raising it makes contact sharper and self-collision tighter without changing what any material is. Lowering it makes thin features merge, because two particles closer than a cell can no longer be told apart by the only machinery that lets them interact.

Everything after the gather is per-particle

Once new_vel and new_C exist, mpm_g2p.comp keeps going for several hundred more lines. Buoyancy for burning particles, air drag, the whole magnetic force family, held-magnet dragging, bulk lift, chain force, surface Maxwell stress, corrugation, bridging. All of them do the same thing: new_vel += force * u_dt.

None of them touch the grid. They cannot — the grid pass has already run and will not run again until the next substep.

The consequence

A force applied in the tail is not seen by any other particle until the next P2G scatter. It arrives one substep late, and it arrives having already moved the particle it acted on. For most forces this is invisible. For anything that should produce a reaction — a magnet pulling ferro toward it should feel the pull back — the reaction simply is not there.

That is not a bug in the tail. It is the reason the magnetic force pipeline is drawn the way it is on the Relationship Map: the field solve is two-way, but the force application is not.

The material-specific speed caps

The gathered velocity is clamped, and the cap is per material. The default is 3 cells per substep; the magnetic family gets far less:

Material	Cap	Why
default	3.00	The CFL-ish safety limit; the same number appears in P2G.
para mist (69)	1.40	A light aerosol should move briskly.
rubber (58)	1.20	Room for bend response.
heavy ferro (67)	0.82	Headroom for cursor-following lift.
permanent ferro (81)	0.72	Heavy and remanent, but must still wick upward.
ferrofluid (31)	0.60	Must chase a moving pole without going ballistic.
diamagnetic (68)	0.50	Moderate cap; it is allowed to drift away.
sticky ferro (70)	0.48	Stays cohesive, so it stays slow.

These are not physical properties. They are an admission that the magnetic forces can exceed what the transfer can carry stably, and the cap is where that is absorbed. When the clamp fires it scales new_C by the same factor, so the affine field never outruns the velocity it was built from — clamping velocity alone would leave a particle spinning faster than it translates.

Fixed point, because there is no atomic float

GLSL has no portable atomic float add on an SSBO. The scatter needs one — nine nodes, many particles, all racing. So the four grid accumulators are declared int and everything crosses scaled:

const float FIXED_SCALE = 1000000.0;
atomicAdd(grid_mass[cell], int(m * FIXED_SCALE));

Property	Value
Resolution	10⁻⁶ of a unit
Range before int32 wraps	±2147.48 units at a single node
Temperature accumulator scale	FIXED_SCALE * 0.001

The temperature channel is deliberately a thousand times coarser. It accumulates mass × temp, and temperatures reach 1500 K, so at the full scale a hot, dense node would be much closer to the wrap point than a cold one. Dropping three decimal places buys back that headroom. It is the same trade every fixed-point format makes, made once, per channel.

The velocity clamp, twice

The same 3-cells-per-substep limit appears on both sides of the grid:

Before the scatter (mpm_p2g.comp), because a particle moving more than its own stencil width per step would deposit its momentum into nodes it has already left. This one also catches NaN explicitly — spd != spd zeroes the velocity rather than propagating a NaN into the grid, where it would poison every particle that shares a node.

After the gather (mpm_g2p.comp), per material, as above.

Neither clamp is a physical model, and both are load-bearing. A single NaN reaching atomicAdd is not recoverable at the particle level: the node is corrupt, and every particle that gathers from it next step inherits the corruption. The scatter-side check is the last place that can be stopped cheaply.

Sixteen storage blocks

mpm_p2g.comp, mpm_g2p.comp and sph_force.comp each declare exactly 16 std430 buffer blocks. That is not a coincidence and it is not a style choice — it is the driver’s cap on shader storage blocks per stage, and all three shaders are sitting on it.

Anything that needs new per-particle state in those shaders cannot simply add a buffer. The escape hatches actually used in this codebase, in rough order of preference:

• A samplerBuffer lookup table on a free texture unit — textures are a separate binding namespace. This is how the authored material, bake, thermal and fracture profile tables get in.
• An imageBuffer, likewise a separate namespace, when writes are needed.
• A dedicated compute pass with its own binding budget — how the damage scatter/gather pair exists at all.
• Bit-packing into an existing lane. material_ids carries the z-occupancy in bits 8–15, which is why every read of it is masked with & 0x7Fu.
• Reusing a lane that a given material does not need — the explosives borrow Jp_data for crack, described on The Explosive Lab.

Two particle index spaces

This is the single most dangerous detail in the codebase, because getting it wrong produces no error of any kind.

Space	Index	Sized	Buffers
Global	p = u_offset + idx	SPH + MPM capacity	positions, velocities, material_ids, masses, colors, temperatures
MPM-local	idx	MPM capacity only	damage, Jp_data, fiber, phase_data, deform_grads, affine_moms

The offset is the SPH capacity, currently 300000. Both index forms appear within a few lines of each other in essentially every material branch:

temperatures[p] = clamp(temp, 0.0, 1400.0);
damage[idx]    = fuse;

What this page does not cover

• The constitutive models themselves — the sixty-odd material branches that produce stress. Each is its own design; several are documented on The Explosive Lab and the Technique Guide.
• The SVD. Every branch that needs principal stretches calls svd2x2(F) and works in sigma; the decomposition itself is a closed-form 2×2 routine and is not interesting.
• The damage scatter/gather pair, which is a second, independent grid round trip layered on top of this one for the AnisoMPM fracture model.
• Everything outside MPM: SPH, the Eulerian air solve, XPBD, and the field solves each have their own step. The frame loop figure places them in order.

Reading the shaders yourself

If you go looking, the order that makes sense is:

• shaders/common/mpm_utils.glsl — twenty lines, and the weights are all of it.
• shaders/physics/mpm_p2g.comp, from the bottom. The last forty lines are the transfer. Everything above them is one enormous else if (mat_type == ...) chain producing stress, and you can skip all of it on a first read.
• shaders/physics/mpm_grid_op.comp — short, and reads top to bottom in the order the grid pass actually happens.
• shaders/physics/mpm_g2p.comp — the gather is the first thirty lines. The rest is the tail.

One stale comment to ignore: mpm_utils.glsl describes an GL_NV_shader_atomic_float path and says the P2G shader declares its grid buffers as float. It does not. They are declared int at mpm_p2g.comp:21-24 and the transfer is fixed-point, as described above.