scirs2-integrate 0.4.2

//! Method of Manufactured Solutions (MMS) toolkit for verification
//!
//! This module provides tools for verifying numerical methods using the Method of
//! Manufactured Solutions. MMS is a powerful technique for code verification where
//! an exact solution is chosen, the governing equation is modified by adding source
//! terms, and the numerical solution is compared against the exact solution.
//!
//! # Examples
//!
//! ## ODE Verification
//! ```
//! use scirs2_integrate::verification::{MMSODEProblem, polynomial_solution};
//! use scirs2_core::ndarray::Array1;
//!
//! // Create a manufactured ODE problem with polynomial solution
//! let exact_solution = polynomial_solution(vec![1.0, 2.0, 3.0]); // 1 + 2t + 3t²
//! let problem = MMSODEProblem::new(exact_solution, [0.0, 1.0]);
//!
//! // The manufactured source term is automatically computed
//! // Solve numerically and verify order of accuracy
//! ```
//!
//! ## PDE Verification  
//! ```
//! use scirs2_integrate::verification::{MMSPDEProblem, trigonometric_solution_2d};
//!
//! // Create manufactured 2D Poisson problem
//! let exact = trigonometric_solution_2d(1.0, 2.0); // sin(x) * cos(2y)
//! let problem = MMSPDEProblem::new_poisson_2d(exact, [0.0, 1.0], [0.0, 1.0]);
//! ```

use crate::common::IntegrateFloat;
use crate::error::{IntegrateError, IntegrateResult};
use scirs2_core::ndarray::{ArrayView1, ArrayView2};
use std::f64::consts::PI;
use std::fmt;

/// Trait for exact solutions in MMS problems
pub trait ExactSolution<F: IntegrateFloat> {
    /// Evaluate the exact solution at given point(s)
    fn evaluate(&self, coordinates: &[F]) -> F;

    /// Evaluate the first derivative with respect to specified variable
    fn derivative(&self, coordinates: &[F], variable: usize) -> F;

    /// Evaluate the second derivative with respect to specified variable
    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F;

    /// Evaluate mixed partial derivatives (for PDEs)
    fn mixed_derivative(&self, coordinates: &[F], var1: usize, var2: usize) -> F {
        // Default implementation using finite differences
        let h = F::from(1e-8).expect("Failed to convert constant to float");
        let mut coords_plus = coordinates.to_vec();
        let mut coords_minus = coordinates.to_vec();

        coords_plus[var2] += h;
        coords_minus[var2] -= h;

        let deriv_plus = self.derivative(&coords_plus, var1);
        let deriv_minus = self.derivative(&coords_minus, var1);

        (deriv_plus - deriv_minus)
            / (F::from(2.0).expect("Failed to convert constant to float") * h)
    }

    /// Get the dimensionality of the problem
    fn dimension(&self) -> usize;
}

/// Polynomial exact solution for ODE problems
#[derive(Debug, Clone)]
pub struct PolynomialSolution<F: IntegrateFloat> {
    coefficients: Vec<F>,
}

impl<F: IntegrateFloat> PolynomialSolution<F> {
    /// Create a new polynomial solution: `sum(coeff[i] * t^i)`
    pub fn new(coefficients: Vec<F>) -> Self {
        Self { coefficients }
    }
}

impl<F: IntegrateFloat> ExactSolution<F> for PolynomialSolution<F> {
    fn evaluate(&self, coordinates: &[F]) -> F {
        let t = coordinates[0];
        let mut result = F::zero();
        let mut t_power = F::one();

        for &coeff in &self.coefficients {
            result += coeff * t_power;
            t_power *= t;
        }

        result
    }

    fn derivative(&self, coordinates: &[F], variable: usize) -> F {
        let t = coordinates[0];
        let mut result = F::zero();
        let mut t_power = F::one();

        for (i, &coeff) in self.coefficients.iter().enumerate().skip(1) {
            result += F::from(i).expect("Failed to convert to float") * coeff * t_power;
            t_power *= t;
        }

        result
    }

    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F {
        let t = coordinates[0];
        let mut result = F::zero();
        let mut t_power = F::one();

        for (i, &coeff) in self.coefficients.iter().enumerate().skip(2) {
            let factor = F::from(i * (i - 1)).expect("Operation failed");
            result += factor * coeff * t_power;
            t_power *= t;
        }

        result
    }

    fn dimension(&self) -> usize {
        1
    }
}

/// Trigonometric exact solution for 2D PDE problems  
#[derive(Debug, Clone)]
pub struct TrigonometricSolution2D<F: IntegrateFloat> {
    freq_x: F,
    freq_y: F,
    phase_x: F,
    phase_y: F,
}

impl<F: IntegrateFloat> TrigonometricSolution2D<F> {
    /// Create sin(freq_x * x + phase_x) * cos(freq_y * y + phase_y)
    pub fn new(freq_x: F, freq_y: F, phase_x: F, phase_y: F) -> Self {
        Self {
            freq_x,
            freq_y,
            phase_x,
            phase_y,
        }
    }

    /// Create sin(freq_x * x) * cos(freq_y * y)
    pub fn simple(freq_x: F, freqy: F) -> Self {
        Self::new(freq_x, freqy, F::zero(), F::zero())
    }
}

impl<F: IntegrateFloat> ExactSolution<F> for TrigonometricSolution2D<F> {
    fn evaluate(&self, coordinates: &[F]) -> F {
        let x = coordinates[0];
        let y = coordinates[1];
        (self.freq_x * x + self.phase_x).sin() * (self.freq_y * y + self.phase_y).cos()
    }

    fn derivative(&self, coordinates: &[F], variable: usize) -> F {
        let x = coordinates[0];
        let y = coordinates[1];

        match variable {
            0 => {
                self.freq_x
                    * (self.freq_x * x + self.phase_x).cos()
                    * (self.freq_y * y + self.phase_y).cos()
            }
            1 => {
                -self.freq_y
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).sin()
            }
            _ => F::zero(),
        }
    }

    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F {
        let x = coordinates[0];
        let y = coordinates[1];

        match variable {
            0 => {
                -self.freq_x
                    * self.freq_x
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
            }
            1 => {
                -self.freq_y
                    * self.freq_y
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
            }
            _ => F::zero(),
        }
    }

    fn dimension(&self) -> usize {
        2
    }
}

/// Manufactured ODE problem for verification
#[derive(Debug)]
pub struct MMSODEProblem<F: IntegrateFloat, S: ExactSolution<F>> {
    exact_solution: S,
    time_span: [F; 2],
    _phantom: std::marker::PhantomData<F>,
}

impl<F: IntegrateFloat, S: ExactSolution<F>> MMSODEProblem<F, S> {
    /// Create a new manufactured ODE problem
    pub fn new(exact_solution: S, timespan: [F; 2]) -> Self {
        Self {
            exact_solution,
            time_span: timespan,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Get the manufactured source term for: y' = f(t, y) + source(t)
    /// where f is the original RHS and source is computed from exact solution
    pub fn source_term(&self, t: F) -> F {
        let coords = [t];
        // For y' = f(t, y), source = y_exact'(t) - f(t, y_exact(t))
        // If f is linear: y' = a*y + b, then source = y_exact'(t) - a*y_exact(t) - b
        // For autonomous f: y' = f(y), source = y_exact'(t) - f(y_exact(t))
        self.exact_solution.derivative(&coords, 0)
    }

    /// Get initial condition from exact solution
    pub fn initial_condition(&self) -> F {
        self.exact_solution.evaluate(&[self.time_span[0]])
    }

    /// Evaluate exact solution at given time
    pub fn exact_at(&self, t: F) -> F {
        self.exact_solution.evaluate(&[t])
    }

    /// Get time span
    pub fn time_span(&self) -> [F; 2] {
        self.time_span
    }
}

/// Manufactured PDE problem for verification  
#[derive(Debug)]
pub struct MMSPDEProblem<F: IntegrateFloat, S: ExactSolution<F>> {
    exact_solution: S,
    domain_x: [F; 2],
    domain_y: [F; 2],
    domain_z: Option<[F; 2]>,
    pde_type: PDEType,
    parameters: PDEParameters<F>,
    _phantom: std::marker::PhantomData<F>,
}

/// Parameters for different PDE types
#[derive(Debug, Clone)]
pub struct PDEParameters<F: IntegrateFloat> {
    /// Diffusion coefficient (for diffusion/heat equation)
    pub diffusion_coeff: F,
    /// Wave speed (for wave equation)
    pub wave_speed: F,
    /// Advection velocity (for advection-diffusion)
    pub advection_velocity: Vec<F>,
    /// Reaction coefficient (for reaction-diffusion)
    pub reaction_coeff: F,
    /// Helmholtz parameter (for Helmholtz equation)
    pub helmholtz_k: F,
}

impl<F: IntegrateFloat> Default for PDEParameters<F> {
    fn default() -> Self {
        Self {
            diffusion_coeff: F::one(),
            wave_speed: F::one(),
            advection_velocity: vec![F::zero()],
            reaction_coeff: F::zero(),
            helmholtz_k: F::one(),
        }
    }
}

#[derive(Debug, Clone)]
pub enum PDEType {
    Poisson2D,
    Poisson3D,
    Diffusion2D,
    Diffusion3D,
    Wave2D,
    Wave3D,
    AdvectionDiffusion2D,
    AdvectionDiffusion3D,
    Helmholtz2D,
    Helmholtz3D,
}

impl<F: IntegrateFloat, S: ExactSolution<F>> MMSPDEProblem<F, S> {
    /// Create manufactured 2D Poisson problem: -∇²u = f
    pub fn new_poisson_2d(exact_solution: S, domain_x: [F; 2], domainy: [F; 2]) -> Self {
        Self {
            exact_solution,
            domain_x,
            domain_y: domainy,
            domain_z: None,
            pde_type: PDEType::Poisson2D,
            parameters: PDEParameters::default(),
            _phantom: std::marker::PhantomData,
        }
    }

    /// Create manufactured 3D Poisson problem: -∇²u = f
    pub fn new_poisson_3d(
        exact_solution: S,
        domain_x: [F; 2],
        domain_y: [F; 2],
        domain_z: [F; 2],
    ) -> Self {
        Self {
            exact_solution,
            domain_x,
            domain_y,
            domain_z: Some(domain_z),
            pde_type: PDEType::Poisson3D,
            parameters: PDEParameters::default(),
            _phantom: std::marker::PhantomData,
        }
    }

    /// Create manufactured 2D diffusion problem: ∂u/∂t = α∇²u + f
    pub fn new_diffusion_2d(
        exact_solution: S,
        domain_x: [F; 2],
        domain_y: [F; 2],
        diffusion_coeff: F,
    ) -> Self {
        let mut params = PDEParameters::default();
        params.diffusion_coeff = diffusion_coeff;
        Self {
            exact_solution,
            domain_x,
            domain_y,
            domain_z: None,
            pde_type: PDEType::Diffusion2D,
            parameters: params,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Create manufactured 2D wave problem: ∂²u/∂t² = c²∇²u + f
    pub fn new_wave_2d(
        exact_solution: S,
        domain_x: [F; 2],
        domain_y: [F; 2],
        wave_speed: F,
    ) -> Self {
        let mut params = PDEParameters::default();
        params.wave_speed = wave_speed;
        Self {
            exact_solution,
            domain_x,
            domain_y,
            domain_z: None,
            pde_type: PDEType::Wave2D,
            parameters: params,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Create manufactured 2D Helmholtz problem: ∇²u + k²u = f
    pub fn new_helmholtz_2d(exact_solution: S, domain_x: [F; 2], domainy: [F; 2], k: F) -> Self {
        let mut params = PDEParameters::default();
        params.helmholtz_k = k;
        Self {
            exact_solution,
            domain_x,
            domain_y: domainy,
            domain_z: None,
            pde_type: PDEType::Helmholtz2D,
            parameters: params,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Get manufactured source term for the PDE
    pub fn source_term(&self, coordinates: &[F]) -> F {
        match self.pde_type {
            PDEType::Poisson2D => {
                // For -∇²u = f, source f = -∇²u_exact
                -(self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1))
            }
            PDEType::Poisson3D => {
                // For -∇²u = f in 3D
                -(self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1)
                    + self.exact_solution.second_derivative(coordinates, 2))
            }
            PDEType::Diffusion2D => {
                // For ∂u/∂t = α∇²u + f
                // source f = ∂u_exact/∂t - α∇²u_exact
                let time_dim = coordinates.len() - 1;
                let spatial_laplacian = self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1);
                self.exact_solution.derivative(coordinates, time_dim)
                    - self.parameters.diffusion_coeff * spatial_laplacian
            }
            PDEType::Wave2D => {
                // For ∂²u/∂t² = c²∇²u + f
                // source f = ∂²u_exact/∂t² - c²∇²u_exact
                // Simplified: assume time coordinate is last
                let spatial_laplacian = self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1);
                // Second time derivative would need higher-order implementation
                F::zero()
                    - self.parameters.wave_speed * self.parameters.wave_speed * spatial_laplacian
            }
            PDEType::Helmholtz2D => {
                // For ∇²u + k²u = f
                // source f = ∇²u_exact + k²u_exact
                let laplacian = self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1);
                laplacian
                    + self.parameters.helmholtz_k
                        * self.parameters.helmholtz_k
                        * self.exact_solution.evaluate(coordinates)
            }
            PDEType::AdvectionDiffusion2D => {
                // For ∂u/∂t + v·∇u = α∇²u + f
                // source f = ∂u_exact/∂t + v·∇u_exact - α∇²u_exact
                let time_dim = coordinates.len() - 1;
                let time_deriv = self.exact_solution.derivative(coordinates, time_dim);
                let spatial_laplacian = self.exact_solution.second_derivative(coordinates, 0)
                    + self.exact_solution.second_derivative(coordinates, 1);

                let mut advection_term = F::zero();
                for (i, &v_i) in self
                    .parameters
                    .advection_velocity
                    .iter()
                    .enumerate()
                    .take(2)
                {
                    advection_term += v_i * self.exact_solution.derivative(coordinates, i);
                }

                time_deriv + advection_term - self.parameters.diffusion_coeff * spatial_laplacian
            }
            _ => F::zero(),
        }
    }

    /// Get boundary condition from exact solution (2D)
    pub fn boundary_condition(&self, x: F, y: F) -> F {
        self.exact_solution.evaluate(&[x, y])
    }

    /// Get boundary condition from exact solution (3D)
    pub fn boundary_condition_3d(&self, x: F, y: F, z: F) -> F {
        self.exact_solution.evaluate(&[x, y, z])
    }

    /// Evaluate exact solution (2D)
    pub fn exact_at(&self, x: F, y: F) -> F {
        self.exact_solution.evaluate(&[x, y])
    }

    /// Evaluate exact solution (3D)
    pub fn exact_at_3d(&self, x: F, y: F, z: F) -> F {
        self.exact_solution.evaluate(&[x, y, z])
    }

    /// Get domain bounds (2D)
    pub fn domain(&self) -> ([F; 2], [F; 2]) {
        (self.domain_x, self.domain_y)
    }

    /// Get domain bounds (3D)
    pub fn domain_3d(&self) -> ([F; 2], [F; 2], [F; 2]) {
        (
            self.domain_x,
            self.domain_y,
            self.domain_z.unwrap_or([F::zero(), F::one()]),
        )
    }

    /// Get PDE parameters
    pub fn parameters(&self) -> &PDEParameters<F> {
        &self.parameters
    }

    /// Check if problem is 3D
    pub fn is_3d(&self) -> bool {
        self.domain_z.is_some()
    }
}

/// Convergence analysis results
#[derive(Debug, Clone)]
pub struct ConvergenceAnalysis<F: IntegrateFloat> {
    /// Grid sizes or time steps used
    pub grid_sizes: Vec<F>,
    /// Error norms for each grid size
    pub errors: Vec<F>,
    /// Estimated order of accuracy
    pub order: F,
    /// Confidence interval for the order
    pub order_confidence: (F, F),
}

impl<F: IntegrateFloat> ConvergenceAnalysis<F> {
    /// Compute order of accuracy from grid sizes and errors
    pub fn compute_order(_gridsizes: Vec<F>, errors: Vec<F>) -> IntegrateResult<Self> {
        if _gridsizes.len() != errors.len() || _gridsizes.len() < 2 {
            return Err(IntegrateError::ValueError(
                "Need at least 2 points for convergence analysis".to_string(),
            ));
        }

        // Use least squares to fit log(error) = log(C) + p*log(h)
        let n = _gridsizes.len();
        let mut sum_log_h = F::zero();
        let mut sum_log_e = F::zero();
        let mut sum_log_h_sq = F::zero();
        let mut sum_log_h_log_e = F::zero();

        for (h, e) in _gridsizes.iter().zip(errors.iter()) {
            if *e <= F::zero() || *h <= F::zero() {
                return Err(IntegrateError::ValueError(
                    "Grid sizes and errors must be positive".to_string(),
                ));
            }

            let log_h = h.ln();
            let log_e = e.ln();

            sum_log_h += log_h;
            sum_log_e += log_e;
            sum_log_h_sq += log_h * log_h;
            sum_log_h_log_e += log_h * log_e;
        }

        let n_f = F::from(n).expect("Failed to convert to float");
        let denominator = n_f * sum_log_h_sq - sum_log_h * sum_log_h;

        if denominator.abs() < F::from(1e-12).expect("Failed to convert constant to float") {
            return Err(IntegrateError::ComputationError(
                "Cannot compute order - insufficient variation in grid sizes".to_string(),
            ));
        }

        let order = (n_f * sum_log_h_log_e - sum_log_h * sum_log_e) / denominator;

        // Simple confidence interval (±0.1 for demonstration)
        let confidence_delta = F::from(0.1).expect("Failed to convert constant to float");
        let order_confidence = (order - confidence_delta, order + confidence_delta);

        Ok(Self {
            grid_sizes: _gridsizes,
            errors,
            order,
            order_confidence,
        })
    }

    /// Check if the computed order matches expected order within tolerance
    pub fn verify_order(&self, expectedorder: F, tolerance: F) -> bool {
        (self.order - expectedorder).abs() <= tolerance
    }
}

impl<F: IntegrateFloat> fmt::Display for ConvergenceAnalysis<F> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "Convergence Analysis Results:")?;
        writeln!(f, "Grid Size    Error")?;
        writeln!(f, "─────────────────────")?;

        for (h, e) in self.grid_sizes.iter().zip(self.errors.iter()) {
            writeln!(f, "{h:9.2e}  {e:9.2e}")?;
        }

        writeln!(f, "─────────────────────")?;
        writeln!(f, "Estimated order: {:.3}", self.order)?;
        writeln!(
            f,
            "95% confidence: ({:.3}, {:.3})",
            self.order_confidence.0, self.order_confidence.1
        )?;

        Ok(())
    }
}

/// Error analysis utilities
pub struct ErrorAnalysis;

impl ErrorAnalysis {
    /// Compute L2 norm of error
    pub fn l2_norm<F: IntegrateFloat>(
        exact: ArrayView1<F>,
        numerical: ArrayView1<F>,
    ) -> IntegrateResult<F> {
        if exact.len() != numerical.len() {
            return Err(IntegrateError::ValueError(
                "Arrays must have same length".to_string(),
            ));
        }

        let mut sum_sq = F::zero();
        for (e, n) in exact.iter().zip(numerical.iter()) {
            let diff = *e - *n;
            sum_sq += diff * diff;
        }

        Ok((sum_sq / F::from(exact.len()).expect("Operation failed")).sqrt())
    }

    /// Compute maximum norm of error
    pub fn max_norm<F: IntegrateFloat>(
        exact: ArrayView1<F>,
        numerical: ArrayView1<F>,
    ) -> IntegrateResult<F> {
        if exact.len() != numerical.len() {
            return Err(IntegrateError::ValueError(
                "Arrays must have same length".to_string(),
            ));
        }

        let mut max_error = F::zero();
        for (e, n) in exact.iter().zip(numerical.iter()) {
            let diff = (*e - *n).abs();
            if diff > max_error {
                max_error = diff;
            }
        }

        Ok(max_error)
    }

    /// Compute L2 norm of error for 2D arrays
    pub fn l2_norm_2d<F: IntegrateFloat>(
        exact: ArrayView2<F>,
        numerical: ArrayView2<F>,
    ) -> IntegrateResult<F> {
        if exact.shape() != numerical.shape() {
            return Err(IntegrateError::ValueError(
                "Arrays must have same shape".to_string(),
            ));
        }

        let mut sum_sq = F::zero();
        let mut count = 0;

        for (e, n) in exact.iter().zip(numerical.iter()) {
            let diff = *e - *n;
            sum_sq += diff * diff;
            count += 1;
        }

        Ok((sum_sq / F::from(count).expect("Failed to convert to float")).sqrt())
    }
}

/// Convenience functions for creating common exact solutions
#[allow(dead_code)]
pub fn polynomial_solution<F: IntegrateFloat>(coefficients: Vec<F>) -> PolynomialSolution<F> {
    PolynomialSolution::new(coefficients)
}

#[allow(dead_code)]
pub fn trigonometric_solution_2d<F: IntegrateFloat>(
    freq_x: F,
    freq_y: F,
) -> TrigonometricSolution2D<F> {
    TrigonometricSolution2D::simple(freq_x, freq_y)
}

/// Exponential exact solution for problems with exponential behavior
#[derive(Debug, Clone)]
pub struct ExponentialSolution<F: IntegrateFloat> {
    amplitude: F,
    decay_rate: F,
    phase: F,
}

impl<F: IntegrateFloat> ExponentialSolution<F> {
    /// Create exp(decay_rate * t + phase) solution
    pub fn new(_amplitude: F, decayrate: F, phase: F) -> Self {
        Self {
            amplitude: _amplitude,
            decay_rate: decayrate,
            phase,
        }
    }

    /// Create simple exp(decay_rate * t) solution
    pub fn simple(amplitude: F, decayrate: F) -> Self {
        Self::new(amplitude, decayrate, F::zero())
    }
}

impl<F: IntegrateFloat> ExactSolution<F> for ExponentialSolution<F> {
    fn evaluate(&self, coordinates: &[F]) -> F {
        let t = coordinates[0];
        self.amplitude * (self.decay_rate * t + self.phase).exp()
    }

    fn derivative(&self, coordinates: &[F], variable: usize) -> F {
        let t = coordinates[0];
        self.amplitude * self.decay_rate * (self.decay_rate * t + self.phase).exp()
    }

    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F {
        let t = coordinates[0];
        self.amplitude
            * self.decay_rate
            * self.decay_rate
            * (self.decay_rate * t + self.phase).exp()
    }

    fn dimension(&self) -> usize {
        1
    }
}

/// Combined solution: polynomial + trigonometric + exponential
#[derive(Debug, Clone)]
pub struct CombinedSolution<F: IntegrateFloat> {
    polynomial: Option<PolynomialSolution<F>>,
    trigonometric: Option<TrigonometricSolution2D<F>>,
    exponential: Option<ExponentialSolution<F>>,
    dimension: usize,
}

impl<F: IntegrateFloat> CombinedSolution<F> {
    /// Create a new combined solution
    pub fn new(dimension: usize) -> Self {
        Self {
            polynomial: None,
            trigonometric: None,
            exponential: None,
            dimension,
        }
    }

    /// Add polynomial component
    pub fn with_polynomial(mut self, poly: PolynomialSolution<F>) -> Self {
        self.polynomial = Some(poly);
        self
    }

    /// Add trigonometric component (for 2D problems)
    pub fn with_trigonometric(mut self, trig: TrigonometricSolution2D<F>) -> Self {
        self.trigonometric = Some(trig);
        self
    }

    /// Add exponential component
    pub fn with_exponential(mut self, exp: ExponentialSolution<F>) -> Self {
        self.exponential = Some(exp);
        self
    }
}

impl<F: IntegrateFloat> ExactSolution<F> for CombinedSolution<F> {
    fn evaluate(&self, coordinates: &[F]) -> F {
        let mut result = F::zero();

        if let Some(ref poly) = self.polynomial {
            result += poly.evaluate(coordinates);
        }

        if let Some(ref trig) = self.trigonometric {
            if coordinates.len() >= 2 {
                result += trig.evaluate(coordinates);
            }
        }

        if let Some(ref exp) = self.exponential {
            result += exp.evaluate(coordinates);
        }

        result
    }

    fn derivative(&self, coordinates: &[F], variable: usize) -> F {
        let mut result = F::zero();

        if let Some(ref poly) = self.polynomial {
            result += poly.derivative(coordinates, variable);
        }

        if let Some(ref trig) = self.trigonometric {
            if coordinates.len() >= 2 {
                result += trig.derivative(coordinates, variable);
            }
        }

        if let Some(ref exp) = self.exponential {
            result += exp.derivative(coordinates, variable);
        }

        result
    }

    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F {
        let mut result = F::zero();

        if let Some(ref poly) = self.polynomial {
            result += poly.second_derivative(coordinates, variable);
        }

        if let Some(ref trig) = self.trigonometric {
            if coordinates.len() >= 2 {
                result += trig.second_derivative(coordinates, variable);
            }
        }

        if let Some(ref exp) = self.exponential {
            result += exp.second_derivative(coordinates, variable);
        }

        result
    }

    fn dimension(&self) -> usize {
        self.dimension
    }
}

/// 3D trigonometric solution for 3D PDE problems
#[derive(Debug, Clone)]
pub struct TrigonometricSolution3D<F: IntegrateFloat> {
    freq_x: F,
    freq_y: F,
    freq_z: F,
    phase_x: F,
    phase_y: F,
    phase_z: F,
}

impl<F: IntegrateFloat> TrigonometricSolution3D<F> {
    /// Create sin(freq_x * x + phase_x) * cos(freq_y * y + phase_y) * sin(freq_z * z + phase_z)
    #[allow(clippy::too_many_arguments)]
    pub fn new(freq_x: F, freq_y: F, freq_z: F, phase_x: F, phase_y: F, phase_z: F) -> Self {
        Self {
            freq_x,
            freq_y,
            freq_z,
            phase_x,
            phase_y,
            phase_z,
        }
    }

    /// Create sin(freq_x * x) * cos(freq_y * y) * sin(freq_z * z)
    pub fn simple(freq_x: F, freq_y: F, freq_z: F) -> Self {
        Self::new(freq_x, freq_y, freq_z, F::zero(), F::zero(), F::zero())
    }
}

impl<F: IntegrateFloat> ExactSolution<F> for TrigonometricSolution3D<F> {
    fn evaluate(&self, coordinates: &[F]) -> F {
        let x = coordinates[0];
        let y = coordinates[1];
        let z = coordinates[2];
        (self.freq_x * x + self.phase_x).sin()
            * (self.freq_y * y + self.phase_y).cos()
            * (self.freq_z * z + self.phase_z).sin()
    }

    fn derivative(&self, coordinates: &[F], variable: usize) -> F {
        let x = coordinates[0];
        let y = coordinates[1];
        let z = coordinates[2];

        match variable {
            0 => {
                self.freq_x
                    * (self.freq_x * x + self.phase_x).cos()
                    * (self.freq_y * y + self.phase_y).cos()
                    * (self.freq_z * z + self.phase_z).sin()
            }
            1 => {
                -self.freq_y
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).sin()
                    * (self.freq_z * z + self.phase_z).sin()
            }
            2 => {
                self.freq_z
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
                    * (self.freq_z * z + self.phase_z).cos()
            }
            _ => F::zero(),
        }
    }

    fn second_derivative(&self, coordinates: &[F], variable: usize) -> F {
        let x = coordinates[0];
        let y = coordinates[1];
        let z = coordinates[2];

        match variable {
            0 => {
                -self.freq_x
                    * self.freq_x
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
                    * (self.freq_z * z + self.phase_z).sin()
            }
            1 => {
                -self.freq_y
                    * self.freq_y
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
                    * (self.freq_z * z + self.phase_z).sin()
            }
            2 => {
                -self.freq_z
                    * self.freq_z
                    * (self.freq_x * x + self.phase_x).sin()
                    * (self.freq_y * y + self.phase_y).cos()
                    * (self.freq_z * z + self.phase_z).sin()
            }
            _ => F::zero(),
        }
    }

    fn dimension(&self) -> usize {
        3
    }
}

/// System of equations verification (simplified without trait objects)
#[derive(Debug, Clone)]
pub struct SystemVerification<F: IntegrateFloat> {
    /// Number of components in the system
    pub system_size: usize,
    /// Component names for identification
    pub component_names: Vec<String>,
    /// Phantom data to maintain type parameter
    _phantom: std::marker::PhantomData<F>,
}

impl<F: IntegrateFloat> SystemVerification<F> {
    /// Create new system verification
    pub fn new(systemsize: usize) -> Self {
        let component_names = (0..systemsize).map(|i| format!("Component {i}")).collect();
        Self {
            system_size: systemsize,
            component_names,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Create with custom component names
    pub fn with_names(componentnames: Vec<String>) -> Self {
        let system_size = componentnames.len();
        Self {
            system_size,
            component_names: componentnames,
            _phantom: std::marker::PhantomData,
        }
    }

    /// Verify system solution using provided exact solutions
    pub fn verify_system<S1, S2>(
        &self,
        exact_solutions: &[S1],
        numerical_solutions: &[S2],
        coordinates: &[F],
    ) -> Vec<F>
    where
        S1: ExactSolution<F>,
        S2: Fn(&[F]) -> F,
    {
        assert_eq!(exact_solutions.len(), self.system_size);
        assert_eq!(numerical_solutions.len(), self.system_size);

        let mut errors = Vec::with_capacity(self.system_size);
        for i in 0..self.system_size {
            let exact = exact_solutions[i].evaluate(coordinates);
            let numerical = numerical_solutions[i](coordinates);
            errors.push((exact - numerical).abs());
        }
        errors
    }
}

/// Automated verification workflow for numerical methods
pub struct VerificationWorkflow<F: IntegrateFloat> {
    /// Test cases to verify
    pub test_cases: Vec<VerificationTestCase<F>>,
}

#[derive(Debug, Clone)]
pub struct VerificationTestCase<F: IntegrateFloat> {
    /// Test name
    pub name: String,
    /// Expected order of accuracy
    pub expected_order: F,
    /// Tolerance for order verification
    pub order_tolerance: F,
    /// Grid sizes to test
    pub grid_sizes: Vec<F>,
    /// Expected errors (if known)
    pub expected_errors: Option<Vec<F>>,
}

impl<F: IntegrateFloat> VerificationWorkflow<F> {
    /// Create new verification workflow
    pub fn new() -> Self {
        Self {
            test_cases: Vec::new(),
        }
    }

    /// Add test case to workflow
    pub fn add_test_case(&mut self, testcase: VerificationTestCase<F>) {
        self.test_cases.push(testcase);
    }

    /// Run all verification tests
    pub fn run_verification<S: Fn(&[F]) -> IntegrateResult<F>>(
        &self,
        solver: S,
    ) -> Vec<VerificationResult<F>> {
        let mut results = Vec::new();

        for test_case in &self.test_cases {
            let mut errors = Vec::new();

            for &grid_size in &test_case.grid_sizes {
                match solver(&[grid_size]) {
                    Ok(error) => errors.push(error),
                    Err(_) => {
                        results.push(VerificationResult {
                            test_name: test_case.name.clone(),
                            passed: false,
                            computed_order: None,
                            error_message: Some("Solver failed".to_string()),
                        });
                        break;
                    }
                }
            }

            if errors.len() == test_case.grid_sizes.len() {
                match ConvergenceAnalysis::compute_order(test_case.grid_sizes.clone(), errors) {
                    Ok(analysis) => {
                        let passed = analysis
                            .verify_order(test_case.expected_order, test_case.order_tolerance);
                        results.push(VerificationResult {
                            test_name: test_case.name.clone(),
                            passed,
                            computed_order: Some(analysis.order),
                            error_message: if passed {
                                None
                            } else {
                                Some("Order verification failed".to_string())
                            },
                        });
                    }
                    Err(e) => {
                        results.push(VerificationResult {
                            test_name: test_case.name.clone(),
                            passed: false,
                            computed_order: None,
                            error_message: Some(format!("Convergence analysis failed: {e}")),
                        });
                    }
                }
            }
        }

        results
    }
}

#[derive(Debug, Clone)]
pub struct VerificationResult<F: IntegrateFloat> {
    /// Test name
    pub test_name: String,
    /// Whether test passed
    pub passed: bool,
    /// Computed order of accuracy
    pub computed_order: Option<F>,
    /// Error message if test failed
    pub error_message: Option<String>,
}

impl<F: IntegrateFloat> Default for VerificationWorkflow<F> {
    fn default() -> Self {
        Self::new()
    }
}

#[allow(dead_code)]
pub fn exponential_solution<F: IntegrateFloat>(
    amplitude: F,
    decay_rate: F,
) -> ExponentialSolution<F> {
    ExponentialSolution::simple(amplitude, decay_rate)
}

#[allow(dead_code)]
pub fn trigonometric_solution_3d<F: IntegrateFloat>(
    freq_x: F,
    freq_y: F,
    freq_z: F,
) -> TrigonometricSolution3D<F> {
    TrigonometricSolution3D::simple(freq_x, freq_y, freq_z)
}

#[allow(dead_code)]
pub fn combined_solution<F: IntegrateFloat>(dimension: usize) -> CombinedSolution<F> {
    CombinedSolution::new(dimension)
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;
    use scirs2_core::ndarray::Array1;

    #[test]
    fn test_polynomial_solution() {
        // Test polynomial: 1 + 2t + 3t²
        let poly = polynomial_solution(vec![1.0, 2.0, 3.0]);

        // Test evaluation
        assert_abs_diff_eq!(poly.evaluate(&[0.0]), 1.0);
        assert_abs_diff_eq!(poly.evaluate(&[1.0]), 6.0); // 1 + 2 + 3

        // Test derivative: 2 + 6t
        assert_abs_diff_eq!(poly.derivative(&[0.0], 0), 2.0);
        assert_abs_diff_eq!(poly.derivative(&[1.0], 0), 8.0);

        // Test second derivative: 6
        assert_abs_diff_eq!(poly.second_derivative(&[0.0], 0), 6.0);
        assert_abs_diff_eq!(poly.second_derivative(&[1.0], 0), 6.0);
    }

    #[test]
    fn test_trigonometric_solution_2d() {
        // Test sin(x) * cos(y)
        let trig = trigonometric_solution_2d(1.0, 1.0);

        // Test evaluation
        assert_abs_diff_eq!(trig.evaluate(&[0.0, 0.0]), 0.0); // sin(0) * cos(0) = 0 * 1 = 0
        assert_abs_diff_eq!(trig.evaluate(&[PI / 2.0, 0.0]), 1.0); // sin(π/2) * cos(0) = 1 * 1 = 1

        // Test derivatives
        // ∂/∂x[sin(x)cos(y)] = cos(x)cos(y)
        assert_abs_diff_eq!(trig.derivative(&[0.0, 0.0], 0), 1.0); // cos(0)*cos(0) = 1

        // ∂/∂y[sin(x)cos(y)] = -sin(x)sin(y)
        assert_abs_diff_eq!(trig.derivative(&[0.0, 0.0], 1), 0.0); // -sin(0)*sin(0) = 0
    }

    #[test]
    fn test_mms_ode_problem() {
        // Test with y = t²
        let poly = polynomial_solution(vec![0.0, 0.0, 1.0]); // t²
        let problem = MMSODEProblem::new(poly, [0.0, 1.0]);

        // Initial condition should be 0
        assert_abs_diff_eq!(problem.initial_condition(), 0.0);

        // Source term should be y' = 2t
        assert_abs_diff_eq!(problem.source_term(0.0), 0.0);
        assert_abs_diff_eq!(problem.source_term(1.0), 2.0);

        // Exact solution at t=0.5 should be 0.25
        assert_abs_diff_eq!(problem.exact_at(0.5), 0.25);
    }

    #[test]
    fn test_convergence_analysis() {
        // Test with theoretical O(h²) convergence
        let grid_sizes = vec![0.1, 0.05, 0.025, 0.0125];
        let errors = vec![0.01, 0.0025, 0.000625, 0.00015625]; // ∝ h²

        let analysis =
            ConvergenceAnalysis::compute_order(grid_sizes, errors).expect("Operation failed");

        // Should detect order ≈ 2
        assert!((analysis.order - 2.0_f64).abs() < 0.1);
        assert!(analysis.verify_order(2.0, 0.2));
    }

    #[test]
    fn test_error_analysis() {
        let exact = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
        let numerical = Array1::from_vec(vec![1.1, 1.9, 3.05, 3.98]);

        let l2_error =
            ErrorAnalysis::l2_norm(exact.view(), numerical.view()).expect("Operation failed");
        let max_error =
            ErrorAnalysis::max_norm(exact.view(), numerical.view()).expect("Operation failed");

        // L2 error should be around sqrt((0.1² + 0.1² + 0.05² + 0.02²)/4)
        assert!(l2_error > 0.0 && l2_error < 0.2);

        // Max error should be 0.1
        assert_abs_diff_eq!(max_error, 0.1, epsilon = 1e-10);
    }

    #[test]
    fn test_exponential_solution() {
        // Test exponential solution: 2 * exp(-3t)
        let exp_sol = exponential_solution(2.0, -3.0);

        // Test evaluation
        assert_abs_diff_eq!(exp_sol.evaluate(&[0.0]), 2.0);
        assert_abs_diff_eq!(exp_sol.evaluate(&[1.0]), 2.0 * (-3.0_f64).exp());

        // Test derivative: d/dt[2 * exp(-3t)] = -6 * exp(-3t)
        assert_abs_diff_eq!(exp_sol.derivative(&[0.0], 0), -6.0);
        assert_abs_diff_eq!(exp_sol.derivative(&[1.0], 0), -6.0 * (-3.0_f64).exp());

        // Test second derivative: d²/dt²[2 * exp(-3t)] = 18 * exp(-3t)
        assert_abs_diff_eq!(exp_sol.second_derivative(&[0.0], 0), 18.0);
        assert_abs_diff_eq!(
            exp_sol.second_derivative(&[1.0], 0),
            18.0 * (-3.0_f64).exp()
        );
    }

    #[test]
    fn test_combined_solution() {
        // Test combined solution: polynomial + exponential
        let poly = polynomial_solution(vec![1.0, 2.0]); // 1 + 2t
        let exp = exponential_solution(1.0, -1.0); // exp(-t)

        let combined = combined_solution(1)
            .with_polynomial(poly)
            .with_exponential(exp);

        // At t=0: (1 + 2*0) + exp(-0) = 1 + 1 = 2
        assert_abs_diff_eq!(combined.evaluate(&[0.0]), 2.0);

        // At t=1: (1 + 2*1) + exp(-1) = 3 + exp(-1)
        let expected = 3.0 + (-1.0_f64).exp();
        assert_abs_diff_eq!(combined.evaluate(&[1.0]), expected, epsilon = 1e-10);

        // Test derivative at t=0: d/dt[1 + 2t + exp(-t)] = 2 - exp(-t)
        // At t=0: 2 - 1 = 1
        assert_abs_diff_eq!(combined.derivative(&[0.0], 0), 1.0);
    }

    #[test]
    fn test_trigonometric_solution_3d() {
        // Test sin(x) * cos(y) * sin(z)
        let trig3d = trigonometric_solution_3d(1.0, 1.0, 1.0);

        // Test evaluation at (π/2, 0, π/2)
        // sin(π/2) * cos(0) * sin(π/2) = 1 * 1 * 1 = 1
        assert_abs_diff_eq!(
            trig3d.evaluate(&[PI / 2.0, 0.0, PI / 2.0]),
            1.0,
            epsilon = 1e-10
        );

        // Test derivative with respect to x at (0, 0, π/2)
        // ∂/∂x[sin(x)cos(y)sin(z)] = cos(x)cos(y)sin(z)
        // cos(0)*cos(0)*sin(π/2) = 1*1*1 = 1
        assert_abs_diff_eq!(
            trig3d.derivative(&[0.0, 0.0, PI / 2.0], 0),
            1.0,
            epsilon = 1e-10
        );

        // Test second derivative with respect to x at (0, 0, π/2)
        // ∂²/∂x²[sin(x)cos(y)sin(z)] = -sin(x)cos(y)sin(z)
        // -sin(0)*cos(0)*sin(π/2) = 0
        assert_abs_diff_eq!(
            trig3d.second_derivative(&[0.0, 0.0, PI / 2.0], 0),
            0.0,
            epsilon = 1e-10
        );
    }

    #[test]
    fn test_3d_poisson_problem() {
        // Test 3D Poisson problem with trigonometric solution
        let exact = trigonometric_solution_3d(PI, PI, PI);
        let problem = MMSPDEProblem::new_poisson_3d(exact, [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]);

        // Test that it's recognized as 3D
        assert!(problem.is_3d());

        // Test domain
        let (dx, dy, dz) = problem.domain_3d();
        assert_eq!(dx, [0.0, 1.0]);
        assert_eq!(dy, [0.0, 1.0]);
        assert_eq!(dz, [0.0, 1.0]);

        // Test source term computation
        // For -∇²u = f where u = sin(πx)cos(πy)sin(πz)
        // ∇²u = -π²sin(πx)cos(πy)sin(πz) - π²sin(πx)cos(πy)sin(πz) - π²sin(πx)cos(πy)sin(πz)
        //     = -3π²sin(πx)cos(πy)sin(πz) = -3π²u
        // So f = -∇²u = 3π²u
        let coords = [0.5, 0.0, 0.5]; // Where cos(πy) = cos(0) = 1
        let u_exact = problem.exact_at_3d(coords[0], coords[1], coords[2]);
        let f_computed = problem.source_term(&coords);
        let f_expected = 3.0 * PI * PI * u_exact;
        assert_abs_diff_eq!(f_computed, f_expected, epsilon = 1e-10);
    }

    #[test]
    fn test_helmholtz_problem() {
        // Test 2D Helmholtz problem: ∇²u + k²u = f
        let exact = trigonometric_solution_2d(PI, PI);
        let k = 2.0;
        let problem = MMSPDEProblem::new_helmholtz_2d(exact, [0.0, 1.0], [0.0, 1.0], k);

        // For u = sin(πx)cos(πy), ∇²u = -2π²sin(πx)cos(πy) = -2π²u
        // So f = ∇²u + k²u = -2π²u + k²u = (k² - 2π²)u
        let coords = [0.5, 0.0]; // Where cos(πy) = cos(0) = 1
        let u_exact = problem.exact_at(coords[0], coords[1]);
        let f_computed = problem.source_term(&coords);
        let f_expected = (k * k - 2.0 * PI * PI) * u_exact;
        assert_abs_diff_eq!(f_computed, f_expected, epsilon = 1e-10);
    }

    #[test]
    fn test_verification_workflow() {
        let mut workflow = VerificationWorkflow::new();

        // Add a test case for second-order method
        let test_case = VerificationTestCase {
            name: "Second-order test".to_string(),
            expected_order: 2.0,
            order_tolerance: 0.1,
            grid_sizes: vec![0.1, 0.05, 0.025],
            expected_errors: None,
        };
        workflow.add_test_case(test_case);

        // Mock solver that produces second-order errors
        let mock_solver = |grid_sizes: &[f64]| -> IntegrateResult<f64> {
            let h = grid_sizes[0];
            Ok(0.1 * h * h) // O(h²) error
        };

        let results = workflow.run_verification(mock_solver);
        assert_eq!(results.len(), 1);
        assert!(results[0].passed);
        assert!(
            results[0]
                .computed_order
                .expect("Test: computed_order missing")
                > 1.8
        );
        assert!(
            results[0]
                .computed_order
                .expect("Test: computed_order missing")
                < 2.2
        );
    }

    #[test]
    fn test_system_verification() {
        // Test system with two components: polynomial and trigonometric
        let poly = polynomial_solution(vec![1.0, 2.0]);
        let trig = trigonometric_solution_2d(1.0, 1.0);

        let system = SystemVerification::<f64>::new(2);
        assert_eq!(system.system_size, 2);

        // Test verification at (0, 0) - handle solutions separately due to different types
        let coords = vec![0.0, 0.0];

        // Test polynomial solution
        let poly_exact = poly.evaluate(&coords);
        let poly_numerical = 1.0 + 2.0 * coords[0]; // approximate polynomial

        // Test trigonometric solution
        let trig_exact = trig.evaluate(&coords);
        let trig_numerical = coords[0] * coords[1]; // approximate trigonometric

        // Calculate errors manually since we can't use mixed-type arrays
        let poly_error = (poly_exact as f64 - poly_numerical).abs() as f64;
        let trig_error = (trig_exact as f64 - trig_numerical).abs() as f64;

        // At (0,0): exact poly = 1, numerical = 1, error = 0
        assert_abs_diff_eq!(poly_error, 0.0);
        // At (0,0): exact trig = 0, numerical = 0, error = 0
        assert_abs_diff_eq!(trig_error, 0.0);

        // Test with custom names
        let named_system: SystemVerification<f64> = SystemVerification::with_names(vec![
            "Polynomial".to_string(),
            "Trigonometric".to_string(),
        ]);
        assert_eq!(named_system.component_names[0], "Polynomial");
        assert_eq!(named_system.component_names[1], "Trigonometric");
    }
}

/// Advanced verification framework with automatic error estimation
pub struct AdvancedVerificationFramework {
    /// Grid refinement strategy
    pub refinement_strategy: RefinementStrategy,
    /// Error estimation method
    pub error_estimation_method: ErrorEstimationMethod,
    /// Convergence criteria
    pub convergence_criteria: ConvergenceCriteria,
    /// Statistical analysis settings
    pub statistical_analysis: bool,
}

/// Grid refinement strategies
#[derive(Debug, Clone, Copy)]
pub enum RefinementStrategy {
    /// Uniform refinement (h → h/2)
    Uniform,
    /// Adaptive refinement based on error indicators
    Adaptive,
    /// Custom refinement ratios
    Custom(f64),
}

/// Error estimation methods
#[derive(Debug, Clone, Copy)]
pub enum ErrorEstimationMethod {
    /// Richardson extrapolation
    Richardson,
    /// Embedded method comparison
    Embedded,
    /// Bootstrap sampling
    Bootstrap,
    /// Cross-validation
    CrossValidation,
}

/// Convergence criteria for verification
#[derive(Debug, Clone)]
pub struct ConvergenceCriteria {
    /// Minimum number of grid levels
    pub min_levels: usize,
    /// Maximum number of grid levels
    pub max_levels: usize,
    /// Required R² for order estimation
    pub min_r_squared: f64,
    /// Tolerance for order estimation
    pub order_tolerance: f64,
    /// Target accuracy
    pub target_accuracy: f64,
}

impl Default for ConvergenceCriteria {
    fn default() -> Self {
        Self {
            min_levels: 3,
            max_levels: 8,
            min_r_squared: 0.95,
            order_tolerance: 0.2,
            target_accuracy: 1e-10,
        }
    }
}

impl Default for AdvancedVerificationFramework {
    fn default() -> Self {
        Self {
            refinement_strategy: RefinementStrategy::Uniform,
            error_estimation_method: ErrorEstimationMethod::Richardson,
            convergence_criteria: ConvergenceCriteria::default(),
            statistical_analysis: true,
        }
    }
}