Technical report · verified differentiable CFD in Rust
One open Rust stack spanning every solver paradigm — each checked against an analytic or reference oracle, differentiable end to end, and deployed self-verifying to WebAssembly.
The Charlot Lab, Institute for Physical AI @ BMI
CC0 1.0 Universal — dedicated to the public domain. Reference implementation: ferromotion-fluid (Rust).Deployed robots do not run computational fluid dynamics in the control loop, and the discipline's own 2026 consensus is that verification — not throughput or model size — is its weak layer. We report a fluid-dynamics stack that treats verification as the primitive: every solver ships with an oracle it must reproduce before it is admitted. The stack spans all five solver paradigms — pressure-projection (MAC), lattice-Boltzmann (D2Q9/D3Q19, CPU and GPU), smoothed-particle hydrodynamics with volume-of-fluid, Fourier-spectral Navier–Stokes, compressible finite-volume (HLL/Sod), and structured and unstructured finite volume up to a coupled Taylor–Hood incompressible Navier–Stokes solver. Each is verified against an analytic or reference solution (Ghia lid-driven cavity, Taylor–Green decay, the exact Riemann problem, the Method of Manufactured Solutions) at the tolerances tabulated in §3. The stack is differentiable by construction: gradients through the lattice, through the global pressure projection (a self-adjoint solve reused backward), through the fluid–structure loop, and through the unstructured sparse assembly all match finite differences to machine precision (§4). Two products follow from verification-first design: a surrogate honesty harness that reads the physics directly and flags a prediction which shares a competitor's mean-squared error but violates incompressibility or the energy spectrum (§5); and a coefficient-model diagnostic that shows exactly where lumped drag breaks under a gust and reconstructs the classical unsteady-force hierarchy from the resolved solver (§6). The entire stack is pure Rust, hardware-unbound (CPU floor, GPU across silicon through a portable compute layer, WebAssembly in the browser), permissively licensed, and — as far as we can establish — the only differentiable CFD in Rust and the only browser-deployable differentiable fluid stack in existence. It runs, and self-verifies, in a web page.
The premise is that no capability in fluid dynamics is enclosable, because every incumbent is bound to something — a hardware platform (OpenCL, CUDA), a license (non-commercial GPU codes), a runtime (Python/HPC), or a maintenance cliff (single-maintainer cultures; abandoned differentiable frameworks). The open architecture is the lever: pure Rust, portable across every silicon, wasm-clean to the browser, permissively licensed, verification-first, differentiable. Three facts fix the strategy.
Fact 1 — deployed robots do not run CFD in the loop, and that is the design point. Real CFD lives in design/certification and in environment precomputation. Control loops everywhere run lumped coefficients plus learned residuals; adaptation of the Neural-Fly class beats CFD-in-the-loop for gust rejection [8]. A stack premised on "robots will run CFD onboard" fights the evidence. Ours builds the seam those systems actually live on (§6).
Fact 2 — the capital is in design surrogates, not robotics fluids; and foundation-models-for-physics are serious research with open limits — the training data for a general one has not been paid for. Marketing that pairs "physical AI" with in-loop CFD surrogates is largely aspirational.
Fact 3 — the trust layer is the acknowledged weak point. Vendor surrogate tooling now ships "don't trust the model, test the physics" checks; vanilla physics-informed networks are consensus-unfit for forward turbulence; neural-operator spectral bias is a known, actively-patched failure mode. Verification is where the field admits it is thin — and therefore where an open stack should plant its flag.
From these follows the empty niche we occupy. Every 2026 browser fluid demo is eye-candy without physics accountability; every verified solver is install-only. The intersection — a solver that renders and self-verifies in the page — was empty. It is now the interactive companion to this report.
The governing equations are the incompressible Navier–Stokes system (and, for the compressible regime, the Euler conservation laws):
The discipline of the stack is that no discretization of this system is admitted until it reproduces a known answer. We use four classes of oracle, in order of strength: (i) analytic solutions — Taylor–Green decay $\mathbf{u}\,e^{-2\nu k^2 t}$, the shear-wave momentum-diffusion mode, the plane-wave advection–diffusion field, the exact Riemann solution of the shock tube; (ii) reference tabulations — the Ghia et al. (1982) lid-driven-cavity data [1]; (iii) the Method of Manufactured Solutions — choose a field, feed the analytic source it implies, and measure the error's convergence order; and (iv) conservation and structural invariants — mass, momentum, and energy to machine precision where the boundary fluxes vanish. Every solver in §3 carries at least one such oracle in its test suite; the numbers reported are those tests' outputs, not claims.
error(h) = C · hᵖ, h = 1/N (grid spacing) observed order p = log(err₁/err₂) / log(h₁/h₂) verified 2nd-order ⇒ p ≈ 2 ⇒ refine 2× (h→h/2) ⇒ error → error / 4
The stack covers all five CFD paradigms. Each row is a solver, its governing oracle, and the verified result.
| Solver | Paradigm | Oracle | Verified result |
|---|---|---|---|
| MAC projection | elliptic/parabolic | Taylor–Green; Ghia cavity [1] | order 1.99; worst dev. 0.0042 |
| LBM D2Q9 (CPU) | kinetic | Taylor–Green (diffusive scaling) | order 3.7; mass drift <1e-12 |
| LBM D2Q9 (GPU) | kinetic | CPU trajectory (ordering-exact) | 1.9e-6; ~3.4 GLUP/s |
| LBM D3Q19 | kinetic (3-D) | shear-wave diffusion; 2-D reduction | order 1.99; 2-D↔3-D 4e-16 |
| SPH (WCSPH) | particle | momentum; hydrostatic profile | drift 3e-18; $dp/dz$ to 2.5e-4 |
| VOF (minmod TVD) | interface | solid-body rotation; boundedness | vol. drift 1.4e-17; $C\in[0,1]$ |
| Spectral operators | spectral | band-limited derivative; Poisson | 5e-15 (exponential conv.) |
| Spectral Navier–Stokes | spectral | Taylor–Green decay (nonlinear ≡ 0) | 4.66e-15 (machine) |
| Euler (HLL) | hyperbolic/compressible | Sod exact Riemann [3] | star $p^*,u^*,\rho^*$ <1% |
| FVM (structured) | finite volume | plane-wave adv.–diff.; conservation | order 2.00; drift 0.0 |
| FVM (unstructured) | finite volume (triangles) | MMS on a jittered mesh | order 2.08 |
| Coupled NS (Taylor–Hood) | incompressible, coupled | MMS (Stokes + Picard NS) | order 3.86; Picard 6 iters |
The lattice-Boltzmann update is the BGK relaxation of the discrete populations $f_k$ toward their equilibrium, followed by streaming:
with $\nu=(\tau-\tfrac12)/3$. The GPU path reproduces the CPU reference's collide-then-stream ordering exactly, so the verified CPU solver — not a fresh statistical comparison — is its oracle. On the coupled incompressible problem we avoid the collocated velocity–pressure instability entirely by using the inf-sup-stable Taylor–Hood mixed element (quadratic velocity, linear pressure), which is stable by construction rather than by a Rhie–Chow patch; the steady Stokes saddle system is $\big[\begin{smallmatrix}\nu K & -B^{\!\top}\\ B & 0\end{smallmatrix}\big]$, and the full Navier–Stokes solution is recovered by a Picard iteration over the Oseen convection operator, converging in six iterations at the Reynolds number tested.
A verified solver becomes an instrument for inverse problems only if one can differentiate through it. The stack is differentiable throughout, and every gradient is checked against a central finite difference at the tolerance shown.
Kinetic — generic-scalar dual. The lattice update is a rational polynomial in the populations, so carrying a dual number through the full simulation yields exact gradients. Identifying the viscosity from an energy-decay trace recovers $\tau$ to relative error $6\times10^{-11}$.
Projection — self-adjoint reuse. A MAC step is an explicit predictor followed by a global pressure-Poisson projection. The projection $(\mathbf{u}^*)\mapsto(\mathbf{u})$ is linear, parameter-free, and self-adjoint in the staggered inner product, so its vector–Jacobian product is the same pre-factored Poisson solve run backward — one solve, no re-derivation. Only the nonlinear predictor needs a hand-transposed stencil. The result reads:
Finite-difference checks: $\partial J/\partial\mathbf{u}_0$ to $1.6\times10^{-8}$, $\partial J/\partial(\text{lid})$ to $6\times10^{-12}$, $\partial J/\partial\nu$ to $1.4\times10^{-8}$; viscosity identification from an observed field converges to floating-point zero.
Fluid–structure — the coupled sensitivity. A self-propelled undulatory swimmer earns thrust from the fluid with no prescribed forward velocity; the exact $\partial(\text{displacement})/\partial(\text{gait})$, threaded through the coupled immersed-boundary loop including the body's own free degree of freedom, matches finite differences to $5.9\times10^{-11}$, and the body learns to swim on that gradient.
Unstructured — the discrete adjoint. For $-\nabla\!\cdot(k\nabla\phi)=S$ on a triangular mesh, the objective gradient with respect to a per-element conductivity field reuses the same sparse Cholesky factor (self-adjoint), giving the entire gradient field in one extra solve; it matches finite differences to $2.6\times10^{-6}$. To our knowledge this is the only differentiable CFD implementation in Rust and the only browser-deployable one in existence; the dead differentiable frameworks it succeeds were bound to platforms it is not.
Verification-first design yields a product the field is asking for: an audit that reads the physics rather than the error metric. The claim a mean-squared error cannot see is that two predictions can share an identical error against ground truth while one is physically valid and the other a fluent lie. The harness computes receipts directly — the incompressibility residual $\lVert\nabla\!\cdot\!\mathbf{u}\rVert$ and a fast (FFT-free) spectral-roughness proxy $\lVert\nabla^2\mathbf{u}\rVert/\lVert\mathbf{u}\rVert$ that stands in for the high-frequency energy a neural operator smears — and grades a prediction against the ground truth's own receipts. On an adversarial pair constructed to share an MSE:
| Cheat mode | Shared MSE | Honest receipt | Cheat receipt | Verdict |
|---|---|---|---|---|
| injected divergence | 0.106 | div 6.9e-15 | div 0.941 | flagged |
| high-frequency noise | 0.071 | rough 84.5 | rough 729 | flagged |
The metric cannot separate the pair; the physics receipt separates it cleanly. A surrogate and its audit ship together. As a concrete surrogate to audit, a Dynamic-Mode-Decomposition operator [6] fit to spectral-NS snapshots — reduced to rank three through the method of snapshots, with the whole computation expressed through the $m\times m$ snapshot Gram matrices — reconstructs the training window to $2.2\times10^{-4}$ and extrapolates beyond it to $2.1\times10^{-3}$; it learned the dynamics, not a lookup.
Because deployed control uses coefficients, the contribution is to say exactly where they break and how the resolved solver repairs them. Drive a body through a gust with the resolved immersed-boundary solver and fit the classical unsteady-force hierarchy to the measured force:
In steady approach the quasi-steady coefficient already fits to a few percent — the coefficient model is not wrong, only incomplete. During the gust it leaves ≈94% of the force unexplained; the added-mass term drops the residual to ≈7%, and the history term halves it again. The resolved solver reconstructs the unsteady-force physics from data — which is precisely how a learned residual [8] is made, measured, and trusted. The same stack closes the loop to the robot: from a divergence-free environment wind field (verified by the harness to $1.6\times10^{-14}$), a Zermelo minimum-time planner routes a vehicle to ride tailwinds and skirt headwind pockets, beating a wind-naive straight line by a margin that grows with the wind.
The disclosed object is a design discipline and its reference implementation: every layer of CFD — solver to surrogate to syllabus — reimplemented open, verified against an analytic or reference oracle, differentiable by construction, and hardware-unbound, in Rust. Every solver paradigm is present and oracle-checked at the tolerances of §3; the stack is differentiable end to end at the tolerances of §4; the honesty harness discriminates physics violations that error metrics cannot (§5); and the coefficient-break diagnostic supplies the seam deployed systems actually use (§6). The entire stack compiles to WebAssembly and self-verifies in a browser page. No constituent method is enclosable — all are public numerical analysis — and no term of the discipline is bound to a platform, license, or runtime.