A solver of one's own

A minimal GPU collision solver

This solver is available on GitHub under MIT License.

This is a collision solver I have designed on and off for the past few months. It's intended for use on Apple Silicon Macs, although it can be adapted to run on CUDA with minimal modifications. On a 12-core M4 Pro, it can run scenes with 5,000 active box-shaped solids at around 50fps, can form stable heaps of thousands of solids, and requires relatively little in terms of memory. It uses an energy-conserving contact model, with dissipation added explicitly on top.

Clip 1 - the scene I started work on this solver for, the initial goal was to animate a rockface crumbling. Color coding optional

Clip 2 - the solver, with the help of its built-in sleep logic, can hold up pretty convincing heaps, here in the expected conical shape assumed by granular materials

Clip 3 - A scree giving way under load. Here solids are placed on a slope and spawned above the camera view.

Earlier this year, I wanted to attempt to build a scene for a set of VJing tools I'm putting together, with the idea behind them being to use physics-driven animations for visual rendering. My starting point was that using physical systems reads as more "real" than preprogrammed or parametric displacement and animations, and I wanted to try and see how far LLM tooling could get me, as reactive animation in OpenGL or a game engine would have been close to hopeless for me a few years ago. I didn't want to have to deal with a full physics engine for simple visuals, and the whole point of the project was to derive movement and visuals from physics using code I could understand, modify and adapt instead of relying on convoluted compositing pipelines in dedicated software. I also believed energy correctness could get me movement that reads more natural, and I liked the idea of trying to integrate "correct" physics into my visuals.

On the topic of LLM coding tools, this is the first time I have let LLM tools "drive" the development process to such a degree, while finding out about the wonderful world of physics solvers along the way. It started with what I can only call hubris in hindsight: I initially believed it wasn't that hard of a problem to simply take the equations of motion in discrete time and run them step by step. This is my attempt at untangling the mess, and also proves the point that real expertise still matters: I have no doubt someone with the relevant expertise would have built something better and faster with the same tools.

In terms of attribution, this is the farthest thing from hand-coded I have built short of throwaway scripts and small utilities for my own use (because let's face it, a working, optimized penalty solver produced by someone with zero prior experience in under six months, part-time would be too good to be true), this is roughly how the work split:

optimization was automated, while I specified the metrics, the invariants to watch, the validation gates and the scenes to measure as well as what to profile. However I made few direct optimization decisions myself, and essentially let the model go wild with a CPU implementation as the reference.
literature search was hybrid, some papers recommended by models, some I found using old-school Google Scholar and ResearchGate
decisions about which invariants mattered, the tradeoffs (stability vs penetration vs correctness vs speed) were mine, as were the contact models and the properties they should have. The choice of penalty over constraint fell out of this.

This roughly maps to how well I understand each bit: I am by no means a performance or GPU computation engineer, and I mostly stayed out of the optimization work, while I felt it was important to understand the maths behind all this, their implications, and more than anything else, the actual physics of the solver: its invariants, conservation laws, symmetries, and the assumptions it makes. At the same time, I will admit I have only a dim idea of the exact memory management schemes, thread dispatch tricks and broadphase lookup methods, which is why I didn't dwell too much on these.

Still, the finished solver leans more "AI generated" than I would have liked, and compared to domains closer to my core domain of expertise (signal and image processing), I had the opportunity to observe first-hand how one can be sucked into an illusion of understanding what comes out of the machine and to claim ownership over it when, in reality, one has only a vague and hazy mental map of how one's code actually functions. In a way, this page is something that exists as much for my own benefit as for its readers, since it was also how I managed to build this understanding.

As for the solver itself, it's something I am now quite happy with, at least enough to share it on my GitHub. I intend to use it for my own interactive visuals, and maybe to expand it a bit in the future to include 3D solids of arbitrary shape, linkages, or tweak it to handle soft body physics. I like the fact I can get natural-looking piles and "rockfalls" out of it, and that tuning it is easier than the contact solvers I did try, since its parameters map to physical values shared by real materials and legible to anyone trained in classical mechanics.

The most important lesson here for me was the fundamental compromise inherent to any collision solver: interpenetration, energy correctness and cost all draw from the same mostly fixed budget. And even running on GPU with a lot of optimization and scientific computation frameworks, my "flagship" piece so far (the one I started work on a physics solver to build in the first place) really taxes my system: a wall eroding with 50 blocks detaching from it every second, and up to 6,000 active bodies would be unreasonable live and is as of now only viable as a demo.

What follows is a set of interactive figures that aim at building some understanding of the constraints behind the solver I built. I will not claim any authority when it comes to collision solvers. I hope these figures can also get you to understand the design space better, and they turned out to be a useful tool for building my own intuitions and helping me tune parameters. Because these figures need to run in browser, they are 2D analogues for specific properties that matter in the final 3D solver. If you are at all familiar with physics engines and rigid body dynamics (I wasn't when I embarked on this), you may skip figures 1 through 3 and jump directly to figure 4, there won't be anything new for you there.

The first choice I had to make was penalty vs constraint. Most videogame engines go for constraint, but I wanted something a bit different, and could identify a few advantages to penalty solvers:

Energy fidelity: a penalty solver works by integrating an elastic force over the duration of a contact. The main "drawback" of these solvers, interpenetration, is really the fundamental thing that allows them to work at all: energy correctness is predicated on the force being derived from a potential, itself differentiable. As the simple disk bounce example in figure 1 shows, this is only really possible if there's a potential to integrate in the first place, and this in turn translates as interpenetration between solids in contact: some measurement derived from the depth, volume or some other derived metric of the interpenetration region is turned into an elastic force.
Performance: for a penalty solver, the work to do only depends on the individual contacts, while constraint solvers have to find a solution for all simultaneous contacts. This means penalty solvers demand the machine do the same operations as many times as there are contact pairs, while constraint solvers are more complex and can fail to converge. The price is paid in a straightforward timestep vs stability tradeoff.
Verifiability: energy conservation (without dissipation) and sink (with dissipation) as a target makes debugging easier, since it offers a simple and powerful validation gate
Contrarianism: I wanted something different, and a more "nerdy" solution that incorporated Hamiltonian dynamics instead of code-based tricks, although the final solver isn't 100% Hamiltonian, since it incorporates maximum velocity safety gates, a friction cap, and sleep logic.

The first figure shows how a generic penalty solver handles a simple bounce of a disk on a surface: the elastic treatment asks the solver to actually resolve the contact, meaning accepting some penetration.

Even the simple bounce example of figure 1 requires one to make a few choices as to how to treat it. Penetration has to be given a value: the solver able to handle this case needs to model the elastic force vs penetration (here treated as a simple normal force monotonic in the overlap area), and the equations of motion for the solid need to be approximated numerically in discrete time using integrators.

Integrators are something I had encountered in the past, working on some gravitational codes (simulating three-body problems, fitting light curves to orbit and the like back when I was working in astronomy). For the unaware, these are the numerical methods with which the equations of motion (or any differential equation) get sampled in the discrete domain. Here I chose to showcase four of them to explain some of the tradeoffs, and also how simply swapping two lines in the code can lead to completely different results, similar to how in image processing, box blur then dilation is a totally different operation from dilation then box blur.

The four integrators examined here split into three energy behaviours: systematic injection, conservation and systematic sink. Energy injection would essentially have to be fought with dissipation at every step and is out for simulating a system of real-life solids where collisions are dissipative. An energy sink is potentially useful, wrapping dissipation into the integration of the equations of motion themselves. This is essentially trading energy fidelity for stability, since unconditional stability (no energy injection for linear dissipative systems) comes with this behaviour. Constraint solvers can employ this family of integrators, however as mentioned above the road I went down was energy conservation (in aggregate), with symplectic Euler and Velocity Verlet remaining two candidates until late in the development process.

One neat trick of velocity Verlet I came across while making these figures was time symmetry. This particular property of the integrator didn't end up having much bearing on the reason I ultimately chose it, but it's an elegant one and I was drawn to the fact it mapped onto the time reversibility of the equations of motion.

Integrating the equations of motion is a sampling problem, where a value that's fundamentally continuous in the time domain gets chopped up into discrete intervals. This interactive figure is an attempt to get you to understand what the tradeoffs are, since energy correctness and solver stability depend on the time resolution at which a contact can be sampled. The upshot is that penetration depth (set by stiffness) fights available computational resources (set by substeps), here exemplified by seriously throttling the JavaScript code driving the animations.

Stable piles are generally the hard test of a collision solver, and penalty solvers in particular face specific challenges in this regard: stability requires limiting stiffness, but stiffness is necessary for the bottom of a pile to hold its top up. This is essentially a compromise, and in this figure I also present a minimal position solver to show how a constraint solver deals with piles as well.

This figure sort of jumps the gun as well, introducing friction, which I'll explain in some more detail in the next figure.

The contact model I used here, replacing SAT in an earlier prototype, was selected specifically for energy conservation. In this figure, I showcase a simplified friction model, the least frozen part of the solver at the moment (I may tweak it to avoid flat-faced solids sitting askew in equilibrium in the future), and how friction can introduce specific stability issues for box-shaped solids, which led me to select velocity Verlet as the integrator.

This is the final figure of the writeup, and should hopefully give you an idea of how the solver functions, and why I chose the specifics in this way.

On what remains to be done on this solver, it's good enough to run some specific pieces, but it deserves a few adjustments and expansions I haven't gotten to yet, such as:

some more optimization, the power draw of my GPU running the solver caps at 17W, I believe it points to some headroom being available still
a better rotation model, the current one doesn't provide physically correct precession nor the "racket flip" in freefall
the ability for friction and the elastic force to provide a righting moment, to avoid flat-faced solids settling askew on another flat surface and provide more physically accurate rolling resistance
an improved wake/sleep logic. The current one is adequate, but needs a lot of hand-tuning depending on the render I want.

Also note this writeup is a (I hope) concise explainer of the solver and of some of the fundamental principles behind it. I have intentionally skipped explaining rotation models, how the code is optimized, narrating my embarrassing early attempts, or detailing the sleep model to keep things as short and digestible as possible, and because future versions may change a few things about the solver itself.