//! Lattice Boltzmann Method simulation of Poiseuille flow
//!
//! Validates LBM implementation against analytical solution for pressure-driven
//! flow between parallel plates. The analytical velocity profile is:
//!
//! u(y) = (dp/dx) * y * (H - y) / (2 * μ)
//!
//! where H is the channel height.

use nalgebra::DVector;
use rtx_cfd::{
    CfdConfig, CfdResult,
    solvers::lbm::{D2Q9Parameters, D2Q9Solver},
};
use std::time::Instant;

/// Poiseuille flow simulation parameters
#[derive(Debug, Clone)]
pub struct PoiseuilleConfig {
    /// Channel length
    pub length: f64,
    /// Channel height
    pub height: f64,
    /// Pressure gradient (Pa/m)
    pub pressure_gradient: f64,
    /// Dynamic viscosity
    pub viscosity: f64,
    /// Density
    pub density: f64,
    /// Grid resolution [nx, ny]
    pub grid_size: [usize; 2],
    /// Lattice spacing
    pub dx: f64,
    /// Time step
    pub dt: f64,
    /// Simulation time
    pub total_time: f64,
    /// Output frequency
    pub output_frequency: usize,
    /// Relaxation time
    pub tau: f64,
}

impl PoiseuilleConfig {
    /// Create standard Poiseuille flow configuration
    pub fn standard() -> Self {
        let height = 1.0;
        let length = 10.0 * height;
        let pressure_gradient = -1000.0; // Pa/m
        let viscosity = 1e-3; // Water
        let density = 1000.0; // Water

        let ny = 64;
        let nx = ny * 10;
        let dx = height / (ny as f64 - 1.0);
        let dt = 1e-4;

        // Calculate relaxation time for LBM
        let cs2 = 1.0 / 3.0; // Sound speed squared in lattice units
        let nu_lattice = viscosity / (density * cs2);
        let tau = 3.0 * nu_lattice + 0.5;

        Self {
            length,
            height,
            pressure_gradient,
            viscosity,
            density,
            grid_size: [nx, ny],
            dx,
            dt,
            total_time: 5.0,
            output_frequency: 1000,
            tau,
        }
    }

    /// Create high Reynolds number case
    pub fn high_reynolds() -> Self {
        let mut config = Self::standard();
        config.pressure_gradient = -10000.0; // Higher pressure gradient
        config.viscosity = 1e-4; // Lower viscosity
        config.tau = 0.6; // Closer to stability limit
        config.total_time = 10.0;
        config
    }

    /// Calculate Reynolds number
    pub fn reynolds_number(&self) -> f64 {
        let u_max = self.max_velocity();
        self.density * u_max * self.height / self.viscosity
    }

    /// Calculate maximum velocity (analytical)
    pub fn max_velocity(&self) -> f64 {
        -self.pressure_gradient * self.height * self.height / (8.0 * self.viscosity)
    }

    /// Calculate analytical velocity profile
    pub fn analytical_velocity(&self, y: f64) -> f64 {
        let h = self.height;
        -self.pressure_gradient * y * (h - y) / (2.0 * self.viscosity)
    }

    /// Calculate volumetric flow rate (analytical)
    pub fn analytical_flow_rate(&self) -> f64 {
        let h = self.height;
        -self.pressure_gradient * h.powi(3) / (12.0 * self.viscosity)
    }
}

/// LBM Poiseuille flow simulation
pub struct LbmPoiseuille {
    /// Configuration
    config: PoiseuilleConfig,
    /// CFD configuration
    cfd_config: CfdConfig,
    /// LBM solver
    solver: D2Q9Solver,
    /// Current time
    current_time: f64,
    /// Velocity history for convergence check
    velocity_history: Vec<f64>,
}

impl LbmPoiseuille {
    /// Create new Poiseuille flow simulation
    pub fn new(config: PoiseuilleConfig) -> CfdResult<Self> {
        // Set up CFD configuration
        let cfd_config = CfdConfig::new()
            .with_density(config.density)
            .with_viscosity(config.viscosity)
            .with_reference_velocity(config.max_velocity())
            .with_reference_length(config.height);

        cfd_config.validate()?;

        println!("Setting up LBM Poiseuille flow simulation:");
        println!("  Channel: {:.2} × {:.2} m", config.length, config.height);
        println!("  Grid: {} × {}", config.grid_size[0], config.grid_size[1]);
        println!("  Pressure gradient: {:.0} Pa/m", config.pressure_gradient);
        println!("  Reynolds number: {:.1}", config.reynolds_number());
        println!(
            "  Max velocity (analytical): {:.4} m/s",
            config.max_velocity()
        );
        println!("  Relaxation time τ: {:.3}", config.tau);

        // Create LBM solver
        let nx = config.grid_size[0];
        let ny = config.grid_size[1];

        let params = D2Q9Parameters::new(config.tau);

        let solver = D2Q9Solver::new(nx, ny, params);

        // Note: Boundary conditions and initialization will be handled in the actual D2Q9Solver implementation

        Ok(Self {
            config,
            cfd_config,
            solver,
            current_time: 0.0,
            velocity_history: Vec::new(),
        })
    }

    // Note: Boundary conditions and initialization are now handled internally by D2Q9Solver
    // These methods are kept as placeholders for documentation but are not used

    /// Run simulation
    pub fn run_simulation(&mut self) -> CfdResult<PoiseuilleResults> {
        println!("\nStarting LBM simulation...");
        println!("  Time step: {:.2e} s", self.config.dt);
        println!("  Total time: {:.2} s", self.config.total_time);

        let start_time = Instant::now();
        let mut time_step = 0;
        let mut velocity_errors = Vec::new();

        while self.current_time < self.config.total_time {
            // Perform LBM time step
            self.solver.step();

            self.current_time += self.config.dt;
            time_step += 1;

            // Check convergence and collect data
            if time_step % self.config.output_frequency == 0 {
                let max_velocity = self.calculate_max_velocity()?;
                self.velocity_history.push(max_velocity);

                // Calculate error against analytical solution
                let velocity_error = self.calculate_velocity_error()?;
                velocity_errors.push(velocity_error);

                println!(
                    "  Time: {:.3}, Step: {}, Max velocity: {:.4}, Error: {:.2e}",
                    self.current_time, time_step, max_velocity, velocity_error
                );

                // Check convergence
                if self.is_converged() {
                    println!("  Simulation converged!");
                    break;
                }
            }
        }

        let elapsed = start_time.elapsed();
        println!(
            "LBM simulation completed in {:.2} seconds",
            elapsed.as_secs_f64()
        );

        // Calculate final results
        let final_velocity_profile = self.extract_velocity_profile()?;
        let final_flow_rate = self.calculate_flow_rate()?;
        let final_error = velocity_errors.last().copied().unwrap_or(1.0);

        Ok(PoiseuilleResults {
            time_steps: time_step,
            final_time: self.current_time,
            elapsed_time: elapsed,
            velocity_profile: final_velocity_profile,
            velocity_history: self.velocity_history.clone(),
            velocity_errors,
            flow_rate: final_flow_rate,
            analytical_flow_rate: self.config.analytical_flow_rate(),
            max_velocity: self.config.max_velocity(),
            final_error,
            converged: self.is_converged(),
        })
    }

    /// Calculate maximum velocity in the channel
    fn calculate_max_velocity(&self) -> CfdResult<f64> {
        let nx = self.config.grid_size[0];
        let ny = self.config.grid_size[1];
        let mut max_vel = 0.0_f64;

        for j in 0..ny {
            for i in 0..nx {
                let vars = self.solver.macroscopic_variables_at(i, j);
                max_vel = max_vel.max(vars.velocity.x.abs());
            }
        }

        Ok(max_vel)
    }

    /// Calculate L2 error against analytical solution
    fn calculate_velocity_error(&self) -> CfdResult<f64> {
        let nx = self.config.grid_size[0];
        let ny = self.config.grid_size[1];
        let mut error_sum = 0.0_f64;
        let mut count = 0;

        // Sample at channel center
        let i = nx / 2;

        for j in 1..ny - 1 {
            // Exclude walls
            let y = j as f64 * self.config.dx;
            let vars = self.solver.macroscopic_variables_at(i, j);
            let u_lbm = vars.velocity.x;
            let u_analytical = self.config.analytical_velocity(y);

            error_sum += (u_lbm - u_analytical).powi(2);
            count += 1;
        }

        Ok((error_sum / count as f64).sqrt())
    }

    /// Extract velocity profile at channel center
    fn extract_velocity_profile(&self) -> CfdResult<Vec<(f64, f64)>> {
        let nx = self.config.grid_size[0];
        let ny = self.config.grid_size[1];
        let mut profile = Vec::new();

        let i = nx / 2; // Channel center

        for j in 0..ny {
            let y = j as f64 * self.config.dx;
            let vars = self.solver.macroscopic_variables_at(i, j);
            profile.push((y, vars.velocity.x));
        }

        Ok(profile)
    }

    /// Calculate volumetric flow rate
    fn calculate_flow_rate(&self) -> CfdResult<f64> {
        let nx = self.config.grid_size[0];
        let ny = self.config.grid_size[1];
        let mut flow_rate = 0.0;

        let i = nx / 2; // Channel center

        for j in 1..ny - 1 {
            // Exclude walls
            let vars = self.solver.macroscopic_variables_at(i, j);
            flow_rate += vars.density * vars.velocity.x * self.config.dx; // Per unit depth
        }

        Ok(flow_rate)
    }

    /// Check if simulation has converged
    fn is_converged(&self) -> bool {
        if self.velocity_history.len() < 10 {
            return false;
        }

        // Check relative change in maximum velocity
        let recent = &self.velocity_history[self.velocity_history.len() - 5..];
        let avg_recent = recent.iter().sum::<f64>() / recent.len() as f64;

        let prev = &self.velocity_history
            [self.velocity_history.len() - 10..self.velocity_history.len() - 5];
        let avg_prev = prev.iter().sum::<f64>() / prev.len() as f64;

        let relative_change = (avg_recent - avg_prev).abs() / avg_prev.max(1e-10);
        relative_change < 1e-6
    }

    /// Get configuration
    pub fn config(&self) -> &PoiseuilleConfig {
        &self.config
    }
}

/// Poiseuille simulation results
#[derive(Debug, Clone)]
pub struct PoiseuilleResults {
    /// Number of time steps
    pub time_steps: usize,
    /// Final simulation time
    pub final_time: f64,
    /// Elapsed wall time
    pub elapsed_time: std::time::Duration,
    /// Final velocity profile
    pub velocity_profile: Vec<(f64, f64)>,
    /// Maximum velocity history
    pub velocity_history: Vec<f64>,
    /// Velocity error history
    pub velocity_errors: Vec<f64>,
    /// Calculated flow rate
    pub flow_rate: f64,
    /// Analytical flow rate
    pub analytical_flow_rate: f64,
    /// Maximum velocity
    pub max_velocity: f64,
    /// Final error
    pub final_error: f64,
    /// Convergence status
    pub converged: bool,
}

impl PoiseuilleResults {
    /// Print summary
    pub fn print_summary(&self) {
        println!("\n=== LBM Poiseuille Flow Results ===");
        println!("Time steps: {}", self.time_steps);
        println!("Final time: {:.3} s", self.final_time);
        println!("Elapsed time: {:.2} s", self.elapsed_time.as_secs_f64());
        println!("Converged: {}", self.converged);

        println!("\n--- Flow Validation ---");
        println!("LBM flow rate: {:.6} m²/s", self.flow_rate);
        println!(
            "Analytical flow rate: {:.6} m²/s",
            self.analytical_flow_rate
        );
        let flow_error =
            (self.flow_rate - self.analytical_flow_rate).abs() / self.analytical_flow_rate;
        println!("Flow rate error: {:.2}%", flow_error * 100.0);

        println!("Final velocity error: {:.2e}", self.final_error);

        if let Some(&max_vel) = self.velocity_history.last() {
            println!("Final max velocity: {:.4} m/s", max_vel);
            println!("Analytical max velocity: {:.4} m/s", self.max_velocity);
            let vel_error = (max_vel - self.max_velocity).abs() / self.max_velocity;
            println!("Max velocity error: {:.2}%", vel_error * 100.0);
        }
    }

    /// Calculate convergence rate
    pub fn convergence_rate(&self) -> Option<f64> {
        if self.velocity_errors.len() < 2 {
            return None;
        }

        // Fit exponential decay to error
        let n = self.velocity_errors.len();
        let dt = 1.0; // Normalized time step

        let mut sum_log_ratio = 0.0;
        let mut count = 0;

        for i in 1..n {
            if self.velocity_errors[i] > 0.0 && self.velocity_errors[i - 1] > 0.0 {
                sum_log_ratio += (self.velocity_errors[i] / self.velocity_errors[i - 1]).ln();
                count += 1;
            }
        }

        if count > 0 {
            Some(-sum_log_ratio / (count as f64 * dt))
        } else {
            None
        }
    }
}

/// Validation against analytical solution
pub fn validate_poiseuille_solution(
    config: &PoiseuilleConfig,
    results: &PoiseuilleResults,
) -> ValidationResults {
    let mut position_errors = Vec::new();
    let mut velocity_errors = Vec::new();

    // Compare velocity profile
    for &(y, u_lbm) in &results.velocity_profile {
        if y > 0.0 && y < config.height {
            // Exclude walls
            let u_analytical = config.analytical_velocity(y);
            let error = (u_lbm - u_analytical).abs();
            let relative_error = error / u_analytical.abs().max(1e-10);

            position_errors.push(y);
            velocity_errors.push(relative_error);
        }
    }

    let max_error = velocity_errors.iter().copied().fold(0.0, f64::max);
    let mean_error = velocity_errors.iter().sum::<f64>() / velocity_errors.len() as f64;
    let rms_error = (velocity_errors.iter().map(|e| e.powi(2)).sum::<f64>()
        / velocity_errors.len() as f64)
        .sqrt();

    ValidationResults {
        max_error,
        mean_error,
        rms_error,
        position_errors,
        velocity_errors,
        flow_rate_error: (results.flow_rate - results.analytical_flow_rate).abs()
            / results.analytical_flow_rate,
    }
}

/// Validation results
#[derive(Debug, Clone)]
pub struct ValidationResults {
    /// Maximum relative error
    pub max_error: f64,
    /// Mean relative error
    pub mean_error: f64,
    /// RMS relative error
    pub rms_error: f64,
    /// Position array
    pub position_errors: Vec<f64>,
    /// Velocity error array
    pub velocity_errors: Vec<f64>,
    /// Flow rate relative error
    pub flow_rate_error: f64,
}

impl ValidationResults {
    /// Print validation summary
    pub fn print_summary(&self) {
        println!("\n=== Validation Against Analytical Solution ===");
        println!("Maximum relative error: {:.2}%", self.max_error * 100.0);
        println!("Mean relative error: {:.2}%", self.mean_error * 100.0);
        println!("RMS relative error: {:.2}%", self.rms_error * 100.0);
        println!(
            "Flow rate relative error: {:.2}%",
            self.flow_rate_error * 100.0
        );

        // Assess accuracy
        if self.max_error < 0.01 {
            println!("Accuracy assessment: EXCELLENT (<1% error)");
        } else if self.max_error < 0.05 {
            println!("Accuracy assessment: GOOD (<5% error)");
        } else if self.max_error < 0.10 {
            println!("Accuracy assessment: ACCEPTABLE (<10% error)");
        } else {
            println!("Accuracy assessment: POOR (>10% error)");
        }
    }
}

fn main() -> CfdResult<()> {
    println!("=== LBM Poiseuille Flow Validation ===");

    // Run standard case
    println!("\n--- Standard Case ---");
    let mut poiseuille_std = LbmPoiseuille::new(PoiseuilleConfig::standard())?;
    let results_std = poiseuille_std.run_simulation()?;
    results_std.print_summary();

    let validation_std = validate_poiseuille_solution(poiseuille_std.config(), &results_std);
    validation_std.print_summary();

    // Run high Reynolds number case
    println!("\n--- High Reynolds Number Case ---");
    let mut poiseuille_high_re = LbmPoiseuille::new(PoiseuilleConfig::high_reynolds())?;
    let results_high_re = poiseuille_high_re.run_simulation()?;
    results_high_re.print_summary();

    let validation_high_re =
        validate_poiseuille_solution(poiseuille_high_re.config(), &results_high_re);
    validation_high_re.print_summary();

    if let Some(rate) = results_high_re.convergence_rate() {
        println!("Convergence rate: {:.3}", rate);
    }

    println!("\n=== LBM validation completed successfully! ===");

    Ok(())
}

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

    #[test]
    fn test_poiseuille_config() {
        let config = PoiseuilleConfig::standard();
        assert_eq!(config.height, 1.0);
        assert!(config.pressure_gradient < 0.0);
        assert!(config.viscosity > 0.0);
        assert!(config.tau > 0.5);
    }

    #[test]
    fn test_analytical_velocity() {
        let config = PoiseuilleConfig::standard();

        // Velocity should be zero at walls
        assert_eq!(config.analytical_velocity(0.0), 0.0);
        assert_eq!(config.analytical_velocity(config.height), 0.0);

        // Maximum velocity at center
        let u_center = config.analytical_velocity(config.height / 2.0);
        let u_max = config.max_velocity();
        assert!((u_center - u_max).abs() < 1e-10);
    }

    #[test]
    fn test_reynolds_number() {
        let config = PoiseuilleConfig::standard();
        let re = config.reynolds_number();
        assert!(re > 0.0);
    }

    #[test]
    fn test_analytical_flow_rate() {
        let config = PoiseuilleConfig::standard();
        let q = config.analytical_flow_rate();
        assert!(q > 0.0); // Flow in positive direction (negative pressure gradient)
    }

    #[test]
    fn test_high_re_config() {
        let config = PoiseuilleConfig::high_reynolds();
        let std_config = PoiseuilleConfig::standard();

        assert!(config.reynolds_number() > std_config.reynolds_number());
        assert!(config.pressure_gradient.abs() > std_config.pressure_gradient.abs());
    }
}