Sequential Impulse Solvers in GLSL Compute

This article provides an overview of how the Projected Gauss-Seidel (PGS) algorithm, widely known in game physics as the Sequential Impulse solver, resolves rigid body contact impulses within GLSL compute shaders. It explores the underlying mathematics of velocity-level constraint resolution, addresses the race conditions inherent to massively parallel GPU architectures, and details the graph-coloring techniques used to execute the solver efficiently on modern graphics hardware.

The Core Algorithm: Projected Gauss-Seidel

The definitive algorithm for resolving rigid body contact impulses is the Projected Gauss-Seidel (PGS) method, adapted by Erin Catto as the Sequential Impulse (SI) solver. In rigid body dynamics, non-penetration and friction constraints form a Linear Complementarity Problem (LCP). Rather than solving this system simultaneously using complex matrix inversions, the Sequential Impulse method iterates through individual constraints one by one, calculating a corrective impulse and immediately projecting (clamping) the accumulated impulse within valid physical limits.

For a contact constraint between body \(A\) and body \(B\), the relative contact velocity \(v_{rel}\) along the contact normal \(n\) must satisfy:

\[J v \ge 0\]

Where \(J\) is the constraint Jacobian. The impulse magnitude \(\Delta \lambda\) required to correct relative penetration velocity is calculated as:

\[\Delta \lambda = -\frac{J v + b}{J M^{-1} J^T}\]

Here, \(M^{-1}\) represents the inverse mass and inertia tensor matrix, and \(b\) is the Baumgarte stabilization term used to resolve positional penetration (penetration slop). The accumulated impulse is clamped to zero to ensure bodies only push away from each other and never pull:

\[\lambda_{accum} = \max(\lambda_{accum} + \Delta \lambda, 0)\]

Friction is resolved sequentially along two orthogonal tangent directions, clamping the frictional impulse between \(-\mu \lambda_{normal}\) and \(\mu \lambda_{normal}\) according to Coulomb’s law.

Adapting Sequential Impulses to GLSL Compute

The traditional Sequential Impulse algorithm is inherently serial: resolving one contact updates the velocities of its bodies immediately, which feeds into the next contact's resolution. In a GLSL compute shader, executing this process naively causes race conditions when multiple threads attempt to write to the velocity buffers of the same rigid body simultaneously.

To resolve this on the GPU, physics engines use one of two primary architectural patterns:

1. Contact Graph Coloring

The contact graph consists of rigid bodies as nodes and contact pairs as edges. A graph-coloring algorithm executed on the CPU (or in a pre-pass compute shader) assigns a "color" to each contact such that no two contacts of the same color share a rigid body.

In GLSL, the solver iterates through each color phase sequentially. Within a single color phase, an individual compute dispatch processes all non-interfering contacts in parallel across GPU workgroups. Because no two threads touch the same body velocity in a given phase, thread contention and write hazards are eliminated without requiring atomic operations.

2. Jacobi Solver with Under-Relaxation

Alternatively, a Jacobi-style formulation computes all constraint impulses simultaneously based strictly on velocities from the start of the iteration. The resulting impulses are accumulated on each rigid body using atomic additions (atomicAdd) or reduction passes. Because standard Jacobi solvers are prone to oscillation and instability, an under-relaxation parameter (typically between 0.2 and 0.5) must be applied to the calculated impulses to ensure convergence.

Compute Shader Execution Pipeline

A typical GLSL compute pipeline implementing the Sequential Impulse solver operates in distinct stages:

  1. Broad-phase and Narrow-phase: Contacts and collision manifolds (contact points, normals, penetrations) are detected and stored in uniform/storage buffers (SSBO).
  2. Setup Pass: Compute threads calculate the constraint mass (\(J M^{-1} J^T\)) and initial warm-starting impulses, applying previous frame impulses directly to the velocity buffers.
  3. Solver Loop: A host loop repeatedly dispatches the impulse resolution shader across the colored constraint batches (typically 4 to 10 iterations):
    • Read current linear and angular velocities from the body SSBO.
    • Calculate normal and friction impulse deltas.
    • Clamp accumulated impulses.
    • Write updated velocities back to the body SSBO.
  4. Integration Pass: A final compute pass updates rigid body positions and orientations using the newly resolved velocities:
// GLSL snippet for linear and angular velocity update
vec3 deltaLinearVelocityA = -impulse * invMassA;
vec3 deltaAngularVelocityA = -invInertiaA * cross(contactArmA, impulse);

velocities[bodyA].linear += deltaLinearVelocityA;
velocities[bodyA].angular += deltaAngularVelocityA;

Through graph-colored Projected Gauss-Seidel or under-relaxed Jacobi methods, the Sequential Impulse framework operates efficiently in GLSL, delivering stable, real-time rigid body dynamics on modern GPUs.