Entropyk/crates/solver/src/initializer.rs

//! Smart initialization heuristic for thermodynamic system solvers.
//!
//! This module provides [`SmartInitializer`], which generates physically
//! reasonable initial guesses for the solver state vector from source and sink
//! temperatures. It uses the Antoine equation to estimate saturation pressures
//! for common refrigerants without requiring an external fluid backend.
//!
//! # Algorithm
//!
//! 1. Estimate evaporator pressure: `P_evap = P_sat(T_source - ΔT_approach)`
//! 2. Estimate condenser pressure: `P_cond = P_sat(T_sink + ΔT_approach)`
//! 3. Clamp `P_evap` to `0.5 * P_critical` if it exceeds the critical pressure
//! 4. Fill the state vector with `[ṁ, P, h_default]` per edge, using circuit topology
//!
//! # Supported Fluids
//!
//! Built-in Antoine coefficients are provided for:
//! - R134a, R410A, R32, R744 (CO2), R290 (Propane)
//!
//! Unknown fluids return an explicit error; no pressure guesses are invented.
//!
//! # No-Allocation Guarantee
//!
//! [`SmartInitializer::populate_state`] writes to a pre-allocated `&mut [f64]`
//! slice and performs no heap allocation.

use entropyk_components::port::FluidId;
use entropyk_core::{Enthalpy, Pressure, Temperature};
use thiserror::Error;

use crate::system::System;
use serde::{Deserialize, Serialize};

// ─────────────────────────────────────────────────────────────────────────────
// Error types
// ─────────────────────────────────────────────────────────────────────────────

/// Errors that can occur during smart initialization.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum InitializerError {
    /// Source or sink temperature exceeds the critical temperature for the fluid.
    ///
    /// Antoine equation is not valid above the critical temperature. The caller
    /// should either use a different fluid or provide a manual initial state.
    #[error("Temperature {temp_celsius:.1}°C exceeds critical temperature for {fluid}")]
    TemperatureAboveCritical {
        /// Temperature that triggered the error (°C).
        temp_celsius: f64,
        /// Fluid identifier string.
        fluid: String,
    },

    /// The provided state slice length does not match the system state vector length.
    #[error("State slice length {actual} does not match system state vector length {expected}")]
    StateLengthMismatch {
        /// Expected length (from `system.state_vector_len()`).
        expected: usize,
        /// Actual length of the provided slice.
        actual: usize,
    },

    /// No Antoine coefficients are available for the configured fluid.
    #[error("No Antoine saturation-pressure coefficients are available for {fluid}")]
    UnsupportedFluid {
        /// Fluid identifier string.
        fluid: String,
    },
}

// ─────────────────────────────────────────────────────────────────────────────
// Antoine coefficients
// ─────────────────────────────────────────────────────────────────────────────

/// Antoine equation coefficients for saturation pressure estimation.
///
/// The Antoine equation (log₁₀ form) is:
///
/// ```text
/// log10(P_sat [Pa]) = A - B / (C + T [°C])
/// ```
///
/// Coefficients are tuned for the −40°C to +80°C range. Accuracy is within 5%
/// of NIST/CoolProp values — sufficient for initialization purposes.
#[derive(Debug, Clone, PartialEq)]
pub struct AntoineCoefficients {
    /// Antoine constant A (dimensionless, log₁₀ scale, Pa units).
    pub a: f64,
    /// Antoine constant B (°C).
    pub b: f64,
    /// Antoine constant C (°C offset).
    pub c: f64,
    /// Critical pressure of the fluid (Pa).
    pub p_critical_pa: f64,
}

impl AntoineCoefficients {
    /// Returns the built-in coefficients for the given fluid identifier string.
    ///
    /// Matching is case-insensitive. Returns `None` for unknown fluids.
    pub fn for_fluid(fluid_str: &str) -> Option<&'static AntoineCoefficients> {
        // Normalize: uppercase, strip dashes/spaces
        let normalized = fluid_str.to_uppercase().replace(['-', ' '], "");
        ANTOINE_TABLE
            .iter()
            .find(|(name, _)| *name == normalized.as_str())
            .map(|(_, coeffs)| coeffs)
    }
}

/// Compute saturation pressure (Pa) from temperature (°C) using Antoine equation.
///
/// `log10(P_sat [Pa]) = A - B / (C + T [°C])`
///
/// This is a pure arithmetic function with no heap allocation.
pub fn antoine_pressure(t_celsius: f64, coeffs: &AntoineCoefficients) -> f64 {
    let log10_p = coeffs.a - coeffs.b / (coeffs.c + t_celsius);
    10f64.powf(log10_p)
}

/// Built-in Antoine coefficient table for common refrigerants.
///
/// Coefficients valid for approximately −40°C to +80°C.
/// Accuracy: within 5% of NIST saturation pressure values.
///
/// Formula: `log10(P_sat [Pa]) = A - B / (C + T [°C])`
///
/// A values are derived from NIST reference saturation pressures:
/// - R134a: P_sat(0°C) = 292,800 Pa  → A = log10(292800) + 1766/243 = 12.739
/// - R410A: P_sat(0°C) = 798,000 Pa  → A = log10(798000) + 1885/243 = 13.659
/// - R32:   P_sat(0°C) = 810,000 Pa  → A = log10(810000) + 1780/243 = 13.233
/// - R744:  P_sat(20°C) = 5,730,000 Pa → A = log10(5730000) + 1347.8/293 = 11.357
/// - R290:  P_sat(0°C) = 474,000 Pa  → A = log10(474000) + 1656/243 = 12.491
///
/// | Fluid  | A (for Pa) | B       | C     | P_critical (Pa) |
/// |--------|------------|---------|-------|-----------------|
/// | R134a  | 12.739     | 1766.0  | 243.0 | 4,059,280       |
/// | R410A  | 13.659     | 1885.0  | 243.0 | 4,901,200       |
/// | R32    | 13.233     | 1780.0  | 243.0 | 5,782,000       |
/// | R744   | 11.357     | 1347.8  | 273.0 | 7,377,300       |
/// | R290   | 12.491     | 1656.0  | 243.0 | 4,247,200       |
static ANTOINE_TABLE: &[(&str, AntoineCoefficients)] = &[
    (
        "R134A",
        AntoineCoefficients {
            a: 12.739,
            b: 1766.0,
            c: 243.0,
            p_critical_pa: 4_059_280.0,
        },
    ),
    (
        "R410A",
        AntoineCoefficients {
            a: 13.659,
            b: 1885.0,
            c: 243.0,
            p_critical_pa: 4_901_200.0,
        },
    ),
    (
        "R32",
        AntoineCoefficients {
            a: 13.233,
            b: 1780.0,
            c: 243.0,
            p_critical_pa: 5_782_000.0,
        },
    ),
    (
        "R744",
        AntoineCoefficients {
            a: 11.357,
            b: 1347.8,
            c: 273.0,
            p_critical_pa: 7_377_300.0,
        },
    ),
    (
        "R290",
        AntoineCoefficients {
            a: 12.491,
            b: 1656.0,
            c: 243.0,
            p_critical_pa: 4_247_200.0,
        },
    ),
];

// ─────────────────────────────────────────────────────────────────────────────
// Initializer configuration
// ─────────────────────────────────────────────────────────────────────────────

/// Configuration for [`SmartInitializer`].
#[derive(Debug, Clone, PartialEq)]
pub struct InitializerConfig {
    /// Fluid identifier used for Antoine coefficient lookup.
    pub fluid: FluidId,

    /// Temperature approach difference for pressure estimation (K).
    ///
    /// - Evaporator: `P_evap = P_sat(T_source - dt_approach)`
    /// - Condenser:  `P_cond = P_sat(T_sink + dt_approach)`
    ///
    /// Default: 5.0 K.
    pub dt_approach: f64,
}

/// Optional start values for one solver edge or auxiliary unknown group.
///
/// These values are numerical guesses only. They must never be interpreted as
/// imposed boundary conditions or component equations.
#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct StartValues {
    /// Pressure start value [Pa].
    pub pressure_pa: Option<f64>,
    /// Enthalpy start value [J/kg].
    pub enthalpy_j_kg: Option<f64>,
    /// Mass-flow start value [kg/s].
    pub mass_flow_kg_s: Option<f64>,
    /// Temperature start value [K], useful for diagnostics and backend conversion.
    pub temperature_k: Option<f64>,
    /// Vapour quality start value [-], when the intended regime is two-phase.
    pub vapor_quality: Option<f64>,
}

/// Regime label used to explain why a start value was assigned.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
#[serde(rename_all = "snake_case")]
pub enum InitializationRegime {
    /// High-pressure superheated vapour, typically compressor discharge.
    HighPressureVapor,
    /// High-pressure liquid, typically condenser outlet.
    HighPressureLiquid,
    /// Low-pressure two-phase mixture, typically EXV outlet.
    LowPressureTwoPhase,
    /// Low-pressure superheated vapour, typically compressor suction.
    LowPressureVapor,
    /// Secondary water/brine/air branch.
    Secondary,
    /// Generic fallback seed.
    Generic,
    /// Boundary condition seed from a source/sink component.
    Boundary,
    /// Control or actuator unknown seed.
    Control,
}

/// One initialization diagnostic entry.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct InitializationSeed {
    /// Human-readable edge/control label.
    pub label: String,
    /// Assigned regime.
    pub regime: InitializationRegime,
    /// Values written to the state vector.
    pub values: StartValues,
}

/// Diagnostics emitted by an initialization pass.
#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct InitializationDiagnostics {
    /// Ordered seed records for edges and control variables.
    pub seeds: Vec<InitializationSeed>,
}

impl InitializationDiagnostics {
    /// Appends a diagnostic seed record.
    pub fn push(
        &mut self,
        label: impl Into<String>,
        regime: InitializationRegime,
        values: StartValues,
    ) {
        self.seeds.push(InitializationSeed {
            label: label.into(),
            regime,
            values,
        });
    }
}

impl Default for InitializerConfig {
    fn default() -> Self {
        Self {
            fluid: FluidId::new("R134a"),
            dt_approach: 5.0,
        }
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// SmartInitializer
// ─────────────────────────────────────────────────────────────────────────────

/// Smart initialization heuristic for thermodynamic solver state vectors.
///
/// Uses the Antoine equation to estimate saturation pressures from source and
/// sink temperatures, then fills a pre-allocated state vector with physically
/// reasonable initial guesses.
///
/// # Example
///
/// ```rust,no_run
/// use entropyk_solver::initializer::{SmartInitializer, InitializerConfig};
/// use entropyk_core::{Temperature, Enthalpy};
///
/// let init = SmartInitializer::new(InitializerConfig::default());
/// let (p_evap, p_cond) = init
///     .estimate_pressures(
///         Temperature::from_celsius(5.0),
///         Temperature::from_celsius(40.0),
///     )
///     .unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct SmartInitializer {
    /// Configuration for this initializer.
    pub config: InitializerConfig,
}

impl SmartInitializer {
    /// Creates a new `SmartInitializer` with the given configuration.
    pub fn new(config: InitializerConfig) -> Self {
        Self { config }
    }

    /// Estimate `(P_evap, P_cond)` from source and sink temperatures.
    ///
    /// Uses the Antoine equation with the configured fluid and approach ΔT:
    /// - `P_evap = P_sat(T_source - ΔT_approach)`, clamped to `0.5 * P_critical`
    /// - `P_cond = P_sat(T_sink + ΔT_approach)`
    ///
    /// # Errors
    ///
    /// Returns [`InitializerError::TemperatureAboveCritical`] if the adjusted
    /// source temperature exceeds the critical temperature for a known fluid.
    /// Returns [`InitializerError::UnsupportedFluid`] if no Antoine coefficients
    /// are available for the configured fluid.
    pub fn estimate_pressures(
        &self,
        t_source: Temperature,
        t_sink: Temperature,
    ) -> Result<(Pressure, Pressure), InitializerError> {
        let fluid_str = self.config.fluid.to_string();

        match AntoineCoefficients::for_fluid(&fluid_str) {
            None => Err(InitializerError::UnsupportedFluid { fluid: fluid_str }),
            Some(coeffs) => {
                let t_source_c = t_source.to_celsius();
                let t_sink_c = t_sink.to_celsius();

                // Evaporator: T_source - ΔT_approach
                let t_evap_c = t_source_c - self.config.dt_approach;
                let p_evap_pa = antoine_pressure(t_evap_c, coeffs);

                // Clamp P_evap to 0.5 * P_critical (AC: #2)
                let p_evap_pa = if p_evap_pa >= coeffs.p_critical_pa {
                    tracing::warn!(
                        fluid = %fluid_str,
                        t_evap_celsius = t_evap_c,
                        p_evap_pa = p_evap_pa,
                        p_critical_pa = coeffs.p_critical_pa,
                        "Estimated P_evap exceeds critical pressure — clamping to 0.5 * P_critical"
                    );
                    0.5 * coeffs.p_critical_pa
                } else {
                    p_evap_pa
                };

                // Condenser: T_sink + ΔT_approach (AC: #3)
                let t_cond_c = t_sink_c + self.config.dt_approach;
                let p_cond_pa = antoine_pressure(t_cond_c, coeffs);

                // Clamp P_cond to 0.5 * P_critical if it exceeds critical
                let p_cond_pa = if p_cond_pa >= coeffs.p_critical_pa {
                    tracing::warn!(
                        fluid = %fluid_str,
                        t_cond_celsius = t_cond_c,
                        p_cond_pa = p_cond_pa,
                        p_critical_pa = coeffs.p_critical_pa,
                        "Estimated P_cond exceeds critical pressure — clamping to 0.5 * P_critical"
                    );
                    0.5 * coeffs.p_critical_pa
                } else {
                    p_cond_pa
                };

                tracing::debug!(
                    fluid = %fluid_str,
                    t_source_celsius = t_source_c,
                    t_sink_celsius = t_sink_c,
                    p_evap_bar = p_evap_pa / 1e5,
                    p_cond_bar = p_cond_pa / 1e5,
                    "SmartInitializer: estimated pressures"
                );

                Ok((
                    Pressure::from_pascals(p_evap_pa),
                    Pressure::from_pascals(p_cond_pa),
                ))
            }
        }
    }

    /// Fill a pre-allocated state vector with smart initial guesses.
    ///
    /// No heap allocation is performed. The `state` slice must have length equal
    /// to `system.state_vector_len()` (i.e., `3 * edge_count` for a system of
    /// refrigerant/hydraulic edges).
    ///
    /// State layout per edge: `[ṁ_edge_i, P_edge_i, h_edge_i]`
    ///
    /// Pressure assignment follows circuit topology:
    /// - Edges in circuit 0 → `p_evap`
    /// - Edges in circuit 1+ → `p_cond`
    /// - Single-circuit systems: all edges use `p_evap`
    ///
    /// # Errors
    ///
    /// Returns [`InitializerError::StateLengthMismatch`] if `state.len()` does
    /// not match `system.full_state_vector_len()` (edges plus any inverse-control
    /// and coupling auxiliary unknowns).
    pub fn populate_state(
        &self,
        system: &System,
        p_evap: Pressure,
        p_cond: Pressure,
        h_default: Enthalpy,
        state: &mut [f64],
    ) -> Result<(), InitializerError> {
        // Size against the FULL state vector (the length Newton/Picard expect):
        // base edge unknowns + inverse-control mappings + coupling residual slots.
        let expected = system.full_state_vector_len();
        if state.len() != expected {
            return Err(InitializerError::StateLengthMismatch {
                expected,
                actual: state.len(),
            });
        }

        let p_evap_pa = p_evap.to_pascals();
        let p_cond_pa = p_cond.to_pascals();
        let h_jkg = h_default.to_joules_per_kg();

        for edge_idx in system.edge_indices() {
            let circuit = system.edge_circuit(edge_idx);
            let p = if circuit.0 == 0 { p_evap_pa } else { p_cond_pa };
            let (m_idx, p_idx, h_idx) = system.edge_state_indices_full(edge_idx);
            // CM1.4: m_idx is BRANCH-shared — multiple edges in the same series
            // branch point to the same slot. Writing the same seed value multiple
            // times is idempotent and stays within state bounds (m_idx < state_len).
            state[m_idx] = crate::system::DEFAULT_MASS_FLOW_SEED_KG_S;
            state[p_idx] = p;
            state[h_idx] = h_jkg;
        }

        // Seed inverse-control unknowns (fan speed, opening, frequency, …) to the
        // midpoint of their bounds so the cold start sits inside the feasible box
        // instead of at zero (often an out-of-bounds, non-physical control value).
        // Coupling auxiliary slots keep their 0.0 default.
        for (_, idx) in system.control_variable_indices() {
            if let Some((min, max)) = system.get_bounds_for_state_index(idx) {
                if min.is_finite() && max.is_finite() && min <= max {
                    state[idx] = 0.5 * (min + max);
                }
            }
        }

        Ok(())
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────

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

    #[test]
    fn test_initialization_diagnostics_records_start_values() {
        let mut diagnostics = InitializationDiagnostics::default();
        diagnostics.push(
            "comp->cond",
            InitializationRegime::HighPressureVapor,
            StartValues {
                pressure_pa: Some(1.2e6),
                enthalpy_j_kg: Some(430_000.0),
                mass_flow_kg_s: Some(0.05),
                temperature_k: None,
                vapor_quality: None,
            },
        );

        assert_eq!(diagnostics.seeds.len(), 1);
        assert_eq!(diagnostics.seeds[0].label, "comp->cond");
        assert_eq!(
            diagnostics.seeds[0].regime,
            InitializationRegime::HighPressureVapor
        );
        assert_eq!(diagnostics.seeds[0].values.pressure_pa, Some(1.2e6));
    }

    // ── Antoine equation unit tests ──────────────────────────────────────────

    /// AC: #1, #5 — R134a at 0°C: P_sat ≈ 2.93 bar (293,000 Pa), within 5%
    #[test]
    fn test_antoine_r134a_at_0c() {
        let coeffs = AntoineCoefficients::for_fluid("R134a").unwrap();
        let p_pa = antoine_pressure(0.0, coeffs);
        // Expected: ~2.93 bar = 293,000 Pa
        assert_relative_eq!(p_pa, 293_000.0, max_relative = 0.05);
    }

    /// AC: #5 — R744 (CO2) at 20°C: P_sat ≈ 57.3 bar (5,730,000 Pa), within 5%
    #[test]
    fn test_antoine_r744_at_20c() {
        let coeffs = AntoineCoefficients::for_fluid("R744").unwrap();
        let p_pa = antoine_pressure(20.0, coeffs);
        // Expected: ~57.3 bar = 5,730,000 Pa
        assert_relative_eq!(p_pa, 5_730_000.0, max_relative = 0.05);
    }

    /// AC: #5 — Case-insensitive fluid lookup
    #[test]
    fn test_fluid_lookup_case_insensitive() {
        assert!(AntoineCoefficients::for_fluid("r134a").is_some());
        assert!(AntoineCoefficients::for_fluid("R134A").is_some());
        assert!(AntoineCoefficients::for_fluid("R134a").is_some());
        assert!(AntoineCoefficients::for_fluid("r744").is_some());
        assert!(AntoineCoefficients::for_fluid("R290").is_some());
    }

    /// AC: #5 — Unknown fluid returns None
    #[test]
    fn test_fluid_lookup_unknown() {
        assert!(AntoineCoefficients::for_fluid("R999").is_none());
        assert!(AntoineCoefficients::for_fluid("").is_none());
    }

    // ── SmartInitializer::estimate_pressures tests ───────────────────────────

    /// AC: #2 — P_evap < P_critical for all built-in fluids at T_source = −40°C
    #[test]
    fn test_p_evap_below_critical_all_fluids() {
        let fluids = ["R134a", "R410A", "R32", "R744", "R290"];
        for fluid in fluids {
            let init = SmartInitializer::new(InitializerConfig {
                fluid: FluidId::new(fluid),
                dt_approach: 5.0,
            });
            let (p_evap, _) = init
                .estimate_pressures(
                    Temperature::from_celsius(-40.0),
                    Temperature::from_celsius(40.0),
                )
                .unwrap();
            let coeffs = AntoineCoefficients::for_fluid(fluid).unwrap();
            assert!(
                p_evap.to_pascals() < coeffs.p_critical_pa,
                "P_evap ({:.0} Pa) should be < P_critical ({:.0} Pa) for {}",
                p_evap.to_pascals(),
                coeffs.p_critical_pa,
                fluid
            );
        }
    }

    /// AC: #3 — P_cond = P_sat(T_sink + 5K) for default ΔT_approach
    #[test]
    fn test_p_cond_approach_default() {
        let init = SmartInitializer::new(InitializerConfig::default()); // R134a, dt=5.0
        let t_sink = Temperature::from_celsius(40.0);
        let (_, p_cond) = init
            .estimate_pressures(Temperature::from_celsius(5.0), t_sink)
            .unwrap();

        // Expected: P_sat(45°C) for R134a
        let coeffs = AntoineCoefficients::for_fluid("R134a").unwrap();
        let expected_pa = antoine_pressure(45.0, coeffs);
        assert_relative_eq!(p_cond.to_pascals(), expected_pa, max_relative = 1e-9);
    }

    /// AC: #6 — Unknown fluid returns an explicit error instead of invented pressures.
    #[test]
    fn test_unknown_fluid_returns_error() {
        let init = SmartInitializer::new(InitializerConfig {
            fluid: FluidId::new("R999-Unknown"),
            dt_approach: 5.0,
        });
        let result = init.estimate_pressures(
            Temperature::from_celsius(5.0),
            Temperature::from_celsius(40.0),
        );
        assert_eq!(
            result,
            Err(InitializerError::UnsupportedFluid {
                fluid: "R999-Unknown".to_string()
            })
        );
    }

    /// AC: #1 — Verify evaporator pressure uses T_source - ΔT_approach
    #[test]
    fn test_p_evap_uses_approach_delta() {
        let dt = 5.0;
        let init = SmartInitializer::new(InitializerConfig {
            fluid: FluidId::new("R134a"),
            dt_approach: dt,
        });
        let t_source = Temperature::from_celsius(10.0);
        let (p_evap, _) = init
            .estimate_pressures(t_source, Temperature::from_celsius(40.0))
            .unwrap();

        let coeffs = AntoineCoefficients::for_fluid("R134a").unwrap();
        let expected_pa = antoine_pressure(10.0 - dt, coeffs); // T_source - ΔT
        assert_relative_eq!(p_evap.to_pascals(), expected_pa, max_relative = 1e-9);
    }

    // ── SmartInitializer::populate_state tests ───────────────────────────────

    /// AC: #4, #7 — populate_state fills state vector correctly for a 2-edge system.
    ///
    /// This test verifies the no-allocation signature: the function takes `&mut [f64]`
    /// and writes in-place without allocating.
    #[test]
    fn test_populate_state_2_edges() {
        use crate::system::System;
        use entropyk_components::{
            Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, StateSlice,
        };

        struct MockComp;
        impl Component for MockComp {
            fn compute_residuals(
                &self,
                _s: &StateSlice,
                r: &mut ResidualVector,
            ) -> Result<(), ComponentError> {
                for v in r.iter_mut() {
                    *v = 0.0;
                }
                Ok(())
            }
            fn jacobian_entries(
                &self,
                _s: &StateSlice,
                j: &mut JacobianBuilder,
            ) -> Result<(), ComponentError> {
                j.add_entry(0, 0, 1.0);
                Ok(())
            }
            fn n_equations(&self) -> usize {
                1
            }
            fn get_ports(&self) -> &[ConnectedPort] {
                &[]
            }
        }

        let mut sys = System::new();
        let n0 = sys.add_component(Box::new(MockComp));
        let n1 = sys.add_component(Box::new(MockComp));
        let n2 = sys.add_component(Box::new(MockComp));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.finalize().unwrap();

        let init = SmartInitializer::new(InitializerConfig::default());
        let p_evap = Pressure::from_bar(3.0);
        let p_cond = Pressure::from_bar(15.0);
        let h_default = Enthalpy::from_joules_per_kg(400_000.0);

        // Pre-allocated slice — no allocation in populate_state
        let mut state = vec![0.0f64; sys.state_vector_len()];
        init.populate_state(&sys, p_evap, p_cond, h_default, &mut state)
            .unwrap();

        // CM1.4: 2-edge linear chain → 1 branch → state_len = 1 + 2×2 = 5
        // Layout: [0:ṁ_branch, 1:P_e0, 2:h_e0, 3:P_e1, 4:h_e1]
        assert_eq!(state.len(), 5);
        // Branch ṁ seeded at DEFAULT_MASS_FLOW_SEED_KG_S
        assert_relative_eq!(
            state[0],
            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
            max_relative = 1e-9
        );
        // P and h for edge 0 (circuit 0 → p_evap)
        assert_relative_eq!(state[1], p_evap.to_pascals(), max_relative = 1e-9);
        assert_relative_eq!(state[2], h_default.to_joules_per_kg(), max_relative = 1e-9);
        // P and h for edge 1 (circuit 0 → p_evap)
        assert_relative_eq!(state[3], p_evap.to_pascals(), max_relative = 1e-9);
        assert_relative_eq!(state[4], h_default.to_joules_per_kg(), max_relative = 1e-9);
    }

    /// AC: #4 — populate_state uses P_cond for circuit 1 edges in multi-circuit system.
    #[test]
    fn test_populate_state_multi_circuit() {
        use crate::system::System;
        use entropyk_components::{
            Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, StateSlice,
        };
        use entropyk_core::CircuitId;

        struct MockComp;
        impl Component for MockComp {
            fn compute_residuals(
                &self,
                _s: &StateSlice,
                r: &mut ResidualVector,
            ) -> Result<(), ComponentError> {
                for v in r.iter_mut() {
                    *v = 0.0;
                }
                Ok(())
            }
            fn jacobian_entries(
                &self,
                _s: &StateSlice,
                j: &mut JacobianBuilder,
            ) -> Result<(), ComponentError> {
                j.add_entry(0, 0, 1.0);
                Ok(())
            }
            fn n_equations(&self) -> usize {
                1
            }
            fn get_ports(&self) -> &[ConnectedPort] {
                &[]
            }
        }

        let mut sys = System::new();
        // Circuit 0: evaporator side
        let n0 = sys
            .add_component_to_circuit(Box::new(MockComp), CircuitId(0))
            .unwrap();
        let n1 = sys
            .add_component_to_circuit(Box::new(MockComp), CircuitId(0))
            .unwrap();
        // Circuit 1: condenser side
        let n2 = sys
            .add_component_to_circuit(Box::new(MockComp), CircuitId(1))
            .unwrap();
        let n3 = sys
            .add_component_to_circuit(Box::new(MockComp), CircuitId(1))
            .unwrap();

        sys.add_edge(n0, n1).unwrap(); // circuit 0 edge
        sys.add_edge(n2, n3).unwrap(); // circuit 1 edge
        sys.finalize().unwrap();

        let init = SmartInitializer::new(InitializerConfig::default());
        let p_evap = Pressure::from_bar(3.0);
        let p_cond = Pressure::from_bar(15.0);
        let h_default = Enthalpy::from_joules_per_kg(400_000.0);

        let mut state = vec![0.0f64; sys.state_vector_len()];
        init.populate_state(&sys, p_evap, p_cond, h_default, &mut state)
            .unwrap();

        // CM1.4: 2 isolated 1-edge chains → 2 branches → state_len = 2 + 2×2 = 6
        // Layout: [0:ṁ_B0, 1:ṁ_B1, 2:P_e0, 3:h_e0, 4:P_e1, 5:h_e1]
        assert_eq!(state.len(), 6);
        // Branch ṁ slots seeded at DEFAULT_MASS_FLOW_SEED_KG_S
        assert_relative_eq!(
            state[0],
            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
            max_relative = 1e-9
        );
        assert_relative_eq!(
            state[1],
            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
            max_relative = 1e-9
        );
        // Verify P and h values are seeded correctly for each edge using actual state indices.
        // Edge 0 is circuit 0 (p_evap), edge 1 is circuit 1 (p_cond).
        for edge_idx in sys.edge_indices() {
            let circuit = sys.edge_circuit(edge_idx);
            let (_m, p, h) = sys.edge_state_indices_full(edge_idx);
            let expected_p = if circuit.0 == 0 {
                p_evap.to_pascals()
            } else {
                p_cond.to_pascals()
            };
            assert_relative_eq!(state[p], expected_p, max_relative = 1e-9);
            assert_relative_eq!(state[h], h_default.to_joules_per_kg(), max_relative = 1e-9);
        }
    }

    /// AC: #7 — populate_state returns error on length mismatch (no panic).
    #[test]
    fn test_populate_state_length_mismatch() {
        use crate::system::System;
        use entropyk_components::{
            Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, StateSlice,
        };

        struct MockComp;
        impl Component for MockComp {
            fn compute_residuals(
                &self,
                _s: &StateSlice,
                r: &mut ResidualVector,
            ) -> Result<(), ComponentError> {
                for v in r.iter_mut() {
                    *v = 0.0;
                }
                Ok(())
            }
            fn jacobian_entries(
                &self,
                _s: &StateSlice,
                j: &mut JacobianBuilder,
            ) -> Result<(), ComponentError> {
                j.add_entry(0, 0, 1.0);
                Ok(())
            }
            fn n_equations(&self) -> usize {
                1
            }
            fn get_ports(&self) -> &[ConnectedPort] {
                &[]
            }
        }

        let mut sys = System::new();
        let n0 = sys.add_component(Box::new(MockComp));
        let n1 = sys.add_component(Box::new(MockComp));
        sys.add_edge(n0, n1).unwrap();
        sys.finalize().unwrap();

        let init = SmartInitializer::new(InitializerConfig::default());
        let p_evap = Pressure::from_bar(3.0);
        let p_cond = Pressure::from_bar(15.0);
        let h_default = Enthalpy::from_joules_per_kg(400_000.0);

        // Wrong length: system has 3 state entries (1 edge × 3), we provide 5
        let mut state = vec![0.0f64; 5];
        let result = init.populate_state(&sys, p_evap, p_cond, h_default, &mut state);
        assert!(matches!(
            result,
            Err(InitializerError::StateLengthMismatch {
                expected: 3,
                actual: 5
            })
        ));
    }

    /// AC: #2 — P_evap is clamped to 0.5 * P_critical when above critical.
    ///
    /// We use R744 (CO2) at a very high source temperature to trigger clamping.
    #[test]
    fn test_p_evap_clamped_above_critical() {
        // R744 critical: 7,377,300 Pa (~73.8 bar), critical T ≈ 31°C
        // At T_source = 40°C, T_evap = 35°C → P_sat > P_critical → should clamp
        let init = SmartInitializer::new(InitializerConfig {
            fluid: FluidId::new("R744"),
            dt_approach: 5.0,
        });
        let (p_evap, _) = init
            .estimate_pressures(
                Temperature::from_celsius(40.0),
                Temperature::from_celsius(50.0),
            )
            .unwrap();

        let coeffs = AntoineCoefficients::for_fluid("R744").unwrap();
        // Must be clamped to 0.5 * P_critical
        assert_relative_eq!(
            p_evap.to_pascals(),
            0.5 * coeffs.p_critical_pa,
            max_relative = 1e-9
        );
    }
}