Writing · Flexcompute Foil
What I Learned
Porting XFOIL to Rust
Living inside Drela’s code for two months, rebuilding XFOIL in Rust and WebAssembly, and then forcing it to justify itself directly against the original.
Harry Smith
I wanted a browser-based airfoil analysis tool with real-time visualisation — adaptive streamlines, surface pressure vectors, smoke particle tracing, GPU-accelerated rendering, the whole thing. But the only honest way to get there was to port the real solver. Not a simplified panel method. Not an “inspired by” clone. The actual XFOIL logic, compiled to WASM, running in a browser tab.
The result is Flexcompute Foil. Try it at foil.flexcompute.com.
Try the Applet
NACA 2412 · α=4° · Re=106 · viscous
What the Validation Actually Looks Like
Not screenshots, not marketing art — the actual cached Plotly comparisons from the parity sweep against XFOIL. The full dashboard has everything in one place.
A caveat on the non-NACA foils: the XFOIL side of the E387 and DAE11 comparisons is incomplete. These sweeps were run using my instrumented fork of XFOIL, which writes detailed JSON debug dumps at every internal checkpoint. That instrumentation adds significant memory and I/O overhead, and for harder foils that require many Newton iterations per operating point, the instrumented binary struggles to complete a full −15° to +25° ASEQ sweep. The convergence failures you see on those foils may be an artefact of the instrumentation rather than a real solver limitation. Stock XFOIL without the debug hooks would almost certainly produce fuller polars. I haven't yet had time to run the comparison through an uninstrumented build. If anyone wants to finish that particular piece, I'd genuinely welcome the help — the sweep scripts, the dataset, and the Plotly embed generator are all in the repo.
Lift Polars Across Three Foils
CL vs α for NACA 0012, 2412 and 4412 from −15° to +25° in 0.5° steps at Re = 10⁶. This is the quickest way to see that the broad aerodynamic behaviour is lining up across symmetric and cambered sections. It reaches well into and beyond stall, so you can see both the clean linear region and where the nonlinear behaviour starts to separate the curves.
Cp: NACA 0012
The surface pressure distribution is where you find out whether the solver is reproducing the actual flow physics rather than just landing on roughly the right integrated coefficients. On a symmetric section the Cp must be antisymmetric at zero alpha and the suction peak must track correctly as incidence increases.
Cp: NACA 2412
Mild camber is a good stress test because small changes in loading and stagnation behaviour show up clearly in Cp. The front-loaded pressure distribution on the upper surface and the aft recovery on the lower surface are both sensitive to how the solver handles the leading-edge velocity peak.
Cp: NACA 4412
This more heavily cambered section is where you start to see whether the agreement still holds once the pressure recovery gets less forgiving. The larger adverse gradient on the upper surface pushes the boundary layer harder and the trailing-edge loading differences become more visible.
Error Envelopes
I do not trust eyeballing overlays on its own. This view tracks the maximum Cp and γ discrepancies across the full cached validation sweep, so you can see directly how the worst-case errors evolve with angle of attack for each foil family rather than cherry-picking a single good-looking comparison.
The validation looks clean.
Here is what it took to get there.
The Safeguard That Froze the Solver
The Safeguard That Froze the Solver
Early on, RustFoil’s drag was 50% too high and lift errors ranged from −8% to +14%. It took days. The root cause was an edge velocity update being clamped by a safety limit we had added: a ±5% cap on how much the edge velocity could change per Newton iteration, relative to the inviscid baseline.
The problem was that the clamp was computed relative to the inviscid value, not the current value. After the first iteration the edge velocity shifted slightly. On iteration two, the correction was computed against the original again, so the clamp effectively prevented any further change. The edge velocity froze after iteration one. XFOIL does not have this safeguard. Drela just lets the Newton solver do its job.
If your solver needs a clamp to not blow up, the clamp is not the solution — understanding why it blows up is.
The Factor of Two
Near the leading edge, RustFoil’s inviscid edge velocities were consistently about twice XFOIL’s. Further downstream, they converged to within a few percent. Too systematic for a random bug. Not constant enough for a simple normalisation error.
The answer: repanelling. XFOIL’s PANE command distributes panel nodes with specific curvature-based refinement near the leading edge. RustFoil’s repanelling was producing slightly different node positions. Near the stagnation point, even tiny differences in panel placement produce large differences in computed velocity. When we fed XFOIL’s own panelling into RustFoil’s solver, every velocity matched to five significant figures.
This also taught me a painful testing lesson. I had a helper, make_naca0012(160), that generates a NACA 0012 with 160 points. It produces different geometry from XFOIL’s own panelled output. For weeks I was comparing one geometry against another and wondering why the numbers disagreed.
The Cascading Transition Disaster
XFOIL predicts laminar-to-turbulent transition using an envelope eN method. RustFoil was triggering transition 23× earlier than XFOIL. At α=4°, we were transitioning at 0.65% chord instead of 14.75%. The entire upper surface was turbulent when it should have had 15% laminar flow.
The error chain was brutal. Each step amplified the previous one. A few percent error in panel placement cascaded into a 2,300% error in transition location. This is exactly why aerodynamics is hard.
The Version Nobody Told You About
RustFoil’s stall prediction was off in a Reynolds-number-dependent way: too early at low Re, too late at high Re. On a cambered NACA 4412 at Re=1M, the error was 78%.
We had ported XFOIL’s 1987 amplification rate function (DAMPL) instead of the 1996 update (DAMPL2). XFOIL’s source contains both. There is no comment explaining that one supersedes the other. The newer version adds an exponential correction for laminar separation bubbles and blends to an Orr-Sommerfeld rate when the kinematic shape factor exceeds 3.5.
The kind of bug that looks like a physics modelling deficiency until you realise you are running the wrong version of the physics model.
Four Bugs in One Equation
After transition, the turbulent boundary layer was blowing up. The kinematic shape factor would collapse from a healthy 2.5 to 1.05, and the shear stress coefficient would hit its maximum clamp within two stations.
The shear-lag equation in xblsys.f had four independent bugs in our port: USA (averaging Hs instead of computing Us), CQA (using Hk where XFOIL uses CQ), CQ itself (not computed at all), and DEA (recomputed from an approximate formula instead of the stored blvar value).
Finding it required reading the Fortran carefully enough to understand what each variable is, not just what its name suggests. USA is not the country.
When the Newton Solver Points the Wrong Way
This one was infuriating. The Newton iteration was producing updates with the wrong sign. XFOIL would compute dθ ≈ −5×10−6, and RustFoil would compute dθ ≈ +15.7×10−6. The 4×4 solve was correct. The closures were validated against 60 test vectors. The Jacobian was within 0.1%. And the update still pointed the wrong direction.
I lost a week on this. The debugging plan had five phases — matrix assembly, sign conventions, XFOIL’s formulation of J·dx = −r versus J·dx = r, the lot. The fix turned out to be in the shape equation Jacobian row. Two missing terms. But finding them took running both solvers on the same input, dumping the full 4×4 system at each station, and diffing element by element.
The Wake Ghost
XFOIL’s stagnation point can move during the viscous iteration. When it does, panels get reclassified: wake rows become airfoil rows, or the other way around. We had a bug where a reclassified row retained its wake state — a wake-gap thickness value (dw) that made the last airfoil row roughly 10× too heavy in the next mass reconstruction.
A related bug: store_to_row() never wrote dw for wake rows and never cleared it for non-wake rows. On a clean attached-flow case, invisible. On a stalled NACA 4412 at 20°, it blew up the entire Newton iteration.
I also found that the transition station index (itran) was stored as a 0-based Rust index when the rest of the state model expected XFOIL’s 1-based numbering. Two lines to fix. A full day to find.
The Sign in the Stagnation Coupling
This one I genuinely couldn’t believe. The lower-side inviscid velocity was being seeded with (-q_inv).abs() instead of the signed velocity component. For attached flow, the magnitudes are close enough that it barely matters. For post-stall flow at 20°, where the stagnation point is migrating rapidly and the lower leading-edge velocities change sign? It matters completely.
The wrong sign propagated through the coupled Newton assembly into the cross-surface forcing term (DULE2), which fed the upper-surface momentum equation, which drove the relaxation factor, which determined whether the stagnation point moved on iteration 9 or iteration 10. One .abs() call. One Newton cycle off.
How 365 Tests and an Instrumented Fortran Binary Got Us Here
The bulk of the Fortran-to-Rust translation was done by AI coding agents. They were remarkably good at the syntactic conversion: translating array indexing, handling Fortran’s implicit typing, converting GOTO-based control flow into structured Rust. For the boring, mechanical parts — porting a 200-line subroutine line by line — they saved me weeks.
What they could not do was understand the intent. When a subroutine applies a seemingly arbitrary correction factor, the agent faithfully ports the arithmetic. But deciding whether that factor is a physical model, a numerical stabilisation trick, or a bug that happened to work — that requires understanding the aerodynamics. Every bug above was found by a human reading the Fortran, not by an AI.
The instrumented Fortran binary is as important as the tests. I forked XFOIL and added JSON debug dumps at every major checkpoint: BLVAR, BLDIF, MRCHUE, SETBL, BLSOLV, UPDATE, QDCALC, TRCHEK2, UESET, STFIND, STMOVE. When a parity test fails, I run both solvers on the same geometry and diff element by element.
Building the instrumentation was one of the few places where the agents were genuinely faster than me. XFOIL’s Fortran has no debug output to speak of. My approach was to add structured JSON dumps at the entry and exit of every major subroutine: inputs going in, outputs coming out, key intermediate values along the way. This is exactly the kind of work agents are good at and that I find soul-destroying: read a 200-line Fortran subroutine, identify every variable that matters, write a WRITE statement that serialises it to JSON, handle the array indexing, make sure the I/O does not change the execution order. I pointed an agent at each subroutine, told it what I needed, and it produced the instrumentation. It did in hours what would have taken me days of tedious Fortran editing.
The real value was not in any single dump — it was in the comparison pipeline. Once both solvers emit the same structured data at the same checkpoints, you can write a harness that runs a test case through both, parses both JSON outputs, and diffs them field by field. When the NACA 4412 at 20° diverges at iteration 9, you do not have to guess where. You can see that SETBL row IV=2 has a Jacobian entry that is 2.3× larger in Rust, trace that back to a BLDIF sensitivity that disagrees, and from there to the specific closure derivative that is wrong. You can actually find things instead of guessing.
The instrumentation is also what made the “four bugs in one equation” discovery possible. Without the ability to compare BLVAR outputs field by field at each station, I would have been staring at a divergent shape factor and guessing which of a dozen closure quantities was wrong. With the dumps, I could see that CQ was zero in Rust and nonzero in XFOIL, that US was just an average of Hs instead of the proper formula, and that DEA was close but not quite right — all in the same comparison run. Four bugs, found in one sitting, because the data was there to look at.
There is one thing the validation sweep still does not do automatically, and it should: split the polar. XFOIL’s viscous solver is path-dependent — each point uses the previous converged solution as its initial guess. The right technique is two legs: 0° → −15° and 0° → +25°, so you always approach from a known-good starting point. Some remaining parity misses near stall are almost certainly sweep-direction artefacts.
The AI did the boring parts. Drela did the brilliant parts forty years ago. I did the bit in the middle — understanding enough of both to connect them, and writing enough tests to prove the connection holds.
Why This Changed My View of XFOIL
I should say: I know Mark — we have been in a few of the same meetings, and I asked his permission to do this over beers rather than over email. He said yes, which was gracious. I have worked with Guppy extensively, and I have always known he was sharp. But there is a difference between using someone’s work and living inside it for two months, trying to reproduce every numerical decision they made.
This project gave me a much deeper appreciation for quite how brilliant the work is. The cleverness is not where you would expect — it is not in the headline algorithms that get written up in papers. It is in the hundreds of small decisions that keep the whole thing from falling apart: the sign convention at the stagnation point, the specific way the wake-to-airfoil handoff clears state, the particular relaxation limits that let the Newton solver run free without blowing up. Every time I found a bug in RustFoil, the fix was always “do what Mark did.”
The result is a solver that, forty years later, remains the standard reference. Every comparison in every paper uses XFOIL as the benchmark. The eN transition model is still the practical standard. The global Newton coupling scheme is still the architecture everyone builds on.
The Rust version is faster in places, dramatically more accessible, and much easier to instrument. But the solver — the part that actually matters — is Mark’s, through and through.
The full RustFoil source will be publicly available shortly. I just need to tidy up the codebase and rename a few variables that acquired their names during late-night debugging sessions and are not suitable for polite company.
Part II: The Geometry Problem
What happens when you deflect a flap and the solver stops converging — and why the fix was porting XFOIL’s GDES FLAP procedure line by line.
Read Part II →