Viscous solver does not converge for Mach>0 (BL march uses M2=0, newton system uses real M2)

[denriquezfirestorm]

Summary
Viscous solves with any non-zero Mach number do not converge. Cases that converge in a few iterations at Mach 0 run to max_iter at Mach 0.25. The residual settles at a constant value and doesn't drop when max_iter is raised. CL and CM also don't change with Mach.

Environment

• flexfoil 1.1.6, built from source at cf303d7 (the same code is in the current main)
• Windows 11 x64, Python 3.13.5 (Anaconda), Rust 1.98.1 (MSVC)

Reproduction

import flexfoil

foil = flexfoil.load("mh60.dat")   # UIUC database file
for mach in (0.0, 0.25):
    for max_iter in (100, 300):
        r = foil.solve(4.0, Re=1.6e6, mach=mach, ncrit=7,
                       max_iter=max_iter, store=False)
        print(mach, max_iter, r.converged, r.iterations,
              f"{r.residual:.2e}", f"{r.cl:.4f}", f"{r.cd:.5f}")

Observed results (α = 4°, Re = 1.6e6, Ncrit = 7)
Airfoil Mach Converged Iterations Residual CL CD
mh60 0.00 Yes 4 4.14e-06 0.5587 0.00595
mh60 0.25 No 100 (and 300) 2.98e-02 (same at both limits) 0.5588 0.00587
e186 0.00 Yes 6 6.84e-05 0.4384 0.00551
e186 0.25 No 100 (and 300) 3.41e-01 (same at both limits) 0.4418 0.00545
NACA 2412 0.25 No 100 / 300 5.01e-01 / 4.81e-01 0.6815 0.00655

In a full sweep (mh60 and e186, α = −5° to 10°, Mach 0.25), 0 of 32 points converged.
Expected
Mach 0.25 should converge about as well as Mach 0. CL should rise by roughly 3% (the Prandtl–Glauert factor 1/√(1−M²)), as it does in XFOIL.
Suspected cause

From reading the code; not yet confirmed by a patched build:

1. Mach² is 0 in the boundary-layer march but not in the Newton system. mrchue and mrchdu hard-code let msq = 0.0_f64; (crates/rustfoil-xfoil/src/march.rs, lines 271 and 778), and setbl doesn't pass Mach to them. The Newton system uses the real value: build_global_system(..., mach * mach, ...) in assembly.rs, and msq = mach * mach in update.rs. Each iteration seems to re-march the boundary layer at M = 0 and then correct toward the M > 0 solution, so the correction stays about the same size every iteration. That would explain a residual that stays constant when max_iter is raised, and why M = 0 converges normally.
2. The boundary-layer closures are only partly compressible. blvar uses the freestream M² instead of the local edge Mach. It also leaves out ∂M²/∂U: let m_u = 0.0; // Simplified for now (crates/rustfoil-bl/src/equations.rs, line 297). blmid has the same simplification: let ma = 0.5 * msq; // Simplified - would need proper M at each station (line 946). There is also no Kármán–Tsien correction of the edge velocity, which XFOIL applies (TKBL in BLKIN).
3. No compressibility correction in forces. compute_panel_forces_from_gamma (crates/rustfoil-xfoil/src/forces.rs) uses cp = 1 - g² and takes no Mach input, so CL and CM don't change with Mach. That matches mh60 above: CL 0.5587 at M = 0 and 0.5588 at M = 0.25.
4. No test coverage. No tests in rustfoil-xfoil use a Mach number above 0.
   Possible fix

• Pass mach * mach from setbl into mrchue and mrchdu instead of using 0.0. This should restore convergence.
• For XFOIL-matching results, also add the Kármán–Tsien edge-velocity correction, local Mach and its velocity derivatives in blvar/blmid, and compressible Cp/CL/CM.
• Add a regression test at Mach 0.25 or higher.
  Workaround
  Run at M = 0 and apply a Prandtl–Glauert or Kármán–Tsien correction to CL and CM afterwards.