feat(components): add ThermoState generators and Eurovent backend demo

crates/solver/Cargo.toml

@@ -0,0 +1,23 @@
+[package]
+name = "entropyk-solver"
+version = "0.1.0"
+edition = "2021"
+authors = ["Sepehr <sepehr@entropyk.com>"]
+description = "System topology and solver engine for Entropyk thermodynamic simulation"
+license = "MIT OR Apache-2.0"
+repository = "https://github.com/entropyk/entropyk"
+
+[dependencies]
+entropyk-components = { path = "../components" }
+entropyk-core = { path = "../core" }
+nalgebra = "0.33"
+petgraph = "0.6"
+thiserror = "1.0"
+tracing = "0.1"
+
+[dev-dependencies]
+approx = "0.5"
+
+[lib]
+name = "entropyk_solver"
+path = "src/lib.rs"

crates/solver/src/coupling.rs

@@ -0,0 +1,435 @@
+//! Thermal coupling between circuits for heat transfer.
+//!
+//! This module provides the infrastructure for modeling heat exchange between
+//! independent fluid circuits. Thermal couplings represent heat exchangers
+//! that transfer heat from a "hot" circuit to a "cold" circuit without
+//! fluid mixing.
+//!
+//! ## Sign Convention
+//!
+//! Heat transfer Q > 0 means heat flows INTO the cold circuit (out of hot circuit).
+//! This follows the convention that the cold circuit receives heat.
+//!
+//! ## Coupling Graph and Circular Dependencies
+//!
+//! Thermal couplings form a directed graph where:
+//! - Nodes are circuits (CircuitId)
+//! - Edges point from hot_circuit to cold_circuit (direction of heat flow)
+//!
+//! Circular dependencies occur when circuits mutually heat each other (A→B and B→A).
+//! Circuits in circular dependencies must be solved simultaneously by the solver.
+
+use entropyk_core::{Temperature, ThermalConductance};
+use petgraph::algo::{is_cyclic_directed, kosaraju_scc};
+use petgraph::graph::{DiGraph, NodeIndex};
+use std::collections::HashMap;
+
+use crate::system::CircuitId;
+
+/// Thermal coupling between two circuits via a heat exchanger.
+///
+/// Heat flows from `hot_circuit` to `cold_circuit` proportional to the
+/// temperature difference and thermal conductance (UA value).
+#[derive(Debug, Clone, PartialEq)]
+pub struct ThermalCoupling {
+    /// Circuit that supplies heat (higher temperature side).
+    pub hot_circuit: CircuitId,
+    /// Circuit that receives heat (lower temperature side).
+    pub cold_circuit: CircuitId,
+    /// Thermal conductance (UA) in W/K. Higher values = more heat transfer.
+    pub ua: ThermalConductance,
+    /// Efficiency factor (0.0 to 1.0). Default is 1.0 (no losses).
+    pub efficiency: f64,
+}
+
+impl ThermalCoupling {
+    /// Creates a new thermal coupling between two circuits.
+    ///
+    /// # Arguments
+    ///
+    /// * `hot_circuit` - Circuit at higher temperature (heat source)
+    /// * `cold_circuit` - Circuit at lower temperature (heat sink)
+    /// * `ua` - Thermal conductance in W/K
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use entropyk_solver::{ThermalCoupling, CircuitId};
+    /// use entropyk_core::ThermalConductance;
+    ///
+    /// let coupling = ThermalCoupling::new(
+    ///     CircuitId(0),
+    ///     CircuitId(1),
+    ///     ThermalConductance::from_watts_per_kelvin(1000.0),
+    /// );
+    /// ```
+    pub fn new(hot_circuit: CircuitId, cold_circuit: CircuitId, ua: ThermalConductance) -> Self {
+        Self {
+            hot_circuit,
+            cold_circuit,
+            ua,
+            efficiency: 1.0,
+        }
+    }
+
+    /// Sets the efficiency factor for the coupling.
+    ///
+    /// Efficiency accounts for heat losses in the heat exchanger.
+    /// A value of 0.9 means 90% of theoretical heat is transferred.
+    pub fn with_efficiency(mut self, efficiency: f64) -> Self {
+        self.efficiency = efficiency.clamp(0.0, 1.0);
+        self
+    }
+}
+
+/// Computes heat transfer for a thermal coupling.
+///
+/// # Formula
+///
+/// Q = η × UA × (T_hot - T_cold)
+///
+/// Where:
+/// - Q is the heat transfer rate (W), positive means heat INTO cold circuit
+/// - η is the efficiency factor
+/// - UA is the thermal conductance (W/K)
+/// - T_hot, T_cold are temperatures (K)
+///
+/// # Sign Convention
+///
+/// - Q > 0: Heat flows from hot to cold (normal operation)
+/// - Q = 0: No temperature difference
+/// - Q < 0: Cold is hotter than hot (reverse flow, unusual)
+///
+/// # Example
+///
+/// ```
+/// use entropyk_solver::{ThermalCoupling, CircuitId, compute_coupling_heat};
+/// use entropyk_core::{Temperature, ThermalConductance};
+///
+/// let coupling = ThermalCoupling::new(
+///     CircuitId(0),
+///     CircuitId(1),
+///     ThermalConductance::from_watts_per_kelvin(1000.0),
+/// );
+///
+/// let t_hot = Temperature::from_kelvin(350.0);
+/// let t_cold = Temperature::from_kelvin(300.0);
+///
+/// let q = compute_coupling_heat(&coupling, t_hot, t_cold);
+/// assert!(q > 0.0, "Heat should flow from hot to cold");
+/// ```
+pub fn compute_coupling_heat(
+    coupling: &ThermalCoupling,
+    t_hot: Temperature,
+    t_cold: Temperature,
+) -> f64 {
+    coupling.efficiency
+        * coupling.ua.to_watts_per_kelvin()
+        * (t_hot.to_kelvin() - t_cold.to_kelvin())
+}
+
+/// Builds a coupling graph for dependency analysis.
+///
+/// Returns a directed graph where:
+/// - Nodes are CircuitIds present in any coupling
+/// - Edges point from hot_circuit to cold_circuit
+fn build_coupling_graph(couplings: &[ThermalCoupling]) -> DiGraph<CircuitId, ()> {
+    let mut graph = DiGraph::new();
+    let mut circuit_to_node: HashMap<CircuitId, NodeIndex> = HashMap::new();
+
+    for coupling in couplings {
+        // Add hot_circuit node if not present
+        let hot_node = *circuit_to_node
+            .entry(coupling.hot_circuit)
+            .or_insert_with(|| graph.add_node(coupling.hot_circuit));
+
+        // Add cold_circuit node if not present
+        let cold_node = *circuit_to_node
+            .entry(coupling.cold_circuit)
+            .or_insert_with(|| graph.add_node(coupling.cold_circuit));
+
+        // Add directed edge: hot -> cold
+        graph.add_edge(hot_node, cold_node, ());
+    }
+
+    graph
+}
+
+/// Checks if the coupling graph contains circular dependencies.
+///
+/// Circular dependencies occur when circuits are mutually thermally coupled
+/// (e.g., A heats B, and B heats A). When circular dependencies exist,
+/// the solver must solve those circuits simultaneously rather than sequentially.
+///
+/// # Example
+///
+/// ```
+/// use entropyk_solver::{ThermalCoupling, CircuitId, has_circular_dependencies};
+/// use entropyk_core::ThermalConductance;
+///
+/// // No circular dependency: A → B → C
+/// let couplings = vec![
+///     ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
+///     ThermalCoupling::new(CircuitId(1), CircuitId(2), ThermalConductance::from_watts_per_kelvin(100.0)),
+/// ];
+/// assert!(!has_circular_dependencies(&couplings));
+///
+/// // Circular dependency: A → B and B → A
+/// let couplings_circular = vec![
+///     ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
+///     ThermalCoupling::new(CircuitId(1), CircuitId(0), ThermalConductance::from_watts_per_kelvin(100.0)),
+/// ];
+/// assert!(has_circular_dependencies(&couplings_circular));
+/// ```
+pub fn has_circular_dependencies(couplings: &[ThermalCoupling]) -> bool {
+    if couplings.is_empty() {
+        return false;
+    }
+    let graph = build_coupling_graph(couplings);
+    is_cyclic_directed(&graph)
+}
+
+/// Returns groups of circuits that must be solved simultaneously.
+///
+/// Groups are computed using strongly connected components (SCC) analysis
+/// of the coupling graph. Circuits in the same SCC have circular thermal
+/// dependencies and must be solved together.
+///
+/// # Returns
+///
+/// A vector of vectors, where each inner vector contains CircuitIds that
+/// must be solved simultaneously. Single-element vectors indicate circuits
+/// that can be solved independently (in topological order).
+///
+/// # Example
+///
+/// ```
+/// use entropyk_solver::{ThermalCoupling, CircuitId, coupling_groups};
+/// use entropyk_core::ThermalConductance;
+///
+/// // A → B, B and C independent
+/// let couplings = vec![
+///     ThermalCoupling::new(CircuitId(0), CircuitId(1), ThermalConductance::from_watts_per_kelvin(100.0)),
+/// ];
+/// let groups = coupling_groups(&couplings);
+/// // Groups will contain individual circuits since there's no cycle
+/// ```
+pub fn coupling_groups(couplings: &[ThermalCoupling]) -> Vec<Vec<CircuitId>> {
+    if couplings.is_empty() {
+        return Vec::new();
+    }
+
+    let graph = build_coupling_graph(couplings);
+    let sccs = kosaraju_scc(&graph);
+
+    sccs.into_iter()
+        .map(|node_indices| node_indices.into_iter().map(|idx| graph[idx]).collect())
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    fn make_coupling(hot: u8, cold: u8, ua_w_per_k: f64) -> ThermalCoupling {
+        ThermalCoupling::new(
+            CircuitId(hot),
+            CircuitId(cold),
+            ThermalConductance::from_watts_per_kelvin(ua_w_per_k),
+        )
+    }
+
+    #[test]
+    fn test_thermal_coupling_creation() {
+        let coupling = ThermalCoupling::new(
+            CircuitId(0),
+            CircuitId(1),
+            ThermalConductance::from_watts_per_kelvin(1000.0),
+        );
+
+        assert_eq!(coupling.hot_circuit, CircuitId(0));
+        assert_eq!(coupling.cold_circuit, CircuitId(1));
+        assert_relative_eq!(coupling.ua.to_watts_per_kelvin(), 1000.0, epsilon = 1e-10);
+        assert_relative_eq!(coupling.efficiency, 1.0, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_thermal_coupling_with_efficiency() {
+        let coupling = ThermalCoupling::new(
+            CircuitId(0),
+            CircuitId(1),
+            ThermalConductance::from_watts_per_kelvin(1000.0),
+        )
+        .with_efficiency(0.85);
+
+        assert_relative_eq!(coupling.efficiency, 0.85, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_efficiency_clamped() {
+        let coupling = make_coupling(0, 1, 100.0).with_efficiency(1.5);
+        assert_relative_eq!(coupling.efficiency, 1.0, epsilon = 1e-10);
+
+        let coupling = make_coupling(0, 1, 100.0).with_efficiency(-0.5);
+        assert_relative_eq!(coupling.efficiency, 0.0, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_compute_coupling_heat_positive() {
+        let coupling = make_coupling(0, 1, 1000.0);
+        let t_hot = Temperature::from_kelvin(350.0);
+        let t_cold = Temperature::from_kelvin(300.0);
+
+        let q = compute_coupling_heat(&coupling, t_hot, t_cold);
+
+        // Q = 1.0 * 1000 * (350 - 300) = 50000 W
+        assert_relative_eq!(q, 50000.0, epsilon = 1e-10);
+        assert!(q > 0.0, "Heat should be positive (into cold circuit)");
+    }
+
+    #[test]
+    fn test_compute_coupling_heat_zero() {
+        let coupling = make_coupling(0, 1, 1000.0);
+        let t_hot = Temperature::from_kelvin(300.0);
+        let t_cold = Temperature::from_kelvin(300.0);
+
+        let q = compute_coupling_heat(&coupling, t_hot, t_cold);
+
+        assert_relative_eq!(q, 0.0, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_compute_coupling_heat_negative() {
+        let coupling = make_coupling(0, 1, 1000.0);
+        let t_hot = Temperature::from_kelvin(280.0);
+        let t_cold = Temperature::from_kelvin(300.0);
+
+        let q = compute_coupling_heat(&coupling, t_hot, t_cold);
+
+        // Q = 1000 * (280 - 300) = -20000 W (reverse flow)
+        assert_relative_eq!(q, -20000.0, epsilon = 1e-10);
+        assert!(q < 0.0, "Heat should be negative (reverse flow)");
+    }
+
+    #[test]
+    fn test_compute_coupling_heat_with_efficiency() {
+        let coupling = make_coupling(0, 1, 1000.0).with_efficiency(0.9);
+        let t_hot = Temperature::from_kelvin(350.0);
+        let t_cold = Temperature::from_kelvin(300.0);
+
+        let q = compute_coupling_heat(&coupling, t_hot, t_cold);
+
+        // Q = 0.9 * 1000 * 50 = 45000 W
+        assert_relative_eq!(q, 45000.0, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_energy_conservation() {
+        // For two circuits coupled, Q_hot = -Q_cold
+        // This means the heat leaving hot circuit equals heat entering cold circuit
+        let coupling = make_coupling(0, 1, 1000.0);
+        let t_hot = Temperature::from_kelvin(350.0);
+        let t_cold = Temperature::from_kelvin(300.0);
+
+        let q_into_cold = compute_coupling_heat(&coupling, t_hot, t_cold);
+        let q_out_of_hot = -q_into_cold; // By convention
+
+        // Heat into cold = - (heat out of hot)
+        assert_relative_eq!(q_into_cold, -q_out_of_hot, epsilon = 1e-10);
+        assert!(q_into_cold > 0.0, "Cold circuit receives heat");
+        assert!(q_out_of_hot < 0.0, "Hot circuit loses heat");
+    }
+
+    #[test]
+    fn test_no_circular_dependency() {
+        // Linear chain: A → B → C
+        let couplings = vec![make_coupling(0, 1, 100.0), make_coupling(1, 2, 100.0)];
+
+        assert!(!has_circular_dependencies(&couplings));
+    }
+
+    #[test]
+    fn test_circular_dependency_detection() {
+        // Mutual: A → B and B → A
+        let couplings = vec![make_coupling(0, 1, 100.0), make_coupling(1, 0, 100.0)];
+
+        assert!(has_circular_dependencies(&couplings));
+    }
+
+    #[test]
+    fn test_circular_dependency_complex() {
+        // Triangle: A → B → C → A
+        let couplings = vec![
+            make_coupling(0, 1, 100.0),
+            make_coupling(1, 2, 100.0),
+            make_coupling(2, 0, 100.0),
+        ];
+
+        assert!(has_circular_dependencies(&couplings));
+    }
+
+    #[test]
+    fn test_empty_couplings_no_cycle() {
+        let couplings: Vec<ThermalCoupling> = vec![];
+        assert!(!has_circular_dependencies(&couplings));
+    }
+
+    #[test]
+    fn test_single_coupling_no_cycle() {
+        let couplings = vec![make_coupling(0, 1, 100.0)];
+        assert!(!has_circular_dependencies(&couplings));
+    }
+
+    #[test]
+    fn test_coupling_groups_no_cycle() {
+        // A → B, C independent
+        let couplings = vec![make_coupling(0, 1, 100.0)];
+
+        let groups = coupling_groups(&couplings);
+
+        // With no cycles, each circuit is its own group
+        assert_eq!(groups.len(), 2);
+
+        // Each group should have exactly one circuit
+        for group in &groups {
+            assert_eq!(group.len(), 1);
+        }
+
+        // Collect all circuit IDs
+        let all_circuits: std::collections::HashSet<CircuitId> =
+            groups.iter().flat_map(|g| g.iter().copied()).collect();
+        assert!(all_circuits.contains(&CircuitId(0)));
+        assert!(all_circuits.contains(&CircuitId(1)));
+    }
+
+    #[test]
+    fn test_coupling_groups_with_cycle() {
+        // A ↔ B (mutual), C → D
+        let couplings = vec![
+            make_coupling(0, 1, 100.0),
+            make_coupling(1, 0, 100.0),
+            make_coupling(2, 3, 100.0),
+        ];
+
+        let groups = coupling_groups(&couplings);
+
+        // Should have 3 groups: [A, B] as one, C as one, D as one
+        assert_eq!(groups.len(), 3);
+
+        // Find the group with 2 circuits (A and B)
+        let large_group: Vec<&Vec<CircuitId>> = groups.iter().filter(|g| g.len() == 2).collect();
+        assert_eq!(large_group.len(), 1);
+
+        let ab_group = large_group[0];
+        assert!(ab_group.contains(&CircuitId(0)));
+        assert!(ab_group.contains(&CircuitId(1)));
+    }
+
+    #[test]
+    fn test_coupling_groups_empty() {
+        let couplings: Vec<ThermalCoupling> = vec![];
+        let groups = coupling_groups(&couplings);
+        assert!(groups.is_empty());
+    }
+}

crates/solver/src/error.rs

@@ -0,0 +1,72 @@
+//! Topology and solver error types.
+
+use thiserror::Error;
+
+/// Errors that can occur during topology validation or system construction.
+#[derive(Error, Debug, Clone, PartialEq)]
+pub enum TopologyError {
+    /// A node has no edges (isolated/dangling node).
+    #[error("Isolated node at index {node_index}: all components must be connected")]
+    IsolatedNode {
+        /// Index of the isolated node in the graph
+        node_index: usize,
+    },
+
+    /// Not all ports are connected (reserved for Story 3.2 port validation).
+    #[error("Unconnected ports: {message}")]
+    #[allow(dead_code)]
+    UnconnectedPorts {
+        /// Description of which ports are unconnected
+        message: String,
+    },
+
+    /// Topology validation failed for another reason.
+    #[error("Invalid topology: {message}")]
+    #[allow(dead_code)]
+    InvalidTopology {
+        /// Description of the validation failure
+        message: String,
+    },
+
+    /// Attempted to connect nodes in different circuits via a flow edge.
+    /// Flow edges must connect nodes within the same circuit. Cross-circuit
+    /// thermal coupling is handled in Story 3.4.
+    #[error("Cross-circuit connection not allowed: source circuit {source_circuit}, target circuit {target_circuit}. Flow edges connect only nodes within the same circuit")]
+    CrossCircuitConnection {
+        /// Circuit ID of the source node
+        source_circuit: u8,
+        /// Circuit ID of the target node
+        target_circuit: u8,
+    },
+
+    /// Too many circuits requested. Maximum is 5 (circuit IDs 0..=4).
+    #[error("Too many circuits: requested {requested}, maximum is 5")]
+    TooManyCircuits {
+        /// The requested circuit ID that exceeded the limit
+        requested: u8,
+    },
+
+    /// Attempted to add thermal coupling with a circuit that doesn't exist.
+    #[error(
+        "Invalid circuit for thermal coupling: circuit {circuit_id} does not exist in the system"
+    )]
+    InvalidCircuitForCoupling {
+        /// The circuit ID that was referenced but doesn't exist
+        circuit_id: u8,
+    },
+}
+
+/// Error when adding an edge with port validation.
+///
+/// Combines port validation errors ([`entropyk_components::ConnectionError`]) and topology errors
+/// ([`TopologyError`]) such as cross-circuit connection attempts.
+#[derive(Error, Debug, Clone, PartialEq)]
+pub enum AddEdgeError {
+    /// Port validation failed (fluid, pressure, enthalpy mismatch).
+    #[error(transparent)]
+    Connection(#[from] entropyk_components::ConnectionError),
+
+    /// Topology validation failed (e.g. cross-circuit connection).
+    #[error(transparent)]
+    Topology(#[from] TopologyError),
+}

crates/solver/src/graph.rs

@@ -0,0 +1,6 @@
+//! Graph building helpers for system topology.
+//!
+//! This module provides utilities for constructing and manipulating
+//! the system graph. The main [`System`](crate::system::System) struct
+//! handles graph operations; this module may be extended with convenience
+//! builders in future stories.

crates/solver/src/initializer.rs

@@ -0,0 +1,675 @@
+//! 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 fall back to sensible defaults (5 bar / 20 bar) with a warning.
+//!
+//! # 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;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// 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,
+    },
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// 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,
+}
+
+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)`
+    ///
+    /// For unknown fluids, returns sensible defaults (5 bar / 20 bar) with a
+    /// `tracing::warn!` log entry.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`InitializerError::TemperatureAboveCritical`] if the adjusted
+    /// source temperature exceeds the critical temperature for a known 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 => {
+                // Unknown fluid: emit warning and return sensible defaults
+                tracing::warn!(
+                    fluid = %fluid_str,
+                    "Unknown fluid for Antoine estimation — using fallback pressures \
+                     (P_evap = 5 bar, P_cond = 20 bar)"
+                );
+                Ok((
+                    Pressure::from_bar(5.0),
+                    Pressure::from_bar(20.0),
+                ))
+            }
+            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., `2 * edge_count`).
+    ///
+    /// State layout per edge: `[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.state_vector_len()`.
+    pub fn populate_state(
+        &self,
+        system: &System,
+        p_evap: Pressure,
+        p_cond: Pressure,
+        h_default: Enthalpy,
+        state: &mut [f64],
+    ) -> Result<(), InitializerError> {
+        let expected = system.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 (i, edge_idx) in system.edge_indices().enumerate() {
+            let circuit = system.edge_circuit(edge_idx);
+            let p = if circuit.0 == 0 { p_evap_pa } else { p_cond_pa };
+            state[2 * i] = p;
+            state[2 * i + 1] = h_jkg;
+        }
+
+        Ok(())
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    // ── 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 fallback (5 bar / 20 bar) without panic
+    #[test]
+    fn test_unknown_fluid_fallback() {
+        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!(result.is_ok(), "Unknown fluid should not return Err");
+        let (p_evap, p_cond) = result.unwrap();
+        assert_relative_eq!(p_evap.to_bar(), 5.0, max_relative = 1e-9);
+        assert_relative_eq!(p_cond.to_bar(), 20.0, max_relative = 1e-9);
+    }
+
+    /// 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, SystemState};
+
+        struct MockComp;
+        impl Component for MockComp {
+            fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
+                for v in r.iter_mut() { *v = 0.0; }
+                Ok(())
+            }
+            fn jacobian_entries(&self, _s: &SystemState, 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();
+
+        // All edges in circuit 0 (single-circuit) → p_evap
+        assert_eq!(state.len(), 4); // 2 edges × 2 entries
+        assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        assert_relative_eq!(state[2], p_evap.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[3], 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::{CircuitId, System};
+        use entropyk_components::{Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState};
+
+        struct MockComp;
+        impl Component for MockComp {
+            fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
+                for v in r.iter_mut() { *v = 0.0; }
+                Ok(())
+            }
+            fn jacobian_entries(&self, _s: &SystemState, 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();
+
+        assert_eq!(state.len(), 4); // 2 edges × 2 entries
+        // Edge 0 (circuit 0) → p_evap
+        assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        // Edge 1 (circuit 1) → p_cond
+        assert_relative_eq!(state[2], p_cond.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[3], 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, SystemState};
+
+        struct MockComp;
+        impl Component for MockComp {
+            fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
+                for v in r.iter_mut() { *v = 0.0; }
+                Ok(())
+            }
+            fn jacobian_entries(&self, _s: &SystemState, 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 2 state entries (1 edge × 2), 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: 2, 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
+        );
+    }
+}

crates/solver/src/jacobian.rs

@@ -67,6 +67,26 @@ impl JacobianMatrix {
         JacobianMatrix(matrix)
     }
 
+    /// Updates an existing Jacobian matrix from sparse entries in-place.
+    ///
+    /// The matrix is first zeroed out, then filled with the new entries.
+    /// This avoids re-allocating memory during iterations, satisfying the
+    /// zero-allocation architecture constraint.
+    ///
+    /// # Arguments
+    ///
+    /// * `entries` - Slice of `(row, col, value)` tuples
+    pub fn update_from_builder(&mut self, entries: &[(usize, usize, f64)]) {
+        self.0.fill(0.0);
+        let n_rows = self.0.nrows();
+        let n_cols = self.0.ncols();
+        for &(row, col, value) in entries {
+            if row < n_rows && col < n_cols {
+                self.0[(row, col)] += value;
+            }
+        }
+    }
+
     /// Creates a zero Jacobian matrix with the given dimensions.
     pub fn zeros(n_rows: usize, n_cols: usize) -> Self {
         JacobianMatrix(DMatrix::zeros(n_rows, n_cols))

crates/solver/src/lib.rs

@@ -0,0 +1,33 @@
+//! # Entropyk Solver
+//!
+//! System topology and solver engine for thermodynamic simulation.
+//!
+//! This crate provides the graph-based representation of thermodynamic systems,
+//! where components are nodes and flow connections are edges. Edges index into
+//! the solver's state vector (P and h per edge).
+
+pub mod coupling;
+pub mod criteria;
+pub mod error;
+pub mod graph;
+pub mod initializer;
+pub mod jacobian;
+pub mod solver;
+pub mod system;
+
+pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
+pub use coupling::{
+    compute_coupling_heat, coupling_groups, has_circular_dependencies, ThermalCoupling,
+};
+pub use entropyk_components::ConnectionError;
+pub use error::{AddEdgeError, TopologyError};
+pub use initializer::{
+    antoine_pressure, AntoineCoefficients, InitializerConfig, InitializerError, SmartInitializer,
+};
+pub use jacobian::JacobianMatrix;
+pub use solver::{
+    ConvergedState, ConvergenceStatus, FallbackConfig, FallbackSolver, JacobianFreezingConfig,
+    NewtonConfig, PicardConfig, Solver, SolverError, SolverStrategy, TimeoutConfig,
+};
+pub use system::{CircuitId, FlowEdge, System};
+

crates/solver/src/solver.rs

@@ -302,6 +302,61 @@ impl Default for TimeoutConfig {
     }
 }
 
+// ─────────────────────────────────────────────────────────────────────────────
+// Jacobian Freezing Configuration (Story 4.8)
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Configuration for Jacobian-freezing optimization.
+///
+/// When enabled, the Newton-Raphson solver reuses the previously computed
+/// Jacobian matrix for up to `max_frozen_iters` consecutive iterations,
+/// provided the residual norm is still decreasing. This avoids expensive
+/// Jacobian assembly and can reduce per-iteration CPU time by up to ~80%.
+///
+/// # Auto-disable on divergence
+///
+/// If the residual norm *increases* while a frozen Jacobian is being used,
+/// the solver immediately forces a fresh Jacobian computation on the next
+/// iteration and resets the frozen-iteration counter.
+///
+/// # Example
+///
+/// ```rust
+/// use entropyk_solver::solver::{NewtonConfig, JacobianFreezingConfig};
+///
+/// let config = NewtonConfig::default()
+///     .with_jacobian_freezing(JacobianFreezingConfig {
+///         max_frozen_iters: 3,
+///         threshold: 0.1,
+///     });
+/// ```
+#[derive(Debug, Clone, PartialEq)]
+pub struct JacobianFreezingConfig {
+    /// Maximum number of consecutive iterations the Jacobian may be reused
+    /// without recomputing.
+    ///
+    /// After this many frozen iterations the solver forces a fresh assembly,
+    /// even if the residual is still decreasing.  Default: 3.
+    pub max_frozen_iters: usize,
+
+    /// Residual-norm ratio threshold below which freezing is considered safe.
+    ///
+    /// Freezing is only attempted when
+    /// `current_norm / previous_norm < (1.0 - threshold)`,
+    /// ensuring that convergence is still progressing sufficiently.
+    /// Default: 0.1 (i.e., at least a 10 % residual decrease per step).
+    pub threshold: f64,
+}
+
+impl Default for JacobianFreezingConfig {
+    fn default() -> Self {
+        Self {
+            max_frozen_iters: 3,
+            threshold: 0.1,
+        }
+    }
+}
+
 // ─────────────────────────────────────────────────────────────────────────────
 // Configuration structs
 // ─────────────────────────────────────────────────────────────────────────────
@@ -393,6 +448,15 @@ pub struct NewtonConfig {
     /// test instead of the raw L2-norm tolerance check. The old `tolerance` field is retained
     /// for backward compatibility and is ignored when this is `Some`.
     pub convergence_criteria: Option<ConvergenceCriteria>,
+
+    /// Jacobian-freezing optimization (Story 4.8).
+    ///
+    /// When `Some`, the solver reuses the previous Jacobian matrix for up to
+    /// `max_frozen_iters` iterations while the residual is decreasing faster than
+    /// the configured threshold.  Auto-disables when the residual increases.
+    ///
+    /// Default: `None` (recompute every iteration — backward-compatible).
+    pub jacobian_freezing: Option<JacobianFreezingConfig>,
 }
 
 impl Default for NewtonConfig {
@@ -410,6 +474,7 @@ impl Default for NewtonConfig {
             previous_state: None,
             initial_state: None,
             convergence_criteria: None,
+            jacobian_freezing: None,
         }
     }
 }
@@ -435,6 +500,17 @@ impl NewtonConfig {
         self
     }
 
+    /// Enables Jacobian-freezing optimization (Story 4.8 — builder pattern).
+    ///
+    /// When set, the solver skips Jacobian re-assembly for iterations where the
+    /// residual is still decreasing, up to `config.max_frozen_iters` consecutive
+    /// frozen steps. Freezing is automatically disabled when the residual
+    /// increases.
+    pub fn with_jacobian_freezing(mut self, config: JacobianFreezingConfig) -> Self {
+        self.jacobian_freezing = Some(config);
+        self
+    }
+
     /// Computes the residual norm (L2 norm of the residual vector).
     fn residual_norm(residuals: &[f64]) -> f64 {
         residuals.iter().map(|r| r * r).sum::<f64>().sqrt()
@@ -658,6 +734,14 @@ impl Solver for NewtonConfig {
         let mut best_state: Vec<f64> = vec![0.0; n_state];
         let mut best_residual: f64;
 
+        // Story 4.8 — Jacobian-freezing tracking state.
+        // `frozen_count` tracks how many consecutive iterations have reused the Jacobian.
+        // `force_recompute` is set when a residual increase is detected.
+        // The Jacobian matrix itself is pre-allocated here (Zero Allocation AC)
+        let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state);
+        let mut frozen_count: usize = 0;
+        let mut force_recompute: bool = true; // Always compute on the very first iteration
+
         // Initial residual computation
         system
             .compute_residuals(&state, &mut residuals)
@@ -728,32 +812,74 @@ impl Solver for NewtonConfig {
                 }
             }
 
-            // Assemble Jacobian (AC: #3)
-            jacobian_builder.clear();
-            let jacobian_matrix = if self.use_numerical_jacobian {
-                // Numerical Jacobian via finite differences
-                let compute_residuals_fn = |s: &[f64], r: &mut [f64]| {
-                    let s_vec = s.to_vec();
-                    let mut r_vec = vec![0.0; r.len()];
-                    let result = system.compute_residuals(&s_vec, &mut r_vec);
-                    r.copy_from_slice(&r_vec);
-                    result.map(|_| ()).map_err(|e| format!("{:?}", e))
-                };
-                JacobianMatrix::numerical(compute_residuals_fn, &state, &residuals, 1e-8).map_err(
-                    |e| SolverError::InvalidSystem {
-                        message: format!("Failed to compute numerical Jacobian: {}", e),
-                    },
-                )?
+            // ── Jacobian Assembly / Freeze Decision (AC: #3, Story 4.8) ──
+            //
+            // Decide whether to recompute or reuse the Jacobian based on the
+            // freezing configuration and convergence behaviour.
+            let should_recompute = if let Some(ref freeze_cfg) = self.jacobian_freezing {
+                if force_recompute {
+                    true
+                } else if frozen_count >= freeze_cfg.max_frozen_iters {
+                    tracing::debug!(
+                        iteration = iteration,
+                        frozen_count = frozen_count,
+                        "Jacobian freeze limit reached — recomputing"
+                    );
+                    true
+                } else {
+                    false
+                }
             } else {
-                // Analytical Jacobian from components
-                system
-                    .assemble_jacobian(&state, &mut jacobian_builder)
-                    .map_err(|e| SolverError::InvalidSystem {
-                        message: format!("Failed to assemble Jacobian: {:?}", e),
-                    })?;
-                JacobianMatrix::from_builder(jacobian_builder.entries(), n_equations, n_state)
+                // No freezing configured — always recompute (backward-compatible)
+                true
             };
 
+            if should_recompute {
+                // Fresh Jacobian assembly (in-place update)
+                jacobian_builder.clear();
+                if self.use_numerical_jacobian {
+                    // Numerical Jacobian via finite differences
+                    let compute_residuals_fn = |s: &[f64], r: &mut [f64]| {
+                        let s_vec = s.to_vec();
+                        let mut r_vec = vec![0.0; r.len()];
+                        let result = system.compute_residuals(&s_vec, &mut r_vec);
+                        r.copy_from_slice(&r_vec);
+                        result.map(|_| ()).map_err(|e| format!("{:?}", e))
+                    };
+                    // Rather than creating a new matrix, compute it and assign
+                    let jm = JacobianMatrix::numerical(compute_residuals_fn, &state, &residuals, 1e-8)
+                        .map_err(|e| SolverError::InvalidSystem {
+                            message: format!("Failed to compute numerical Jacobian: {}", e),
+                        })?;
+                    // Deep copy elements to existing matrix (DMatrix::copy_from does not reallocate)
+                    jacobian_matrix.as_matrix_mut().copy_from(jm.as_matrix());
+                } else {
+                    // Analytical Jacobian from components
+                    system
+                        .assemble_jacobian(&state, &mut jacobian_builder)
+                        .map_err(|e| SolverError::InvalidSystem {
+                            message: format!("Failed to assemble Jacobian: {:?}", e),
+                        })?;
+                    jacobian_matrix.update_from_builder(jacobian_builder.entries());
+                };
+
+                frozen_count = 0;
+                force_recompute = false;
+
+                tracing::debug!(
+                    iteration = iteration,
+                    "Fresh Jacobian computed"
+                );
+            } else {
+                // Reuse the frozen Jacobian (Story 4.8 — AC: #2)
+                frozen_count += 1;
+                tracing::debug!(
+                    iteration = iteration,
+                    frozen_count = frozen_count,
+                    "Reusing frozen Jacobian"
+                );
+            }
+
             // Solve linear system J·Δx = -r (AC: #1)
             let delta = match jacobian_matrix.solve(&residuals) {
                 Some(d) => d,
@@ -811,6 +937,29 @@ impl Solver for NewtonConfig {
                 );
             }
 
+            // ── Story 4.8 — Jacobian-freeze feedback ──
+            //
+            // If the residual norm increased or did not decrease enough
+            // (below the threshold), force a fresh Jacobian on the next
+            // iteration and reset the frozen counter.
+            if let Some(ref freeze_cfg) = self.jacobian_freezing {
+                if previous_norm > 0.0
+                    && current_norm / previous_norm >= (1.0 - freeze_cfg.threshold)
+                {
+                    if frozen_count > 0 || !force_recompute {
+                        tracing::debug!(
+                            iteration = iteration,
+                            current_norm = current_norm,
+                            previous_norm = previous_norm,
+                            ratio = current_norm / previous_norm,
+                            "Residual not decreasing fast enough — unfreezing Jacobian"
+                        );
+                    }
+                    force_recompute = true;
+                    frozen_count = 0;
+                }
+            }
+
             tracing::debug!(
                 iteration = iteration,
                 residual_norm = current_norm,
@@ -1694,10 +1843,12 @@ impl FallbackSolver {
                                 tracing::debug!(
                                     final_residual = final_residual,
                                     threshold = self.config.return_to_newton_threshold,
-                                    "Picard not yet stabilized, continuing with Picard"
+                                    "Picard not yet stabilized, aborting"
                                 );
-                                // Continue with Picard - no allocation overhead
-                                continue;
+                                return Err(SolverError::NonConvergence {
+                                    iterations,
+                                    final_residual,
+                                });
                             }
                         }
                     }
@@ -1958,6 +2109,7 @@ mod tests {
             previous_state: None,
             initial_state: None,
             convergence_criteria: None,
+            jacobian_freezing: None,
         }
         .with_timeout(Duration::from_millis(200));
 

crates/solver/tests/fallback_solver.rs

@@ -0,0 +1,672 @@
+//! Integration tests for Story 4.4: Intelligent Fallback Strategy
+//!
+//! Tests the FallbackSolver behavior:
+//! - Newton diverges → Picard converges
+//! - Newton diverges → Picard stabilizes → Newton returns
+//! - Oscillation prevention (max switches reached)
+//! - Fallback disabled (pure Newton behavior)
+//! - Timeout applies across switches
+//! - No heap allocation during switches
+
+use entropyk_components::{
+    Component, ComponentError, JacobianBuilder, ResidualVector, SystemState,
+};
+use entropyk_solver::solver::{
+    ConvergenceStatus, FallbackConfig, FallbackSolver, NewtonConfig, PicardConfig, Solver,
+    SolverError, SolverStrategy,
+};
+use entropyk_solver::system::System;
+use std::time::Duration;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Mock Components for Testing
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A simple linear system: r = A * x - b
+/// Converges in one Newton step, but can be made to diverge.
+struct LinearSystem {
+    /// System matrix (n x n)
+    a: Vec<Vec<f64>>,
+    /// Right-hand side
+    b: Vec<f64>,
+    /// Number of equations
+    n: usize,
+}
+
+impl LinearSystem {
+    fn new(a: Vec<Vec<f64>>, b: Vec<f64>) -> Self {
+        let n = b.len();
+        Self { a, b, n }
+    }
+
+    /// Creates a well-conditioned 2x2 system that converges easily.
+    fn well_conditioned() -> Self {
+        // A = [[2, 1], [1, 2]], b = [3, 3]
+        // Solution: x = [1, 1]
+        Self::new(vec![vec![2.0, 1.0], vec![1.0, 2.0]], vec![3.0, 3.0])
+    }
+}
+
+impl Component for LinearSystem {
+    fn compute_residuals(
+        &self,
+        state: &SystemState,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        // r = A * x - b
+        for i in 0..self.n {
+            let mut ax_i = 0.0;
+            for j in 0..self.n {
+                ax_i += self.a[i][j] * state[j];
+            }
+            residuals[i] = ax_i - self.b[i];
+        }
+        Ok(())
+    }
+
+    fn jacobian_entries(
+        &self,
+        _state: &SystemState,
+        jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        // J = A (constant Jacobian)
+        for i in 0..self.n {
+            for j in 0..self.n {
+                jacobian.add_entry(i, j, self.a[i][j]);
+            }
+        }
+        Ok(())
+    }
+
+    fn n_equations(&self) -> usize {
+        self.n
+    }
+
+    fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+        &[]
+    }
+}
+
+/// A non-linear system that causes Newton to diverge but Picard to converge.
+/// Uses a highly non-linear residual function.
+struct StiffNonlinearSystem {
+    /// Non-linearity factor (higher = more stiff)
+    alpha: f64,
+    /// Number of equations
+    n: usize,
+}
+
+impl StiffNonlinearSystem {
+    fn new(alpha: f64, n: usize) -> Self {
+        Self { alpha, n }
+    }
+}
+
+impl Component for StiffNonlinearSystem {
+    fn compute_residuals(
+        &self,
+        state: &SystemState,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        // Non-linear residual: r_i = x_i^3 - alpha * x_i - 1
+        // This creates a cubic equation that can have multiple roots
+        for i in 0..self.n {
+            let x = state[i];
+            residuals[i] = x * x * x - self.alpha * x - 1.0;
+        }
+        Ok(())
+    }
+
+    fn jacobian_entries(
+        &self,
+        state: &SystemState,
+        jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        // J_ii = 3 * x_i^2 - alpha
+        for i in 0..self.n {
+            let x = state[i];
+            jacobian.add_entry(i, i, 3.0 * x * x - self.alpha);
+        }
+        Ok(())
+    }
+
+    fn n_equations(&self) -> usize {
+        self.n
+    }
+
+    fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+        &[]
+    }
+}
+
+/// A system that converges slowly with Picard but diverges with Newton
+/// from certain initial conditions.
+struct SlowConvergingSystem {
+    /// Convergence rate (0 < rate < 1)
+    rate: f64,
+    /// Target value
+    target: f64,
+}
+
+impl SlowConvergingSystem {
+    fn new(rate: f64, target: f64) -> Self {
+        Self { rate, target }
+    }
+}
+
+impl Component for SlowConvergingSystem {
+    fn compute_residuals(
+        &self,
+        state: &SystemState,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        // r = x - target (simple, but Newton can overshoot)
+        residuals[0] = state[0] - self.target;
+        Ok(())
+    }
+
+    fn jacobian_entries(
+        &self,
+        _state: &SystemState,
+        jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        jacobian.add_entry(0, 0, 1.0);
+        Ok(())
+    }
+
+    fn n_equations(&self) -> usize {
+        1
+    }
+
+    fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+        &[]
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Helper Functions
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Creates a minimal system with a single component for testing.
+fn create_test_system(component: Box<dyn Component>) -> System {
+    let mut system = System::new();
+    let n0 = system.add_component(component);
+    // Add a self-loop edge to satisfy topology requirements
+    system.add_edge(n0, n0).unwrap();
+    system.finalize().unwrap();
+    system
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Integration Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that FallbackSolver converges on a well-conditioned linear system.
+#[test]
+fn test_fallback_solver_converges_linear_system() {
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+    let mut solver = FallbackSolver::default_solver();
+
+    let result = solver.solve(&mut system);
+    assert!(result.is_ok(), "Should converge on well-conditioned system");
+
+    let converged = result.unwrap();
+    assert!(converged.is_converged());
+    assert!(converged.final_residual < 1e-6);
+}
+
+/// Test that FallbackSolver with fallback disabled behaves like pure Newton.
+#[test]
+fn test_fallback_disabled_pure_newton() {
+    let config = FallbackConfig {
+        fallback_enabled: false,
+        ..Default::default()
+    };
+    let mut solver = FallbackSolver::new(config);
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    let result = solver.solve(&mut system);
+    assert!(
+        result.is_ok(),
+        "Should converge with Newton on well-conditioned system"
+    );
+}
+
+/// Test that FallbackSolver handles empty system correctly.
+#[test]
+fn test_fallback_solver_empty_system() {
+    let mut system = System::new();
+    system.finalize().unwrap();
+
+    let mut solver = FallbackSolver::default_solver();
+    let result = solver.solve(&mut system);
+
+    assert!(result.is_err());
+    match result {
+        Err(SolverError::InvalidSystem { ref message }) => {
+            assert!(message.contains("Empty") || message.contains("no state"));
+        }
+        other => panic!("Expected InvalidSystem, got {:?}", other),
+    }
+}
+
+/// Test timeout enforcement across solver switches.
+#[test]
+fn test_fallback_solver_timeout() {
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    // Very short timeout that should trigger
+    let timeout = Duration::from_micros(1);
+    let mut solver = FallbackSolver::default_solver()
+        .with_timeout(timeout)
+        .with_newton_config(NewtonConfig {
+            max_iterations: 10000,
+            ..Default::default()
+        });
+
+    // The system should either converge very quickly or timeout
+    // Given the simple linear system, it will likely converge before timeout
+    let result = solver.solve(&mut system);
+    // Either convergence or timeout is acceptable
+    match result {
+        Ok(_) => {}                            // Converged before timeout
+        Err(SolverError::Timeout { .. }) => {} // Timed out as expected
+        Err(other) => panic!("Unexpected error: {:?}", other),
+    }
+}
+
+/// Test that FallbackSolver can be used as a trait object.
+#[test]
+fn test_fallback_solver_as_trait_object() {
+    let mut boxed: Box<dyn Solver> = Box::new(FallbackSolver::default_solver());
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    let result = boxed.solve(&mut system);
+    assert!(result.is_ok());
+}
+
+/// Test FallbackConfig customization.
+#[test]
+fn test_fallback_config_customization() {
+    let config = FallbackConfig {
+        fallback_enabled: true,
+        return_to_newton_threshold: 5e-4,
+        max_fallback_switches: 3,
+    };
+
+    let solver = FallbackSolver::new(config.clone());
+    assert_eq!(solver.config, config);
+    assert_eq!(solver.config.return_to_newton_threshold, 5e-4);
+    assert_eq!(solver.config.max_fallback_switches, 3);
+}
+
+/// Test that FallbackSolver with custom Newton config uses that config.
+#[test]
+fn test_fallback_solver_custom_newton_config() {
+    let newton_config = NewtonConfig {
+        max_iterations: 50,
+        tolerance: 1e-8,
+        ..Default::default()
+    };
+
+    let solver = FallbackSolver::default_solver().with_newton_config(newton_config.clone());
+    assert_eq!(solver.newton_config.max_iterations, 50);
+    assert!((solver.newton_config.tolerance - 1e-8).abs() < 1e-15);
+}
+
+/// Test that FallbackSolver with custom Picard config uses that config.
+#[test]
+fn test_fallback_solver_custom_picard_config() {
+    let picard_config = PicardConfig {
+        relaxation_factor: 0.3,
+        max_iterations: 200,
+        ..Default::default()
+    };
+
+    let solver = FallbackSolver::default_solver().with_picard_config(picard_config.clone());
+    assert!((solver.picard_config.relaxation_factor - 0.3).abs() < 1e-15);
+    assert_eq!(solver.picard_config.max_iterations, 200);
+}
+
+/// Test that max_fallback_switches = 0 prevents any switching.
+#[test]
+fn test_fallback_zero_switches() {
+    let config = FallbackConfig {
+        fallback_enabled: true,
+        max_fallback_switches: 0,
+        ..Default::default()
+    };
+
+    let solver = FallbackSolver::new(config);
+    // With 0 switches, Newton should be the only solver used
+    assert_eq!(solver.config.max_fallback_switches, 0);
+}
+
+/// Test that FallbackSolver converges on a simple system with both solvers.
+#[test]
+fn test_fallback_both_solvers_can_converge() {
+    // Create a system that both Newton and Picard can solve
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    // Test with Newton directly
+    let mut newton = NewtonConfig::default();
+    let newton_result = newton.solve(&mut system);
+    assert!(newton_result.is_ok(), "Newton should converge");
+
+    // Reset system
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    // Test with Picard directly
+    let mut picard = PicardConfig::default();
+    let picard_result = picard.solve(&mut system);
+    assert!(picard_result.is_ok(), "Picard should converge");
+
+    // Reset system
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    // Test with FallbackSolver
+    let mut fallback = FallbackSolver::default_solver();
+    let fallback_result = fallback.solve(&mut system);
+    assert!(fallback_result.is_ok(), "FallbackSolver should converge");
+}
+
+/// Test return_to_newton_threshold configuration.
+#[test]
+fn test_return_to_newton_threshold() {
+    let config = FallbackConfig {
+        return_to_newton_threshold: 1e-2, // Higher threshold
+        ..Default::default()
+    };
+
+    let solver = FallbackSolver::new(config);
+    // Higher threshold means Newton return happens earlier
+    assert!((solver.config.return_to_newton_threshold - 1e-2).abs() < 1e-15);
+}
+
+/// Test that FallbackSolver handles a stiff non-linear system with graceful degradation.
+#[test]
+fn test_fallback_stiff_nonlinear() {
+    // Create a stiff non-linear system that challenges both solvers
+    let mut system = create_test_system(Box::new(StiffNonlinearSystem::new(10.0, 2)));
+
+    let mut solver = FallbackSolver::default_solver()
+        .with_newton_config(NewtonConfig {
+            max_iterations: 50,
+            tolerance: 1e-6,
+            ..Default::default()
+        })
+        .with_picard_config(PicardConfig {
+            relaxation_factor: 0.3,
+            max_iterations: 200,
+            tolerance: 1e-6,
+            ..Default::default()
+        });
+
+    let result = solver.solve(&mut system);
+
+    // Verify expected behavior:
+    // 1. Should converge (fallback strategy succeeds)
+    // 2. Or should fail with NonConvergence (didn't converge within iterations)
+    // 3. Or should fail with Divergence (solver diverged)
+    // Should NEVER panic or infinite loop
+    match result {
+        Ok(converged) => {
+            // SUCCESS CASE: Fallback strategy worked
+            // Verify convergence is actually valid
+            assert!(
+                converged.final_residual < 1.0,
+                "Converged residual {} should be reasonable (< 1.0)",
+                converged.final_residual
+            );
+            if converged.is_converged() {
+                assert!(
+                    converged.final_residual < 1e-6,
+                    "Converged state should have residual below tolerance"
+                );
+            }
+        }
+        Err(SolverError::NonConvergence {
+            iterations,
+            final_residual,
+        }) => {
+            // EXPECTED FAILURE: Hit iteration limit without converging
+            // Verify we actually tried to solve (not an immediate failure)
+            assert!(
+                iterations > 0,
+                "NonConvergence should occur after some iterations, not immediately"
+            );
+            // Verify residual is finite (didn't explode)
+            assert!(
+                final_residual.is_finite(),
+                "Non-converged residual should be finite, got {}",
+                final_residual
+            );
+        }
+        Err(SolverError::Divergence { reason }) => {
+            // EXPECTED FAILURE: Solver detected divergence
+            // Verify we have a meaningful reason
+            assert!(!reason.is_empty(), "Divergence error should have a reason");
+            assert!(
+                reason.contains("diverg")
+                    || reason.contains("exceed")
+                    || reason.contains("increas"),
+                "Divergence reason should explain what happened: {}",
+                reason
+            );
+        }
+        Err(other) => {
+            // UNEXPECTED: Any other error type is a problem
+            panic!("Unexpected error type for stiff system: {:?}", other);
+        }
+    }
+}
+
+/// Test that timeout is enforced across solver switches.
+#[test]
+fn test_timeout_across_switches() {
+    let mut system = create_test_system(Box::new(StiffNonlinearSystem::new(5.0, 2)));
+
+    // Very short timeout
+    let timeout = Duration::from_millis(10);
+    let mut solver = FallbackSolver::default_solver()
+        .with_timeout(timeout)
+        .with_newton_config(NewtonConfig {
+            max_iterations: 1000,
+            ..Default::default()
+        })
+        .with_picard_config(PicardConfig {
+            max_iterations: 1000,
+            ..Default::default()
+        });
+
+    let result = solver.solve(&mut system);
+    // Should either converge quickly or timeout
+    match result {
+        Ok(_) => {}                                   // Converged
+        Err(SolverError::Timeout { .. }) => {}        // Timed out
+        Err(SolverError::NonConvergence { .. }) => {} // Didn't converge in time
+        Err(SolverError::Divergence { .. }) => {}     // Diverged
+        Err(other) => panic!("Unexpected error: {:?}", other),
+    }
+}
+
+/// Test that max_fallback_switches config value is respected.
+#[test]
+fn test_max_fallback_switches_config() {
+    let config = FallbackConfig {
+        fallback_enabled: true,
+        max_fallback_switches: 1, // Only one switch allowed
+        ..Default::default()
+    };
+
+    let solver = FallbackSolver::new(config);
+    // With max 1 switch, oscillation is prevented
+    assert_eq!(solver.config.max_fallback_switches, 1);
+}
+
+/// Test oscillation prevention - Newton diverges, switches to Picard, stays on Picard.
+#[test]
+fn test_oscillation_prevention_newton_to_picard_stays() {
+    use entropyk_solver::solver::{
+        FallbackConfig, FallbackSolver, NewtonConfig, PicardConfig, Solver,
+    };
+
+    // Create a system where Newton diverges but Picard converges
+    // Use StiffNonlinearSystem with high alpha to cause Newton divergence
+    let mut system = create_test_system(Box::new(StiffNonlinearSystem::new(100.0, 2)));
+
+    // Configure with max 1 switch - Newton diverges → Picard, should stay on Picard
+    let config = FallbackConfig {
+        fallback_enabled: true,
+        max_fallback_switches: 1,
+        return_to_newton_threshold: 1e-6, // Very low threshold so Newton return won't trigger easily
+        ..Default::default()
+    };
+
+    let mut solver = FallbackSolver::new(config)
+        .with_newton_config(NewtonConfig {
+            max_iterations: 20,
+            tolerance: 1e-6,
+            ..Default::default()
+        })
+        .with_picard_config(PicardConfig {
+            relaxation_factor: 0.2,
+            max_iterations: 500,
+            ..Default::default()
+        });
+
+    // Should either converge (Picard succeeds) or non-converge (but NOT oscillate)
+    let result = solver.solve(&mut system);
+
+    match result {
+        Ok(converged) => {
+            // Success - Picard converged after Newton divergence
+            assert!(converged.is_converged() || converged.final_residual < 1.0);
+        }
+        Err(SolverError::NonConvergence { .. }) => {
+            // Acceptable - didn't converge, but shouldn't have oscillated
+        }
+        Err(SolverError::Divergence { .. }) => {
+            // Picard diverged - acceptable for stiff system
+        }
+        Err(other) => panic!("Unexpected error type: {:?}", other),
+    }
+}
+
+/// Test that Newton re-divergence causes permanent commit to Picard.
+#[test]
+fn test_newton_redivergence_commits_to_picard() {
+    // Create a system that's borderline - Newton might diverge, Picard converges slowly
+    let mut system = create_test_system(Box::new(StiffNonlinearSystem::new(50.0, 2)));
+
+    let config = FallbackConfig {
+        fallback_enabled: true,
+        max_fallback_switches: 3, // Allow multiple switches to test re-divergence
+        return_to_newton_threshold: 1e-2, // Relatively high threshold for return
+        ..Default::default()
+    };
+
+    let mut solver = FallbackSolver::new(config)
+        .with_newton_config(NewtonConfig {
+            max_iterations: 30,
+            tolerance: 1e-8,
+            ..Default::default()
+        })
+        .with_picard_config(PicardConfig {
+            relaxation_factor: 0.25,
+            max_iterations: 300,
+            ..Default::default()
+        });
+
+    let result = solver.solve(&mut system);
+
+    // Should complete without infinite oscillation
+    match result {
+        Ok(converged) => {
+            assert!(converged.final_residual < 1.0 || converged.is_converged());
+        }
+        Err(SolverError::NonConvergence {
+            iterations,
+            final_residual,
+        }) => {
+            // Verify we didn't iterate forever (oscillation would cause excessive iterations)
+            assert!(
+                iterations < 1000,
+                "Too many iterations - possible oscillation"
+            );
+            assert!(final_residual < 1e10, "Residual diverged excessively");
+        }
+        Err(SolverError::Divergence { .. }) => {
+            // Acceptable - system is stiff
+        }
+        Err(other) => panic!("Unexpected error: {:?}", other),
+    }
+}
+
+/// Test that FallbackSolver works with SolverStrategy pattern.
+#[test]
+fn test_fallback_solver_integration() {
+    // Verify FallbackSolver can be used alongside other solvers
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+
+    // Test with SolverStrategy::NewtonRaphson
+    let mut strategy = SolverStrategy::default();
+    let result1 = strategy.solve(&mut system);
+    assert!(result1.is_ok());
+
+    // Reset and test with FallbackSolver
+    let mut system = create_test_system(Box::new(LinearSystem::well_conditioned()));
+    let mut fallback = FallbackSolver::default_solver();
+    let result2 = fallback.solve(&mut system);
+    assert!(result2.is_ok());
+
+    // Both should converge to similar residuals
+    let r1 = result1.unwrap();
+    let r2 = result2.unwrap();
+    assert!((r1.final_residual - r2.final_residual).abs() < 1e-6);
+}
+
+/// Test that FallbackSolver handles convergence at initial state.
+#[test]
+fn test_fallback_already_converged() {
+    // Create a system that's already at solution
+    struct ZeroResidualComponent;
+
+    impl Component for ZeroResidualComponent {
+        fn compute_residuals(
+            &self,
+            _state: &SystemState,
+            residuals: &mut ResidualVector,
+        ) -> Result<(), ComponentError> {
+            residuals[0] = 0.0; // Already zero
+            Ok(())
+        }
+
+        fn jacobian_entries(
+            &self,
+            _state: &SystemState,
+            jacobian: &mut JacobianBuilder,
+        ) -> Result<(), ComponentError> {
+            jacobian.add_entry(0, 0, 1.0);
+            Ok(())
+        }
+
+        fn n_equations(&self) -> usize {
+            1
+        }
+
+        fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+            &[]
+        }
+    }
+
+    let mut system = create_test_system(Box::new(ZeroResidualComponent));
+    let mut solver = FallbackSolver::default_solver();
+
+    let result = solver.solve(&mut system);
+    assert!(result.is_ok());
+
+    let converged = result.unwrap();
+    assert_eq!(converged.iterations, 0); // Should converge immediately
+    assert!(converged.is_converged());
+}

crates/solver/tests/multi_circuit.rs

@@ -0,0 +1,239 @@
+//! Integration tests for multi-circuit machine definition (Story 3.3, FR9).
+//!
+//! Verifies multi-circuit heat pump topology (refrigerant + water) without thermal coupling.
+//! Tests circuits from 2 up to the maximum of 5 circuits (circuit IDs 0-4).
+
+use entropyk_components::{
+    Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
+};
+use entropyk_solver::{CircuitId, System, ThermalCoupling, TopologyError};
+use entropyk_core::ThermalConductance;
+
+/// Mock refrigerant component (e.g. compressor, condenser refrigerant side).
+struct RefrigerantMock {
+    n_equations: usize,
+}
+
+impl Component for RefrigerantMock {
+    fn compute_residuals(
+        &self,
+        _state: &SystemState,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        for r in residuals.iter_mut().take(self.n_equations) {
+            *r = 0.0;
+        }
+        Ok(())
+    }
+
+    fn jacobian_entries(
+        &self,
+        _state: &SystemState,
+        _jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        Ok(())
+    }
+
+    fn n_equations(&self) -> usize {
+        self.n_equations
+    }
+
+    fn get_ports(&self) -> &[ConnectedPort] {
+        &[]
+    }
+}
+
+#[test]
+fn test_two_circuit_heat_pump_topology() {
+    let mut sys = System::new();
+
+    // Circuit 0: refrigerant (compressor -> condenser -> valve -> evaporator)
+    let comp = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 2 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    let cond = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 2 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    let valve = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 2 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    let evap = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 2 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+
+    sys.add_edge(comp, cond).unwrap();
+    sys.add_edge(cond, valve).unwrap();
+    sys.add_edge(valve, evap).unwrap();
+    sys.add_edge(evap, comp).unwrap();
+
+    // Circuit 1: water (pump -> condenser water side -> evaporator water side)
+    let pump = sys
+        .add_component_to_circuit(Box::new(RefrigerantMock { n_equations: 1 }), CircuitId(1))
+        .unwrap();
+    let cond_w = sys
+        .add_component_to_circuit(Box::new(RefrigerantMock { n_equations: 1 }), CircuitId(1))
+        .unwrap();
+    let evap_w = sys
+        .add_component_to_circuit(Box::new(RefrigerantMock { n_equations: 1 }), CircuitId(1))
+        .unwrap();
+
+    sys.add_edge(pump, cond_w).unwrap();
+    sys.add_edge(cond_w, evap_w).unwrap();
+    sys.add_edge(evap_w, pump).unwrap();
+
+    assert_eq!(sys.circuit_count(), 2);
+    assert_eq!(sys.circuit_nodes(CircuitId::ZERO).count(), 4);
+    assert_eq!(sys.circuit_nodes(CircuitId(1)).count(), 3);
+    assert_eq!(sys.circuit_edges(CircuitId::ZERO).count(), 4);
+    assert_eq!(sys.circuit_edges(CircuitId(1)).count(), 3);
+
+    let result = sys.finalize();
+    assert!(
+        result.is_ok(),
+        "finalize should succeed: {:?}",
+        result.err()
+    );
+}
+
+#[test]
+fn test_cross_circuit_rejected_integration() {
+    let mut sys = System::new();
+    let n0 = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 0 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    let n1 = sys
+        .add_component_to_circuit(Box::new(RefrigerantMock { n_equations: 0 }), CircuitId(1))
+        .unwrap();
+
+    let result = sys.add_edge(n0, n1);
+    assert!(result.is_err());
+    assert!(matches!(
+        result,
+        Err(TopologyError::CrossCircuitConnection { .. })
+    ));
+}
+
+#[test]
+fn test_maximum_five_circuits_integration() {
+    // Integration test: Verify maximum of 5 circuits (IDs 0-4) is supported
+    let mut sys = System::new();
+
+    // Create 5 separate circuits, each with 2 nodes forming a cycle
+    for circuit_id in 0..=4 {
+        let n0 = sys
+            .add_component_to_circuit(
+                Box::new(RefrigerantMock { n_equations: 1 }),
+                CircuitId(circuit_id),
+            )
+            .unwrap();
+        let n1 = sys
+            .add_component_to_circuit(
+                Box::new(RefrigerantMock { n_equations: 1 }),
+                CircuitId(circuit_id),
+            )
+            .unwrap();
+        sys.add_edge(n0, n1).unwrap();
+        sys.add_edge(n1, n0).unwrap();
+    }
+
+    assert_eq!(sys.circuit_count(), 5, "should have exactly 5 circuits");
+
+    // Verify each circuit has its own nodes and edges
+    for circuit_id in 0..=4 {
+        assert_eq!(
+            sys.circuit_nodes(CircuitId(circuit_id)).count(),
+            2,
+            "circuit {} should have 2 nodes",
+            circuit_id
+        );
+        assert_eq!(
+            sys.circuit_edges(CircuitId(circuit_id)).count(),
+            2,
+            "circuit {} should have 2 edges",
+            circuit_id
+        );
+    }
+
+    // Verify 6th circuit is rejected
+    let result =
+        sys.add_component_to_circuit(Box::new(RefrigerantMock { n_equations: 1 }), CircuitId(5));
+    assert!(
+        result.is_err(),
+        "circuit 5 should be rejected (exceeds max of 4)"
+    );
+    assert!(matches!(
+        result,
+        Err(TopologyError::TooManyCircuits { requested: 5 })
+    ));
+
+    // Verify system can still be finalized with 5 circuits
+    sys.finalize().unwrap();
+}
+
+#[test]
+fn test_coupling_residuals_basic() {
+    // Two circuits with one thermal coupling; verify coupling_residual_count and coupling_residuals.
+    let mut sys = System::new();
+    let n0 = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 1 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    let n1 = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 1 }),
+            CircuitId::ZERO,
+        )
+        .unwrap();
+    sys.add_edge(n0, n1).unwrap();
+    sys.add_edge(n1, n0).unwrap();
+
+    let n2 = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 1 }),
+            CircuitId(1),
+        )
+        .unwrap();
+    let n3 = sys
+        .add_component_to_circuit(
+            Box::new(RefrigerantMock { n_equations: 1 }),
+            CircuitId(1),
+        )
+        .unwrap();
+    sys.add_edge(n2, n3).unwrap();
+    sys.add_edge(n3, n2).unwrap();
+
+    let coupling = ThermalCoupling::new(
+        CircuitId::ZERO,
+        CircuitId(1),
+        ThermalConductance::from_watts_per_kelvin(1000.0),
+    );
+    sys.add_thermal_coupling(coupling).unwrap();
+
+    sys.finalize().unwrap();
+
+    assert_eq!(sys.coupling_residual_count(), 1);
+
+    let temperatures = [(350.0_f64, 300.0_f64)]; // T_hot, T_cold in K
+    let mut out = [0.0_f64; 4];
+    sys.coupling_residuals(&temperatures, &mut out);
+    // Q = UA * (T_hot - T_cold) = 1000 * 50 = 50000 W into cold circuit
+    assert!(out[0] > 0.0);
+    assert!((out[0] - 50000.0).abs() < 1.0);
+}

crates/solver/tests/newton_convergence.rs

@@ -0,0 +1,480 @@
+//! Comprehensive integration tests for Newton-Raphson solver (Story 4.2).
+//!
+//! Tests cover all Acceptance Criteria:
+//! - AC #1: Quadratic convergence near solution
+//! - AC #2: Line search prevents overshooting
+//! - AC #3: Analytical and numerical Jacobian support
+//! - AC #4: Timeout enforcement
+//! - AC #5: Divergence detection
+//! - AC #6: Pre-allocated buffers
+
+use entropyk_solver::{ConvergenceStatus, JacobianMatrix, NewtonConfig, Solver, SolverError, System};
+use approx::assert_relative_eq;
+use std::time::Duration;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #1: Quadratic Convergence Near Solution
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that Newton-Raphson exhibits quadratic convergence on a simple system.
+/// 
+/// For a well-conditioned system near the solution, the residual norm should
+/// decrease quadratically (roughly square each iteration).
+#[test]
+fn test_quadratic_convergence_simple_system() {
+    // We'll test the Jacobian solve directly since we need a mock system
+    // For J = [[2, 0], [0, 3]] and r = [2, 3], solution is x = [-1, -1]
+    
+    let entries = vec![(0, 0, 2.0), (1, 1, 3.0)];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let residuals = vec![2.0, 3.0];
+    let delta = jacobian.solve(&residuals).expect("non-singular");
+    
+    // J·Δx = -r => Δx = -J^{-1}·r
+    assert_relative_eq!(delta[0], -1.0, epsilon = 1e-10);
+    assert_relative_eq!(delta[1], -1.0, epsilon = 1e-10);
+}
+
+/// Test convergence on a 2x2 linear system.
+#[test]
+fn test_solve_2x2_linear_system() {
+    // J = [[4, 1], [1, 3]], r = [1, 2]
+    // Solution: Δx = -J^{-1}·r
+    let entries = vec![(0, 0, 4.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 3.0)];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let residuals = vec![1.0, 2.0];
+    let delta = jacobian.solve(&residuals).expect("non-singular");
+    
+    // Verify: J·Δx = -r
+    let j00 = 4.0;
+    let j01 = 1.0;
+    let j10 = 1.0;
+    let j11 = 3.0;
+    
+    let computed_r0 = j00 * delta[0] + j01 * delta[1];
+    let computed_r1 = j10 * delta[0] + j11 * delta[1];
+    
+    assert_relative_eq!(computed_r0, -1.0, epsilon = 1e-10);
+    assert_relative_eq!(computed_r1, -2.0, epsilon = 1e-10);
+}
+
+/// Test that a diagonal system converges in one Newton iteration.
+#[test]
+fn test_diagonal_system_one_iteration() {
+    // For a diagonal Jacobian, Newton should converge in 1 iteration
+    // J = [[a, 0], [0, b]], r = [c, d]
+    // Δx = [-c/a, -d/b]
+    
+    let entries = vec![(0, 0, 5.0), (1, 1, 7.0)];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let residuals = vec![10.0, 21.0];
+    let delta = jacobian.solve(&residuals).expect("non-singular");
+    
+    assert_relative_eq!(delta[0], -2.0, epsilon = 1e-10);
+    assert_relative_eq!(delta[1], -3.0, epsilon = 1e-10);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #2: Line Search Prevents Overshooting
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that line search is configured correctly.
+#[test]
+fn test_line_search_configuration() {
+    let cfg = NewtonConfig {
+        line_search: true,
+        line_search_armijo_c: 1e-4,
+        line_search_max_backtracks: 20,
+        ..Default::default()
+    };
+    
+    assert!(cfg.line_search);
+    assert_relative_eq!(cfg.line_search_armijo_c, 1e-4);
+    assert_eq!(cfg.line_search_max_backtracks, 20);
+}
+
+/// Test that line search can be disabled.
+#[test]
+fn test_line_search_disabled_by_default() {
+    let cfg = NewtonConfig::default();
+    assert!(!cfg.line_search);
+}
+
+/// Test Armijo condition constants are sensible.
+#[test]
+fn test_armijo_constant_range() {
+    let cfg = NewtonConfig::default();
+    
+    // Armijo constant should be in (0, 0.5) for typical line search
+    assert!(cfg.line_search_armijo_c > 0.0);
+    assert!(cfg.line_search_armijo_c < 0.5);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #3: Analytical and Numerical Jacobian Support
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that numerical Jacobian can be enabled.
+#[test]
+fn test_numerical_jacobian_configuration() {
+    let cfg = NewtonConfig {
+        use_numerical_jacobian: true,
+        ..Default::default()
+    };
+    
+    assert!(cfg.use_numerical_jacobian);
+}
+
+/// Test that analytical Jacobian is the default.
+#[test]
+fn test_analytical_jacobian_default() {
+    let cfg = NewtonConfig::default();
+    assert!(!cfg.use_numerical_jacobian);
+}
+
+/// Test numerical Jacobian computation matches analytical for linear function.
+#[test]
+fn test_numerical_jacobian_linear_function() {
+    // r[0] = 2*x0 + 3*x1
+    // r[1] = x0 - 2*x1
+    // J = [[2, 3], [1, -2]]
+    
+    let state = vec![1.0, 2.0];
+    let residuals = vec![2.0 * state[0] + 3.0 * state[1], state[0] - 2.0 * state[1]];
+    
+    let compute_residuals = |s: &[f64], r: &mut [f64]| {
+        r[0] = 2.0 * s[0] + 3.0 * s[1];
+        r[1] = s[0] - 2.0 * s[1];
+        Ok(())
+    };
+    
+    let j_num = JacobianMatrix::numerical(compute_residuals, &state, &residuals, 1e-8).unwrap();
+    
+    // Check against analytical Jacobian
+    assert_relative_eq!(j_num.get(0, 0).unwrap(), 2.0, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(0, 1).unwrap(), 3.0, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(1, 0).unwrap(), 1.0, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(1, 1).unwrap(), -2.0, epsilon = 1e-5);
+}
+
+/// Test numerical Jacobian for non-linear function.
+#[test]
+fn test_numerical_jacobian_nonlinear_function() {
+    // r[0] = x0^2 + x1
+    // r[1] = sin(x0) + cos(x1)
+    // J = [[2*x0, 1], [cos(x0), -sin(x1)]]
+    
+    let state = vec![0.5_f64, 1.0_f64];
+    let residuals = vec![state[0].powi(2) + state[1], state[0].sin() + state[1].cos()];
+    
+    let compute_residuals = |s: &[f64], r: &mut [f64]| {
+        r[0] = s[0].powi(2) + s[1];
+        r[1] = s[0].sin() + s[1].cos();
+        Ok(())
+    };
+    
+    let j_num = JacobianMatrix::numerical(compute_residuals, &state, &residuals, 1e-8).unwrap();
+    
+    // Analytical values
+    let j00 = 2.0 * state[0]; // 1.0
+    let j01 = 1.0;
+    let j10 = state[0].cos();
+    let j11 = -state[1].sin();
+    
+    assert_relative_eq!(j_num.get(0, 0).unwrap(), j00, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(0, 1).unwrap(), j01, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(1, 0).unwrap(), j10, epsilon = 1e-5);
+    assert_relative_eq!(j_num.get(1, 1).unwrap(), j11, epsilon = 1e-5);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #4: Timeout Enforcement
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test timeout configuration.
+#[test]
+fn test_timeout_configuration() {
+    let timeout = Duration::from_millis(500);
+    let cfg = NewtonConfig::default().with_timeout(timeout);
+    
+    assert_eq!(cfg.timeout, Some(timeout));
+}
+
+/// Test timeout is None by default.
+#[test]
+fn test_no_timeout_by_default() {
+    let cfg = NewtonConfig::default();
+    assert!(cfg.timeout.is_none());
+}
+
+/// Test timeout error contains correct duration.
+#[test]
+fn test_timeout_error_contains_duration() {
+    let err = SolverError::Timeout { timeout_ms: 1234 };
+    let msg = err.to_string();
+    
+    assert!(msg.contains("1234"));
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #5: Divergence Detection
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test divergence threshold configuration.
+#[test]
+fn test_divergence_threshold_configuration() {
+    let cfg = NewtonConfig {
+        divergence_threshold: 1e8,
+        ..Default::default()
+    };
+    
+    assert_relative_eq!(cfg.divergence_threshold, 1e8);
+}
+
+/// Test default divergence threshold.
+#[test]
+fn test_default_divergence_threshold() {
+    let cfg = NewtonConfig::default();
+    assert_relative_eq!(cfg.divergence_threshold, 1e10);
+}
+
+/// Test divergence error contains reason.
+#[test]
+fn test_divergence_error_contains_reason() {
+    let err = SolverError::Divergence {
+        reason: "Residual increased for 3 consecutive iterations".to_string(),
+    };
+    let msg = err.to_string();
+    
+    assert!(msg.contains("Residual increased"));
+    assert!(msg.contains("3 consecutive"));
+}
+
+/// Test divergence error for threshold exceeded.
+#[test]
+fn test_divergence_error_threshold_exceeded() {
+    let err = SolverError::Divergence {
+        reason: "Residual norm 1e12 exceeds threshold 1e10".to_string(),
+    };
+    let msg = err.to_string();
+    
+    assert!(msg.contains("exceeds threshold"));
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #6: Pre-Allocated Buffers
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that solver handles empty system gracefully (pre-allocated buffers work).
+#[test]
+fn test_preallocated_buffers_empty_system() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = NewtonConfig::default();
+    let result = solver.solve(&mut sys);
+    
+    // Should return error without panic
+    assert!(result.is_err());
+}
+
+/// Test that solver handles configuration variations without panic.
+#[test]
+fn test_preallocated_buffers_all_configs() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    // Test with all features enabled
+    let mut solver = NewtonConfig {
+        max_iterations: 50,
+        tolerance: 1e-8,
+        line_search: true,
+        timeout: Some(Duration::from_millis(100)),
+        use_numerical_jacobian: true,
+        line_search_armijo_c: 1e-3,
+        line_search_max_backtracks: 10,
+        divergence_threshold: 1e8,
+        ..Default::default()
+    };
+    
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err()); // Empty system, but no panic
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Jacobian Matrix Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test singular Jacobian returns None.
+#[test]
+fn test_singular_jacobian_returns_none() {
+    // Singular matrix: [[1, 1], [1, 1]]
+    let entries = vec![(0, 0, 1.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 1.0)];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let residuals = vec![1.0, 2.0];
+    let result = jacobian.solve(&residuals);
+    
+    assert!(result.is_none(), "Singular matrix should return None");
+}
+
+/// Test zero Jacobian returns None.
+#[test]
+fn test_zero_jacobian_returns_none() {
+    let jacobian = JacobianMatrix::zeros(2, 2);
+    
+    let residuals = vec![1.0, 2.0];
+    let result = jacobian.solve(&residuals);
+    
+    assert!(result.is_none(), "Zero matrix should return None");
+}
+
+/// Test Jacobian condition number for well-conditioned matrix.
+#[test]
+fn test_jacobian_condition_number_well_conditioned() {
+    let entries = vec![(0, 0, 1.0), (1, 1, 1.0)];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let cond = jacobian.condition_number().unwrap();
+    assert_relative_eq!(cond, 1.0, epsilon = 1e-10);
+}
+
+/// Test Jacobian condition number for ill-conditioned matrix.
+#[test]
+fn test_jacobian_condition_number_ill_conditioned() {
+    // Nearly singular matrix
+    let entries = vec![
+        (0, 0, 1.0),
+        (0, 1, 1.0),
+        (1, 0, 1.0),
+        (1, 1, 1.0 + 1e-12),
+    ];
+    let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
+    
+    let cond = jacobian.condition_number();
+    assert!(cond.unwrap() > 1e10, "Should be ill-conditioned");
+}
+
+/// Test Jacobian for non-square (overdetermined) system uses least-squares.
+#[test]
+fn test_jacobian_non_square_overdetermined() {
+    // 3 equations, 2 unknowns (overdetermined)
+    let entries = vec![
+        (0, 0, 1.0),
+        (0, 1, 1.0),
+        (1, 0, 1.0),
+        (1, 1, 2.0),
+        (2, 0, 1.0),
+        (2, 1, 3.0),
+    ];
+    let jacobian = JacobianMatrix::from_builder(&entries, 3, 2);
+    
+    let residuals = vec![1.0, 2.0, 3.0];
+    let result = jacobian.solve(&residuals);
+    
+    // Should return a least-squares solution
+    assert!(result.is_some(), "Non-square system should return least-squares solution");
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// ConvergenceStatus Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test ConvergenceStatus::Converged.
+#[test]
+fn test_convergence_status_converged() {
+    use entropyk_solver::ConvergedState;
+    
+    let state = ConvergedState::new(
+        vec![1.0, 2.0],
+        10,
+        1e-8,
+        ConvergenceStatus::Converged,
+    );
+    
+    assert!(state.is_converged());
+    assert_eq!(state.status, ConvergenceStatus::Converged);
+}
+
+/// Test ConvergenceStatus::TimedOutWithBestState.
+#[test]
+fn test_convergence_status_timed_out() {
+    use entropyk_solver::ConvergedState;
+    
+    let state = ConvergedState::new(
+        vec![1.0],
+        50,
+        1e-3,
+        ConvergenceStatus::TimedOutWithBestState,
+    );
+    
+    assert!(!state.is_converged());
+    assert_eq!(state.status, ConvergenceStatus::TimedOutWithBestState);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Error Display Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test NonConvergence error display.
+#[test]
+fn test_non_convergence_display() {
+    let err = SolverError::NonConvergence {
+        iterations: 100,
+        final_residual: 1.23e-4,
+    };
+    let msg = err.to_string();
+    
+    assert!(msg.contains("100"));
+    assert!(msg.contains("1.23"));
+}
+
+/// Test InvalidSystem error display.
+#[test]
+fn test_invalid_system_display() {
+    let err = SolverError::InvalidSystem {
+        message: "Empty system has no equations".to_string(),
+    };
+    let msg = err.to_string();
+    
+    assert!(msg.contains("Empty system"));
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Configuration Validation Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Test that max_iterations must be positive.
+#[test]
+fn test_max_iterations_positive() {
+    let cfg = NewtonConfig::default();
+    assert!(cfg.max_iterations > 0);
+}
+
+/// Test that tolerance must be positive.
+#[test]
+fn test_tolerance_positive() {
+    let cfg = NewtonConfig::default();
+    assert!(cfg.tolerance > 0.0);
+}
+
+/// Test that relaxation factor for Picard is in valid range.
+#[test]
+fn test_picard_relaxation_factor_range() {
+    use entropyk_solver::PicardConfig;
+    
+    let cfg = PicardConfig::default();
+    assert!(cfg.relaxation_factor > 0.0);
+    assert!(cfg.relaxation_factor <= 1.0);
+}
+
+/// Test line search max backtracks is reasonable.
+#[test]
+fn test_line_search_max_backtracks_reasonable() {
+    let cfg = NewtonConfig::default();
+    assert!(cfg.line_search_max_backtracks > 0);
+    assert!(cfg.line_search_max_backtracks <= 100);
+}

crates/solver/tests/newton_raphson.rs

@@ -0,0 +1,254 @@
+//! Integration tests for Newton-Raphson solver (Story 4.2).
+//!
+//! Tests cover:
+//! - AC #1: Solver trait and strategy dispatch
+//! - AC #2: Configuration options
+//! - AC #3: Timeout enforcement
+//! - AC #4: Error handling for empty/invalid systems
+//! - AC #5: Pre-allocated buffers (no panic)
+
+use entropyk_solver::{NewtonConfig, Solver, SolverError, System};
+use approx::assert_relative_eq;
+use std::time::Duration;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #1: Solver Trait and Strategy Dispatch
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_newton_config_default() {
+    let cfg = NewtonConfig::default();
+    
+    assert_eq!(cfg.max_iterations, 100);
+    assert_relative_eq!(cfg.tolerance, 1e-6);
+    assert!(!cfg.line_search);
+    assert!(cfg.timeout.is_none());
+    assert!(!cfg.use_numerical_jacobian);
+    assert_relative_eq!(cfg.line_search_armijo_c, 1e-4);
+    assert_eq!(cfg.line_search_max_backtracks, 20);
+    assert_relative_eq!(cfg.divergence_threshold, 1e10);
+}
+
+#[test]
+fn test_newton_config_with_timeout() {
+    let timeout = Duration::from_millis(500);
+    let cfg = NewtonConfig::default().with_timeout(timeout);
+    
+    assert_eq!(cfg.timeout, Some(timeout));
+}
+
+#[test]
+fn test_newton_config_custom_values() {
+    let cfg = NewtonConfig {
+        max_iterations: 50,
+        tolerance: 1e-8,
+        line_search: true,
+        timeout: Some(Duration::from_millis(500)),
+        use_numerical_jacobian: true,
+        line_search_armijo_c: 1e-3,
+        line_search_max_backtracks: 10,
+        divergence_threshold: 1e8,
+        ..Default::default()
+    };
+    
+    assert_eq!(cfg.max_iterations, 50);
+    assert_relative_eq!(cfg.tolerance, 1e-8);
+    assert!(cfg.line_search);
+    assert_eq!(cfg.timeout, Some(Duration::from_millis(500)));
+    assert!(cfg.use_numerical_jacobian);
+    assert_relative_eq!(cfg.line_search_armijo_c, 1e-3);
+    assert_eq!(cfg.line_search_max_backtracks, 10);
+    assert_relative_eq!(cfg.divergence_threshold, 1e8);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #2: Empty System Handling
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_empty_system_returns_invalid() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = NewtonConfig::default();
+    let result = solver.solve(&mut sys);
+    
+    assert!(result.is_err());
+    match result {
+        Err(SolverError::InvalidSystem { message }) => {
+            assert!(message.contains("Empty") || message.contains("no state"));
+        }
+        other => panic!("Expected InvalidSystem, got {:?}", other),
+    }
+}
+
+#[test]
+#[should_panic(expected = "finalize")]
+fn test_empty_system_without_finalize_panics() {
+    // System panics if solve() is called without finalize()
+    // This is expected behavior - the solver requires a finalized system
+    let mut sys = System::new();
+    // Don't call finalize
+
+    let mut solver = NewtonConfig::default();
+    let _ = solver.solve(&mut sys);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #3: Timeout Enforcement
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_timeout_value_in_error() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let timeout_ms = 10u64;
+    let mut solver = NewtonConfig {
+        timeout: Some(Duration::from_millis(timeout_ms)),
+        ..Default::default()
+    };
+
+    let result = solver.solve(&mut sys);
+    
+    // Empty system returns InvalidSystem immediately (before timeout check)
+    assert!(result.is_err());
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #4: Error Types
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_error_display_non_convergence() {
+    let err = SolverError::NonConvergence {
+        iterations: 42,
+        final_residual: 1.23e-3,
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("42"));
+    assert!(msg.contains("1.23"));
+}
+
+#[test]
+fn test_error_display_timeout() {
+    let err = SolverError::Timeout { timeout_ms: 500 };
+    let msg = err.to_string();
+    assert!(msg.contains("500"));
+}
+
+#[test]
+fn test_error_display_divergence() {
+    let err = SolverError::Divergence {
+        reason: "test reason".to_string(),
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("test reason"));
+}
+
+#[test]
+fn test_error_display_invalid_system() {
+    let err = SolverError::InvalidSystem {
+        message: "test message".to_string(),
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("test message"));
+}
+
+#[test]
+fn test_error_equality() {
+    let e1 = SolverError::NonConvergence {
+        iterations: 10,
+        final_residual: 1e-3,
+    };
+    let e2 = SolverError::NonConvergence {
+        iterations: 10,
+        final_residual: 1e-3,
+    };
+    assert_eq!(e1, e2);
+    
+    let e3 = SolverError::Timeout { timeout_ms: 100 };
+    assert_ne!(e1, e3);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #5: Pre-Allocated Buffers (No Panic)
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_solver_does_not_panic_on_empty_system() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = NewtonConfig::default();
+    
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_does_not_panic_with_line_search() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = NewtonConfig {
+        line_search: true,
+        ..Default::default()
+    };
+    
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_does_not_panic_with_numerical_jacobian() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = NewtonConfig {
+        use_numerical_jacobian: true,
+        ..Default::default()
+    };
+    
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #6: ConvergedState
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_converged_state_is_converged() {
+    use entropyk_solver::ConvergenceStatus;
+    use entropyk_solver::ConvergedState;
+    
+    let state = ConvergedState::new(
+        vec![1.0, 2.0, 3.0],
+        10,
+        1e-8,
+        ConvergenceStatus::Converged,
+    );
+    
+    assert!(state.is_converged());
+    assert_eq!(state.iterations, 10);
+    assert_eq!(state.state, vec![1.0, 2.0, 3.0]);
+}
+
+#[test]
+fn test_converged_state_timed_out() {
+    use entropyk_solver::ConvergenceStatus;
+    use entropyk_solver::ConvergedState;
+    
+    let state = ConvergedState::new(
+        vec![1.0],
+        50,
+        1e-3,
+        ConvergenceStatus::TimedOutWithBestState,
+    );
+    
+    assert!(!state.is_converged());
+}

crates/solver/tests/picard_sequential.rs

@@ -0,0 +1,410 @@
+//! Integration tests for Sequential Substitution (Picard) solver (Story 4.3).
+//!
+//! Tests cover:
+//! - AC #1: Reliable convergence when Newton diverges
+//! - AC #2: Sequential variable update
+//! - AC #3: Configurable relaxation factors
+//! - AC #4: Timeout enforcement
+//! - AC #5: Divergence detection
+//! - AC #6: Pre-allocated buffers
+
+use entropyk_solver::{PicardConfig, Solver, SolverError, System};
+use approx::assert_relative_eq;
+use std::time::Duration;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #1: Solver Trait and Configuration
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_picard_config_default() {
+    let cfg = PicardConfig::default();
+
+    assert_eq!(cfg.max_iterations, 100);
+    assert_relative_eq!(cfg.tolerance, 1e-6);
+    assert_relative_eq!(cfg.relaxation_factor, 0.5);
+    assert!(cfg.timeout.is_none());
+    assert_relative_eq!(cfg.divergence_threshold, 1e10);
+    assert_eq!(cfg.divergence_patience, 5);
+}
+
+#[test]
+fn test_picard_config_with_timeout() {
+    let timeout = Duration::from_millis(500);
+    let cfg = PicardConfig::default().with_timeout(timeout);
+
+    assert_eq!(cfg.timeout, Some(timeout));
+}
+
+#[test]
+fn test_picard_config_custom_values() {
+    let cfg = PicardConfig {
+        max_iterations: 200,
+        tolerance: 1e-8,
+        relaxation_factor: 0.3,
+        timeout: Some(Duration::from_millis(1000)),
+        divergence_threshold: 1e8,
+        divergence_patience: 7,
+        ..Default::default()
+    };
+
+    assert_eq!(cfg.max_iterations, 200);
+    assert_relative_eq!(cfg.tolerance, 1e-8);
+    assert_relative_eq!(cfg.relaxation_factor, 0.3);
+    assert_eq!(cfg.timeout, Some(Duration::from_millis(1000)));
+    assert_relative_eq!(cfg.divergence_threshold, 1e8);
+    assert_eq!(cfg.divergence_patience, 7);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #2: Empty System Handling
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_empty_system_returns_invalid() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig::default();
+    let result = solver.solve(&mut sys);
+
+    assert!(result.is_err());
+    match result {
+        Err(SolverError::InvalidSystem { message }) => {
+            assert!(
+                message.contains("Empty") || message.contains("no state"),
+                "Expected empty system message, got: {}",
+                message
+            );
+        }
+        other => panic!("Expected InvalidSystem, got {:?}", other),
+    }
+}
+
+#[test]
+#[should_panic(expected = "finalize")]
+fn test_picard_empty_system_without_finalize_panics() {
+    // System panics if solve() is called without finalize()
+    // This is expected behavior - the solver requires a finalized system
+    let mut sys = System::new();
+    // Don't call finalize
+
+    let mut solver = PicardConfig::default();
+    let _ = solver.solve(&mut sys);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #3: Relaxation Factor Configuration
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_relaxation_factor_default() {
+    let cfg = PicardConfig::default();
+    assert_relative_eq!(cfg.relaxation_factor, 0.5);
+}
+
+#[test]
+fn test_relaxation_factor_full_update() {
+    // omega = 1.0: Full update (fastest, may oscillate)
+    let cfg = PicardConfig {
+        relaxation_factor: 1.0,
+        ..Default::default()
+    };
+    assert_relative_eq!(cfg.relaxation_factor, 1.0);
+}
+
+#[test]
+fn test_relaxation_factor_heavy_damping() {
+    // omega = 0.1: Heavy damping (slow but very stable)
+    let cfg = PicardConfig {
+        relaxation_factor: 0.1,
+        ..Default::default()
+    };
+    assert_relative_eq!(cfg.relaxation_factor, 0.1);
+}
+
+#[test]
+fn test_relaxation_factor_moderate() {
+    // omega = 0.5: Moderate damping (default, good balance)
+    let cfg = PicardConfig {
+        relaxation_factor: 0.5,
+        ..Default::default()
+    };
+    assert_relative_eq!(cfg.relaxation_factor, 0.5);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #4: Timeout Enforcement
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_timeout_value_stored() {
+    let timeout = Duration::from_millis(250);
+    let cfg = PicardConfig::default().with_timeout(timeout);
+
+    assert_eq!(cfg.timeout, Some(timeout));
+}
+
+#[test]
+fn test_timeout_preserves_other_fields() {
+    let cfg = PicardConfig {
+        max_iterations: 150,
+        tolerance: 1e-7,
+        relaxation_factor: 0.25,
+        timeout: None,
+        divergence_threshold: 1e9,
+        divergence_patience: 8,
+        ..Default::default()
+    }
+    .with_timeout(Duration::from_millis(300));
+
+    assert_eq!(cfg.max_iterations, 150);
+    assert_relative_eq!(cfg.tolerance, 1e-7);
+    assert_relative_eq!(cfg.relaxation_factor, 0.25);
+    assert_eq!(cfg.timeout, Some(Duration::from_millis(300)));
+    assert_relative_eq!(cfg.divergence_threshold, 1e9);
+    assert_eq!(cfg.divergence_patience, 8);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #5: Divergence Detection Configuration
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_divergence_threshold_default() {
+    let cfg = PicardConfig::default();
+    assert_relative_eq!(cfg.divergence_threshold, 1e10);
+}
+
+#[test]
+fn test_divergence_patience_default() {
+    let cfg = PicardConfig::default();
+    assert_eq!(cfg.divergence_patience, 5);
+}
+
+#[test]
+fn test_divergence_patience_higher_than_newton() {
+    // Newton uses hardcoded patience of 3
+    // Picard should be more tolerant (5 by default)
+    let cfg = PicardConfig::default();
+    assert!(
+        cfg.divergence_patience >= 5,
+        "Picard divergence_patience ({}) should be >= 5 (more tolerant than Newton's 3)",
+        cfg.divergence_patience
+    );
+}
+
+#[test]
+fn test_divergence_threshold_custom() {
+    let cfg = PicardConfig {
+        divergence_threshold: 1e6,
+        ..Default::default()
+    };
+    assert_relative_eq!(cfg.divergence_threshold, 1e6);
+}
+
+#[test]
+fn test_divergence_patience_custom() {
+    let cfg = PicardConfig {
+        divergence_patience: 10,
+        ..Default::default()
+    };
+    assert_eq!(cfg.divergence_patience, 10);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #6: Pre-Allocated Buffers (No Panic)
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_solver_does_not_panic_on_empty_system() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig::default();
+
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_does_not_panic_with_small_relaxation() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig {
+        relaxation_factor: 0.1,
+        ..Default::default()
+    };
+
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_does_not_panic_with_full_relaxation() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig {
+        relaxation_factor: 1.0,
+        ..Default::default()
+    };
+
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_does_not_panic_with_timeout() {
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig {
+        timeout: Some(Duration::from_millis(10)),
+        ..Default::default()
+    };
+
+    // Should complete without panic
+    let result = solver.solve(&mut sys);
+    assert!(result.is_err());
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Error Types
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_error_display_non_convergence() {
+    let err = SolverError::NonConvergence {
+        iterations: 100,
+        final_residual: 5.67e-4,
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("100"));
+    assert!(msg.contains("5.67"));
+}
+
+#[test]
+fn test_error_display_timeout() {
+    let err = SolverError::Timeout { timeout_ms: 250 };
+    let msg = err.to_string();
+    assert!(msg.contains("250"));
+}
+
+#[test]
+fn test_error_display_divergence() {
+    let err = SolverError::Divergence {
+        reason: "residual increased for 5 consecutive iterations".to_string(),
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("residual increased"));
+}
+
+#[test]
+fn test_error_display_invalid_system() {
+    let err = SolverError::InvalidSystem {
+        message: "State dimension does not match equation count".to_string(),
+    };
+    let msg = err.to_string();
+    assert!(msg.contains("State dimension"));
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// ConvergedState
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_converged_state_is_converged() {
+    use entropyk_solver::{ConvergedState, ConvergenceStatus};
+
+    let state = ConvergedState::new(
+        vec![1.0, 2.0, 3.0],
+        25,
+        1e-7,
+        ConvergenceStatus::Converged,
+    );
+
+    assert!(state.is_converged());
+    assert_eq!(state.iterations, 25);
+    assert_eq!(state.state, vec![1.0, 2.0, 3.0]);
+    assert_relative_eq!(state.final_residual, 1e-7);
+}
+
+#[test]
+fn test_converged_state_timed_out() {
+    use entropyk_solver::{ConvergedState, ConvergenceStatus};
+
+    let state = ConvergedState::new(
+        vec![0.5],
+        75,
+        1e-2,
+        ConvergenceStatus::TimedOutWithBestState,
+    );
+
+    assert!(!state.is_converged());
+    assert_eq!(state.status, ConvergenceStatus::TimedOutWithBestState);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SolverStrategy Integration
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_solver_strategy_picard_dispatch() {
+    use entropyk_solver::SolverStrategy;
+
+    let mut strategy = SolverStrategy::SequentialSubstitution(PicardConfig::default());
+    let mut system = System::new();
+    system.finalize().unwrap();
+
+    let result = strategy.solve(&mut system);
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_solver_strategy_picard_with_timeout() {
+    use entropyk_solver::SolverStrategy;
+
+    let strategy =
+        SolverStrategy::SequentialSubstitution(PicardConfig::default())
+            .with_timeout(Duration::from_millis(100));
+
+    match strategy {
+        SolverStrategy::SequentialSubstitution(cfg) => {
+            assert_eq!(cfg.timeout, Some(Duration::from_millis(100)));
+        }
+        other => panic!("Expected SequentialSubstitution, got {:?}", other),
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Dimension Mismatch Handling
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[test]
+fn test_picard_dimension_mismatch_returns_error() {
+    // Picard requires state dimension == equation count
+    // This is validated in solve() before iteration begins
+    let mut sys = System::new();
+    sys.finalize().unwrap();
+
+    let mut solver = PicardConfig::default();
+    let result = solver.solve(&mut sys);
+
+    // Empty system should return InvalidSystem
+    assert!(result.is_err());
+    match result {
+        Err(SolverError::InvalidSystem { message }) => {
+            assert!(
+                message.contains("Empty") || message.contains("no state"),
+                "Expected empty system message, got: {}",
+                message
+            );
+        }
+        other => panic!("Expected InvalidSystem, got {:?}", other),
+    }
+}

crates/solver/tests/smart_initializer.rs

@@ -0,0 +1,267 @@
+//! Integration tests for Story 4.6: Smart Initialization Heuristic (AC: #8)
+//!
+//! Tests cover:
+//! - AC #8: Integration with FallbackSolver via `with_initial_state`
+//! - Cold-start convergence: SmartInitializer → FallbackSolver
+//! - `initial_state` respected by NewtonConfig and PicardConfig
+//! - `with_initial_state` builder on FallbackSolver delegates to both sub-solvers
+
+use entropyk_components::{Component, ComponentError, JacobianBuilder, ResidualVector, SystemState};
+use entropyk_core::{Enthalpy, Pressure, Temperature};
+use entropyk_solver::{
+    solver::{FallbackSolver, NewtonConfig, PicardConfig, Solver},
+    InitializerConfig, SmartInitializer, System,
+};
+use approx::assert_relative_eq;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Mock Components for Testing
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A simple linear component whose residual is r_i = x_i - target_i.
+/// The solution is x = target. Used to verify initial_state is copied correctly.
+struct LinearTargetSystem {
+    /// Target values (solution)
+    targets: Vec<f64>,
+}
+
+impl LinearTargetSystem {
+    fn new(targets: Vec<f64>) -> Self {
+        Self { targets }
+    }
+}
+
+impl Component for LinearTargetSystem {
+    fn compute_residuals(
+        &self,
+        state: &SystemState,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        for (i, &t) in self.targets.iter().enumerate() {
+            residuals[i] = state[i] - t;
+        }
+        Ok(())
+    }
+
+    fn jacobian_entries(
+        &self,
+        _state: &SystemState,
+        jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        for i in 0..self.targets.len() {
+            jacobian.add_entry(i, i, 1.0);
+        }
+        Ok(())
+    }
+
+    fn n_equations(&self) -> usize {
+        self.targets.len()
+    }
+
+    fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+        &[]
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Helpers
+// ─────────────────────────────────────────────────────────────────────────────
+
+fn build_system_with_targets(targets: Vec<f64>) -> System {
+    let mut sys = System::new();
+    let n0 = sys.add_component(Box::new(LinearTargetSystem::new(targets)));
+    sys.add_edge(n0, n0).unwrap();
+    sys.finalize().unwrap();
+    sys
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// AC #8: Integration with Solver — initial_state accepted via builders
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// AC #8 — `NewtonConfig::with_initial_state` starts from provided state.
+///
+/// We build a 2-entry system where target = [3e5, 4e5].
+/// Starting from zeros → needs to close the gap.
+/// Starting from the exact solution → should converge in 0 additional iterations
+/// (already converged at initial check).
+#[test]
+fn test_newton_with_initial_state_converges_at_target() {
+    // 2-entry state (1 edge × 2 entries: P, h)
+    let targets = vec![300_000.0, 400_000.0];
+    let mut sys = build_system_with_targets(targets.clone());
+
+    let mut solver = NewtonConfig::default().with_initial_state(targets.clone());
+    let result = solver.solve(&mut sys);
+
+    assert!(result.is_ok(), "Should converge: {:?}", result.err());
+    let converged = result.unwrap();
+    // Started exactly at solution → 0 iterations needed
+    assert_eq!(converged.iterations, 0, "Should converge at initial state (0 iterations)");
+    assert!(converged.final_residual < 1e-6);
+}
+
+/// AC #8 — `PicardConfig::with_initial_state` starts from provided state.
+#[test]
+fn test_picard_with_initial_state_converges_at_target() {
+    let targets = vec![300_000.0, 400_000.0];
+    let mut sys = build_system_with_targets(targets.clone());
+
+    let mut solver = PicardConfig::default().with_initial_state(targets.clone());
+    let result = solver.solve(&mut sys);
+
+    assert!(result.is_ok(), "Should converge: {:?}", result.err());
+    let converged = result.unwrap();
+    assert_eq!(converged.iterations, 0, "Should converge at initial state (0 iterations)");
+    assert!(converged.final_residual < 1e-6);
+}
+
+/// AC #8 — `FallbackSolver::with_initial_state` delegates to both newton and picard.
+#[test]
+fn test_fallback_solver_with_initial_state_delegates() {
+    let state = vec![300_000.0, 400_000.0];
+
+    let solver = FallbackSolver::default_solver().with_initial_state(state.clone());
+
+    // Verify both sub-solvers received the initial state
+    assert_eq!(
+        solver.newton_config.initial_state.as_deref(),
+        Some(state.as_slice()),
+        "NewtonConfig should have the initial state"
+    );
+    assert_eq!(
+        solver.picard_config.initial_state.as_deref(),
+        Some(state.as_slice()),
+        "PicardConfig should have the initial state"
+    );
+}
+
+/// AC #8 — `FallbackSolver::with_initial_state` causes early convergence at exact solution.
+#[test]
+fn test_fallback_solver_with_initial_state_at_solution() {
+    let targets = vec![300_000.0, 400_000.0];
+    let mut sys = build_system_with_targets(targets.clone());
+
+    let mut solver = FallbackSolver::default_solver().with_initial_state(targets.clone());
+    let result = solver.solve(&mut sys);
+
+    assert!(result.is_ok(), "Should converge: {:?}", result.err());
+    let converged = result.unwrap();
+    assert_eq!(converged.iterations, 0, "Should converge immediately at initial state");
+}
+
+/// AC #8 — Smart initial state reduces iterations vs. zero initial state.
+///
+/// We use a system where the solution is far from zero (large P, h values).
+/// Newton from zero must close a large gap; Newton from SmartInitializer's output
+/// starts close and should converge in fewer iterations.
+#[test]
+fn test_smart_initializer_reduces_iterations_vs_zero_start() {
+    // System solution: P = 300_000, h = 400_000
+    let targets = vec![300_000.0_f64, 400_000.0_f64];
+
+    // Run 1: from zeros
+    let mut sys_zero = build_system_with_targets(targets.clone());
+    let mut solver_zero = NewtonConfig::default();
+    let result_zero = solver_zero.solve(&mut sys_zero).expect("zero-start should converge");
+
+    // Run 2: from smart initial state (we directly provide the values as an approximation)
+    // Use 95% of target as "smart" initial — simulating a near-correct heuristic
+    let smart_state: Vec<f64> = targets.iter().map(|&t| t * 0.95).collect();
+    let mut sys_smart = build_system_with_targets(targets.clone());
+    let mut solver_smart = NewtonConfig::default().with_initial_state(smart_state);
+    let result_smart = solver_smart.solve(&mut sys_smart).expect("smart-start should converge");
+
+    // Smart start should converge at least as fast (same or fewer iterations)
+    // For a linear system, Newton always converges in 1 step regardless of start,
+    // so both should use ≤ 1 iteration and achieve tolerance
+    assert!(result_zero.final_residual < 1e-6, "Zero start should converge to tolerance");
+    assert!(result_smart.final_residual < 1e-6, "Smart start should converge to tolerance");
+    assert!(
+        result_smart.iterations <= result_zero.iterations,
+        "Smart start ({} iters) should not need more iterations than zero start ({} iters)",
+        result_smart.iterations,
+        result_zero.iterations
+    );
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SmartInitializer API — cold-start pressure estimation
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// AC #8 — SmartInitializer produces pressures and populate_state works end-to-end.
+///
+/// Full integration: estimate pressures → populate state → verify no allocation.
+#[test]
+fn test_cold_start_estimate_then_populate() {
+    let init = SmartInitializer::new(InitializerConfig {
+        fluid: entropyk_components::port::FluidId::new("R134a"),
+        dt_approach: 5.0,
+    });
+
+    let t_source = Temperature::from_celsius(5.0);
+    let t_sink = Temperature::from_celsius(40.0);
+
+    let (p_evap, p_cond) = init
+        .estimate_pressures(t_source, t_sink)
+        .expect("R134a estimation should succeed");
+
+    // Both pressures should be physically reasonable
+    assert!(p_evap.to_bar() > 0.5, "P_evap should be > 0.5 bar");
+    assert!(p_cond.to_bar() > p_evap.to_bar(), "P_cond should exceed P_evap");
+    assert!(p_cond.to_bar() < 50.0, "P_cond should be < 50 bar (not supercritical)");
+
+    // Build a 2-edge system and populate state
+    let mut sys = System::new();
+    let n0 = sys.add_component(Box::new(LinearTargetSystem::new(vec![1.0, 1.0])));
+    let n1 = sys.add_component(Box::new(LinearTargetSystem::new(vec![1.0, 1.0])));
+    let n2 = sys.add_component(Box::new(LinearTargetSystem::new(vec![1.0, 1.0])));
+    sys.add_edge(n0, n1).unwrap();
+    sys.add_edge(n1, n2).unwrap();
+    sys.finalize().unwrap();
+
+    let h_default = Enthalpy::from_joules_per_kg(420_000.0);
+    let mut state = vec![0.0f64; sys.state_vector_len()]; // pre-allocated, no allocation in populate_state
+
+    init.populate_state(&sys, p_evap, p_cond, h_default, &mut state)
+        .expect("populate_state should succeed");
+
+    assert_eq!(state.len(), 4); // 2 edges × [P, h]
+
+    // All edges in single circuit → P_evap used for all
+    assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
+    assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
+    assert_relative_eq!(state[2], p_evap.to_pascals(), max_relative = 1e-9);
+    assert_relative_eq!(state[3], h_default.to_joules_per_kg(), max_relative = 1e-9);
+}
+
+/// AC #8 — Verify initial_state length mismatch falls back gracefully (doesn't panic).
+///
+/// In release mode the solver silently falls back to zeros; in debug mode
+/// debug_assert fires but we can't test that here (it would abort). We verify
+/// the release-mode behavior: a mismatched initial_state causes fallback to zeros
+/// and the solver still converges.
+#[test]
+fn test_initial_state_length_mismatch_fallback() {
+    // System has 2 state entries (1 edge × 2)
+    let targets = vec![300_000.0, 400_000.0];
+    let mut sys = build_system_with_targets(targets.clone());
+
+    // Provide wrong-length initial state (3 instead of 2)
+    // In release mode: solver falls back to zeros, still converges
+    // In debug mode: debug_assert panics — we skip this test in debug
+    #[cfg(not(debug_assertions))]
+    {
+        let wrong_state = vec![1.0, 2.0, 3.0]; // length 3, system needs 2
+        let mut solver = NewtonConfig::default().with_initial_state(wrong_state);
+        let result = solver.solve(&mut sys);
+        // Should still converge (fell back to zeros)
+        assert!(result.is_ok(), "Should converge even with mismatched initial_state in release mode");
+    }
+
+    #[cfg(debug_assertions)]
+    {
+        // In debug mode, skip this test (debug_assert would abort)
+        let _ = (sys, targets); // suppress unused variable warnings
+    }
+}