I drew OpenAI's math chart with a 2D Navier-Stokes solver. Here's the script.

Set one constant in this Python script and the smoke stops at that height. It's a real 2D Navier-Stokes solver in about 60 lines. A bigger version drew this:

That's OpenAI's chart of selected open math problems. Silver is GPT-6 Astra at its highest plotted compute, about 17%. Blue is OpenAI's internal model at its lowest, about 26%. Blue used about 12× less compute, so it gets 1 thread of smoke against 12. They're smoke threads, not CPU threads.

The chart came out the same day as the OpenAI Navier-Stokes proof. So I drew it with Navier-Stokes. I had Claude write the solver.

The script

# A 2D Boussinesq plume, pseudo-spectral, rising into a stable layer. Writes plume.png. import numpy as np from PIL import Image from scipy.fft import fftfreq, irfft2, rfft2, rfftfreq from scipy.ndimage import gaussian_filter CEILING = 110 # where the stable layer starts, in grid cells above the floor. Try 60 to 170. N2 = 0.05 # how stable that layer is. Set it to 0 and there's no ceiling. NX, NY, FLOOR, NOZZLE, DT, STEPS = 200, 250, 30, 36, 0.125, 6000 NU, KAPPA, DRAG, SOURCE = 0.05, 0.025, 0.004, 0.022 rng = np.random.default_rng(17) y, x = np.arange(NY, dtype=float)[:, None], np.arange(NX, dtype=float)[None, :] kx, ky = 2 * np.pi * rfftfreq(NX)[None, :], 2 * np.pi * fftfreq(NY)[:, None] k2 = kx**2 + ky**2 k2_inv = np.where(k2 > 0, 1 / np.where(k2 > 0, k2, 1), 0) dealias = np.exp(-36 * (k2 / k2.max()) ** 4) ifft = lambda f: irfft2(f, s=(NY, NX)) ramp = lambda t: np.clip(t, 0, 1) ** 2 * (3 - 2 * np.clip(t, 0, 1)) stable = N2 * ramp((y - FLOOR - CEILING) / 8) # stable air, not a wall sponge = 0.12 * (ramp((FLOOR - 22 - y) / 18) + ramp((y - NY + 28) / 20)) # soaks up waves, top and bottom nozzle = np.exp(-(((x - NX / 2) / (NOZZLE / 2)) ** 4) - ((y - FLOOR) / 5) ** 2) def velocity(wh): # vorticity to velocity, through the streamfunction  psi = wh * k2_inv return ifft(1j * ky * psi), ifft(-1j * kx * psi) def rhs(wh, bh): # vorticity and buoyancy; the nonlinear products are taken on the grid  u, v = velocity(wh) wx, wy, bx, by = ifft(1j * kx * wh), ifft(1j * ky * wh), ifft(1j * kx * bh), ifft(1j * ky * bh) return ( rfft2(-(u * wx + v * wy)) + 1j * kx * bh - (NU * k2 + DRAG) * wh, rfft2(-(u * bx + v * by) - stable * v) - KAPPA * k2 * bh, ) def at(f, px, py): # bilinear lookup at particle positions  i, j, a, b = np.floor(px).astype(int), np.floor(py).astype(int), px % 1, py % 1 i0, i1, j0, j1 = i % NX, (i + 1) % NX, j % NY, (j + 1) % NY return (f[j0, i0] * (1 - a) + f[j0, i1] * a) * (1 - b) + (f[j1, i0] * (1 - a) + f[j1, i1] * a) * b wh = np.zeros((NY, NX // 2 + 1), dtype=complex) bh, px, py = wh.copy(), np.empty(0), np.empty(0) for t in range(STEPS): k1w, k1b = rhs(wh, bh) k2w, k2b = rhs(wh + DT * k1w, bh + DT * k1b) wh, bh = (wh + DT / 2 * (k1w + k2w)) * dealias, (bh + DT / 2 * (k1b + k2b)) * dealias if not np.isfinite(wh).all(): raise SystemExit(f"Solver blew up at step {t}. That's the numerics, not the equations. Lower DT or SOURCE.") w, b = ifft(wh), ifft(bh) b += DT * 0.5 * nozzle * (SOURCE * (1 + 0.25 * rng.standard_normal()) - b) # the flicker seeds turbulence  w += DT * 0.002 * nozzle * rng.standard_normal((NY, NX)) wh, bh = rfft2(w * (1 - DT * sponge)), rfft2(b * (1 - DT * sponge)) u, v = velocity(wh) px, py = np.r_[px, NX / 2 + (rng.random(6) - 0.5) * NOZZLE], np.r_[py, FLOOR + rng.random(6) * 3] px, py = px + DT * at(u, px, py), py + DT * at(v, px, py) px, py = px[py > FLOOR - 15][-24000:], py[py > FLOOR - 15][-24000:] # the oldest smoke drops out  smoke = np.histogram2d(py, px % NX, bins=(NY * 3, NX * 3), range=((0, NY), (0, NX)))[0] pixels = (255 * (s := gaussian_filter(smoke, 1.2) ** 0.6) / s.max()).astype("uint8") pixels[3 * (FLOOR + CEILING), ::8] = 200 # a dotted line where the stable layer starts Image.fromarray(pixels[::-1]).save("plume.png") print(f"Stable layer starts at {CEILING}. 99% of the smoke is below {np.percentile(py, 99) - FLOOR:.0f}.") 

Save it as plume.py, then:

uv run --no-project --with numpy --with scipy --with pillow python plume.py 

Run this and you get plume.png in under a minute on an M1 Pro, and this line:

Stable layer starts at 110. 99% of the smoke is below 116. 

The image shows a column of smoke rising from the floor. At the dotted line it spreads sideways in curls. The seed is fixed, but another machine may differ slightly.

Four edits to try

The last row is the guard working. The script checks for non-finite values every step, so it stops instead of writing a garbage image.

What the code solves

These are the equations in rhs, where ω is vorticity, b is buoyancy, and μ is DRAG:

∂ω/∂t + u·∇ω = ∂b/∂x + ν∇²ω − μω ∂b/∂t + u·∇b = −N²(y) v + κ∇²b ∇²ψ = −(ω − ⟨ω⟩), u = ∂ψ/∂y, v = −∂ψ/∂x 

That's 2D incompressible Navier-Stokes in vorticity form, with Boussinesq buoyancy. Small density differences drive the flow, and density is otherwise treated as constant.

• velocity turns vorticity into velocity through the streamfunction. Division by k² is one line in Fourier space.
• rhs is pseudo-spectral. Derivatives happen in Fourier space. Nonlinear products happen on the grid.
• dealias is an exponential filter. It damps the highest frequencies. That limits aliasing, but doesn't fully remove it.
• The time step is Heun's method, so rhs runs twice per step.
• stable is the ceiling. N² is zero below the line and ramps up over 8 cells. Rising fluid loses buoyancy there. It's stable air and not a wall, so a little overshoots.
• sponge damps vorticity and buoyancy near the top and bottom. The FFT makes the domain periodic, so otherwise waves would wrap around.
• nozzle feeds in hot fluid with a random flicker. The flicker seeds the turbulence.
• The particles are forward Euler with a bilinear lookup. They're the smoke. They don't push back on the flow.

What the full version adds

The solver behind the GIF has the same equations on a 540×675 grid, plus an added upward wind on the single thread. Four hand-set choices shape what you see.

The heights. Each plume gets a stable layer at its number. A strong enough stable layer can hold back rising air, and this one is set strong enough. The renderer also fades anything drawn above that line.

The threads. 12 against 1 is the compute gap, drawn. Each thread gets the same number of particles.

The wind. The single thread rises in a slow upward co-flow, like a smoke wire in a wind tunnel. Without it the thread breaks into a cloud. That wind is added by hand, and it isn't strictly incompressible. It's a choice made for the picture. The 12 threads rise on buoyancy alone.

The tidying. It has drag and sponge layers too. Near its ceiling, the single thread's particles also age out fast. That's why it ends in a curl.

So the numbers come from the chart. The simulation doesn't produce them, and the picture isn't evidence. Change one constant and the smoke stops somewhere else. That's why the picture can't check the chart.

A solver blowing up proves nothing about Navier-Stokes

The full solver diverged three times while it was being tuned. That's numerical instability, not mathematical blowup.

A timestep past the method's CFL limit can do it. So can aliasing or too little resolution. Non-finite numbers alone don't tell you which. The official problem description warns about interpreting apparent numerical blowup in Euler.

2D can't show the claimed blowup either, because that claim is in 3D. The same description says classical 2D flow, with smooth data and positive viscosity, stays regular. This solver adds buoyancy, so that statement doesn't automatically cover it. Either way, a 2D picture can't demonstrate a 3D singularity.

The chart readings

OpenAI published the chart on September 8. @chetaslua shared an annotated version on November 3. The readings here are checked against OpenAI's original. They're pixel reads, not benchmark data.

The sum, so you can check it. In the 1199×858 image, the two markers sit near x=379 and x=841. One decade on the log axis is about 430 pixels.

compute gap 10^((841 − 379) / 430) ≈ 11.9, so about 12× pass rate gap 0.263 − 0.171 ≈ 0.09, so about 9 percentage points relative rise 0.263 / 0.171 − 1 ≈ 54% 

The axis has no numbers on it, so the 12× assumes the repeating ticks mark decades. The chart, the announcement and the proof repository give no problem count, problem list or pass definition. Lowest and highest are the ends of each plotted line, and the chart shows no presets. It's one evaluation. No price, latency or hardware claim follows.

The proof itself

OpenAI has published a claimed blowup proof for 3D Navier-Stokes with smooth external forcing, with a Lean formalization. Independent formal checks and mathematical analysis have followed. So has criticism of how the written argument matches the Lean. Clay's notice says the problem has "apparently been settled", while its evaluation continues.

OpenAI reports about 10,000 concurrent agents and a result after 88 hours. That's its own account, and it's separate from the chart's compute sweep. In a third-party transcript, Sam Altman says "the model that is going to be released soon, the next GPT". That transcript isn't audio-checked, so I call it the internal model.

The maintained version, with the full source list, is on argbe.tech.