//! Turbulence modeling for CFD
//!
//! This module provides various turbulence models for simulating turbulent flows:
//!
//! - **RANS Models**: k-ε (standard, RNG, realizable), k-ω, k-ω SST
//! - **LES Models**: Smagorinsky, Dynamic Smagorinsky, WALE
//! - **Wall Functions**: Log-law, enhanced wall treatment
//! - **Transition Models**: γ-Reθ, k-kL-ω

pub mod k_epsilon;
/// GPU-accelerated turbulence models
#[cfg(feature = "cuda")]
// PERF-2 (2026-09-16): the legacy k-epsilon GPU model no longer matches
// `KEpsilonConstants` and has not compiled since the initial commit; kept out
// of the `cuda` build until someone needs it.
// pub mod k_epsilon_gpu;
pub mod smagorinsky;
#[cfg(feature = "cuda")]
pub mod smagorinsky_gpu;
pub mod transition;
pub mod wall_functions;

pub use k_epsilon::{KEpsilonConstants, KEpsilonModel, KEpsilonVariant};
#[cfg(feature = "cuda")]
// pub use k_epsilon_gpu::KEpsilonGpuModel;
pub use smagorinsky::{SmagorinskyConstants, SmagorinskyModel};
pub use wall_functions::{EnhancedWallTreatment, LogLawWallFunction, WallFunction};

use crate::error::CfdResult;
use nalgebra::{DVector, Vector3};

/// Trait for turbulence models
pub trait TurbulenceModel {
    /// Calculate turbulent viscosity
    fn turbulent_viscosity(&self, state: &TurbulenceState) -> CfdResult<DVector<f64>>;

    /// Calculate production terms
    fn production_terms(&self, state: &TurbulenceState) -> CfdResult<TurbulenceProduction>;

    /// Update turbulence quantities
    fn update(&mut self, state: &TurbulenceState, dt: f64) -> CfdResult<()>;

    /// Get model name
    fn name(&self) -> &str;

    /// Check if model is RANS-based
    fn is_rans(&self) -> bool;

    /// Check if model is LES-based
    fn is_les(&self) -> bool;
}

/// Turbulence state containing flow variables
#[derive(Debug, Clone)]
pub struct TurbulenceState {
    /// Velocity field [u, v, w]
    pub velocity: Vec<Vector3<f64>>,
    /// Velocity gradients ∇u
    pub velocity_gradients: Vec<[[f64; 3]; 3]>,
    /// Pressure field
    pub pressure: DVector<f64>,
    /// Turbulent kinetic energy (for RANS models)
    pub turbulent_ke: Option<DVector<f64>>,
    /// Turbulent dissipation rate (for k-ε models)
    pub epsilon: Option<DVector<f64>>,
    /// Specific dissipation rate (for k-ω models)
    pub omega: Option<DVector<f64>>,
    /// Distance to wall
    pub wall_distance: DVector<f64>,
    /// Cell volumes
    pub cell_volumes: DVector<f64>,
    /// Dynamic viscosity
    pub molecular_viscosity: f64,
    /// Density
    pub density: f64,
}

impl TurbulenceState {
    /// Create new turbulence state
    #[must_use]
    pub fn new(n_cells: usize) -> Self {
        Self {
            velocity: vec![Vector3::zeros(); n_cells],
            velocity_gradients: vec![[[0.0; 3]; 3]; n_cells],
            pressure: DVector::zeros(n_cells),
            turbulent_ke: None,
            epsilon: None,
            omega: None,
            wall_distance: DVector::zeros(n_cells),
            cell_volumes: DVector::from_element(n_cells, 1.0),
            molecular_viscosity: 1e-5,
            density: 1.0,
        }
    }

    /// Initialize k-ε model variables
    pub fn initialize_k_epsilon(&mut self, k_init: f64, epsilon_init: f64) {
        let n_cells = self.velocity.len();
        self.turbulent_ke = Some(DVector::from_element(n_cells, k_init));
        self.epsilon = Some(DVector::from_element(n_cells, epsilon_init));
    }

    /// Initialize k-ω model variables
    pub fn initialize_k_omega(&mut self, k_init: f64, omega_init: f64) {
        let n_cells = self.velocity.len();
        self.turbulent_ke = Some(DVector::from_element(n_cells, k_init));
        self.omega = Some(DVector::from_element(n_cells, omega_init));
    }

    /// Calculate strain rate magnitude
    #[must_use]
    pub fn strain_rate_magnitude(&self, cell_idx: usize) -> f64 {
        if cell_idx >= self.velocity_gradients.len() {
            return 0.0;
        }

        let grad = &self.velocity_gradients[cell_idx];
        let mut s_mag = 0.0;

        // S_ij = 0.5 * (∂u_i/∂x_j + ∂u_j/∂x_i)
        // For strain rate magnitude: |S| = sqrt(2 * S_ij * S_ij)
        // Sum over all components (symmetric tensor)
        for i in 0..3 {
            for j in 0..3 {
                let s_ij = 0.5 * (grad[i][j] + grad[j][i]);
                s_mag += s_ij * s_ij;
            }
        }

        // The factor of 2 accounts for the definition |S| = sqrt(2*Sij*Sij)
        (2.0 * s_mag).sqrt()
    }

    /// Calculate vorticity magnitude
    #[must_use]
    pub fn vorticity_magnitude(&self, cell_idx: usize) -> f64 {
        if cell_idx >= self.velocity_gradients.len() {
            return 0.0;
        }

        let grad = &self.velocity_gradients[cell_idx];
        let mut omega_mag = 0.0;

        // Ω_ij = 0.5 * (∂u_i/∂x_j - ∂u_j/∂x_i)
        // For antisymmetric tensor, sum over all components
        for i in 0..3 {
            for j in 0..3 {
                let omega_ij = 0.5 * (grad[i][j] - grad[j][i]);
                omega_mag += omega_ij * omega_ij;
            }
        }

        // The factor of 2 accounts for the definition |Ω| = sqrt(2*Ωij*Ωij)
        (2.0 * omega_mag).sqrt()
    }
}

/// Production terms for turbulence models
#[derive(Debug, Clone)]
pub struct TurbulenceProduction {
    /// Production of turbulent kinetic energy
    pub pk: DVector<f64>,
    /// Production of dissipation
    pub pe: Option<DVector<f64>>,
    /// Production of specific dissipation
    pub pw: Option<DVector<f64>>,
}

impl TurbulenceProduction {
    /// Create new production terms
    #[must_use]
    pub fn new(n_cells: usize) -> Self {
        Self {
            pk: DVector::zeros(n_cells),
            pe: None,
            pw: None,
        }
    }

    /// Initialize for k-ε model
    #[must_use]
    pub fn for_k_epsilon(n_cells: usize) -> Self {
        Self {
            pk: DVector::zeros(n_cells),
            pe: Some(DVector::zeros(n_cells)),
            pw: None,
        }
    }

    /// Initialize for k-ω model
    #[must_use]
    pub fn for_k_omega(n_cells: usize) -> Self {
        Self {
            pk: DVector::zeros(n_cells),
            pe: None,
            pw: Some(DVector::zeros(n_cells)),
        }
    }
}

/// Turbulence intensity calculation utilities
pub struct TurbulenceIntensity;

impl TurbulenceIntensity {
    /// Calculate turbulence intensity from velocity fluctuations
    #[must_use]
    pub fn from_fluctuations(u_rms: f64, v_rms: f64, w_rms: f64, u_mean: f64) -> f64 {
        let turbulent_ke = 0.5 * (u_rms * u_rms + v_rms * v_rms + w_rms * w_rms);
        turbulent_ke.sqrt() / u_mean.abs().max(1e-10)
    }

    /// Estimate turbulence intensity for external flows
    #[must_use]
    pub fn external_flow_estimate(reynolds_number: f64) -> f64 {
        // Empirical correlation for external flows
        0.16 * reynolds_number.powf(-1.0 / 8.0)
    }

    /// Estimate turbulence intensity for internal flows
    #[must_use]
    pub fn internal_flow_estimate(reynolds_number: f64) -> f64 {
        // Empirical correlation for pipe flows
        0.16 * reynolds_number.powf(-1.0 / 8.0).min(0.1)
    }

    /// Calculate turbulent length scale
    #[must_use]
    pub fn turbulent_length_scale(characteristic_length: f64, turbulence_intensity: f64) -> f64 {
        // Empirical estimate
        0.07 * characteristic_length * (1.0 + 10.0 * turbulence_intensity)
    }
}

/// Reynolds number utilities for turbulence modeling
pub struct ReynoldsNumber;

impl ReynoldsNumber {
    /// Calculate turbulent Reynolds number
    #[must_use]
    pub fn turbulent(k: f64, epsilon: f64, nu: f64) -> f64 {
        k * k / (epsilon * nu).max(1e-15)
    }

    /// Calculate wall Reynolds number y+
    #[must_use]
    pub fn y_plus(y: f64, u_tau: f64, nu: f64) -> f64 {
        y * u_tau / nu
    }

    /// Calculate friction velocity
    #[must_use]
    pub fn friction_velocity(wall_shear_stress: f64, density: f64) -> f64 {
        (wall_shear_stress / density).sqrt()
    }

    /// Calculate wall shear stress from velocity gradient
    #[must_use]
    pub fn wall_shear_stress(du_dy: f64, mu: f64) -> f64 {
        mu * du_dy
    }
}

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

    #[test]
    fn test_turbulence_state_creation() {
        let state = TurbulenceState::new(10);
        assert_eq!(state.velocity.len(), 10);
        assert_eq!(state.pressure.len(), 10);
        assert!(state.turbulent_ke.is_none());
        assert!(state.epsilon.is_none());
    }

    #[test]
    fn test_turbulence_state_k_epsilon_init() {
        let mut state = TurbulenceState::new(5);
        state.initialize_k_epsilon(1.0, 0.1);

        assert!(state.turbulent_ke.is_some());
        assert!(state.epsilon.is_some());
        assert_eq!(state.turbulent_ke.as_ref().unwrap().len(), 5);
        assert_relative_eq!(
            state.turbulent_ke.as_ref().unwrap()[0],
            1.0,
            epsilon = 1e-10
        );
    }

    #[test]
    fn test_turbulence_state_k_omega_init() {
        let mut state = TurbulenceState::new(5);
        state.initialize_k_omega(1.0, 10.0);

        assert!(state.turbulent_ke.is_some());
        assert!(state.omega.is_some());
        assert_eq!(state.omega.as_ref().unwrap().len(), 5);
        assert_relative_eq!(state.omega.as_ref().unwrap()[0], 10.0, epsilon = 1e-10);
    }

    #[test]
    fn test_strain_rate_magnitude() {
        let mut state = TurbulenceState::new(1);
        // Set up simple shear flow: du/dy = 1, others = 0
        state.velocity_gradients[0][0][1] = 1.0; // du/dy = 1

        let s_mag = state.strain_rate_magnitude(0);
        // S_01 = S_10 = 0.5 * (1 + 0) = 0.5, all others = 0
        // Sum of S_ij^2 = 2 * 0.5^2 = 0.5
        // |S| = sqrt(2 * 0.5) = 1.0
        assert_relative_eq!(s_mag, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_vorticity_magnitude() {
        let mut state = TurbulenceState::new(1);
        // Set up rotation: du/dy = 1, dv/dx = -1
        state.velocity_gradients[0][0][1] = 1.0; // du/dy = 1
        state.velocity_gradients[0][1][0] = -1.0; // dv/dx = -1

        let omega_mag = state.vorticity_magnitude(0);
        // Ω_01 = 0.5 * (1 - (-1)) = 1, Ω_10 = 0.5 * ((-1) - 1) = -1
        // Sum of Ω_ij^2 = 1^2 + (-1)^2 = 2
        // |Ω| = sqrt(2 * 2) = 2
        assert_relative_eq!(omega_mag, 2.0, epsilon = 1e-10);
    }

    #[test]
    fn test_turbulence_production_creation() {
        let prod = TurbulenceProduction::new(5);
        assert_eq!(prod.pk.len(), 5);
        assert!(prod.pe.is_none());
        assert!(prod.pw.is_none());
    }

    #[test]
    fn test_turbulence_production_k_epsilon() {
        let prod = TurbulenceProduction::for_k_epsilon(5);
        assert_eq!(prod.pk.len(), 5);
        assert!(prod.pe.is_some());
        assert!(prod.pw.is_none());
        assert_eq!(prod.pe.as_ref().unwrap().len(), 5);
    }

    #[test]
    fn test_turbulence_production_k_omega() {
        let prod = TurbulenceProduction::for_k_omega(5);
        assert_eq!(prod.pk.len(), 5);
        assert!(prod.pe.is_none());
        assert!(prod.pw.is_some());
        assert_eq!(prod.pw.as_ref().unwrap().len(), 5);
    }

    #[test]
    fn test_turbulence_intensity_from_fluctuations() {
        let ti = TurbulenceIntensity::from_fluctuations(1.0, 1.0, 1.0, 10.0);
        // k = 0.5 * (1 + 1 + 1) = 1.5
        // TI = sqrt(k) / U = sqrt(1.5) / 10
        assert_relative_eq!(ti, 1.5_f64.sqrt() / 10.0, epsilon = 1e-10);
    }

    #[test]
    fn test_turbulence_intensity_estimates() {
        let re = 1e6;
        let ti_ext = TurbulenceIntensity::external_flow_estimate(re);
        let ti_int = TurbulenceIntensity::internal_flow_estimate(re);

        assert!(ti_ext > 0.0);
        assert!(ti_int > 0.0);
        assert!(ti_ext < 1.0);
        assert!(ti_int < 1.0);
    }

    #[test]
    fn test_turbulent_length_scale() {
        let l_scale = TurbulenceIntensity::turbulent_length_scale(1.0, 0.05);
        assert!(l_scale > 0.0);
    }

    #[test]
    fn test_reynolds_number_turbulent() {
        let re_t = ReynoldsNumber::turbulent(1.0, 1.0, 1e-6);
        assert_relative_eq!(re_t, 1e6, epsilon = 1e-10);
    }

    #[test]
    fn test_reynolds_number_y_plus() {
        let y_plus = ReynoldsNumber::y_plus(1e-5, 0.1, 1e-6);
        assert_relative_eq!(y_plus, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_friction_velocity() {
        let u_tau = ReynoldsNumber::friction_velocity(0.1, 1.0);
        assert_relative_eq!(u_tau, 0.1_f64.sqrt(), epsilon = 1e-10);
    }

    #[test]
    fn test_wall_shear_stress() {
        let tau_w = ReynoldsNumber::wall_shear_stress(100.0, 1e-3);
        assert_relative_eq!(tau_w, 0.1, epsilon = 1e-10);
    }
}