fix(bem): fix windmill/helicopter branch discontinuity at lambda_climb boundary

Cargo.toml

@@ -4,7 +4,7 @@ members = ["dynbem", "dynbem_rs", "validation_rs"]
 
 # Fields below are inherited by member crates via `field.workspace = true`.
 [workspace.package]
-version = "0.5.0"
+version = "0.6.0"
 edition = "2021"
 license = "MIT"
 authors = ["Kristof"]

docs/VPM_DESIGN.md

@@ -677,13 +677,13 @@ a reasonable wake size (N=5,000) evaluated with the Barnes-Hut tree, parallel.
 
 | Model | ms/step | x Oye |
 |---|---:|---:|
-| Oye (2-stage annular filter) | 0.105 | 1x |
-| Pitt-Peters (3-state L-matrix) | 0.107 | 1x |
-| Quasi-static BEM | 10.95 | 104x |
-| VPM (Barnes-Hut, parallel, N=5,000) | 54 | 514x |
+| Oye (2-stage annular filter) | 0.016 | 1x |
+| Pitt-Peters (3-state L-matrix) | 0.016 | 1x |
+| Quasi-static BEM | 2.49 | 156x |
+| VPM (Barnes-Hut, parallel, N=5,000) | 4.1 | 256x |
 
 Oye and Pitt-Peters are algebraic inflow relations (near-free). The BEM is
-~100x slower from its per-station Brent root-find. The VPM is not an inflow
+~150x slower from its per-station Brent root-find. The VPM is not an inflow
 model but a wake-resolving tool: its cost is dominated by the Biot-Savart
 evaluation and scales with the particle count N.
 
@@ -692,18 +692,18 @@ evaluation and scales with the particle count N.
 
 | N | direct seq | direct par | BH seq | BH par |
 |---:|---:|---:|---:|---:|
-| 2,000 | 21 | 11 | 39 | 15 |
-| 5,000 | 198 | 59 | 113 | 53 |
-| 10,000 | 778 | 227 | 197 | 108 |
-| 16,000 | 2,633 | 428 | 708 | 139 |
-| 32,000 | 8,750 | 2,265 | 1,058 | 337 |
+| 2,000 | 4 | 1 | 4 | 2 |
+| 5,000 | 25 | 4 | 12 | 4 |
+| 10,000 | 97 | 14 | 32 | 7 |
+| 16,000 | 297 | 33 | 59 | 12 |
+| 32,000 | 1,075 | 134 | 121 | 30 |
 
 The direct sum is O(N^2) (cost quadruples per doubling of N); the Barnes-Hut
 tree (`VpmRotorConfig::barnes_hut`, off by default) lumps distant particle
 clusters into single equivalent vortices and grows O(N log N), so it overtakes
 the direct sum as the wake grows -- comparable at N=5k, ~3x faster at N=16k,
-~7x at N=32k. Rayon parallelism (default) gives ~2-4x over sequential for the
-direct sum and ~2x for the tree (`--seq`/`--par` isolate it). Absolute ms are
+~4x at N=32k. Rayon parallelism (default) gives ~4-8x over sequential for the
+direct sum and ~2-4x for the tree (`--seq`/`--par` isolate it). Absolute ms are
 machine-dependent; the ratios are the point.
 
 ### 7.1 Barnes-Hut tree (O(N log N))
@@ -744,9 +744,9 @@ vs direct ($<5\%$ of peak at $\theta = 0.5$).
 
 | Validation (`validation_rs`) | Reference | What it checks | Result |
 |---|---|---|---|
-| `blade_element_hover` | Combined BEMT, hover (Leishman ch. 3) | Hover thrust coefficient vs closed form | within ~15-25% |
+| `blade_element_hover` | Combined BEMT, hover (Leishman ch. 3) | Hover thrust coefficient vs closed form | within ~14-18% |
 | `climb_momentum` | Axial-climb momentum theory | `C_T` ~ 2 lam_i (lam_i + lam_c); loads fall with climb | PASS (monotone) |
-| `glauert_forward_inflow` | Glauert forward-flight inflow | Disk inflow + wake-skew angle | skew <1.2%; inflow <26% |
+| `glauert_forward_inflow` | Glauert forward-flight inflow | Disk inflow + wake-skew angle | skew <1.2%; inflow ~26% |
 | `wake_skew` | Wake-skew geometry | Skew grows with mu; covariant under X/Y rotation | PASS |
 | `prandtl_tip_loss` | Directional | Tip-loss flag reduces global loads | PASS |
 | `autorotation` | Directional | Negative-torque branch reached in descent + edgewise | PASS |
@@ -756,7 +756,7 @@ vs direct ($<5\%$ of peak at $\theta = 0.5$).
 | `servo_flap` | Directional | Kaman servo-flap feathering (zero / collective / cyclic) | PASS |
 | `cyclic_phase_servo` | Directional | Direct-mech pitching `My` vs servo-flap rolling `Mx` | PASS |
 
-Hover thrust tracks BEMT to within ~15-25%; coning `a0` and longitudinal
+Hover thrust tracks BEMT to within ~14-18%; coning `a0` and longitudinal
 flapping `a1` match the closed forms to ~14% (the disk tilts back by the right
 amount). Lateral flapping `b1` is currently under-predicted -- see the TODO.
 

dynbem/pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "dynbem"
-version = "0.5.0"
+version = "0.6.0"
 description = "Rotor aerodynamics: BEM, Pitt-Peters, Oye dynamic inflow, and VPM free-wake -- all in one library"
 requires-python = ">=3.10"
 license = "MIT"

dynbem/python/dynbem/__init__.py

@@ -6,10 +6,7 @@
 """
 
 from ._dynbem import (  # noqa: F401
-    vrs_lambda1,
     cyclic_coeffs as _cyclic_coeffs_rust,
-    prandtl_tip_loss,
-    prandtl_hub_loss,
     LinearPolar,
     TabulatedPolar,
     QuasiStaticRotorState,
@@ -109,8 +106,7 @@ def VpmRotor(defn, polar=None, **config):  # noqa: N802
 
 __all__ = [
     # functions
-    "vrs_lambda1", "cyclic_coeffs",
-    "prandtl_tip_loss", "prandtl_hub_loss",
+    "cyclic_coeffs",
     "create_aero", "build_polar", "load_tabulated_polar",
     "solve_trim_cyclic", "relax_inflow",
     # types

dynbem/python/dynbem/bem.py

@@ -1,17 +1,11 @@
 """dynbem.bem submodule (compat shim).
 
-Re-exports QuasiStaticBEM (legacy alias BEMModel), the Prandtl loss
-helpers, and the per-annulus solve_bem_element (with BEMElementResult)
-so legacy dotted-path imports continue to work.
+Re-exports QuasiStaticBEM (legacy alias BEMModel) so legacy dotted-path
+imports continue to work.
 """
-from . import BEMModel, QuasiStaticBEM, prandtl_hub_loss, prandtl_tip_loss
-from ._dynbem import BEMElementResult, solve_bem_element
+from . import BEMModel, QuasiStaticBEM
 
 __all__ = [
     "QuasiStaticBEM",
     "BEMModel",
-    "BEMElementResult",
-    "prandtl_hub_loss",
-    "prandtl_tip_loss",
-    "solve_bem_element",
 ]

dynbem/python/dynbem/pitt_peters.py

@@ -1,4 +1,4 @@
 """dynbem.pitt_peters submodule (compat shim)."""
-from . import PittPetersModel, vrs_lambda1
+from . import PittPetersModel
 
-__all__ = ["PittPetersModel", "vrs_lambda1"]
+__all__ = ["PittPetersModel"]

dynbem/src/lib.rs

@@ -14,11 +14,6 @@ mod wrappers;
 use trim_py::{relax_inflow_py, solve_trim_cyclic_py, PyTrimResult};
 use wrappers::*;
 
-#[pyfunction]
-fn vrs_lambda1(lambda2: f64) -> f64 {
-    dynbem_rs::common::vrs_lambda1(lambda2)
-}
-
 #[pyfunction]
 #[pyo3(signature = (tilt_lon, tilt_lat, control = None))]
 fn cyclic_coeffs(tilt_lon: f64, tilt_lat: f64, control: Option<PyControlProperties>) -> (f64, f64) {
@@ -35,122 +30,9 @@ fn cyclic_coeffs(tilt_lon: f64, tilt_lat: f64, control: Option<PyControlProperti
     dynbem_rs::cyclic::cyclic_coeffs(tilt_lon, tilt_lat, gains)
 }
 
-#[pyfunction]
-fn prandtl_tip_loss(n_blades: usize, x: f64, phi_rad: f64) -> f64 {
-    dynbem_rs::quasi_static_bem::prandtl_tip_loss(n_blades, x, phi_rad)
-}
-
-#[pyfunction]
-fn prandtl_hub_loss(n_blades: usize, x: f64, x_hub: f64, phi_rad: f64) -> f64 {
-    dynbem_rs::quasi_static_bem::prandtl_hub_loss(n_blades, x, x_hub, phi_rad)
-}
-
-// ---------------------------------------------------------------------------
-// solve_bem_element: per-annulus BEM solver, exposed for diagnostics +
-// the spanwise-CL verification scripts. Mirrors the legacy Python
-// dynbem.bem.solve_bem_element NamedTuple API.
-// ---------------------------------------------------------------------------
-
-#[pyclass(name = "BEMElementResult", module = "dynbem._dynbem")]
-#[derive(Clone, Debug)]
-#[allow(non_snake_case)]
-pub struct PyBEMElementResult {
-    #[pyo3(get)]
-    pub lambda_r: f64,
-    #[pyo3(get)]
-    pub a_prime: f64,
-    #[pyo3(get)]
-    pub dT: f64,
-    #[pyo3(get)]
-    pub dQ: f64,
-    #[pyo3(get)]
-    pub momentum_residual: f64,
-}
-
-#[pyfunction]
-#[pyo3(signature = (
-    r, dr, chord, twist_rad, collective_rad, omega, v_climb, rho,
-    n_blades, radius_m, polar, use_tip_loss,
-    v_t_extra = 0.0, root_cutout_m = 0.0,
-))]
-fn solve_bem_element(
-    r: f64,
-    dr: f64,
-    chord: f64,
-    twist_rad: f64,
-    collective_rad: f64,
-    omega: f64,
-    v_climb: f64,
-    rho: f64,
-    n_blades: usize,
-    radius_m: f64,
-    polar: &Bound<'_, PyAny>,
-    use_tip_loss: bool,
-    v_t_extra: f64,
-    root_cutout_m: f64,
-) -> PyResult<PyBEMElementResult> {
-    let polar = wrappers::extract_polar(polar)?;
-    let res = match polar {
-        wrappers::ResolvedPolar::Linear(p) => {
-            let geom = dynbem_rs::quasi_static_bem::BEMElementGeometry::new(
-                r,
-                dr,
-                chord,
-                twist_rad,
-                omega,
-                rho,
-                n_blades,
-                radius_m,
-                &p,
-                use_tip_loss,
-                root_cutout_m,
-            );
-            dynbem_rs::quasi_static_bem::solve_bem_element(
-                &geom,
-                collective_rad,
-                v_climb,
-                v_t_extra,
-            )
-        }
-        wrappers::ResolvedPolar::Tabulated(p) => {
-            let geom = dynbem_rs::quasi_static_bem::BEMElementGeometry::new(
-                r,
-                dr,
-                chord,
-                twist_rad,
-                omega,
-                rho,
-                n_blades,
-                radius_m,
-                &p,
-                use_tip_loss,
-                root_cutout_m,
-            );
-            dynbem_rs::quasi_static_bem::solve_bem_element(
-                &geom,
-                collective_rad,
-                v_climb,
-                v_t_extra,
-            )
-        }
-    };
-    Ok(PyBEMElementResult {
-        lambda_r: res.lambda_r,
-        a_prime: res.a_prime,
-        dT: res.d_t,
-        dQ: res.d_q,
-        momentum_residual: res.momentum_residual,
-    })
-}
-
 #[pymodule]
 fn _dynbem(m: &Bound<'_, PyModule>) -> PyResult<()> {
-    m.add_function(wrap_pyfunction!(vrs_lambda1, m)?)?;
     m.add_function(wrap_pyfunction!(cyclic_coeffs, m)?)?;
-    m.add_function(wrap_pyfunction!(prandtl_tip_loss, m)?)?;
-    m.add_function(wrap_pyfunction!(prandtl_hub_loss, m)?)?;
-    m.add_function(wrap_pyfunction!(solve_bem_element, m)?)?;
-    m.add_class::<PyBEMElementResult>()?;
     m.add_class::<PyLinearPolar>()?;
     m.add_class::<PyTabulatedPolar>()?;
     m.add_class::<PyBladeGeometry>()?;

dynbem_rs/src/common.rs

@@ -58,3 +58,41 @@ pub fn vrs_lambda1(k: f64) -> f64 {
     let k4 = k2 * k2;
     1.0 + VRS_C[0] * k + VRS_C[1] * k2 + VRS_C[2] * k3 + VRS_C[3] * k4
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// At lambda2 = 0 (hover), induced velocity ratio must equal 1.0 exactly
+    /// (polynomial anchor point).
+    #[test]
+    fn test_vrs_lambda1_hover_boundary() {
+        assert!(
+            (vrs_lambda1(0.0) - 1.0).abs() < 1e-10,
+            "vrs_lambda1(0) = {:.10}, expected 1.0",
+            vrs_lambda1(0.0)
+        );
+    }
+
+    /// At lambda2 = 2.0 (wind-turbine-brake / WBS entry), the ratio returns
+    /// close to 1.0 (+/-5%), marking the upper end of the VRS region.
+    #[test]
+    fn test_vrs_lambda1_wbs_boundary() {
+        let v = vrs_lambda1(2.0);
+        assert!(
+            (0.95..=1.10).contains(&v),
+            "vrs_lambda1(2.0) = {v:.4}, expected 0.95..1.10"
+        );
+    }
+
+    /// In mid-VRS (lambda2 ~ 1), induced velocity exceeds the hover value --
+    /// this is the VRS peak.
+    #[test]
+    fn test_vrs_lambda1_peak_above_hover() {
+        assert!(
+            vrs_lambda1(1.0) > 1.2,
+            "vrs_lambda1(1.0) = {:.4}, expected > 1.2",
+            vrs_lambda1(1.0)
+        );
+    }
+}