//! Incompressible flow solvers
//!
//! This module implements pressure-velocity coupling algorithms for incompressible flows:
//! - SIMPLE (Semi-Implicit Method for Pressure Linked Equations)
//! - PISO (Pressure-Implicit with Splitting of Operators)
//! - SIMPLER (SIMPLE Revised)

use crate::{CfdConfig, CfdError, CfdResult};
// use nalgebra::{DMatrix, DVector};
// use std::collections::HashMap;

/// ALE solver on a moving tensor-product staggered grid
pub mod ale;
/// Boundary conditions
pub mod boundary_conditions;
/// Collocated PISO on a structured curvilinear patch (the overset patch)
pub mod curvilinear;
/// PISO on the fixed grid with an embedded body
pub mod embedded;
pub mod embedded3;
/// Embedded-body geometry, classification and loads
pub mod embedded_body;
/// Flow field data structures
pub mod flow_field;
/// The overset hybrid: curvilinear patch over the fixed background
pub mod overset;
/// PISO algorithm implementation
pub mod piso;
/// GPU-accelerated PISO algorithm implementation
#[cfg(feature = "cuda")]
pub mod piso_gpu;
/// Five-point Poisson problems and the multigrid-preconditioned CG solver
pub mod poisson;
pub mod polygon_sdf;
/// SIMPLE algorithm implementation
pub mod simple;
/// GPU-accelerated SIMPLE algorithm implementation
#[cfg(feature = "cuda")]
pub mod simple_gpu;
/// CSR matrix + Jacobi-BiCGSTAB for the curvilinear pressure equation
pub mod sparse_bicgstab;
pub mod three_d;

// Re-export main types
pub use ale::{
    AleBoundaries, AleField, AleParameters, AlePisoSolver, AleResult, SideBoundary, SweptFaceRule,
};
pub use boundary_conditions::{
    BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
};
pub use curvilinear::{
    CurvilinearParameters, CurvilinearPisoSolver, CurvilinearResult, CurvilinearSolverState,
    NormalDiffusion, Operators, PatchBalance, PatchBoundaries, PatchConvection, PatchField,
    PatchLoad, RobinWall, SideBc, StepGeometry,
};
pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
pub use embedded_body::{
    EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_interface_velocity,
    polygon_signed_distance,
};
pub use flow_field::FlowField;
pub use overset::{
    CellClass, MomentumResidual, OverlapMap, OversetField, OversetParameters, OversetPisoSolver,
    OversetResult, OversetSolverState, ResidualBucket, StepTimers,
};
pub use piso::{PisoParameters, PisoResult, PisoSolver};
#[cfg(feature = "cuda")]
pub use piso_gpu::PisoGpuSolver;
pub use poisson::{
    LevelExport, MgPrecision, MgSmoother, MultigridParameters, PcgCache, PoissonProblem,
    PoissonSolution, PoissonSolverKind, configure_threads, export_hierarchy, plane_counters,
    plane_streams, set_plane_lane, set_plane_streams, solve_multigrid_pcg,
    solve_multigrid_pcg_cached, vcycle_f32_reference, vcycle_f32_work,
};
pub use polygon_sdf::PolygonSdf;
pub use simple::{ConvectionScheme, SimpleParameters, SimpleResult, SimpleSolver};
#[cfg(feature = "cuda")]
pub use simple_gpu::SimpleGpuSolver;

/// Common solver parameters
#[derive(Debug, Clone)]
pub struct SolverParameters {
    /// Maximum number of iterations
    pub max_iterations: usize,
    /// Convergence tolerance
    pub tolerance: f64,
    /// Time step size
    pub time_step: f64,
    /// Under-relaxation factors
    pub relaxation: RelaxationFactors,
}

/// Under-relaxation factors for stability
#[derive(Debug, Clone)]
pub struct RelaxationFactors {
    /// Pressure relaxation factor (typically 0.2-0.8)
    pub pressure: f64,
    /// Velocity relaxation factor (typically 0.5-0.8)
    pub velocity: f64,
    /// Turbulence relaxation factor
    pub turbulence: f64,
}

impl Default for SolverParameters {
    fn default() -> Self {
        Self {
            max_iterations: 1000,
            tolerance: 1e-6,
            time_step: 0.001,
            relaxation: RelaxationFactors::default(),
        }
    }
}

impl Default for RelaxationFactors {
    fn default() -> Self {
        Self {
            pressure: 0.3,
            velocity: 0.7,
            turbulence: 0.5,
        }
    }
}

/// Common solver result information
#[derive(Debug, Clone)]
pub struct SolverResult {
    /// Whether the solver converged
    pub converged: bool,
    /// Number of iterations performed
    pub iterations: usize,
    /// Final residual norm
    pub final_residual: f64,
    /// Residual history
    pub residual_history: Vec<f64>,
    /// Computational time
    pub solve_time: std::time::Duration,
}

/// Pressure-velocity coupling algorithms
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum CouplingAlgorithm {
    /// Semi-Implicit Method for Pressure Linked Equations
    Simple,
    /// Pressure-Implicit with Splitting of Operators
    Piso,
    /// SIMPLE Revised
    Simpler,
}

/// Common trait for incompressible solvers
#[async_trait::async_trait]
pub trait IncompressibleSolver {
    /// Solver-specific parameters
    type Parameters;
    /// Solver-specific result
    type Result;

    /// Create new solver instance
    fn new(config: CfdConfig, params: Self::Parameters) -> CfdResult<Self>
    where
        Self: Sized;

    /// Solve one time step
    async fn solve_time_step(
        &mut self,
        flow_field: &mut FlowField,
        boundary_conditions: &BoundaryConditions,
        dt: f64,
    ) -> CfdResult<Self::Result>;

    /// Solve to steady state
    async fn solve(
        &mut self,
        flow_field: &mut FlowField,
        boundary_conditions: &BoundaryConditions,
    ) -> CfdResult<Self::Result>;

    /// Get solver configuration
    fn config(&self) -> &CfdConfig;

    /// Get solver parameters
    fn parameters(&self) -> &Self::Parameters;
}

/// Utility functions for incompressible solvers
pub mod utils {
    use super::{CfdError, CfdResult};

    /// Compute Courant number
    #[must_use]
    pub fn compute_courant_number(u_max: f64, v_max: f64, dx: f64, dy: f64, dt: f64) -> f64 {
        let u_cfl = u_max * dt / dx;
        let v_cfl = v_max * dt / dy;
        (u_cfl * u_cfl + v_cfl * v_cfl).sqrt()
    }

    /// Compute viscous CFL number
    #[must_use]
    pub fn compute_viscous_cfl(nu: f64, dx: f64, dy: f64, dt: f64) -> f64 {
        nu * dt * (1.0 / (dx * dx) + 1.0 / (dy * dy))
    }

    /// Check stability criteria
    pub fn check_stability(courant: f64, viscous_cfl: f64) -> CfdResult<()> {
        if courant > 1.0 {
            return Err(CfdError::physics(format!(
                "Convective CFL condition violated: CFL = {courant:.3} > 1.0"
            )));
        }

        if viscous_cfl > 0.5 {
            return Err(CfdError::physics(format!(
                "Viscous CFL condition violated: CFL_visc = {viscous_cfl:.3} > 0.5"
            )));
        }

        Ok(())
    }

    /// Compute Reynolds number based on flow conditions
    #[must_use]
    pub fn compute_reynolds_number(
        u_characteristic: f64,
        length_characteristic: f64,
        kinematic_viscosity: f64,
    ) -> f64 {
        u_characteristic * length_characteristic / kinematic_viscosity
    }
}