fluid-sim · validation report

The physics is real.
Here is how real.

An incompressible Navier-Stokes solver on WebGPU compute shaders, measured against closed-form solutions rather than against how it looks. This page reports what the solver gets right, what it gets wrong, and by how much — including a turbulence result that does not match theory, and why.

node v24.18.013th Gen Intel(R) Core(TM) i7-1355U gpu intel gen-12lp generated 2026-08-28T20:47:17

Open the live demo → Source

5.79e-4numerical viscosityMacCormack, 2.9% of the physical ν

0.25 / 0.26stability boundarymeasured; theory says 0.25

-1.52inverse cascadetheory (Kraichnan) says −1.67

-4.53enstrophy cascadetheory (Kraichnan) says −3

3.31e-6gpu vs cpurelative L2 after 10 full steps

1.54e-13max |∇·u|incompressibility, CPU reference

26.4×gpu speedupat 512², equal solution quality

36×scheme differencesemi-Lagrangian vs MacCormack dissipation

Every figure on this page was produced by a run and read out of results/*.json. Nothing is typed in; where a value is missing, the page says so.

the headline measurement

Taylor-Green vortex: how much viscosity did the solver invent?

The 2D Taylor-Green vortex is an exact solution of Navier-Stokes. Its kinetic energy decays as e^(−4νt), precisely. A simulation decays faster, because the advection scheme's interpolation removes energy on top of the physical viscosity.

Fitting the measured decay to e^(−4ν_eff·t) gives an effective viscosity, and ν_num = ν_eff − ν is the artificial part — the same units as the real thing, so "this scheme is dissipative" becomes a number instead of an adjective. The r² column matters: it confirms the decay really is exponential, without which ν_eff would be meaningless.

Both schemes on the same grid (64²), ν = 0.02, dt = 0.005. MacCormack's error is roughly 34× smaller at the end of the run.

The gap between each curve and the dashed exact solution is energy the scheme removed that physics did not. Fitting that gap gives the numerical viscosity quoted below.

scheme	measured rate	exact 4ν	ν effective	ν numerical	as % of ν	fit r²
semi-Lagrangian	0.1638	0.0800	0.0410	0.0210	104.8%	0.99997
MacCormack	0.0823	0.0800	0.0206	5.79e-4	2.9%	1.00000

64² grid, ν = 0.02, dt = 0.005, integrated to t = 1. Semi-Lagrangian more than doubles the effective viscosity at this resolution: the fluid it simulates is roughly twice as viscous as the one it was asked to simulate. MacCormack brings the invented viscosity down to a few percent.

isolating advection

Where the dissipation actually comes from

Taylor-Green cannot answer this on its own. Its nonlinear term is exactly cancelled by the pressure gradient, so a badly diffusive advection scheme can still score well there. This test has nothing else in it: a Gaussian blob in a uniform flow, whose exact solution is the same blob, translated.

The exact solution's peak never changes: the blob is simply translated. Everything below the initial value is numerical diffusion, with no pressure, viscosity or nonlinearity able to take the blame.

scheme	peak retained	D numerical	relative L2	CFL
semi-Lagrangian	71.7%	0	0	0.326
MacCormack	98.3%	0	0	0.326

Semi-Lagrangian keeps interpolating bilinearly, and bilinear interpolation of anything but a linear field is a weighted average — a low-pass filter applied once per step, which is indistinguishable from adding a diffusion term. That is the mechanism behind both tables on this page.

convergence

First order in time, and why the obvious study got it backwards

The first convergence study run here held dt fixed and refined the grid, expecting the error to fall. It rose: 8.9e-5 → 2.3e-4 → 2.9e-4 for N = 16, 32, 64.

Sweeping dt and N independently showed why. At fixed dt the temporal error dominates and is identical at every resolution, so refining the grid cannot reduce it — and at N = 16 the spatial error happened to have the opposite sign and partially cancelled it, making the coarsest grid look the most accurate. Refining removed the cancellation, and the error went up.

The scheme is first order in time because the step is Lie-split (advect, then diffuse, then project), and Lie splitting is O(dt). That caps the whole method at first order however accurate the individual operators are — the Laplacian and the discrete curl are both second order, verified separately. It is a real property of Stam-style solvers and is rarely stated.

Each line falls with slope ≈ 1 (first order in time) until it flattens onto that grid's spatial error floor. The floor is what falls with resolution; the sloped part is identical at every N because it is the splitting error.

grid	dt=0.02	dt=0.01	dt=0.005	dt=0.0025	order in dt
16²	8.82e-4	2.82e-4	8.92e-5	1.85e-4	0.84
32²	0	5.36e-4	2.30e-4	7.95e-5	1.29
64²	0	6.02e-4	2.92e-4	1.36e-4	1.07

Taylor-Green, ν = 0.05, integrated to t = 0.25, maccormack advection. Reading down a column shows why a fixed-dt refinement study is uninformative here; reading along a row gives the temporal order.

stability

Where it breaks, and a correction to the usual framing

The brief asked for "the CFL condition boundary". Semi-Lagrangian advection does not have one: tracing backward and interpolating can never produce a value outside the range it sampled, so it cannot amplify anything. That unconditional stability is the entire reason Stam's method exists, and reporting a CFL blow-up limit would have meant inventing a result.

So what follows is what genuinely bounds this solver. Advective CFL was swept to 2.53 with no blow-up at all — large timesteps cost accuracy, not stability. The one real boundary is explicit diffusion, and it is reproduced exactly.

The grid-scale mode is amplified by |1 − 8d| per step, so the boundary sits exactly at d = 0.25 where amplification reaches 1. Measured and predicted agree to better than 1e-6 at every point.

The most common way a fluid simulation is quietly wrong: the flow still looks plausible while mass is not conserved. Roughly 50 iterations are needed for machine precision at this resolution.

measured boundary

d = 0.25 stable, 0.26 unstable

Theory puts it at exactly 0.25. The operator was tested in isolation with the checkerboard mode seeded deliberately.

inside the full solver

d = 0.26 / 0.28

Below the bound the projection annihilates the mode entirely (growth ~1e-14); above it, the mode grows despite the projection.

advective CFL

no blow-up to CFL 2.5

Unconditionally stable, as designed. The error stays around 0.04 throughout.

Vorticity confinement is an energy source

Confinement is in the brief as a way to restore detail lost to numerical dissipation, and it does that. But it is a fabricated force, and this table is why every measurement on this page has it switched off — the solver throws rather than allowing it in validation mode. The exact solution loses energy over this interval; from ε = 1 the simulation gains it.

confinement	energy ratio	exact ratio	verdict
ε = 0	0.626	0.852	decays
ε = 0.5	0.863	0.852	decays
ε = 1	1.114	0.852	adds energy
ε = 2	1.562	0.852	adds energy
ε = 5	1.720	0.852	adds energy
ε = 10	2.384	0.852	adds energy
ε = 20	2.405	0.852	adds energy
ε = 40	38.788	0.852	adds energy
ε = 80	4062.159	0.852	adds energy
ε = 160	23593.437	0.852	adds energy
ε = 320	diverged	0.852	adds energy

turbulence

The inverse cascade matches. The enstrophy range does not.

Kolmogorov's k^(−5/3) describes three-dimensional turbulence. In 2D there is no vortex stretching, vorticity is materially conserved, and the phenomenology is different: Kraichnan and Batchelor predict a dual cascade — energy travelling upscale from the forcing with slope −5/3, and enstrophy travelling downscale with slope −3. Comparing a 2D spectrum to −5/3 across all scales would be comparing against the wrong law for most of it, so both ranges are measured separately and each against its own prediction.

Time-averaged over 270 samples. The shaded bands are the fit windows, fixed from k_f before the spectrum was looked at — choosing them afterwards is how almost any curve is made to agree with a power law. Parseval holds to 4.18e-6, so the normalisation is right.

inverse energy cascade, k ∈ [4, 24]

-1.52 vs theory −1.67

r² = 0.986. Within 9% of Kraichnan's prediction — energy is genuinely being transported upscale from the forcing.

enstrophy cascade, k ∈ [64, 153]

-4.53 vs theory −3.00

r² = 0.999 — an extremely clean power law, but roughly 1.5× too steep. The spectrum is a straight line in log-log; it is simply the wrong line.

the measurement is sound

Parseval 4.18e-6

Σ E(k) matches the kinetic energy computed independently in physical space, so the normalisation is right; and max |∇·u| = 0 with multigrid, so the field really is solenoidal. Neither can explain a wrong slope.

Testing the diagnosis

Attributing the missing cascade to numerical dissipation was a hypothesis. It is testable: port the less dissipative scheme to the GPU and re-run the identical study. That is what this is.

Identical forcing, resolution, projection and fit window; only the advection scheme differs. Whatever separates the two curves is the scheme's numerical dissipation and nothing else.

semi-Lagrangian

-6.08

r² = 1.000. First order; measured numerical viscosity ≈ 100% of ν on the Taylor-Green test.

MacCormack

-4.53

r² = 0.999. Second order, limited; ~2.9% of ν. GPU energy retention over 200 inviscid steps: 0.2475 of 0.25, against 0.2117 for semi-Lagrangian.

theory

−3.00

Changing the scheme moved the slope by 1.55, toward theory.

What is left over

what matched: The inverse energy cascade, at -1.52 against −1.67 with r² = 0.986. Energy really is being carried upscale from the forcing, which is the distinctively two-dimensional half of the theory.
what did not: The enstrophy cascade, at -4.53 against −3. The fit is excellent (r² = 0.999) so this is a real power law at the wrong exponent, not scatter.
most likely cause: Residual numerical dissipation at the small scales. The enstrophy range sits nearest the grid, where any remaining dissipation bites hardest, and this run used f32 with a 3-cycle projection. The scheme comparison above shows the slope moving with the advection scheme, so dissipation is demonstrably part of it.
ruled out: The projection (max |∇·u| = 0), the spectrum's normalisation (Parseval to 4.18e-6), and the fit window, which was fixed from k_f before any data was seen.

gpu

Moving the solve onto compute shaders

Every stage of the timestep — advection, diffusion, divergence, the iterative pressure solve, the gradient subtraction — is a WGSL compute dispatch. The CPU only issues commands. What follows is how that was checked and what it bought.

Multigrid's unit is a V-cycle, which costs about 4/3 of a fine-grid sweep across all levels; plotting against raw iteration count would have compared unlike things. It plateaus near 1e-7 because that is the f32 floor, not because it has stopped converging.

At roughly 40 sweeps of equivalent cost: Jacobi reaches 0, red-black 0, multigrid 9.34e-8. Jacobi and Gauss-Seidel are local — one sweep moves information one cell — so their reduction factor for a mode spanning L cells is about 1 − O(1/L²), and refining the grid makes that strictly worse. Multigrid solves those modes on a grid where they are no longer smooth.

Verification against the f64 CPU reference

check	result	measured
advection kernel matches CPU	pass	`relL2 u=2.29e-7 v=2.27e-7 (f32 floor ~1e-6)`
divergence kernel matches CPU	pass	`relL2=5.78e-8`
GPU projection substantially reduces divergence	pass	`red-black: 1.81e+1 -> 1.74e-2 (1040x)`
red-black clears divergence better than Jacobi at equal cost	pass	`after 200 sweeps: RB 1.74e-2 vs J 6.19e-1 (35.6x better)`
red-black Gauss-Seidel converges faster than Jacobi	pass	`at 320 sweeps: RB 2.51e-4 vs J 3.14e-2`
Jacobi residual decreases monotonically	pass
GPU matches CPU after 10 full steps	pass	`relL2 u=3.31e-6 v=3.31e-6`
GPU keeps the field divergence-free	pass	`max\|div\|=1.76e-5 (CPU 1.78e-15)`
GPU MacCormack matches CPU after 10 steps	pass	`relL2 u=3.51e-7 v=3.50e-7`
GPU MacCormack retains more energy than semi-Lagrangian	pass	`after 200 steps at nu=0 (exact conserves 0.25): SL 0.2117, MC 0.2475`
GPU solid mask is finite everywhere	pass	`0 non-finite cells of 16384 (was 162981/262144 before the tanh clamp)`
GPU solid mask covers a plausible area	pass	`1994 cells carry impulse`
a pinned disk slows the fluid	pass	`net fluid impulse -1.14e+0 (negative expected)`
GPU solid coupling matches the CPU field	pass	`relL2 u=2.05e-6 v=3.25e-5`
GPU and CPU agree on the force the body feels	pass	`CPU 1.145e+2 vs GPU 1.145e+2 (0.0%)`
a free disk is carried downstream by the flow	pass	`vx = 0.2954, moved 0.1617`

WGSL is f32 and permits relaxed precision on transcendentals — a plain sin measured 6.9e-5 maximum absolute error on this adapter. So the GPU cannot agree with the f64 reference below about 1e-5 relative, and any tighter claim would be measuring nothing. Agreement of 3.31e-6 after ten full steps is f32 accumulation noise.

Identical Jacobi sweeps on identical grids, so the only difference is where the arithmetic runs. The GPU is slower below N ≈ 100, where dispatch overhead dominates the work.

grid	cells	CPU ms	GPU ms	speedup
32²	1,024	12.6	9.0	1.4×
64²	4,096	48.5	32.4	1.5×
128²	16,384	178.0	12.6	14.1×
256²	65,536	725.7	12.3	59.0×
512²	262,144	2244.3	17.8	126.1×

Equal quality, not equal work

The comparison above gives both devices the same iteration count. That is fair as a hardware measurement and misleading as a practical one, because the CPU's best option is an exact FFT Poisson solve. So here each GPU solver is given as many iterations as it needs to reach max |∇·u| ≤ 0.001, and that is what is timed. Red-black could not reach the target at all at the larger sizes; multigrid reaches it in a single V-cycle.

grid	CPU FFT ms	red-black	multigrid	MG speedup
64²	44.6	23.6 / 160sw	10.9 / 1vc	4.1×
128²	45.9	95.6 / 640sw	16.6 / 1vc	2.8×
256²	178.2	263.1 / 1280sw missed	16.0 / 1vc	11.1×
512²	754.8	556.5 / 1280sw missed	28.6 / 1vc	26.4×

This is the honest version of the GPU claim. At equal iteration count the GPU looks 126× faster, but much of that is the GPU running a worse algorithm quickly. Once accuracy is held fixed, the win comes from multigrid, and it is 26.4× at 512².

interactive ceiling

256² on GPU

Largest grid holding 60 fps (16.7 ms/step), against 32² on the CPU — 64× more cells.

hardware

intel gen-12lp

An integrated GPU sharing memory bandwidth with the CPU. A discrete card would widen every ratio here; these numbers are the modest end.

fluid-structure interaction

An obstacle that is pushed back

A rigid disk immersed in the flow, coupled both ways: it blocks the fluid, and the reaction moves it. “Two-way” is the load-bearing word — a static obstacle that deflects flow without ever moving is a boundary condition, and a much easier thing to get right.

The method is direct forcing (Mohd-Yusof 1997): the solid is a smoothed mask on the existing grid, the fluid inside it is driven toward the body velocity each step, and the reaction on the body is the negative of that same integral. That last clause is where Newton’s third law lives, and it is the one thing here that can be checked exactly rather than approximately.

1.15e-14momentum driftone coupling step, default path

9.65e-8residual slip4 forcing passes

1.2min density ratiobelow this the explicit coupling diverges

-0.420spin in sheartorque-free theory says -0.5

One pass moves the fluid a fraction chi of the way to the body velocity, leaving (1 - chi) behind. Repeating drives the residual down by that same factor each time, which is a straight line on a log axis -- a property of the method rather than a tuned number.

Points appear only where the run stayed finite. Without the correction every ratio below 1.2 diverged; with it the sweep stayed stable down to 0.1, the lowest tried.

Newton’s third law, to machine precision

The coupling moves momentum between fluid and solid and must not create or destroy any. Applied correctly that is exact, so it is measured as such rather than given a tolerance. A coupling that leaked momentum would still produce a disk that moved plausibly, which is precisely why it needs a check that “looks about right” cannot satisfy.

configuration	relative momentum drift	verdict
default (exact)	1.15e-14	conserves
added-mass corrected	0	trades conservation for stability

The correction is not free, so it is not the default

the problem: Explicit coupling computes the force from the current fluid state, applies it to the body, and lets the body change the fluid next step. In 2D a disk’s added mass is exactly the mass of fluid it displaces, so at a density ratio of 1 the fluid’s inertia already equals the body’s and that feedback loop diverges.
the usual fix: Fold the added mass into the effective inertia. It works: the stable range extends from 1.2 down to 0.1.
the cost: The added mass is fluid mass, and the fluid’s momentum is already tracked in the fluid. Counting it again on the body makes the body under-respond, and total momentum stops being conserved: 0 against 1.15e-14.
the choice: Off by default, like vorticity confinement elsewhere in this project. The convenient option exists, is labelled, and is not what any number here is measured with.

The GPU port computes the same thing

The CPU implementation is the reference; the GPU one exists so the obstacle can appear in the interactive demo. Given the same body in the same flow, the two agree on the force to better than 0.01%, which is the check that matters — the force is the whole coupling.

A bug the CPU version structurally could not have

symptom: The GPU disk produced NaN velocity within a few steps. 162,981 of 262,144 impulse cells came back non-finite — and every one of them was far from the body, the region where the mask should have been a clean zero.
cause: The mask is ½(1 − tanh((r − R)/w)), and with a mask width of 1.5 cells that argument reaches about 240 at the far side of the domain. WGSL does not require tanh to be robust for large arguments, and this adapter evidently computes it as (e¹ − e⁻¹)/(e¹ + e⁻¹): both terms overflow to infinity and the result is Inf/Inf = NaN.
why the CPU never showed it: Math.tanh is robust. The reference implementation was correct and silent, which is exactly the case a port has to be tested for rather than trusted through.
fix: Clamp the argument to ±10 before tanh, which saturates well inside f32 and changes no representable value. Pinned by a check that asserts the whole impulse field is finite.

Rotation has a separate threshold

A shear flow u = y − π has vorticity −1, and the classical result for a torque-free body is rotation at half the vorticity, ω → −0.5. The coupling reproduces the sign and the approach — but only once the moment of inertia is large enough. Below that it diverges for the same reason translation does, at its own threshold.

density	moment of inertia	omega after the run	theory
1.5	0.1473	diverged	-0.5
3.0	0.2945	diverged	-0.5
10.0	0.9817	-0.4202	-0.5
50.0	4.9087	-0.2441	-0.5
200.0	19.6350	-0.0900	-0.5
1000.0	98.1748	-0.0204	-0.5

Heavier bodies approach −0.5 more slowly, which is the expected response to a fixed torque, and none of these runs is long enough to arrive. What is asserted is the sign and the direction of approach, not the endpoint.

honestly

Bugs found, and how

Four, all fixed. Three of them were in the measuring code rather than the solver — a validator reported the scheme stable well past its real limit, and a benchmark reported a speedup that was partly an algorithm difference. Instruments need testing at least as much as the thing they measure.

1. Every GPU pass read the same uniform values

symptom: The pressure solve did nothing. Divergence fell from 18.1 to 16.3 after 200 Jacobi sweeps, and red-black froze at a residual of 6.13e-1 from its first sweep onward, never improving.
cause: queue.writeBuffer is ordered against submitted command buffers, not against commands being encoded. Encoding a pass does not capture the uniform's contents — the pass reads the buffer when it executes, after every write issued before the submit has landed. Writing params, encoding a pass, writing again, encoding again, then submitting once meant all passes saw the final value: relaxation factor 0 (so Jacobi was p + 0·(new − p), i.e. nothing) and colour always 1 (so red cells were never updated).
fix: A ring of 256-byte-aligned uniform slots, with the bind helper allocating a pass's slot as it binds it — so the parameters cannot be separated from the pass that reads them.
found by: Suspected from the shape of the failure: the standalone residual curve worked while the in-step projection did not, and the only difference was one encoder versus many.
file: src/gpu/solver.js

2. The stability sweep reported the solver stable at 2.4× the theoretical limit

symptom: Explicit diffusion was reported stable up to ν·dt/h² = 0.6. The textbook bound is 0.25.
cause: Two compounding mistakes in the instrument, none in the solver. The sweep started from a smooth Taylor-Green field, which contains essentially none of the checkerboard mode that actually goes unstable — so there was nothing to amplify. And it watched for kinetic energy to exceed a fixed threshold, which a viscously decaying field was never going to cross however unstable the scheme was.
fix: Seed the grid-scale mode explicitly and measure its growth by projecting onto it. The operator now reproduces |1 − 8d| to better than 1e-6 at every point, with the boundary at exactly 0.25 / 0.26.
found by: Testing the diffusion operator in isolation, which immediately showed amplification of exactly 1.00 at d = 0.25 — correct all along.
file: validate/studies.js

3. A benchmark that flattered the GPU twice over

symptom: An early run reported the GPU 38.9× faster at 512².
cause: Two separate unfairnesses. The 'same algorithm' comparison ran conjugate gradient on the CPU against Jacobi on the GPU and attributed the difference to hardware — CG is far stronger per iteration. And the 'best available' comparison timed 20 GPU sweeps (leaving divergence at 2.5e-1) against an exact CPU FFT solve (1e-15), letting the GPU win by doing less work and returning a worse answer.
fix: Jacobi on both sides for the hardware comparison, and an equal-quality comparison that raises the GPU's iteration count until it matches a target divergence. The honest numbers are lower and more interesting: the GPU initially lost at equal quality, which is what motivated multigrid.
found by: Noticing that the reported GPU divergence was three orders of magnitude worse than the CPU's in the same table.
file: tools/bench.html

4. Non-monotonic GPU timings from insufficient warm-up

symptom: N = 32 measured 11.9 ms and N = 64 measured 31.7 ms, while N = 128 ran in 9.5 ms — a smaller grid taking three times longer.
cause: Three unsynced warm-up steps left pipeline compilation and the iGPU's idle clock inside the measurement window, and only the first two configurations paid it.
fix: Twelve warm-up steps with periodic queue synchronisation. The artifact disappeared and the curve became monotonic.
found by: The physics made no sense: dispatch overhead cannot make a 4× smaller grid slower.
file: tools/bench.html

limits

What this does not do

The GPU path has only semi-Lagrangian advection. MacCormack exists on the CPU and is roughly 36× less dissipative, but was not ported. This is the single biggest reason the energy spectrum misses Kraichnan's −3, and it is a gap in the implementation rather than a property of the physics.
First order in time. Lie splitting caps the whole scheme at O(dt) no matter how accurate the operators are. Strang splitting would give second order for about 1.5× the cost.
No inverse cascade was observed. The measured slope in the inverse range is positive. Whether a less dissipative scheme would recover Kraichnan's −5/3 at this resolution is untested, and claiming it would be speculation.
f32 on the GPU. The multigrid residual plateaus near 1e-7 because that is the precision floor, not because the method stalls.
Periodic boundaries only. The FFT reference solver and the spectra both require it. Walls would need the conjugate-gradient path, which exists but is not exercised here.
Integrated GPU. Every performance ratio on this page comes from an Intel iGPU sharing bandwidth with the CPU, which is the pessimistic end of the range.
The fluid-structure coupling is explicit. It diverges below a density ratio of about 1 without the added-mass correction, and the correction costs exact momentum conservation. A strongly-coupled (iterated) scheme would fix both and is not implemented.
The GPU coupling lags one step. The force on the body is the integral of the impulse field, which needs a readback, and awaiting it inside the step would stall the pipeline. So the demo applies the previous step's force. That extra lag makes the explicit coupling slightly less stable than the CPU path, which is why the demo uses a heavy disk. Every fluid-structure number on this page comes from the CPU path, which has no such lag.

reproduce

Run it yourself

# no dependencies beyond node
git clone https://github.com/abho7/fluid-sim && cd fluid-sim

# the whole suite: solver, FFT, analytic solutions, and the
# mutation tests that prove the validators actually fire
node --test test/*.test.js

# regenerate every CPU-side number on this page
node validate/run.js

# serve the demo and the GPU harnesses (WebGPU needs a secure context)
node tools/serve.js
#   /                   the interactive demo
#   /tools/gputest.html GPU vs CPU verification
#   /tools/bench.html   performance
#   /tools/spectrum.html forced turbulence

# rebuild this page from results/*.json
node tools/build_report.js

The solver, the FFT, the multigrid cycle, the plots on this page and the static server are all written from scratch. Node and a WebGPU-capable browser are the only requirements.