fix: reference-matched native drag decomposition

crates/rustfoil-solver/src/viscous/forces.rs

@@ -115,8 +115,8 @@ pub fn compute_forces_from_canonical_state(
 ) -> AeroForces {
     let (cl, cm) = compute_panel_forces_from_gamma(panel_x, panel_y, state.panel_gamma(), alpha_rad);
 
-    let cd_friction_upper = compute_row_friction_drag(&state.upper_rows);
-    let cd_friction_lower = compute_row_friction_drag(&state.lower_rows);
+    let cd_friction_upper = compute_row_friction_drag(&state.upper_rows, alpha_rad);
+    let cd_friction_lower = compute_row_friction_drag(&state.lower_rows, alpha_rad);
     let cd_friction = cd_friction_upper + cd_friction_lower;
 
     let lower_wake: Vec<&CanonicalBlRow> = state.lower_rows.iter().filter(|row| row.is_wake).collect();
@@ -621,44 +621,18 @@ pub fn compute_friction_drag(stations: &[BlStation]) -> f64 {
     cdf
 }
 
-fn compute_row_friction_drag(rows: &[CanonicalBlRow]) -> f64 {
-    if rows.len() < 4 {
-        return 0.0;
-    }
-
-    let max_physical_cf = 0.02;
-    let r_theta_min = 80.0;
-
-    rows[3..]
-        .windows(2)
+fn compute_row_friction_drag(rows: &[CanonicalBlRow], alpha_rad: f64) -> f64 {
+    let (sa, ca) = alpha_rad.sin_cos();
+    // XFOIL CDCALC: IBL=3..IBLTE integrates real rows 1→2 onward, excluding the stagnation anchor.
+    rows.get(1..).unwrap_or_default().windows(2)
+        .filter(|pair| !pair[0].is_wake && !pair[1].is_wake)
         .map(|pair| {
-            if pair[0].is_wake || pair[1].is_wake {
-                return 0.0;
-            }
-            if pair[0].theta < 1.0e-10 || pair[1].theta < 1.0e-10 {
-                return 0.0;
-            }
-            if pair[0].r_theta < r_theta_min || pair[1].r_theta < r_theta_min {
-                return 0.0;
-            }
-
-            let dx = (pair[1].x_coord - pair[0].x_coord).abs();
-            let dx = if dx < 1.0e-12 {
-                (pair[1].x - pair[0].x).abs()
-            } else {
-                dx
-            };
-
-            let cf1 = pair[0].cf.clamp(0.0, max_physical_cf);
-            let cf2 = pair[1].cf.clamp(0.0, max_physical_cf);
-            let cf_avg = 0.5 * (cf1 + cf2);
-            let ue1 = pair[0].uedg.abs().max(0.01);
-            let ue2 = pair[1].uedg.abs().max(0.01);
-            let ue_sq_avg = 0.5 * (ue1 * ue1 + ue2 * ue2);
-
-            ue_sq_avg * cf_avg * dx
-        })
-        .sum()
+            let dx = (pair[1].x_coord - pair[0].x_coord) * ca
+                + (pair[1].y_coord - pair[0].y_coord) * sa;
+            let twice_tau1 = pair[0].cf * pair[0].uedg.powi(2);
+            let twice_tau2 = pair[1].cf * pair[1].uedg.powi(2);
+            0.5 * (twice_tau1 + twice_tau2) * dx
+        }).sum()
 }
 
 fn compute_cd_from_te_rows(upper_rows: &[CanonicalBlRow], lower_rows: &[CanonicalBlRow]) -> f64 {
@@ -1015,6 +989,27 @@ pub fn transition_drag_penalty(stations: &[BlStation], x_transition: f64) -> f64
 mod tests {
     use super::*;
 
+    #[test]
+    fn canonical_friction_projects_wall_shear_and_keeps_low_reynolds_rows() {
+        let rows = vec![
+            CanonicalBlRow { cf: 1000.0, ..Default::default() },
+            CanonicalBlRow { x_coord: 0.1, y_coord: 0.0, cf: 0.01, uedg: 2.0, r_theta: 20.0, ..Default::default() },
+            CanonicalBlRow { x_coord: 0.4, y_coord: 0.2, cf: 0.02, uedg: 1.0, r_theta: 30.0, ..Default::default() },
+            CanonicalBlRow { x_coord: 1.0, cf: 1000.0, uedg: 1.0, is_wake: true, ..Default::default() },
+        ];
+        assert!((compute_row_friction_drag(&rows, std::f64::consts::FRAC_PI_2) - 0.006).abs() < 1e-12);
+    }
+
+    #[test]
+    fn canonical_friction_retains_signed_shear_and_projection() {
+        let mut rows = vec![CanonicalBlRow::default(),
+            CanonicalBlRow { x_coord: 0.1, cf: -0.01, uedg: 1.0, ..Default::default() },
+            CanonicalBlRow { x_coord: 0.6, cf: -0.01, uedg: 1.0, ..Default::default() }];
+        assert!((compute_row_friction_drag(&rows, 0.0) + 0.005).abs() < 1e-12);
+        rows[2].x_coord = -0.4;
+        assert!((compute_row_friction_drag(&rows, 0.0) - 0.005).abs() < 1e-12);
+    }
+
     fn create_test_stations(n: usize) -> Vec<BlStation> {
         (0..n)
             .map(|i| {

crates/rustfoil-xfoil/src/forces.rs

@@ -119,18 +119,18 @@ fn canonical_force_state(
     let mut canonical = XfoilLikeViscousState::new(state.n_panel_nodes());
     canonical.gam = state.gam.clone();
     canonical.gam_a = state.gam_a.clone();
-    canonical.upper_rows = upper_stations.iter().map(|s| canonical_row_from_station(s, mach)).collect();
-    canonical.lower_rows = lower_stations.iter().map(|s| canonical_row_from_station(s, mach)).collect();
+    canonical.upper_rows = upper_stations.iter().zip(&state.upper_rows).map(|(s, row)| canonical_row_from_station(s, mach, row.y_coord)).collect();
+    canonical.lower_rows = lower_stations.iter().zip(&state.lower_rows).map(|(s, row)| canonical_row_from_station(s, mach, row.y_coord)).collect();
     canonical
 }
 
-fn canonical_row_from_station(station: &BlStation, mach: f64) -> CanonicalBlRow {
+fn canonical_row_from_station(station: &BlStation, mach: f64, y_coord: f64) -> CanonicalBlRow {
     let u = rustfoil_bl::compressibility::edge_velocity(station.u, mach * mach).0;
     let density = rustfoil_bl::compressibility::edge_thermodynamics(u, station.theta, mach * mach, 1.0).density;
     CanonicalBlRow {
         x: station.x,
         x_coord: station.x_coord,
-        y_coord: 0.0,
+        y_coord,
         panel_idx: station.panel_idx,
         uedg: u,
         uinv: station.u,

docs-site/docs/internals/compressible-forces.mdx

@@ -7,7 +7,7 @@ import Inspector from '@site/src/components/Compressibility/Inspector';
 
 ## What it is
 
-The faithful solver's pressure/force family converts normalized incompressible panel gamma to pressure coefficient at the requested freestream Mach. It integrates that pressure on the existing contour for CL and quarter-chord CM. Total wake drag uses physical edge speed and density. The separate friction/pressure drag decomposition retains a legacy integrator and is not reference matched. Python/CLI coefficient results, trial operating corrections and the shipping WASM Cp array use the same conventions. Coordinates and thicknesses are normalized by chord; gamma and U are normalized by freestream speed. Angles passed to force integration are radians.
+The faithful solver's pressure/force family converts normalized incompressible panel gamma to pressure coefficient at the requested freestream Mach. It integrates that pressure on the existing contour for CL and quarter-chord CM. Total wake drag uses physical edge speed and density. The separate CDf and CD−CDf drag fields use XFOIL CDCALC wall-shear projection and have bounded reference coverage. Python/CLI coefficient results, trial operating corrections and the shipping WASM Cp array use the same conventions. Coordinates and thicknesses are normalized by chord; gamma and U are normalized by freestream speed. Angles passed to force integration are radians.
 
 ## Theory
 
@@ -59,4 +59,22 @@ The change retains existing contour integration and canonical wake/friction rout
 
 <SourceCode path="docs-site/src/components/Compressibility/math.js" name="pressure" />
 
-The gate validates **total CD**, not the separate `cd_friction` / `cd_pressure` fields. The native `cd_pressure` field is total CD minus integrated skin friction. XFOIL's iteration report uses that same CD−CDf decomposition, but its polar-file CDp column is a distinct surface-pressure integral from CLCALC; those outputs must not be compared as if they were the same quantity. The legacy friction integrator excludes low momentum-thickness Reynolds-number rows and differs in wall-shear projection. Correcting and validating the matching decomposition needs separate reference-backed tests; total-CD agreement alone does not establish drag-component agreement.
+The reference gate now checks total CD, `cd_friction`, and `cd_pressure` on all 48 cases. `cd_pressure` means **CD−CDf**, matching XFOIL's iteration report. The polar-file CDp column is a distinct surface-pressure integral from CLCALC and is not this field. No new surface-pressure-drag output is added.
+
+CDCALC integrates local wall shear along the freestream direction, beginning with the first pair of real surface stations (IBL=3 in Fortran), excluding the synthetic stagnation row and wake:
+
+$$
+C_{Df}=\sum_{IS=1}^{2}\sum_{IBL=3}^{IBLTE(IS)}\frac{(\tau_{IBL} + \tau_{IBL-1})}{\rho_\infty U_\infty^2}\left[\Delta x\cos\alpha+\Delta y\sin\alpha\right],\qquad \tau=\frac12\rho_e U_e^2 C_{f,e}\label{eq:compressible-friction}
+$$
+
+Equation $\eqref{eq:compressible-friction}$ keeps the sign of wall shear and projection and includes low-Reynolds-number surface stations. The canonical force bridge carries physical speed, density-weighted Cf and actual surface y coordinates. Legacy clamp/low-Rtheta exclusions and absolute-x projection are removed from this canonical route; the separate legacy station-based engine remains unvalidated.
+
+<SourceCode path="crates/rustfoil-solver/src/viscous/forces.rs" name="compute_row_friction_drag" />
+<SourceCode path="crates/rustfoil-xfoil/src/forces.rs" name="canonical_force_state" />
+
+Measured maxima across the 48 cases are CDf 4.91e-6 and CD−CDf 4.99e-6 versus XFOIL's five-decimal iteration output, within the 1e-5 acceptance tolerance. Total CD and the other previously gated outputs keep their existing tolerances. This is bounded numerical agreement, not a new physical-validity or global-derivative claim.
+
+## Interactive wall-shear projection
+
+<Inspector mode="friction" />
+<SourceCode path="docs-site/src/components/Compressibility/Inspector.jsx" name="CompressibilityInspector" />

docs-site/docs/internals/solver-reference-gate.mdx

@@ -69,3 +69,5 @@ The documented nose/mode/platform ceilings would be raised with dedicated matche
 <SourceCode path="scripts/solver_gate.py" name="compare_distributions" />
 
 <SourceCode path="docs-site/src/components/SolverGate/Inspector.jsx" name="GateInspector" />
+
+The matching drag decomposition is read from the final converged XFOIL iteration report, not the polar-file CDp column. Both CDf and CD−CDf use absolute tolerance 1e-5. Their measured maximum differences on the 48-case matrix are 4.91e-6 and 4.99e-6.

docs-site/docs/python-api.mdx

@@ -412,3 +412,5 @@ The [`examples/`](https://github.com/flexcompute/flexfoil/tree/main/packages/fle
 | `10_flap_study.py` | Sweep flap deflections (0-20 deg) |
 | `11_flap_hinge_sweep.py` | Find optimal hinge position for L/D |
 | `12_matrix_sweep.py` | Alpha x Re matrix for a flapped airfoil |
+
+Low-level native `cd_friction` is signed freestream-projected wall-shear drag; `cd_pressure` is total CD minus that friction contribution, matching XFOIL’s iteration-report decomposition. It is not the separate surface-pressure integral labelled CDp in an XFOIL polar file. See [force reporting](./internals/compressible-forces) for the formula and bounded reference tests.

docs-site/src/components/Compressibility/Inspector.jsx

@@ -5,15 +5,22 @@ export default function CompressibilityInspector({ mode = 'edge' }) {
   const [mach, setMach] = useState(.25);
   const [ui, setUi] = useState(1.2);
   const [itran, setItran] = useState(80);
+  const [alpha, setAlpha] = useState(4);
+  const [cf, setCf] = useState(.005);
   const state = edgeState(ui, mach, .001, 1e6, 2.6);
   if (mode === 'transition') return <fieldset><legend>Surface versus wake index</legend>
     <label>ITRAN (one based): {itran} <input type="range" min="1" max="84" value={itran} onChange={e => setItran(Number(e.target.value))} /></label>
     <p>Last surface index IBLTE (zero based): 80</p>
     <output aria-live="polite">{surfaceTransition(itran, 80) ? 'Report the interpolated surface location, including the final surface interval' : 'Transition lies in the wake; report 1.0 chord'}</output>
   </fieldset>;
-  return <fieldset><legend>{mode === 'pressure' ? 'Karman–Tsien pressure' : 'Compressible station'}</legend>
+  return <fieldset><legend>{mode === 'pressure' ? 'Karman–Tsien pressure' : mode === 'friction' ? 'Signed wall-shear projection' : 'Compressible station'}</legend>
     <label style={{ display: 'block' }}>Freestream Mach: {mach.toFixed(2)} <input type="range" min="0" max="0.5" step="0.01" value={mach} onChange={e => setMach(Number(e.target.value))} /></label>
     <label style={{ display: 'block' }}>Incompressible speed / freestream: {ui.toFixed(2)} <input type="range" min="0.5" max="1.5" step="0.01" value={ui} onChange={e => setUi(Number(e.target.value))} /></label>
+    {mode === 'friction' && <>
+      <label style={{display:'block'}}>Angle (degrees): {alpha} <input type="range" min="-90" max="90" value={alpha} onChange={e=>setAlpha(Number(e.target.value))} /></label>
+      <label style={{display:'block'}}>Local Cf: {cf.toFixed(3)} <input type="range" min="-0.02" max="0.02" step="0.001" value={cf} onChange={e=>setCf(Number(e.target.value))} /></label>
+      <p>Equal shear at both endpoints; Δx/c = Δy/c = 0.02. Drag contribution = {(state.density * state.u ** 2 * cf * .02 * (Math.cos(alpha*Math.PI/180)+Math.sin(alpha*Math.PI/180))).toFixed(7)}. Shear and projection keep their signs.</p>
+    </>}
     <output aria-live="polite">{mode === 'pressure' ? `Cp = ${pressure(ui, mach).toFixed(6)} (incompressible Cp = ${(1-ui*ui).toFixed(6)})` : `Physical U = ${state.u.toFixed(6)}, local Mach = ${Math.sqrt(state.localMsq).toFixed(6)}, Hk = ${state.hk.toFixed(6)}, Rtheta = ${state.rtheta.toFixed(6)}`}</output>
     <p>Illustrative algebra only. The station uses theta = 0.001 chord, Re = 1,000,000 and H = 2.6. It does not run a coupled flow solution.</p>
   </fieldset>;

flexfoil-ui/package.json

@@ -2,7 +2,7 @@
   "name": "flexfoil-ui",
   "private": true,
   "license": "MIT",
-  "version": "1.1.15",
+  "version": "1.1.16",
   "type": "module",
   "scripts": {
     "dev": "vite",

flexfoil-ui/src/lib/version.ts

@@ -28,6 +28,9 @@ export interface ChangelogEntry {
 }
 
 export const CHANGELOG: ChangelogEntry[] = [
+  { version: '1.1.16', date: '2026-10-14', items: [
+    { category: 'fixed', text: 'Match native friction and total-minus-friction drag reporting to the reference wall-shear integration' },
+  ] },
   { version: '1.1.15', date: '2026-10-14', items: [
     { category: 'fixed', text: 'Preserve solver attempts and reuse only converged results from matching solver artifacts and operating settings' },
   ] },

packages/flexfoil-python/pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "flexfoil"
-version = "1.1.15"
+version = "1.1.16"
 description = "Airfoil analysis in Python — XFOIL-faithful solver with a local web UI"
 readme = "README.md"
 license = { text = "MIT" }

scripts/solver_gate.py

@@ -11,7 +11,7 @@
 ROOT = Path(__file__).resolve().parents[1]
 # XFOIL polar print half-units: CL/CM/trips 5e-5, CD 5e-6; both solvers stop at 1e-4 RMS.
 # These are bounded regression tolerances for this matrix, not physical-accuracy claims.
-TOLERANCES = dict(cl=2e-4, cd=1e-5, cm=1e-4, x_tr_upper=2e-4, x_tr_lower=2e-4)
+TOLERANCES = dict(cl=2e-4, cd=1e-5, cd_pressure=1e-5, cd_friction=1e-5, cm=1e-4, x_tr_upper=2e-4, x_tr_lower=2e-4)
 
 
 def require_close(actual, expected, tolerance, label):
@@ -48,6 +48,11 @@ def run_reference(binary, geometry, alpha, reynolds, mach, ncrit, trips, directo
     row = values[0]
     require_close(row[0], alpha, 1e-8, 'reference alpha')
     result = dict(zip(['alpha', 'cl', 'cd', 'cdp', 'cm', 'x_tr_upper', 'x_tr_lower'], row))
+    # XFOIL polar CDp is the surface-pressure integral, not its reported CD-CDf component.
+    components = re.findall(r'CDf =\s*([+\-\d.eE]+)\s+CDp =\s*([+\-\d.eE]+)', process.stdout)
+    if not components:
+        raise AssertionError('Reference drag components missing')
+    result['cd_friction'], result['cd_pressure'] = map(float, components[-1])
     distribution = [list(map(float, line.split())) for line in
                     (directory / 'bl.txt').read_text().splitlines()[1:] if line.strip()]
     points = [list(map(float, line.split())) for line in geometry.read_text().splitlines()[1:]]

scripts/test_solver_gate.py

@@ -33,6 +33,22 @@ def test_reference_process_failure_is_fatal(self):
                 with self.assertRaisesRegex(AssertionError, 'XFOIL exit 2'):
                     run_reference(Path('unused'), geometry, 4., 1e6, 0., 9., (1., 1.), root)
 
+    def test_drag_components_use_final_iteration_not_polar_pressure(self):
+        with tempfile.TemporaryDirectory() as name:
+            root = Path(name)
+            geometry = root / 'input.dat'
+            geometry.write_text('test\n1 0\n0 0\n1 -.1\n')
+            (root / 'polar.txt').write_text('1 1 Reynolds number fixed\nxtrf = 1.000 (top) 1.000 (bottom)\nMach = 0.000 Re = 1.000 e 6 Ncrit = 9.000 9.000\nalpha CL CD CDp CM XtrT XtrB\n------ ----- -----\n4 .5 .00716 -.00010 -.05 .3 1 0 0\n')
+            (root / 'bl.txt').write_text('s x y ue ds th cf hk\n0 1 0 1 .001 .0001 .002 2.5\n1 0 0 1 .001 .0001 .002 2.5\n2 1 -.1 1 .001 .0001 .002 2.5\n')
+            with patch('solver_gate.subprocess.run', return_value=subprocess.CompletedProcess([], 0, 'CDf = 0.00600 CDp = 0.00116\nCDf = 0.00465 CDp = 0.00250', '')):
+                result, _ = run_reference(Path('unused'), geometry, 4., 1e6, 0., 9., (1., 1.), root)
+            self.assertEqual(result['cdp'], -.00010)
+            self.assertEqual(result['cd_friction'], .00465)
+            self.assertEqual(result['cd_pressure'], .00250)
+            with patch('solver_gate.subprocess.run', return_value=subprocess.CompletedProcess([], 0, '', '')):
+                with self.assertRaisesRegex(AssertionError, 'drag components missing'):
+                    run_reference(Path('unused'), geometry, 4., 1e6, 0., 9., (1., 1.), root)
+
 
 if __name__ == '__main__':
     unittest.main()