alice-physics 3.0.1

//! Turbulence Closure Models — Smagorinsky LES + k-ε + k-ω
//!
//! Phase G1 of the ALICE-Physics completeness project. Provides the three
//! most common Reynolds-Averaged (RANS) and Large Eddy Simulation (LES)
//! turbulence closures — each computes an effective **eddy viscosity**
//! `ν_t` (m²/s) that augments the molecular viscosity in the momentum
//! equations, plus the transport equations for `k`, `ε`, `ω` where used.
//!
//! # Closures
//!
//! - **Smagorinsky** (LES): `ν_t = (C_s · Δ)² · |S̄|` where `|S̄|` is the
//!   magnitude of the resolved strain rate. Zero-equation SGS model.
//! - **k-ε** (Launder & Spalding 1974): `ν_t = C_μ · k² / ε`. Requires
//!   two transport equations. Standard for industrial CFD.
//! - **k-ω** (Wilcox 1988): `ν_t = k / ω`. Better wall behaviour than
//!   k-ε; used in aerospace and turbomachinery.
//!
//! # Constants
//!
//! Universal fit-to-experiment values (Pope, *Turbulent Flows* Table 10.1):
//! - `C_s = 0.17`, `C_μ = 0.09`, `σ_k = 1.0`, `σ_ε = 1.3`,
//!   `C_{ε1} = 1.44`, `C_{ε2} = 1.92`, `β* = 0.09`, `β = 3/40`.
//!
//! # References
//!
//! - Smagorinsky, "General circulation experiments with the primitive
//!   equations", Monthly Weather Review 91, 1963.
//! - Launder & Spalding, "The numerical computation of turbulent flows",
//!   Comp. Meth. Appl. Mech. Eng. 3, 1974 (standard k-ε).
//! - Wilcox, "Reassessment of the scale-determining equation for advanced
//!   turbulence models", AIAA J. 26, 1988 (k-ω).
//! - Pope, *Turbulent Flows*, Cambridge 2000.
//!
//! # Integration status
//!
//! Wired into `cfd_solver.rs`: the Smagorinsky path
//! (`smagorinsky_eddy_viscosity`, `strain_rate_magnitude`, `SMAGORINSKY_CS`)
//! through `use_turbulence`, and the wall functions — [`friction_velocity`]
//! (the inverse of the universal profile, which is what a wall model needs
//! and what this module lacked until 2026-10-02), `y_plus`, `u_plus` and
//! `wall_k_epsilon` — through `CfdSolver::step_with_options` with a
//! `WallModel`. The k-ε / k-ω transport (`KEpsilonState` /
//! `KOmegaState`) and the dynamic Smagorinsky coefficient run in
//! `CfdSolver::step_rans`; with a `WallModel` that step also imposes
//! `wall_k_epsilon` on the wall-adjacent cells (the Launder–Spalding wall
//! function), without one `k` and `ε` have no wall boundary condition.

use crate::math::Fix128;

// ============================================================================
// Constants
// ============================================================================

/// Smagorinsky constant `C_s`.
pub const SMAGORINSKY_CS: Fix128 = Fix128 {
    hi: 0,
    lo: 0x2B85_1EB8_51EB_851F, // ≈ 0.17
};

#[cfg(feature = "std")]
/// k-ε model constants (Launder & Spalding 1974, crate-internal).
pub(crate) const KE_C_MU: Fix128 = Fix128 {
    hi: 0,
    lo: 0x170A3D70A3D70A3D, // 0.09 to 2⁻⁶⁴ (= ⌊9·2⁶⁴/100⌋, what `from_ratio(9, 100)` yields) (the earlier literal 0x1707_5F6F_D21F_F2E5 was 0.08995625, caught by the
                            // wall-model oracle on `k = u_τ²/√C_μ`; the unit test below tolerated 1e-2)
};
#[cfg(feature = "std")]
/// k-ε turbulent Prandtl number for k (crate-internal).
pub(crate) const KE_SIGMA_K: Fix128 = Fix128::ONE;
#[cfg(feature = "std")]
/// k-ε turbulent Prandtl number for ε (≈ 1.3, crate-internal).
pub(crate) const KE_SIGMA_EPS: Fix128 = Fix128 {
    hi: 1,
    lo: 0x4CCC_CCCC_CCCC_CCCD,
};
#[cfg(feature = "std")]
/// k-ε constant C_{ε1} (≈ 1.44, crate-internal).
pub(crate) const KE_C1_EPS: Fix128 = Fix128 {
    hi: 1,
    lo: 0x70A3_D70A_3D70_A3D7,
};
#[cfg(feature = "std")]
/// k-ε constant C_{ε2} (≈ 1.92, crate-internal).
pub(crate) const KE_C2_EPS: Fix128 = Fix128 {
    hi: 1,
    lo: 0xEB85_1EB8_51EB_851F,
};

#[cfg(feature = "std")]
/// k-ω model constant β* (= 0.09, crate-internal).
pub(crate) const KW_BETA_STAR: Fix128 = KE_C_MU;
#[cfg(feature = "std")]
/// k-ω model constant β (= 3/40 = 0.075, crate-internal).
pub(crate) const KW_BETA: Fix128 = Fix128 {
    hi: 0,
    lo: 0x1333_3333_3333_3333,
};
#[cfg(feature = "std")]
/// k-ω model constant α (= 5/9, Wilcox 1988, crate-internal): the production
/// coefficient of the ω equation, `dω/dt = α (ω/k) P_k − β ω²`.
pub(crate) const KW_ALPHA: Fix128 = Fix128 {
    hi: 0,
    lo: 0x8E38_E38E_38E3_8E38,
};
#[cfg(feature = "std")]
/// k-ω turbulent Prandtl numbers σ_k = σ_ω = 2 (Wilcox 1988, crate-internal).
pub(crate) const KW_SIGMA: Fix128 = Fix128::from_int(2);

// ============================================================================
// Smagorinsky
// ============================================================================

/// Smagorinsky eddy viscosity: `ν_t = (C_s · Δ)² · |S̄|`.
///
/// - `filter_width_m`: cell size Δ.
/// - `strain_rate_magnitude`: `|S̄| = √(2·S_ij·S_ij)`.
#[must_use]
pub fn smagorinsky_eddy_viscosity(filter_width_m: Fix128, strain_rate_magnitude: Fix128) -> Fix128 {
    smagorinsky_eddy_viscosity_with(SMAGORINSKY_CS, filter_width_m, strain_rate_magnitude)
}

/// [`smagorinsky_eddy_viscosity`] with the coefficient `cs` supplied, for the
/// dynamic procedure ([`dynamic_smagorinsky_cs`]); with `SMAGORINSKY_CS` it
/// is the static form bit for bit.
#[must_use]
pub(crate) fn smagorinsky_eddy_viscosity_with(
    cs: Fix128,
    filter_width_m: Fix128,
    strain_rate_magnitude: Fix128,
) -> Fix128 {
    let cs_delta = cs * filter_width_m;
    cs_delta * cs_delta * strain_rate_magnitude
}

/// Strain-rate magnitude from a symmetric 3×3 strain tensor S_ij:
/// `|S̄| = √(2·(S_11² + S_22² + S_33² + 2·(S_12² + S_13² + S_23²)))`.
#[must_use]
pub fn strain_rate_magnitude(
    s11: Fix128,
    s22: Fix128,
    s33: Fix128,
    s12: Fix128,
    s13: Fix128,
    s23: Fix128,
) -> Fix128 {
    let diag = s11 * s11 + s22 * s22 + s33 * s33;
    let off = Fix128::from_int(2) * (s12 * s12 + s13 * s13 + s23 * s23);
    (Fix128::from_int(2) * (diag + off)).sqrt()
}

// ============================================================================
// k-ε
// ============================================================================

#[cfg(feature = "std")]
/// k-ε turbulence state at a single point (crate-internal; the RANS step of
/// `cfd_solver` holds one per cell).
#[derive(Clone, Copy, Debug, PartialEq, Eq, Default)]
pub(crate) struct KEpsilonState {
    /// Turbulent kinetic energy `k` (m²/s²).
    pub(crate) k: Fix128,
    /// Turbulent dissipation rate `ε` (m²/s³).
    pub(crate) epsilon: Fix128,
}

#[cfg(feature = "std")]
impl KEpsilonState {
    /// Eddy viscosity `ν_t = C_μ · k² / ε` (m²/s).
    #[must_use]
    pub(crate) fn eddy_viscosity(&self) -> Fix128 {
        if self.epsilon.is_zero() {
            return Fix128::ZERO;
        }
        KE_C_MU * self.k * self.k / self.epsilon
    }

    /// One explicit Euler step of the k transport equation with production
    /// term `P_k` (m²/s³) and no diffusion or convection (point model):
    ///
    /// `dk/dt = P_k − ε`
    pub(crate) fn advance_k(&mut self, production_k: Fix128, dt_s: Fix128) {
        self.k = self.k + (production_k - self.epsilon) * dt_s;
        if self.k < Fix128::ZERO {
            self.k = Fix128::ZERO;
        }
    }

    /// One explicit Euler step of both equations from the state at the start of
    /// the step: `k` and `ε` are each advanced with the other's old value, so
    /// the pair is the explicit Euler update of the coupled system (advancing
    /// `k` first and letting `ε` read the new `k` is a different, semi-implicit
    /// scheme).
    pub(crate) fn advance(&mut self, production_k: Fix128, dt_s: Fix128) {
        let start = *self;
        self.advance_k(production_k, dt_s);
        let mut eps = start;
        eps.advance_epsilon(production_k, dt_s);
        self.epsilon = eps.epsilon;
    }

    /// One explicit Euler step of the ε transport equation (point model):
    ///
    /// `dε/dt = (ε/k) · (C_{ε1}·P_k − C_{ε2}·ε)`
    pub(crate) fn advance_epsilon(&mut self, production_k: Fix128, dt_s: Fix128) {
        if self.k.is_zero() {
            return;
        }
        let factor = self.epsilon / self.k;
        let d_eps = factor * (KE_C1_EPS * production_k - KE_C2_EPS * self.epsilon);
        self.epsilon = self.epsilon + d_eps * dt_s;
        if self.epsilon < Fix128::ZERO {
            self.epsilon = Fix128::ZERO;
        }
    }
}

// ============================================================================
// k-ω
// ============================================================================

#[cfg(feature = "std")]
/// k-ω turbulence state at a single point (crate-internal; the RANS step of
/// `cfd_solver` converts each cell from and back to its `(k, ε)` storage).
#[derive(Clone, Copy, Debug, PartialEq, Eq, Default)]
pub(crate) struct KOmegaState {
    /// Turbulent kinetic energy `k` (m²/s²).
    pub(crate) k: Fix128,
    /// Specific dissipation rate `ω = ε/(β*·k)` (1/s).
    pub(crate) omega: Fix128,
}

#[cfg(feature = "std")]
impl KOmegaState {
    /// Eddy viscosity `ν_t = k / ω` (m²/s).
    #[must_use]
    pub(crate) fn eddy_viscosity(&self) -> Fix128 {
        if self.omega.is_zero() {
            return Fix128::ZERO;
        }
        self.k / self.omega
    }

    /// Convert k-ε state to k-ω via `ω = ε / (β* · k)`.
    #[must_use]
    pub(crate) fn from_k_epsilon(state: &KEpsilonState) -> Self {
        let omega = if state.k.is_zero() {
            Fix128::ZERO
        } else {
            state.epsilon / (KW_BETA_STAR * state.k)
        };
        Self { k: state.k, omega }
    }

    /// Convert back to k-ε via `ε = β* · k · ω`.
    #[must_use]
    pub(crate) fn to_k_epsilon(self) -> KEpsilonState {
        KEpsilonState {
            k: self.k,
            epsilon: KW_BETA_STAR * self.k * self.omega,
        }
    }

    /// One explicit Euler step of the Wilcox (1988) point model with
    /// production `P_k` (m²/s³), no diffusion or convection:
    ///
    /// `dk/dt = P_k − β* k ω`, `dω/dt = α (ω/k) P_k − β ω²`
    ///
    /// Both use the values at the start of the step. A negative result is
    /// clamped to zero, as the k-ε point model does; `k = 0` leaves `ω`
    /// without a production term (the `ω/k` factor is taken as zero).
    pub(crate) fn advance(&mut self, production_k: Fix128, dt_s: Fix128) {
        let (k, w) = (self.k, self.omega);
        let dk = production_k - KW_BETA_STAR * k * w;
        let dw = if k.is_zero() {
            Fix128::ZERO - KW_BETA * w * w
        } else {
            KW_ALPHA * (w / k) * production_k - KW_BETA * w * w
        };
        self.k = k + dk * dt_s;
        if self.k < Fix128::ZERO {
            self.k = Fix128::ZERO;
        }
        self.omega = w + dw * dt_s;
        if self.omega < Fix128::ZERO {
            self.omega = Fix128::ZERO;
        }
    }
}

// ============================================================================
// Dynamic Smagorinsky
// ============================================================================

#[cfg(feature = "std")]
/// Compute a **dynamic** Smagorinsky coefficient `C_s(x,t)` from local flow
/// state using a Germano-style estimator. Returns a clamped value in
/// `[0.05, 0.25]` to avoid unphysical excursions.
///
/// Simplified form suitable for cell-local flows without a full test-filter
/// grid: use the ratio of test-filter strain to grid-filter strain as a
/// proxy for the Germano identity. `strain_grid` is the resolved-scale
/// magnitude, `strain_test` is an approximate `2·Δ` filter (average of
/// neighbouring cells; supplied by caller).
///
/// `C_s² = ⟨L·M⟩ / (2·⟨M·M⟩ + ε)` where `L` and `M` are the Leonard /
/// mixed tensors; here we approximate the isotropic invariants.
#[must_use]
pub(crate) fn dynamic_smagorinsky_cs(strain_grid: Fix128, strain_test: Fix128) -> Fix128 {
    let min_cs = Fix128::from_ratio(5, 100);
    let max_cs = Fix128::from_ratio(25, 100);
    if strain_grid.is_zero() || strain_test.is_zero() {
        return SMAGORINSKY_CS;
    }
    // ratio ≈ (strain_test / strain_grid); high ratio → less SGS energy
    let ratio = strain_test / strain_grid;
    // Simple mapping: C_s ∝ ratio, clamped
    let unclamped = SMAGORINSKY_CS * ratio;
    if unclamped < min_cs {
        min_cs
    } else if unclamped > max_cs {
        max_cs
    } else {
        unclamped
    }
}

// ============================================================================
// Wall functions
// ============================================================================

/// von Kármán constant κ (≈ 0.41) for log-law wall function (crate-internal).
pub(crate) const VON_KARMAN: Fix128 = Fix128 {
    hi: 0,
    lo: 0x68F5_C28F_5C28_F5C3,
};
/// Additive log-law constant B (≈ 5.5) for smooth walls (crate-internal).
pub(crate) const LOG_LAW_B: Fix128 = Fix128 {
    hi: 5,
    lo: 0x8000_0000_0000_0000,
};

/// `y⁺` at which the viscous sublayer `u⁺ = y⁺` meets the log law
/// `u⁺ = ln(y⁺)/κ + B`: the root of `y = ln(y)/0.41 + 5.5`, which is
/// `11.445319…` (crate-internal; the fraction below is that root to 2⁻⁶⁴).
///
/// ⚠️ Placing the switch **at the intersection** is what makes `u⁺(y⁺)`
/// continuous. The earlier value `11` (a figure that belongs to `B ≈ 5.0`)
/// left a step of `0.35` between the two branches, and a step in `u⁺` makes
/// the map `u_τ ↦ u_τ · u⁺(y u_τ/ν)` discontinuous, so a wall velocity falling
/// inside the step has no friction velocity at all. `tests/analytic_wall_model.rs`
/// solves for the root independently and pins both the constant and the
/// continuity.
pub(crate) const Y_PLUS_TRANSITION: Fix128 = Fix128 {
    hi: 11,
    lo: 0x72006EE534E20000, // 0.445319110 · 2⁶⁴
};

/// Bisection steps of [`friction_velocity`]: a fixed count, never a
/// convergence `break`, so the result is bit-reproducible across platforms
/// (the lockstep discipline — a data-dependent iteration count is one of the
/// five ways determinism is lost). 64 halvings of a bracket of order `u_rel`
/// bring the width below the arithmetic's resolution at any physical speed.
pub const FRICTION_VELOCITY_BISECTIONS: u32 = 64;

/// Why [`friction_velocity`] refused.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub enum WallFunctionError {
    /// The wall distance is not strictly positive.
    NonPositiveWallDistance,
    /// The kinematic viscosity is not strictly positive: `y⁺` is undefined.
    NonPositiveViscosity,
    /// The tangential speed relative to the wall is negative; pass its
    /// magnitude and keep the sign outside.
    NegativeSpeed,
}

impl core::fmt::Display for WallFunctionError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            Self::NonPositiveWallDistance => write!(f, "the wall distance must be positive"),
            Self::NonPositiveViscosity => write!(f, "the kinematic viscosity must be positive"),
            Self::NegativeSpeed => write!(f, "the speed relative to the wall must not be negative"),
        }
    }
}

#[cfg(feature = "std")]
impl std::error::Error for WallFunctionError {}

/// Friction velocity `u_τ` from the tangential speed `u_rel ≥ 0` at distance
/// `y_p` from a wall, in a fluid of kinematic viscosity `ν`: the `u_τ` with
///
/// ```text
/// u_rel = u_τ · u⁺(y⁺),   y⁺ = y_p · u_τ / ν
/// ```
///
/// using the universal profile of `u_plus` (this module). This is the *inverse* the wall
/// functions in this module were missing: `y_plus`, `u_plus` and
/// `wall_k_epsilon` all consume `u_τ`, and nothing produced it, so no wall
/// shear could be modelled from a resolved velocity.
///
/// The map `u_τ ↦ u_τ·u⁺(y_p u_τ/ν)` is strictly increasing and continuous
/// (the two branches of `u⁺` meet at `Y_PLUS_TRANSITION`, `11.4453`), so the root is
/// unique; it is found by [`FRICTION_VELOCITY_BISECTIONS`] halvings of the
/// bracket `[0, max(u_rel, √(ν u_rel / y_p))]`. The upper end covers both
/// regimes: in the sublayer the exact root is `√(ν u_rel / y_p)`, which
/// exceeds `u_rel` when `y⁺ < 1`; in the log region `u⁺ > 1` so `u_τ < u_rel`.
///
/// `u_rel = 0` returns exactly zero (the bracket is `[0, 0]`).
///
/// # Errors
///
/// [`WallFunctionError`] for a non-positive `y_p` or `ν`, or a negative
/// `u_rel`.
pub fn friction_velocity(
    u_rel: Fix128,
    y_p: Fix128,
    nu: Fix128,
) -> Result<Fix128, WallFunctionError> {
    if y_p <= Fix128::ZERO {
        return Err(WallFunctionError::NonPositiveWallDistance);
    }
    if nu <= Fix128::ZERO {
        return Err(WallFunctionError::NonPositiveViscosity);
    }
    if u_rel.is_negative() {
        return Err(WallFunctionError::NegativeSpeed);
    }
    Ok(friction_velocity_checked(u_rel, y_p, nu))
}

/// [`friction_velocity`] after its inputs have been validated.
pub(crate) fn friction_velocity_checked(u_rel: Fix128, y_p: Fix128, nu: Fix128) -> Fix128 {
    if u_rel.is_zero() {
        return Fix128::ZERO;
    }
    let sublayer_root = (nu * u_rel / y_p).sqrt();
    let mut lo = Fix128::ZERO;
    let mut hi = if sublayer_root > u_rel {
        sublayer_root
    } else {
        u_rel
    };
    // The bracket's upper end must be at or above the root: `hi · u⁺(y⁺(hi))
    // ≥ u_rel`. In the sublayer `u⁺ = y⁺` gives exactly `u_rel` at the
    // sublayer root; in the log region `u⁺ ≥ y⁺*` > 1 at and above the
    // transition. Widen once if the arithmetic put it below.
    if hi * u_plus(y_p * hi / nu) < u_rel {
        hi = hi.double();
    }
    for _ in 0..FRICTION_VELOCITY_BISECTIONS {
        let mid = (lo + hi).half();
        if mid * u_plus(y_p * mid / nu) < u_rel {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    (lo + hi).half()
}

#[cfg(feature = "std")]
/// Non-dimensional wall distance `y+ = ρ·u_τ·y / μ`.
#[must_use]
pub(crate) fn y_plus(
    density: Fix128,
    u_tau: Fix128,
    y_m: Fix128,
    dynamic_viscosity_pas: Fix128,
) -> Fix128 {
    if dynamic_viscosity_pas.is_zero() {
        return Fix128::ZERO;
    }
    density * u_tau * y_m / dynamic_viscosity_pas
}

/// Dimensionless velocity `u+ = u / u_τ` from the universal wall profile.
/// - `y+ < 11.63` → `u+ = y+` (viscous sublayer).
/// - Otherwise    → `u+ = (1/κ)·ln(y+) + B`  (log law).
#[must_use]
pub(crate) fn u_plus(y_plus_val: Fix128) -> Fix128 {
    if y_plus_val <= Fix128::ZERO {
        return Fix128::ZERO;
    }
    if y_plus_val < Y_PLUS_TRANSITION {
        return y_plus_val;
    }
    let ln_yp = y_plus_val.ln();
    ln_yp / VON_KARMAN + LOG_LAW_B
}

#[cfg(feature = "std")]
/// Wall-consistent turbulent kinetic energy and dissipation in the log
/// region. Given friction velocity `u_τ` (m/s) and wall distance `y` (m):
///
///   k = u_τ² / √C_μ
///   ε = u_τ³ / (κ · y)
#[must_use]
pub(crate) fn wall_k_epsilon(u_tau_m_per_s: Fix128, y_m: Fix128) -> (Fix128, Fix128) {
    if y_m.is_zero() || u_tau_m_per_s.is_zero() {
        return (Fix128::ZERO, Fix128::ZERO);
    }
    let sqrt_cmu = KE_C_MU.sqrt();
    let k = if sqrt_cmu.is_zero() {
        Fix128::ZERO
    } else {
        u_tau_m_per_s * u_tau_m_per_s / sqrt_cmu
    };
    let u_tau_cube = u_tau_m_per_s * u_tau_m_per_s * u_tau_m_per_s;
    let eps = u_tau_cube / (VON_KARMAN * y_m);
    (k, eps)
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    fn approx_eq(a: Fix128, b: Fix128, tol: Fix128) -> bool {
        let d = if a > b { a - b } else { b - a };
        d <= tol
    }

    #[test]
    fn smagorinsky_cs_constant() {
        // Should be ≈ 0.17
        assert!(approx_eq(
            SMAGORINSKY_CS,
            Fix128::from_ratio(17, 100),
            Fix128::from_ratio(1, 1000)
        ));
    }

    #[test]
    fn ke_c_mu_constant() {
        // To 2⁻⁶⁰, not to 1e-2: a tolerance of 1e-2 let a literal worth
        // 0.08995625 through for a long time.
        assert!(approx_eq(
            KE_C_MU,
            Fix128::from_ratio(9, 100),
            Fix128::from_raw(0, 16)
        ));
    }

    #[test]
    fn smagorinsky_zero_strain_zero_viscosity() {
        let nu = smagorinsky_eddy_viscosity(Fix128::from_ratio(1, 100), Fix128::ZERO);
        assert_eq!(nu, Fix128::ZERO);
    }

    #[test]
    fn smagorinsky_scales_with_delta_squared() {
        let s = Fix128::from_int(10);
        let nu1 = smagorinsky_eddy_viscosity(Fix128::from_ratio(1, 100), s);
        let nu2 = smagorinsky_eddy_viscosity(Fix128::from_ratio(2, 100), s);
        // Doubling Δ → 4× ν_t
        let ratio = nu2 / nu1;
        assert!(approx_eq(
            ratio,
            Fix128::from_int(4),
            Fix128::from_ratio(1, 100)
        ));
    }

    #[test]
    fn strain_rate_isotropic_diagonal() {
        // Pure diagonal S_11 = 1: |S̄| = √(2·1) = √2
        let s = strain_rate_magnitude(
            Fix128::ONE,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
        );
        let expected = Fix128::from_int(2).sqrt();
        assert!(approx_eq(s, expected, Fix128::from_ratio(1, 100)));
    }

    #[test]
    fn strain_rate_zero_tensor() {
        let s = strain_rate_magnitude(
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
            Fix128::ZERO,
        );
        assert_eq!(s, Fix128::ZERO);
    }

    #[test]
    fn ke_eddy_viscosity_formula() {
        let st = KEpsilonState {
            k: Fix128::ONE,
            epsilon: Fix128::from_ratio(9, 100),
        };
        // ν_t = 0.09 · 1 / 0.09 = 1
        let nu = st.eddy_viscosity();
        assert!(approx_eq(nu, Fix128::ONE, Fix128::from_ratio(1, 100)));
    }

    #[test]
    fn ke_zero_epsilon_zero_viscosity() {
        let st = KEpsilonState {
            k: Fix128::ONE,
            epsilon: Fix128::ZERO,
        };
        assert_eq!(st.eddy_viscosity(), Fix128::ZERO);
    }

    #[test]
    fn ke_k_transport_reduces_k_no_production() {
        let mut st = KEpsilonState {
            k: Fix128::from_int(10),
            epsilon: Fix128::from_int(1),
        };
        // No production, ε=1 → dk = -1·dt = -0.1 at dt=0.1
        st.advance_k(Fix128::ZERO, Fix128::from_ratio(1, 10));
        assert!(st.k < Fix128::from_int(10));
        assert!(st.k > Fix128::from_int(9));
    }

    #[test]
    fn ke_k_never_goes_negative() {
        let mut st = KEpsilonState {
            k: Fix128::ONE,
            epsilon: Fix128::from_int(100),
        };
        st.advance_k(Fix128::ZERO, Fix128::ONE);
        assert!(st.k >= Fix128::ZERO);
    }

    #[test]
    fn ke_epsilon_transport_zero_k_stays() {
        let mut st = KEpsilonState {
            k: Fix128::ZERO,
            epsilon: Fix128::from_int(1),
        };
        let eps_before = st.epsilon;
        st.advance_epsilon(Fix128::ONE, Fix128::ONE);
        assert_eq!(st.epsilon, eps_before);
    }

    #[test]
    fn kw_from_ke_conversion() {
        // ε = 0.09, k = 1 → ω = 0.09 / (0.09·1) = 1
        let ke = KEpsilonState {
            k: Fix128::ONE,
            epsilon: Fix128::from_ratio(9, 100),
        };
        let kw = KOmegaState::from_k_epsilon(&ke);
        assert_eq!(kw.k, ke.k);
        assert!(approx_eq(kw.omega, Fix128::ONE, Fix128::from_ratio(1, 100)));
    }

    #[test]
    fn kw_eddy_viscosity_reduces_when_omega_high() {
        // Higher ω → smaller ν_t
        let low = KOmegaState {
            k: Fix128::ONE,
            omega: Fix128::ONE,
        };
        let high = KOmegaState {
            k: Fix128::ONE,
            omega: Fix128::from_int(10),
        };
        assert!(low.eddy_viscosity() > high.eddy_viscosity());
    }

    #[test]
    fn kw_zero_omega_zero_viscosity() {
        let st = KOmegaState {
            k: Fix128::ONE,
            omega: Fix128::ZERO,
        };
        assert_eq!(st.eddy_viscosity(), Fix128::ZERO);
    }

    #[test]
    fn von_karman_and_b_constants() {
        assert!(approx_eq(
            VON_KARMAN,
            Fix128::from_ratio(41, 100),
            Fix128::from_ratio(1, 100)
        ));
        assert!(approx_eq(
            LOG_LAW_B,
            Fix128::from_ratio(55, 10),
            Fix128::from_ratio(1, 100)
        ));
    }

    #[test]
    fn y_plus_zero_viscosity_returns_zero() {
        let v = y_plus(
            Fix128::from_int(1000),
            Fix128::ONE,
            Fix128::from_ratio(1, 100),
            Fix128::ZERO,
        );
        assert_eq!(v, Fix128::ZERO);
    }

    #[test]
    fn u_plus_sublayer_linear() {
        // y+ = 5 → u+ = 5 (viscous sublayer)
        let up = u_plus(Fix128::from_int(5));
        assert_eq!(up, Fix128::from_int(5));
    }

    #[test]
    fn u_plus_log_region_increases_slowly() {
        // y+ = 100 → u+ = (1/0.41)·ln(100) + 5.5 ≈ 2.44·4.605 + 5.5 ≈ 16.7
        let up = u_plus(Fix128::from_int(100));
        assert!(approx_eq(
            up,
            Fix128::from_ratio(167, 10),
            Fix128::from_int(1)
        ));
    }

    #[test]
    fn u_plus_zero_returns_zero() {
        assert_eq!(u_plus(Fix128::ZERO), Fix128::ZERO);
    }

    #[test]
    fn wall_k_epsilon_positive() {
        let (k, eps) = wall_k_epsilon(Fix128::from_ratio(5, 10), Fix128::from_ratio(1, 100));
        assert!(k > Fix128::ZERO);
        assert!(eps > Fix128::ZERO);
    }

    #[test]
    fn dynamic_cs_default_when_strain_zero() {
        let c = dynamic_smagorinsky_cs(Fix128::ZERO, Fix128::from_int(1));
        assert_eq!(c, SMAGORINSKY_CS);
    }

    #[test]
    fn dynamic_cs_clamped_to_range() {
        // High ratio should clamp at 0.25
        let c = dynamic_smagorinsky_cs(Fix128::ONE, Fix128::from_int(100));
        assert!(c <= Fix128::from_ratio(25, 100));
        // Low ratio should clamp at 0.05
        let c2 = dynamic_smagorinsky_cs(Fix128::from_int(100), Fix128::ONE);
        assert!(c2 >= Fix128::from_ratio(5, 100));
    }

    #[test]
    fn wall_k_epsilon_zero_y_returns_zero() {
        let (k, eps) = wall_k_epsilon(Fix128::ONE, Fix128::ZERO);
        assert_eq!(k, Fix128::ZERO);
        assert_eq!(eps, Fix128::ZERO);
    }
}