//! D3Q19 Lattice Boltzmann Method implementation
//!
//! This module implements the D3Q19 (3D, 19 velocities) lattice Boltzmann method
//! for simulating incompressible fluid flows in three dimensions.

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

/// D3Q19 lattice velocities (in lattice units)
const D3Q19_VELOCITIES: [Vector3<i32>; 19] = [
    Vector3::new(0, 0, 0),   // 0: Rest particle
    Vector3::new(1, 0, 0),   // 1: +X
    Vector3::new(-1, 0, 0),  // 2: -X
    Vector3::new(0, 1, 0),   // 3: +Y
    Vector3::new(0, -1, 0),  // 4: -Y
    Vector3::new(0, 0, 1),   // 5: +Z
    Vector3::new(0, 0, -1),  // 6: -Z
    Vector3::new(1, 1, 0),   // 7: +X+Y
    Vector3::new(-1, -1, 0), // 8: -X-Y
    Vector3::new(1, -1, 0),  // 9: +X-Y
    Vector3::new(-1, 1, 0),  // 10: -X+Y
    Vector3::new(1, 0, 1),   // 11: +X+Z
    Vector3::new(-1, 0, -1), // 12: -X-Z
    Vector3::new(1, 0, -1),  // 13: +X-Z
    Vector3::new(-1, 0, 1),  // 14: -X+Z
    Vector3::new(0, 1, 1),   // 15: +Y+Z
    Vector3::new(0, -1, -1), // 16: -Y-Z
    Vector3::new(0, 1, -1),  // 17: +Y-Z
    Vector3::new(0, -1, 1),  // 18: -Y+Z
];

/// D3Q19 lattice weights
const D3Q19_WEIGHTS: [f64; 19] = [
    1.0 / 3.0,  // 0: Rest particle
    1.0 / 18.0, // 1-6: Face neighbors
    1.0 / 18.0,
    1.0 / 18.0,
    1.0 / 18.0,
    1.0 / 18.0,
    1.0 / 18.0,
    1.0 / 36.0, // 7-18: Edge neighbors
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
    1.0 / 36.0,
];

/// Macroscopic variables for 3D flows
#[derive(Debug, Clone, PartialEq)]
pub struct MacroscopicVariables3D {
    /// Fluid density
    pub density: f64,
    /// Velocity vector
    pub velocity: Vector3<f64>,
}

impl MacroscopicVariables3D {
    /// Create new macroscopic variables
    #[must_use]
    pub fn new(density: f64, velocity: Vector3<f64>) -> Self {
        Self { density, velocity }
    }

    /// Zero velocity state
    #[must_use]
    pub fn at_rest(density: f64) -> Self {
        Self {
            density,
            velocity: Vector3::zeros(),
        }
    }
}

/// Parameters for D3Q19 lattice Boltzmann method
#[derive(Debug, Clone)]
pub struct D3Q19Parameters {
    /// Relaxation time for collision operator
    pub tau: f64,
    /// Kinematic viscosity
    pub nu: f64,
}

impl D3Q19Parameters {
    /// Create new D3Q19 parameters
    #[must_use]
    pub fn new(tau: f64) -> Self {
        let nu = (tau - 0.5) / 3.0;
        Self { tau, nu }
    }

    /// Validate parameters
    pub fn validate(&self) -> CfdResult<()> {
        if self.tau <= 0.5 {
            return Err(CfdError::invalid_parameter(
                "Relaxation time must be greater than 0.5 for stability",
            ));
        }
        if self.tau >= 2.0 {
            return Err(CfdError::invalid_parameter(
                "Relaxation time should be less than 2.0 for efficiency",
            ));
        }
        if self.nu <= 0.0 {
            return Err(CfdError::invalid_parameter(
                "Kinematic viscosity must be positive",
            ));
        }
        Ok(())
    }

    /// Create parameters from kinematic viscosity
    #[must_use]
    pub fn from_viscosity(nu: f64) -> Self {
        let tau = 3.0 * nu + 0.5;
        Self { tau, nu }
    }
}

impl Default for D3Q19Parameters {
    fn default() -> Self {
        Self::new(0.6)
    }
}

/// D3Q19 Lattice Boltzmann Method solver
pub struct D3Q19Solver {
    /// Grid dimensions
    nx: usize,
    ny: usize,
    nz: usize,

    /// Distribution functions f[x][y][z][i] where i is velocity direction
    f: Vec<Vec<Vec<Vec<f64>>>>,

    /// Temporary storage for streaming step
    f_temp: Vec<Vec<Vec<Vec<f64>>>>,

    /// Solver parameters
    params: D3Q19Parameters,
}

impl D3Q19Solver {
    /// Create new D3Q19 solver
    #[must_use]
    pub fn new(nx: usize, ny: usize, nz: usize, params: D3Q19Parameters) -> Self {
        params.validate().expect("Invalid D3Q19 parameters");

        let f = vec![vec![vec![vec![0.0; 19]; nz]; ny]; nx];
        let f_temp = vec![vec![vec![vec![0.0; 19]; nz]; ny]; nx];

        Self {
            nx,
            ny,
            nz,
            f,
            f_temp,
            params,
        }
    }

    /// Get lattice velocities
    #[must_use]
    pub fn lattice_velocities(&self) -> Vec<Vector3<i32>> {
        D3Q19_VELOCITIES.to_vec()
    }

    /// Get lattice weights
    #[must_use]
    pub fn weights(&self) -> Vec<f64> {
        D3Q19_WEIGHTS.to_vec()
    }

    /// Calculate equilibrium distribution function
    #[must_use]
    pub fn equilibrium_distribution(&self, density: f64, velocity: &Vector3<f64>) -> Vec<f64> {
        let mut f_eq = vec![0.0; 19];
        let u_sqr = velocity.norm_squared();

        for i in 0..19 {
            let e_i = Vector3::new(
                f64::from(D3Q19_VELOCITIES[i].x),
                f64::from(D3Q19_VELOCITIES[i].y),
                f64::from(D3Q19_VELOCITIES[i].z),
            );
            let e_dot_u = e_i.dot(velocity);

            // Equilibrium distribution: f_i^eq = w_i * rho * (1 + 3*e_i·u + 9/2*(e_i·u)^2 - 3/2*u^2)
            f_eq[i] = D3Q19_WEIGHTS[i]
                * density
                * (1.0 + 3.0 * e_dot_u + 4.5 * e_dot_u * e_dot_u - 1.5 * u_sqr);
        }

        f_eq
    }

    /// Set distribution function at a specific grid point
    pub fn set_distribution_at(&mut self, x: usize, y: usize, z: usize, f_values: &[f64]) {
        assert_eq!(f_values.len(), 19);
        self.f[x][y][z].copy_from_slice(f_values);
    }

    /// Get distribution function at a specific grid point
    #[must_use]
    pub fn distribution_at(&self, x: usize, y: usize, z: usize) -> Vec<f64> {
        self.f[x][y][z].clone()
    }

    /// Extract macroscopic variables (density and velocity) from distribution functions
    #[must_use]
    pub fn macroscopic_variables_at(&self, x: usize, y: usize, z: usize) -> MacroscopicVariables3D {
        let f_local = &self.f[x][y][z];

        // Density: sum of all distribution functions
        let density: f64 = f_local.iter().sum();

        // Momentum: sum of f_i * e_i
        let mut momentum = Vector3::zeros();
        for i in 0..19 {
            let e_i = Vector3::new(
                f64::from(D3Q19_VELOCITIES[i].x),
                f64::from(D3Q19_VELOCITIES[i].y),
                f64::from(D3Q19_VELOCITIES[i].z),
            );
            momentum += f_local[i] * e_i;
        }

        // Velocity: momentum / density
        let velocity = if density > 1e-15 {
            momentum / density
        } else {
            Vector3::zeros()
        };

        MacroscopicVariables3D::new(density, velocity)
    }

    /// BGK collision step (simplified, can be extended to MRT)
    pub fn collision_step(&mut self) {
        let omega = 1.0 / self.params.tau; // Collision frequency

        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    let vars = self.macroscopic_variables_at(x, y, z);
                    let f_eq = self.equilibrium_distribution(vars.density, &vars.velocity);

                    // BGK collision: f_i^new = f_i - omega * (f_i - f_i^eq)
                    for i in 0..19 {
                        self.f[x][y][z][i] -= omega * (self.f[x][y][z][i] - f_eq[i]);
                    }
                }
            }
        }
    }

    /// Streaming step (propagation)
    pub fn streaming_step(&mut self) {
        // Copy current state to temporary storage
        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    self.f_temp[x][y][z].copy_from_slice(&self.f[x][y][z]);
                }
            }
        }

        // Stream particles according to their velocities
        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    for i in 0..19 {
                        let e_i = D3Q19_VELOCITIES[i];
                        let x_src = (x as i32 - e_i.x).rem_euclid(self.nx as i32) as usize;
                        let y_src = (y as i32 - e_i.y).rem_euclid(self.ny as i32) as usize;
                        let z_src = (z as i32 - e_i.z).rem_euclid(self.nz as i32) as usize;

                        self.f[x][y][z][i] = self.f_temp[x_src][y_src][z_src][i];
                    }
                }
            }
        }
    }

    /// Complete LBM time step (collision + streaming)
    pub fn step(&mut self) {
        self.collision_step();
        self.streaming_step();
    }

    /// Complete LBM time step without boundary conditions (for periodic domains)
    pub fn step_periodic(&mut self) {
        self.collision_step();
        self.streaming_step();
    }

    /// Initialize uniform flow field
    pub fn initialize_uniform(&mut self, density: f64, velocity: Vector3<f64>) {
        let f_eq = self.equilibrium_distribution(density, &velocity);

        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    self.f[x][y][z].copy_from_slice(&f_eq);
                }
            }
        }
    }

    /// Calculate total mass in the domain
    #[must_use]
    pub fn total_mass(&self) -> f64 {
        let mut total = 0.0;
        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    let vars = self.macroscopic_variables_at(x, y, z);
                    total += vars.density;
                }
            }
        }
        total
    }

    /// Get grid dimensions
    #[must_use]
    pub fn dimensions(&self) -> (usize, usize, usize) {
        (self.nx, self.ny, self.nz)
    }

    /// Get solver parameters
    #[must_use]
    pub fn parameters(&self) -> &D3Q19Parameters {
        &self.params
    }

    /// Calculate kinetic energy in the domain
    #[must_use]
    pub fn kinetic_energy(&self) -> f64 {
        let mut total_ke = 0.0;
        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    let vars = self.macroscopic_variables_at(x, y, z);
                    total_ke += 0.5 * vars.density * vars.velocity.norm_squared();
                }
            }
        }
        total_ke
    }

    /// Calculate maximum velocity in the domain
    #[must_use]
    pub fn max_velocity(&self) -> f64 {
        let mut max_vel = 0.0f64;
        for x in 0..self.nx {
            for y in 0..self.ny {
                for z in 0..self.nz {
                    let vars = self.macroscopic_variables_at(x, y, z);
                    max_vel = max_vel.max(vars.velocity.norm());
                }
            }
        }
        max_vel
    }

    /// Check CFL condition for numerical stability
    #[must_use]
    pub fn check_cfl_condition(&self) -> bool {
        let max_vel = self.max_velocity();
        // CFL condition: max_velocity * dt / dx < 1
        // In LBM, dt = dx = 1 in lattice units, so we need max_vel < 1
        max_vel < 0.1 // Conservative limit
    }
}

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

    #[test]
    fn test_d3q19_parameters() {
        let params = D3Q19Parameters::new(0.6);
        assert_relative_eq!(params.tau, 0.6);
        assert_relative_eq!(params.nu, (0.6 - 0.5) / 3.0);
    }

    #[test]
    fn test_d3q19_from_viscosity() {
        let nu = 0.1;
        let params = D3Q19Parameters::from_viscosity(nu);
        assert_relative_eq!(params.nu, nu);
        assert_relative_eq!(params.tau, 3.0 * nu + 0.5);
    }

    #[test]
    fn test_d3q19_grid_creation() {
        let solver = D3Q19Solver::new(10, 8, 6, D3Q19Parameters::default());
        assert_eq!(solver.dimensions(), (10, 8, 6));
    }

    #[test]
    fn test_equilibrium_mass_conservation() {
        let solver = D3Q19Solver::new(5, 5, 5, D3Q19Parameters::default());
        let density = 1.5;
        let velocity = Vector3::new(0.1, -0.05, 0.03);

        let f_eq = solver.equilibrium_distribution(density, &velocity);
        let sum: f64 = f_eq.iter().sum();

        assert_relative_eq!(sum, density, epsilon = 1e-15);
    }

    #[test]
    fn test_equilibrium_momentum_conservation() {
        let solver = D3Q19Solver::new(5, 5, 5, D3Q19Parameters::default());
        let density = 1.0;
        let velocity = Vector3::new(0.1, -0.05, 0.03);

        let f_eq = solver.equilibrium_distribution(density, &velocity);

        // Calculate momentum from equilibrium distribution
        let mut momentum = Vector3::zeros();
        for i in 0..19 {
            let e_i = Vector3::new(
                D3Q19_VELOCITIES[i].x as f64,
                D3Q19_VELOCITIES[i].y as f64,
                D3Q19_VELOCITIES[i].z as f64,
            );
            momentum += f_eq[i] * e_i;
        }

        let expected_momentum = density * velocity;
        assert_relative_eq!(momentum.x, expected_momentum.x, epsilon = 1e-15);
        assert_relative_eq!(momentum.y, expected_momentum.y, epsilon = 1e-15);
        assert_relative_eq!(momentum.z, expected_momentum.z, epsilon = 1e-15);
    }
}