GPU and Massively Parallel Simulation

A GPU has thousands of lanes. A physics step wants to use all of them. There are two ways to get there, and they are not equally easy.

Two kinds of parallelism

The first is batched environments: run thousands of independent simulations side by side, one per GPU lane, sharing nothing. This is embarrassingly parallel, and it is how end-to-end GPU reinforcement learning collapses training time. A single workstation GPU can step tens of thousands of copies at once, buying two to three orders of magnitude over a CPU farm.¹ Nothing about the physics changes; you just stamp out many worlds.

The second is parallelism inside one solve. This is the hard one. The constraint solvers that engines rely on are variants of Gauss-Seidel, and Gauss-Seidel is sequential by construction: each constraint reads the results of the ones before it. To parallelize a single large world we have to break that dependency without giving up the fast convergence that made Gauss-Seidel worth using.

The matrix-splitting view

Both classical iterations come from splitting the system matrix. For $A\mathbf{x}=\mathbf{b}$, write $A$ as its diagonal, strict lower, and strict upper parts:

$$A = D + L + U.$$

Jacobi keeps only $D$ on the left and pushes everything else to the right, evaluated at the old iterate:

$$\mathbf{x}^{k+1} = D^{-1}\bigl(\mathbf{b} - (L+U)\,\mathbf{x}^{k}\bigr).$$

Per component this is $x_i^{k+1} = \bigl(b_i - \sum_{j\ne i} A_{ij}\,x_j^{k}\bigr)/A_{ii}$. Every unknown depends only on old values, so all of them update at once. That is perfect for a GPU, but Jacobi converges slowly.

Gauss-Seidel moves $L$ to the left as well, so updates use the fresh values already computed this sweep:

$$\mathbf{x}^{k+1} = (D+L)^{-1}\bigl(\mathbf{b} - U\,\mathbf{x}^{k}\bigr), \qquad x_i^{k+1} = \frac{1}{A_{ii}}\Bigl(b_i - \sum_{j<i} A_{ij}\,x_j^{k+1} - \sum_{j>i} A_{ij}\,x_j^{k}\Bigr).$$

The in-place substitution roughly doubles the convergence rate, but it forces $x_2$ to wait for $x_1$, and $x_3$ for $x_2$. Sequential.

A worked step

Take a small tridiagonal system:

$$A=\begin{bmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{bmatrix},\qquad \mathbf{b}=\begin{bmatrix}2\\4\\10\end{bmatrix},\qquad \mathbf{x}^{*}=\begin{bmatrix}1\\2\\3\end{bmatrix}.$$

Here $D=\operatorname{diag}(4,4,4)$ and $L+U$ holds the $-1$ off-diagonals. Start from $\mathbf{x}^{0}=\mathbf 0$. The first Jacobi step is just $D^{-1}\mathbf b$, since $(L+U)\mathbf x^{0}=\mathbf 0$:

$$\mathbf{x}^{1} = D^{-1}\mathbf{b} = \begin{bmatrix}2/4\\4/4\\10/4\end{bmatrix} = \begin{bmatrix}0.5\\1.0\\2.5\end{bmatrix}.$$

Now one more Jacobi step, feeding $\mathbf x^{1}$ into $x_i^{2}=\bigl(b_i - \sum_{j\ne i} A_{ij}\,x_j^{1}\bigr)/A_{ii}$:

$$\begin{aligned} x_1^{2} &= \tfrac{1}{4}\bigl(2 - (-1)(1.0)\bigr) = \tfrac{3.0}{4} = 0.75,\\[2pt] x_2^{2} &= \tfrac{1}{4}\bigl(4 - (-1)(0.5) - (-1)(2.5)\bigr) = \tfrac{7.0}{4} = 1.75,\\[2pt] x_3^{2} &= \tfrac{1}{4}\bigl(10 - (-1)(1.0)\bigr) = \tfrac{11.0}{4} = 2.75. \end{aligned}$$

So $\mathbf x^{2}=(0.75,\,1.75,\,2.75)$, crawling toward $(1,2,3)$. For contrast, a single Gauss-Seidel sweep from $\mathbf 0$ reuses each value immediately: $x_1=0.5$, then $x_2=(4+0.5)/4=1.125$, then $x_3=(10+1.125)/4=2.781$. One sweep, and it is already closer. That head start is exactly what we do not want to throw away when we go parallel.

Graph coloring: Gauss-Seidel in parallel

The fix is to notice why Gauss-Seidel is sequential. Two constraints conflict only when they share a body, where one writes a value the other is about to read. Constraints that touch disjoint bodies never conflict and can be applied simultaneously with no loss of freshness.

So build the constraint graph, one node per body, one edge per constraint, and color the edges so that no two edges meeting at a node share a color. Each color class is then a set of mutually independent constraints. Sweep the colors in turn; within a color, fire every constraint in parallel. With $K$ colors you get $K$ sequential parallel passes per iteration, and the result matches Gauss-Seidel because no body is ever touched twice inside one pass.²³

Coloring a lattice of distance constraintsinteractive

width W = 4height H = 3

constraints

17

colors K

4

parallel passes

4

◀ prevnext ▶all

Color 1 of 4: these 6 constraints touch disjoint bodies and solve in parallel.

The trade-off is now visible. More colors means more sequential passes, so a good coloring is one with few colors and large, balanced color classes. This is the structure behind GPU-friendly Gauss-Seidel solvers for cloth, soft bodies, and particle systems, where a lattice colors with a handful of colors and each pass saturates the device.²³

Batched worlds give parallelism for free; parallelism inside a solve has to be earned by coloring. With the solver mapped onto the hardware, the last part asks what happens when we want to differentiate straight through it.

Footnotes

1. Makoviychuk et al., Isaac Gym: High Performance GPU-Based Physics Simulation For Robot Learning (2021). ↩

2. Fratarcangeli, Tibaldo & Pellacini, Vivace: a Practical Gauss-Seidel Method for Stable Soft Body Dynamics (SIGGRAPH Asia 2016). ↩ ↩²

3. Macklin, Müller, Chentanez & Kim, Unified Particle Physics for Real-Time Applications (SIGGRAPH 2014). ↩ ↩²