Snapshot WIP: solver HP epic progress, BPHX/HX physics, BMAD skill refresh.

Capture uncommitted solver robustness work (regularization, domain errors, linear solver lifecycle, tube DP/MSH), web workbench updates, and synced BMAD skills across IDE agent folders before starting BPHX pressure-drop.

Co-authored-by: Cursor <cursoragent@cursor.com>

crates/solver/src/lib.rs

@@ -6,6 +6,25 @@
 //! where components are nodes and flow connections are edges. Edges index into
 //! the solver's state vector (P and h per edge).
 
+// Pre-existing lint debt accumulated before CI was introduced. Allowed so the
+// new CI skeleton can enforce warnings-as-errors on new code. Cleanup is tracked
+// in deferred-work.md.
+#![allow(clippy::redundant_field_names)]
+#![allow(clippy::let_and_return)]
+#![allow(clippy::unnecessary_fallible_conversions)]
+#![allow(clippy::while_let_loop)]
+#![allow(clippy::needless_range_loop)]
+#![allow(clippy::len_zero)]
+#![allow(clippy::useless_conversion)]
+#![allow(clippy::unnecessary_cast)]
+#![allow(clippy::derivable_impls)]
+#![allow(clippy::ptr_arg)]
+#![allow(clippy::manual_map)]
+#![allow(clippy::map_flatten)]
+#![allow(clippy::unnecessary_map_or)]
+#![allow(clippy::useless_asref)]
+#![allow(clippy::too_many_arguments)]
+
 pub mod coupling;
 pub mod criteria;
 pub mod dof;
@@ -14,6 +33,7 @@ pub mod graph;
 pub mod initializer;
 pub mod inverse;
 pub mod jacobian;
+pub mod linear;
 pub mod macro_component;
 pub mod metadata;
 pub mod scaling;
@@ -27,13 +47,16 @@ pub mod topology;
 pub use coupling::{
     compute_coupling_heat, coupling_groups, has_circular_dependencies, ThermalCoupling,
 };
+pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
 pub use dof::{
     align_roles, unspecified_roles, ComponentEquationBlock, DofReport, EquationRole,
     SystemDofBalance, SystemDofError, UnknownKind,
 };
-pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
 pub use entropyk_components::ConnectionError;
 pub use entropyk_core::CircuitId;
+pub use entropyk_solver_core::{
+    ConvergenceReason, DomainViolation, LinearSolver, SolveOutcome, SolverCoreError,
+};
 pub use error::{AddEdgeError, ThermoError, TopologyError};
 pub use initializer::{
     antoine_pressure, AntoineCoefficients, InitializationDiagnostics, InitializationRegime,
@@ -41,6 +64,7 @@ pub use initializer::{
 };
 pub use inverse::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
 pub use jacobian::JacobianMatrix;
+pub use linear::NalgebraLuSolver;
 pub use macro_component::{MacroComponent, MacroComponentSnapshot, PortMapping};
 pub use metadata::SimulationMetadata;
 pub use scaling::{equilibrate, unscale_dx};

crates/solver/src/linear.rs

@@ -0,0 +1,358 @@
+//! `NalgebraLuSolver`: the existing dense nalgebra LU path re-implemented
+//! behind the [`LinearSolver`] lifecycle (FR3, Story 1.4).
+//!
+//! The numeric pipeline is **bit-identical** to [`crate::jacobian::JacobianMatrix::solve`]:
+//! Ruiz equilibration ([`crate::scaling::equilibrate`]) → scaled LU with
+//! partial pivoting → scaled solve `J·x = b` → undo column scaling
+//! ([`crate::scaling::unscale_dx`]) → reject non-finite steps. The difference
+//! is structural: assembly+factorization happen once in `set_linearisation`
+//! and the factorization (with the cached scaling factors) is reused across
+//! `solve_in_place` calls — the typed hook FR4's frozen-Jacobian reuse needs.
+//!
+//! ## Worked example (lifecycle, diffsol-style)
+//!
+//! ```rust
+//! use entropyk_solver::linear::NalgebraLuSolver;
+//! use entropyk_solver_core::LinearSolver;
+//!
+//! let mut lu = NalgebraLuSolver::new();
+//! lu.set_problem(2, 2).unwrap();                       // 1. structure/setup
+//! lu.set_linearisation(&[(0, 0, 2.0), (1, 1, 1.0)])    // 2. assemble + factor
+//!     .unwrap();
+//! let mut rhs = vec![4.0, 3.0];                        // J·x = b
+//! lu.solve_in_place(&mut rhs).unwrap();                // 3. solve in place
+//! assert!((rhs[0] - 2.0).abs() < 1e-10);
+//! assert!((rhs[1] - 3.0).abs() < 1e-10);
+//! ```
+
+use entropyk_solver_core::{LinearSolver, SolverCoreError};
+use nalgebra::DMatrix;
+
+/// Stored factorization plus the equilibration scalings used to build it.
+struct Factorized {
+    d_r: Vec<f64>,
+    d_c: Vec<f64>,
+    lu: nalgebra::LU<f64, nalgebra::Dyn, nalgebra::Dyn>,
+}
+
+/// Dense LU backend (nalgebra) behind the object-safe [`LinearSolver`] trait.
+///
+/// Default backend for the Newton/PTC strategies. Send by construction
+/// (nalgebra dense types are `Send`). Scratch buffers are allocated once in
+/// `set_problem` so `solve_in_place` performs **zero heap allocation**.
+pub struct NalgebraLuSolver {
+    n_rows: usize,
+    n_cols: usize,
+    /// Assembled (unscaled) matrix, kept so `set_matrix` can reuse storage.
+    matrix: Option<DMatrix<f64>>,
+    factor: Option<Factorized>,
+    /// Scratch: scaled rhs / unscaled step (allocated in `set_problem`).
+    delta: Vec<f64>,
+}
+
+impl Default for NalgebraLuSolver {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl NalgebraLuSolver {
+    /// Creates an unconfigured backend (call [`LinearSolver::set_problem`]).
+    pub fn new() -> Self {
+        Self {
+            n_rows: 0,
+            n_cols: 0,
+            matrix: None,
+            factor: None,
+            delta: Vec::new(),
+        }
+    }
+
+    /// Inherent fast path: copy an already-assembled dense matrix and factor
+    /// it (one `copy_from`, no triplet conversion). Strategies in this crate
+    /// use it so the assembly container stays [`crate::jacobian::JacobianMatrix`];
+    /// the trait's triplet-based `set_linearisation` shares the same
+    /// factorization step.
+    pub fn set_matrix(&mut self, matrix: &DMatrix<f64>) -> Result<(), SolverCoreError> {
+        let stored = self.matrix.as_mut().ok_or_else(|| SolverCoreError::Usage {
+            message: "set_problem must be called before set_matrix".to_string(),
+        })?;
+        if stored.nrows() != matrix.nrows() || stored.ncols() != matrix.ncols() {
+            return Err(SolverCoreError::Usage {
+                message: format!(
+                    "matrix shape {}x{} does not match problem {}x{}",
+                    matrix.nrows(),
+                    matrix.ncols(),
+                    stored.nrows(),
+                    stored.ncols()
+                ),
+            });
+        }
+        stored.copy_from(matrix);
+        self.factorize()
+    }
+
+    /// Equilibrate + scale + factorize the stored matrix (square path).
+    fn factorize(&mut self) -> Result<(), SolverCoreError> {
+        let Some(matrix) = self.matrix.as_ref() else {
+            // Unreachable via the public API (guarded by set_problem), but the
+            // zero-panic policy forbids `expect` in library code.
+            self.factor = None;
+            return Err(SolverCoreError::Usage {
+                message: "set_problem must be called before factorizing".to_string(),
+            });
+        };
+        if matrix.nrows() != matrix.ncols() {
+            // Non-square problems have no LU path; solve_in_place reports Usage.
+            self.factor = None;
+            return Ok(());
+        }
+        let n = matrix.nrows();
+        let (d_r, d_c) = crate::scaling::equilibrate(matrix);
+        let mut scaled = matrix.clone();
+        for i in 0..n {
+            for j in 0..n {
+                scaled[(i, j)] *= d_r[i] * d_c[j];
+            }
+        }
+        self.factor = Some(Factorized {
+            d_r,
+            d_c,
+            lu: scaled.lu(),
+        });
+        Ok(())
+    }
+}
+
+impl LinearSolver for NalgebraLuSolver {
+    fn set_problem(
+        &mut self,
+        n_equations: usize,
+        n_unknowns: usize,
+    ) -> Result<(), SolverCoreError> {
+        if n_equations == 0 || n_unknowns == 0 {
+            return Err(SolverCoreError::Usage {
+                message: format!("zero-sized problem {n_equations}x{n_unknowns}"),
+            });
+        }
+        self.n_rows = n_equations;
+        self.n_cols = n_unknowns;
+        self.matrix = Some(DMatrix::zeros(n_equations, n_unknowns));
+        self.factor = None;
+        // Scratch buffer for the zero-allocation solve path.
+        self.delta = vec![0.0; n_unknowns];
+        Ok(())
+    }
+
+    fn set_linearisation(
+        &mut self,
+        entries: &[(usize, usize, f64)],
+    ) -> Result<(), SolverCoreError> {
+        let (n_rows, n_cols) = (self.n_rows, self.n_cols);
+        if self.matrix.is_none() {
+            return Err(SolverCoreError::Usage {
+                message: "set_problem must be called before set_linearisation".to_string(),
+            });
+        }
+        // Validate bounds BEFORE touching the matrix: an out-of-bounds entry
+        // must leave neither a partially-filled matrix nor a stale
+        // factorization behind (audit finding, High).
+        for &(row, col, _) in entries {
+            if row >= n_rows || col >= n_cols {
+                self.factor = None;
+                return Err(SolverCoreError::Usage {
+                    message: format!("entry ({row},{col}) out of bounds {n_rows}x{n_cols}"),
+                });
+            }
+        }
+        let matrix = self.matrix.as_mut().ok_or_else(|| SolverCoreError::Usage {
+            message: "set_problem must be called before set_linearisation".to_string(),
+        })?;
+        matrix.fill(0.0);
+        for &(row, col, value) in entries {
+            matrix[(row, col)] += value;
+        }
+        self.factorize()
+    }
+
+    fn solve_in_place(&mut self, rhs: &mut [f64]) -> Result<(), SolverCoreError> {
+        if rhs.len() != self.n_rows || self.n_rows != self.n_cols {
+            return Err(SolverCoreError::Usage {
+                message: format!(
+                    "solve_in_place needs a square rhs of len {}, got {} (problem {}x{})",
+                    self.n_rows,
+                    rhs.len(),
+                    self.n_rows,
+                    self.n_cols
+                ),
+            });
+        }
+        let n = self.n_rows;
+        // Disjoint field borrows (zero-panic, zero-clone): factorization +
+        // scratch buffers allocated in set_problem. The solve runs fully in
+        // place via nalgebra's `solve_mut` over vector views — **zero heap
+        // allocation** (Story 1.4 audit).
+        let Self { factor, delta, .. } = self;
+        let factor = factor.as_ref().ok_or_else(|| SolverCoreError::Usage {
+            message: "set_linearisation must be called before solve_in_place".to_string(),
+        })?;
+        // Scaled right-hand side: D_r · b (into the pre-allocated scratch).
+        for i in 0..n {
+            delta[i] = rhs[i] * factor.d_r[i];
+        }
+        let mut b_view = nalgebra::DVectorViewMut::from_slice(&mut delta[..n], n);
+        if factor.lu.solve_mut(&mut b_view) {
+            // Undo the column scaling in place: x_i = D_c,i · y_i.
+            for i in 0..n {
+                delta[i] *= factor.d_c[i];
+            }
+            if delta[..n].iter().all(|v| v.is_finite()) {
+                rhs.copy_from_slice(&delta[..n]);
+                Ok(())
+            } else {
+                tracing::warn!(
+                    "LU solve produced a non-finite step - Jacobian may contain NaN/Inf"
+                );
+                Err(SolverCoreError::InvalidSystem {
+                    message: "linear solve produced a non-finite step".to_string(),
+                })
+            }
+        } else {
+            tracing::warn!("LU solve failed - Jacobian may be singular");
+            Err(SolverCoreError::InvalidSystem {
+                message: "linear solve failed (singular factorization)".to_string(),
+            })
+        }
+    }
+
+    fn n(&self) -> usize {
+        self.n_rows
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::jacobian::JacobianMatrix;
+
+    /// Deterministic dense system (well-conditioned, diagonally dominant).
+    fn dense_entries(n: usize) -> Vec<(usize, usize, f64)> {
+        let mut entries = Vec::new();
+        for i in 0..n {
+            for j in 0..n {
+                let v = if i == j {
+                    4.0 + 0.1 * (i as f64)
+                } else {
+                    1.0 / (1.0 + (i as f64 - j as f64).abs())
+                };
+                entries.push((i, j, v));
+            }
+        }
+        entries
+    }
+
+    #[test]
+    fn parity_with_jacobian_matrix_solve() {
+        for n in [1, 2, 5, 12] {
+            let entries = dense_entries(n);
+            let jm = JacobianMatrix::from_builder(&entries, n, n);
+            let residuals: Vec<f64> = (0..n)
+                .map(|i| (i as f64 * 0.7 - 1.3).sin() * 100.0)
+                .collect();
+
+            let legacy = jm.solve(&residuals).expect("legacy solve");
+
+            let mut backend = NalgebraLuSolver::new();
+            backend.set_problem(n, n).unwrap();
+            backend.set_linearisation(&entries).unwrap();
+            let mut rhs: Vec<f64> = residuals.iter().map(|r| -*r).collect();
+            backend.solve_in_place(&mut rhs).unwrap();
+
+            for i in 0..n {
+                assert!(
+                    (legacy[i] - rhs[i]).abs() <= 1e-12 * legacy[i].abs().max(1.0),
+                    "n={n} i={i}: legacy={} backend={}",
+                    legacy[i],
+                    rhs[i]
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn parity_via_set_matrix_and_factorization_reuse() {
+        let n = 6;
+        let entries = dense_entries(n);
+        let jm = JacobianMatrix::from_builder(&entries, n, n);
+        let mut backend = NalgebraLuSolver::new();
+        backend.set_problem(n, n).unwrap();
+        backend.set_matrix(jm.as_matrix()).unwrap();
+
+        // Two different RHS against the SAME stored factorization (the FR4
+        // reuse hook): both must match the legacy path bit-for-bit.
+        for k in 0..2 {
+            let residuals: Vec<f64> = (0..n).map(|i| (i + k) as f64 * 3.1 - 5.0).collect();
+            let legacy = jm.solve(&residuals).expect("legacy solve");
+            let mut rhs: Vec<f64> = residuals.iter().map(|r| -*r).collect();
+            backend.solve_in_place(&mut rhs).unwrap();
+            for i in 0..n {
+                assert!(
+                    (legacy[i] - rhs[i]).abs() <= 1e-12 * legacy[i].abs().max(1.0),
+                    "reuse k={k} i={i}"
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn singular_and_error_paths() {
+        let mut backend = NalgebraLuSolver::new();
+        assert!(backend.set_problem(0, 2).is_err());
+        backend.set_problem(2, 2).unwrap();
+        // Solve before set_linearisation → Usage.
+        assert!(matches!(
+            backend.solve_in_place(&mut [1.0, 2.0]),
+            Err(SolverCoreError::Usage { .. })
+        ));
+        // Out-of-bounds entry → Usage.
+        assert!(backend.set_linearisation(&[(2, 0, 1.0)]).is_err());
+        // Singular matrix → same outcome as the legacy path (whatever nalgebra
+        // flags after equilibration; parity, not an idealized Err).
+        let singular = [(0, 0, 1.0), (0, 1, 2.0), (1, 0, 1.0), (1, 1, 2.0)];
+        let jm = JacobianMatrix::from_builder(&singular, 2, 2);
+        let legacy = jm.solve(&[1.0, 2.0]);
+        backend.set_linearisation(&singular).unwrap();
+        let mut rhs = vec![-1.0, -2.0];
+        match (legacy, backend.solve_in_place(&mut rhs)) {
+            (None, Err(SolverCoreError::InvalidSystem { .. })) => {}
+            (Some(legacy_delta), Ok(())) => {
+                for i in 0..2 {
+                    assert!(
+                        (legacy_delta[i] - rhs[i]).abs() <= 1e-9 * legacy_delta[i].abs().max(1.0),
+                        "parity on singular matrix: legacy={} backend={}",
+                        legacy_delta[i],
+                        rhs[i]
+                    );
+                }
+            }
+            (legacy, backend) => {
+                panic!("divergent behavior: legacy={legacy:?} backend={backend:?}")
+            }
+        }
+        assert_eq!(backend.n(), 2);
+    }
+
+    #[test]
+    fn non_finite_matrix_entries_are_rejected() {
+        let mut backend = NalgebraLuSolver::new();
+        backend.set_problem(2, 2).unwrap();
+        backend
+            .set_linearisation(&[(0, 0, f64::NAN), (1, 1, 1.0)])
+            .unwrap();
+        let mut rhs = vec![1.0, 2.0];
+        assert!(matches!(
+            backend.solve_in_place(&mut rhs),
+            Err(SolverCoreError::InvalidSystem { .. })
+        ));
+    }
+}

crates/solver/src/solver.rs

@@ -7,6 +7,8 @@ use serde::{Deserialize, Serialize};
 use std::time::Duration;
 use thiserror::Error;
 
+use entropyk_solver_core::{ConvergenceReason, SolveOutcome};
+
 use crate::criteria::ConvergenceReport;
 use crate::metadata::SimulationMetadata;
 use crate::system::System;
@@ -81,6 +83,18 @@ pub enum SolverError {
         energy_error: f64,
     },
 
+    /// A recoverable component domain violation (KINSOL `> 0` convention)
+    /// could not be absorbed by step reduction.
+    ///
+    /// Strategies convert recoverable violations into smaller trial steps and
+    /// retry; this error is raised only when that retry budget is exhausted
+    /// (or when no retry mechanism exists at the evaluation site). Recoverable
+    /// exhaustion is an *outcome*, not a crash: it classifies as
+    /// [`ConvergenceReason::DomainViolation`] and surfaces through
+    /// [`Solver::solve_outcome`] as `Ok(SolveOutcome { .. })`.
+    #[error("{0}")]
+    DomainViolation(entropyk_solver_core::DomainViolation),
+
     /// Solver failure with attached post-mortem convergence diagnostics.
     ///
     /// This wrapper keeps the existing concrete error variants unchanged, so
@@ -94,6 +108,18 @@ pub enum SolverError {
     },
 }
 
+/// Converts a solver-core domain violation into the typed solver error.
+///
+/// Manual impl (instead of thiserror's `#[from]`) because thiserror 1.x
+/// treats `#[from]` fields as error sources, which would require
+/// `DomainViolation` to implement `std::error::Error` — it is an alloc-gated
+/// `no_std` outcome carrier with only `Display`, by design.
+impl From<entropyk_solver_core::DomainViolation> for SolverError {
+    fn from(violation: entropyk_solver_core::DomainViolation) -> Self {
+        Self::DomainViolation(violation)
+    }
+}
+
 impl SolverError {
     /// Attach diagnostics to an error without changing its display text.
     pub fn with_diagnostics(self, diagnostics: ConvergenceDiagnostics) -> Self {
@@ -139,6 +165,28 @@ impl SolverError {
             _ => self,
         }
     }
+
+    /// Classifies this error as a typed convergence outcome, if it is one (FR8).
+    ///
+    /// Outcome-style terminations map to their [`ConvergenceReason`]:
+    /// - `NonConvergence` → [`ConvergenceReason::MaxIters`]
+    /// - `Timeout` → [`ConvergenceReason::TimedOut`]
+    /// - `Divergence` → [`ConvergenceReason::Stalled`]
+    /// - `DomainViolation` → [`ConvergenceReason::DomainViolation`]
+    /// - `WithDiagnostics` delegates to the wrapped error.
+    ///
+    /// Hard errors (`InvalidSystem`, `Validation`) return `None`: they remain
+    /// `Err` and never become outcome data.
+    pub fn convergence_reason(&self) -> Option<ConvergenceReason> {
+        match self {
+            SolverError::NonConvergence { .. } => Some(ConvergenceReason::MaxIters),
+            SolverError::Timeout { .. } => Some(ConvergenceReason::TimedOut),
+            SolverError::Divergence { .. } => Some(ConvergenceReason::Stalled),
+            SolverError::DomainViolation(_) => Some(ConvergenceReason::DomainViolation),
+            SolverError::WithDiagnostics { error, .. } => error.convergence_reason(),
+            SolverError::InvalidSystem { .. } | SolverError::Validation { .. } => None,
+        }
+    }
 }
 
 // ─────────────────────────────────────────────────────────────────────────────
@@ -311,6 +359,57 @@ pub trait Solver {
     /// - [`SolverError::InvalidSystem`] — system is ill-formed.
     fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError>;
 
+    /// Solve the system and report the termination as typed outcome data (FR8).
+    ///
+    /// Unlike [`Solver::solve`], a non-converged termination (max iterations,
+    /// stall, timeout) is **not** a hard `Err`: it returns
+    /// `Ok(SolveOutcome { reason, .. })` with the matching
+    /// [`ConvergenceReason`]. Hard errors (invalid system, validation) remain
+    /// `Err`. "No solution" is an outcome, not an exception.
+    ///
+    /// The default implementation delegates to [`Solver::solve`] and
+    /// classifies the result; strategies may override it to report richer
+    /// outcomes (e.g. [`ConvergenceReason::LikelyNoSolution`]).
+    ///
+    /// When the underlying termination carried no iteration/residual payload,
+    /// `iterations` is `0` and `final_residual` is `f64::NAN`.
+    fn solve_outcome(&mut self, system: &mut System) -> Result<SolveOutcome, SolverError> {
+        match self.solve(system) {
+            Ok(state) => {
+                let reason = match state.status {
+                    ConvergenceStatus::TimedOutWithBestState => ConvergenceReason::TimedOut,
+                    ConvergenceStatus::Converged | ConvergenceStatus::ControlSaturation => {
+                        ConvergenceReason::Converged
+                    }
+                };
+                Ok(SolveOutcome {
+                    reason,
+                    iterations: state.iterations,
+                    final_residual: state.final_residual,
+                    state: Some(state.state),
+                })
+            }
+            Err(err) => match err.convergence_reason() {
+                Some(reason) => {
+                    let (iterations, final_residual) = match err.base_error() {
+                        SolverError::NonConvergence {
+                            iterations,
+                            final_residual,
+                        } => (*iterations, *final_residual),
+                        _ => (0, f64::NAN),
+                    };
+                    Ok(SolveOutcome {
+                        reason,
+                        iterations,
+                        final_residual,
+                        state: None,
+                    })
+                }
+                None => Err(err),
+            },
+        }
+    }
+
     /// Configure a time budget for this solver (builder pattern).
     ///
     /// If the solver exceeds `timeout`, it stops and returns
@@ -604,6 +703,15 @@ pub struct IterationDiagnostics {
     /// ($\max_i |r_i|$, the $\ell_\infty$ residual norm).
     #[serde(default)]
     pub max_residual: f64,
+
+    /// Number of recoverable component domain violations (KINSOL `> 0`
+    /// events) absorbed by step reduction during this iteration.
+    ///
+    /// Populated by the Newton line search; `0` for strategies without a
+    /// recoverable-aware retry path or when no trial evaluation left the
+    /// component's physical domain.
+    #[serde(default)]
+    pub recoverable_events: u32,
 }
 
 /// Identifies the dominant (largest-magnitude) entry of a residual vector.

crates/solver/src/strategies/fallback.rs

@@ -31,7 +31,7 @@ use crate::solver::{
 };
 use crate::system::System;
 
-use super::{HomotopyConfig, NewtonConfig, PicardConfig};
+use super::{HomotopyConfig, NewtonConfig, PicardConfig, PtcConfig, TrustRegionConfig};
 
 /// Configuration for the intelligent fallback solver.
 ///
@@ -196,6 +196,25 @@ pub struct FallbackSolver {
     /// `iPRVS`): cheap solvers first, robust continuation only when needed.
     /// `None` (default) preserves the original Newton↔Picard-only behaviour.
     pub homotopy_config: Option<HomotopyConfig>,
+    /// Optional pseudo-transient continuation (PTC) globalization stage.
+    ///
+    /// Invoked when the primary Newton/Picard stages fail, before the homotopy
+    /// stage: PTC regularizes the Newton step with a diagonal δ⁻¹ shift and
+    /// grows the timestep as the residual shrinks (Kelley-Keyes 1998). It
+    /// targets exactly the line-search stagnation observed on fully-coupled
+    /// systems with two-phase pressure drops at high EXV opening. `Some` by
+    /// default (it only runs on primary failure); `None` disables the stage.
+    pub ptc_config: Option<PtcConfig>,
+    /// Optional Powell dogleg trust-region globalization stage.
+    ///
+    /// Invoked when the primary Newton/Picard stages *and* PTC have failed,
+    /// before the homotopy stage. Trust-region escapes the local-minimum-of-‖F‖
+    /// trap by blending Newton and Cauchy steps inside a trust ball of radius
+    /// δ, accepting steps based on the gain ratio ρ of actual vs predicted
+    /// ‖r‖ reduction. It targets the exact failure mode where PTC's timestep
+    /// collapses (Newton step almost orthogonal to ∇‖F‖, but still close to a
+    /// root). `Some` by default; `None` disables the stage.
+    pub trust_region_config: Option<TrustRegionConfig>,
 }
 
 impl FallbackSolver {
@@ -211,6 +230,8 @@ impl FallbackSolver {
             newton_config: NewtonConfig::default(),
             picard_config: PicardConfig::default().with_anderson(3),
             homotopy_config: None,
+            ptc_config: Some(PtcConfig::default()),
+            trust_region_config: Some(TrustRegionConfig::default()),
         }
     }
 
@@ -242,6 +263,20 @@ impl FallbackSolver {
         self
     }
 
+    /// Sets a custom pseudo-transient continuation configuration, or disables
+    /// the stage with `None`. See [`PtcConfig`] for the stage semantics.
+    pub fn with_ptc(mut self, config: Option<PtcConfig>) -> Self {
+        self.ptc_config = config;
+        self
+    }
+
+    /// Sets a custom Powell dogleg trust-region configuration, or disables the
+    /// stage with `None`. See [`TrustRegionConfig`] for the stage semantics.
+    pub fn with_trust_region(mut self, config: Option<TrustRegionConfig>) -> Self {
+        self.trust_region_config = config;
+        self
+    }
+
     /// Sets the initial state for cold-start solving (Story 4.6 — builder pattern).
     ///
     /// Delegates to `newton_config`, `picard_config`, and the optional homotopy
@@ -249,6 +284,12 @@ impl FallbackSolver {
     pub fn with_initial_state(mut self, state: Vec<f64>) -> Self {
         self.newton_config.initial_state = Some(state.clone());
         self.picard_config.initial_state = Some(state.clone());
+        if let Some(ref mut p) = self.ptc_config {
+            p.initial_state = Some(state.clone());
+        }
+        if let Some(ref mut tr) = self.trust_region_config {
+            tr.initial_state = Some(state.clone());
+        }
         if let Some(ref mut h) = self.homotopy_config {
             h.initial_state = Some(state);
         }
@@ -390,7 +431,18 @@ impl FallbackSolver {
                             return Err(SolverError::Timeout { timeout_ms }
                                 .with_optional_diagnostics(child_diagnostics));
                         }
-                        SolverError::Divergence { reason } => {
+                        divergence_like @ (SolverError::Divergence { .. }
+                        | SolverError::DomainViolation(_)) => {
+                            // Recoverable-exhaustion failures (DomainViolation)
+                            // are recovery-eligible, like divergence: switch
+                            // strategy or attempt the PTC/trust-region/homotopy
+                            // cascade instead of propagating immediately.
+                            let reason =
+                                if let SolverError::Divergence { reason } = &divergence_like {
+                                    reason.clone()
+                                } else {
+                                    divergence_like.to_string()
+                                };
                             // Handle divergence based on current solver and state
                             if !self.config.fallback_enabled {
                                 tracing::info!(
@@ -401,8 +453,9 @@ impl FallbackSolver {
                                     reason = reason,
                                     "Divergence detected, fallback disabled"
                                 );
-                                return Err(SolverError::Divergence { reason }
-                                    .with_optional_diagnostics(child_diagnostics));
+                                return Err(
+                                    divergence_like.with_optional_diagnostics(child_diagnostics)
+                                );
                             }
 
                             match state.current_solver {
@@ -497,7 +550,7 @@ impl FallbackSolver {
                                             reason = reason,
                                             "Picard diverged, no more fallbacks available"
                                         );
-                                        return Err(SolverError::Divergence { reason }
+                                        return Err(divergence_like
                                             .with_optional_diagnostics(child_diagnostics));
                                     }
                                 }
@@ -668,6 +721,123 @@ impl FallbackSolver {
         }
     }
 
+    /// Attempts pseudo-transient continuation after the primary Newton/Picard
+    /// stages have failed, before falling back to homotopy continuation.
+    ///
+    /// Returns the PTC result on success; otherwise returns the *primary* error
+    /// (the PTC failure is logged but not surfaced, like the homotopy stage).
+    /// Structural `InvalidSystem` errors are never retried.
+    fn try_ptc_recovery(
+        &self,
+        system: &mut System,
+        primary_err: SolverError,
+        remaining: Option<Duration>,
+    ) -> Result<ConvergedState, SolverError> {
+        if matches!(primary_err.base_error(), SolverError::InvalidSystem { .. }) {
+            return Err(primary_err);
+        }
+
+        let Some(mut ptc) = self.ptc_config.clone() else {
+            return Err(primary_err);
+        };
+
+        // Share the cold start with the primary solvers unless explicitly set.
+        if ptc.initial_state.is_none() {
+            ptc.initial_state = self.newton_config.initial_state.clone();
+        }
+        // Clamp the stage budget to the remaining global time budget: the
+        // cascade shares ONE time budget across all stages — the tighter of
+        // the stage's own timeout and what is left of the global budget wins.
+        ptc.timeout = match (ptc.timeout, remaining) {
+            (Some(t), Some(rem)) => Some(t.min(rem)),
+            (None, rem) => rem,
+            (t, None) => t,
+        };
+
+        tracing::info!(
+            error = %primary_err,
+            "Primary solvers failed; attempting pseudo-transient continuation"
+        );
+
+        match ptc.solve(system) {
+            Ok(converged) => {
+                tracing::info!(
+                    iterations = converged.iterations,
+                    final_residual = converged.final_residual,
+                    "PTC recovered convergence"
+                );
+                Ok(converged)
+            }
+            Err(ptc_err) => {
+                tracing::warn!(
+                    error = %ptc_err,
+                    "PTC also failed; returning primary error"
+                );
+                Err(primary_err)
+            }
+        }
+    }
+
+    /// Attempts the Powell dogleg trust-region stage after PTC has failed,
+    /// before the homotopy stage. Trust-region is more effective than PTC when
+    /// Newton stagnates *near* a root (PTC's timestep collapses there because
+    /// the regularized step is too short to make progress).
+    ///
+    /// Returns the trust-region result on success; otherwise returns the
+    /// *primary* error (the trust-region failure is logged but not surfaced,
+    /// like the homotopy stage). Structural `InvalidSystem` errors are never
+    /// retried.
+    fn try_trust_region_recovery(
+        &self,
+        system: &mut System,
+        primary_err: SolverError,
+        remaining: Option<Duration>,
+    ) -> Result<ConvergedState, SolverError> {
+        if matches!(primary_err.base_error(), SolverError::InvalidSystem { .. }) {
+            return Err(primary_err);
+        }
+
+        let Some(mut tr) = self.trust_region_config.clone() else {
+            return Err(primary_err);
+        };
+
+        // Share the cold start with the primary solvers unless explicitly set.
+        if tr.initial_state.is_none() {
+            tr.initial_state = self.newton_config.initial_state.clone();
+        }
+        // Clamp the stage budget to the remaining global time budget: the
+        // cascade shares ONE time budget across all stages — the tighter of
+        // the stage's own timeout and what is left of the global budget wins.
+        tr.timeout = match (tr.timeout, remaining) {
+            (Some(t), Some(rem)) => Some(t.min(rem)),
+            (None, rem) => rem,
+            (t, None) => t,
+        };
+
+        tracing::info!(
+            error = %primary_err,
+            "PTC failed; attempting Powell dogleg trust-region continuation"
+        );
+
+        match tr.solve(system) {
+            Ok(converged) => {
+                tracing::info!(
+                    iterations = converged.iterations,
+                    final_residual = converged.final_residual,
+                    "Trust-region continuation recovered convergence"
+                );
+                Ok(converged)
+            }
+            Err(tr_err) => {
+                tracing::warn!(
+                    error = %tr_err,
+                    "Trust-region continuation also failed; returning primary error"
+                );
+                Err(primary_err)
+            }
+        }
+    }
+
     /// Attempts Newton-homotopy continuation as a last-resort recovery after the
     /// primary Newton/Picard stages have failed.
     ///
@@ -694,10 +864,14 @@ impl FallbackSolver {
         if homotopy.initial_state.is_none() {
             homotopy.initial_state = self.newton_config.initial_state.clone();
         }
-        // Inherit the remaining global time budget if the stage has none.
-        if homotopy.timeout.is_none() {
-            homotopy.timeout = remaining;
-        }
+        // Clamp the stage budget to the remaining global time budget: the
+        // cascade shares ONE time budget across all stages — the tighter of
+        // the stage's own timeout and what is left of the global budget wins.
+        homotopy.timeout = match (homotopy.timeout, remaining) {
+            (Some(t), Some(rem)) => Some(t.min(rem)),
+            (None, rem) => rem,
+            (t, None) => t,
+        };
 
         tracing::info!(
             error = %primary_err,
@@ -764,7 +938,20 @@ impl Solver for FallbackSolver {
             Ok(converged) => Ok(converged),
             Err(primary_err) => {
                 let remaining = timeout.map(|t| t.saturating_sub(start_time.elapsed()));
-                self.try_homotopy_recovery(system, primary_err, remaining)
+                match self.try_ptc_recovery(system, primary_err.clone(), remaining) {
+                    Ok(converged) => Ok(converged),
+                    Err(ptc_err) => {
+                        let remaining = timeout.map(|t| t.saturating_sub(start_time.elapsed()));
+                        match self.try_trust_region_recovery(system, ptc_err.clone(), remaining) {
+                            Ok(converged) => Ok(converged),
+                            Err(tr_err) => {
+                                let remaining =
+                                    timeout.map(|t| t.saturating_sub(start_time.elapsed()));
+                                self.try_homotopy_recovery(system, tr_err, remaining)
+                            }
+                        }
+                    }
+                }
             }
         }
     }
@@ -772,6 +959,12 @@ impl Solver for FallbackSolver {
     fn with_timeout(mut self, timeout: Duration) -> Self {
         self.newton_config.timeout = Some(timeout);
         self.picard_config.timeout = Some(timeout);
+        if let Some(ref mut p) = self.ptc_config {
+            p.timeout = Some(timeout);
+        }
+        if let Some(ref mut tr) = self.trust_region_config {
+            tr.timeout = Some(timeout);
+        }
         if let Some(ref mut h) = self.homotopy_config {
             h.timeout = Some(timeout);
         }

crates/solver/src/strategies/homotopy.rs

@@ -157,6 +157,7 @@ impl HomotopyConfig {
             jacobian_condition: None,
             max_residual_index,
             max_residual,
+            recoverable_events: 0,
         });
         Some(diagnostics)
     }
@@ -179,13 +180,22 @@ impl HomotopyConfig {
         residuals_h: &mut Vec<f64>,
         jacobian: &mut JacobianMatrix,
         jacobian_builder: &mut JacobianBuilder,
-    ) -> Result<usize, ()> {
+    ) -> Result<usize, Option<SolverError>> {
         let offset = 1.0 - lambda;
 
         for k in 0..self.inner_max_iterations {
             // Evaluate F(x) and form the homotopy residual H = F − (1 − λ)·r0.
-            if system.compute_residuals(state, residuals).is_err() {
-                return Err(());
+            // Recoverable domain violations shrink λ (KINSOL convention, Story
+            // 1.3); fatal errors abort the whole homotopy immediately instead
+            // of burning the λ budget on a bug.
+            if let Err(e) = system.compute_residuals(state, residuals) {
+                return Err(if e.is_recoverable() {
+                    None
+                } else {
+                    Some(SolverError::InvalidSystem {
+                        message: format!("Failed to compute residuals: {:?}", e),
+                    })
+                });
             }
             for i in 0..residuals.len() {
                 residuals_h[i] = residuals[i] - offset * r0[i];
@@ -196,7 +206,7 @@ impl HomotopyConfig {
                 return Ok(k);
             }
             if !norm.is_finite() || norm > self.divergence_threshold {
-                return Err(());
+                return Err(None);
             }
 
             // ∂H/∂x = ∂F/∂x, so the Jacobian of F is used unchanged.
@@ -211,12 +221,18 @@ impl HomotopyConfig {
                 };
                 match JacobianMatrix::numerical(compute, state, residuals.as_slice(), eps) {
                     Ok(jm) => jacobian.as_matrix_mut().copy_from(jm.as_matrix()),
-                    Err(_) => return Err(()),
+                    Err(_) => return Err(None),
                 }
             } else {
                 jacobian_builder.clear();
-                if system.assemble_jacobian(state, jacobian_builder).is_err() {
-                    return Err(());
+                if let Err(e) = system.assemble_jacobian(state, jacobian_builder) {
+                    return Err(if e.is_recoverable() {
+                        None
+                    } else {
+                        Some(SolverError::InvalidSystem {
+                            message: format!("Failed to assemble Jacobian: {:?}", e),
+                        })
+                    });
                 }
                 jacobian.update_from_builder(jacobian_builder.entries());
             }
@@ -224,15 +240,21 @@ impl HomotopyConfig {
             // Solve J·Δx = −H (the solve routine negates the supplied residual).
             let delta = match jacobian.solve(residuals_h.as_slice()) {
                 Some(d) => d,
-                None => return Err(()),
+                None => return Err(None),
             };
 
             apply_newton_step(state, &delta, clipping_mask, 1.0);
         }
 
         // Final convergence check after the last step.
-        if system.compute_residuals(state, residuals).is_err() {
-            return Err(());
+        if let Err(e) = system.compute_residuals(state, residuals) {
+            return Err(if e.is_recoverable() {
+                None
+            } else {
+                Some(SolverError::InvalidSystem {
+                    message: format!("Failed to compute residuals: {:?}", e),
+                })
+            });
         }
         for i in 0..residuals.len() {
             residuals_h[i] = residuals[i] - offset * r0[i];
@@ -240,7 +262,7 @@ impl HomotopyConfig {
         if Self::residual_norm(residuals_h.as_slice()) < self.inner_tolerance {
             Ok(self.inner_max_iterations)
         } else {
-            Err(())
+            Err(None)
         }
     }
 }
@@ -373,13 +395,11 @@ impl Solver for HomotopyConfig {
                     // Step succeeded: gently grow the increment for the next step.
                     dlambda = (dlambda * 1.5).min(max_step);
                 }
-                Err(()) => {
-                    // Step failed: restore and halve the increment, then retry.
+                Err(None) => {
+                    // Step failed (recoverable): restore and halve the increment, then retry.
                     state.copy_from_slice(&state_saved);
                     dlambda *= 0.5;
                     if dlambda < min_lambda_step {
-                        // Report the residual at the restored (last-good) state so
-                        // final_residual matches the state we actually return from.
                         let compute_ok = system.compute_residuals(&state, &mut residuals).is_ok();
                         let final_residual = if compute_ok {
                             Self::residual_norm(&residuals)
@@ -399,6 +419,10 @@ impl Solver for HomotopyConfig {
                         .with_optional_diagnostics(diagnostics));
                     }
                 }
+                Err(Some(fatal)) => {
+                    // Fatal evaluation error: do not burn the λ budget on it.
+                    return Err(fatal);
+                }
             }
         }
 

crates/solver/src/strategies/mod.rs

@@ -23,12 +23,16 @@
 mod fallback;
 mod homotopy;
 mod newton_raphson;
+mod pseudo_transient;
 mod sequential_substitution;
+mod trust_region;
 
 pub use fallback::{FallbackConfig, FallbackSolver};
 pub use homotopy::HomotopyConfig;
 pub use newton_raphson::NewtonConfig;
+pub use pseudo_transient::PtcConfig;
 pub use sequential_substitution::PicardConfig;
+pub use trust_region::TrustRegionConfig;
 
 use crate::solver::{ConvergedState, Solver, SolverError};
 use crate::system::System;

crates/solver/src/strategies/newton_raphson.rs

@@ -14,7 +14,8 @@ use crate::solver::{
     SolverType, TimeoutConfig, VerboseConfig,
 };
 use crate::system::System;
-use entropyk_components::JacobianBuilder;
+use entropyk_components::{ComponentError, JacobianBuilder};
+use entropyk_solver_core::{DomainViolation, LinearSolver};
 
 /// Configuration for the Newton-Raphson solver.
 ///
@@ -180,6 +181,11 @@ impl NewtonConfig {
 
     /// Performs Armijo line search. Returns Some(alpha) if valid step found.
     /// hot path. `state_copy` and `new_residuals` must have appropriate lengths.
+    ///
+    /// Recoverable component domain violations (KINSOL `> 0`, see
+    /// [`ComponentError::is_recoverable`]) shrink the step and are counted in
+    /// `recoverable_events` for the iteration diagnostics; a fatal evaluation
+    /// error aborts the search immediately instead of burning backtracks.
     #[allow(clippy::too_many_arguments)]
     fn line_search(
         &self,
@@ -191,7 +197,8 @@ impl NewtonConfig {
         state_copy: &mut [f64],
         new_residuals: &mut Vec<f64>,
         clipping_mask: &[Option<(f64, f64)>],
-    ) -> Option<f64> {
+        recoverable_events: &mut u32,
+    ) -> Result<Option<f64>, SolverError> {
         let mut alpha: f64 = 1.0;
         state_copy.copy_from_slice(state);
         let gradient_dot_delta = -current_norm;
@@ -199,10 +206,18 @@ impl NewtonConfig {
         for _backtrack in 0..self.line_search_max_backtracks {
             apply_newton_step(state, delta, clipping_mask, alpha);
 
-            if system.compute_residuals(state, new_residuals).is_err() {
+            if let Err(e) = system.compute_residuals(state, new_residuals) {
                 state.copy_from_slice(state_copy);
-                alpha *= 0.5;
-                continue;
+                if e.is_recoverable() {
+                    // Recoverable domain violation (KINSOL > 0): shrink and retry.
+                    alpha *= 0.5;
+                    *recoverable_events += 1;
+                    continue;
+                }
+                // Fatal evaluation error: abort without burning backtracks.
+                return Err(SolverError::InvalidSystem {
+                    message: format!("Failed to compute residuals: {:?}", e),
+                });
             }
 
             let new_norm = Self::residual_norm(new_residuals);
@@ -213,7 +228,7 @@ impl NewtonConfig {
                     new_norm,
                     "Line search accepted"
                 );
-                return Some(alpha);
+                return Ok(Some(alpha));
             }
 
             state.copy_from_slice(state_copy);
@@ -224,7 +239,7 @@ impl NewtonConfig {
             "Line search failed after {} backtracks",
             self.line_search_max_backtracks
         );
-        None
+        Ok(None)
     }
 
     fn finalize_failure_diagnostics(
@@ -330,6 +345,17 @@ impl Solver for NewtonConfig {
         let mut frozen_count: usize = 0;
         let mut force_recompute: bool = true;
 
+        // LinearSolver backend (FR3): factorizes at each fresh assembly and
+        // reuses the factorization (with cached scalings) across frozen
+        // iterations and repeated solves — bit-identical to JacobianMatrix::solve.
+        let mut linear_backend = crate::linear::NalgebraLuSolver::new();
+        linear_backend
+            .set_problem(n_equations, n_state)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to initialize linear backend: {e}"),
+            })?;
+        let mut delta_buf: Vec<f64> = vec![0.0; n_state];
+
         // Cached condition number (for verbose mode when Jacobian frozen)
         let mut cached_condition: Option<f64> = None;
 
@@ -338,17 +364,47 @@ impl Solver for NewtonConfig {
             .map(|i| system.get_solver_bounds_for_state_index(i))
             .collect();
 
-        // Initial residual computation
+        // Initial residual computation. A recoverable domain violation at the
+        // cold start is a typed outcome; fatal errors stay InvalidSystem.
         system
             .compute_residuals(&state, &mut residuals)
-            .map_err(|e| SolverError::InvalidSystem {
-                message: format!("Failed to compute initial residuals: {:?}", e),
+            .map_err(|e| match e {
+                ComponentError::DomainViolation(violation) => {
+                    SolverError::DomainViolation(violation)
+                }
+                fatal => SolverError::InvalidSystem {
+                    message: format!("Failed to compute initial residuals: {:?}", fatal),
+                },
             })?;
 
         let mut current_norm = Self::residual_norm(&residuals);
         best_state.copy_from_slice(&state);
         best_residual = current_norm;
 
+        // Optional per-equation residual dump (debugging aid, controlled by env
+        // var so production runs are unaffected). Prints the 8 largest-magnitude
+        // residuals at the cold start so we can pinpoint which component / state
+        // slot dominates the residual norm — exactly the NLPD "bottleneck
+        // equation" diagnostic but at the solver entry point.
+        if std::env::var("ENTROPYK_DUMP_RESIDUALS").is_ok() {
+            let mut indexed: Vec<(usize, f64)> = residuals
+                .iter()
+                .enumerate()
+                .map(|(i, &r)| (i, r.abs()))
+                .collect();
+            indexed.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+            tracing::info!(
+                residual_norm = current_norm,
+                top_residuals = ?indexed
+                    .iter()
+                    .take(8)
+                    .map(|(i, v)| format!("[{i}]={v:.4e}"))
+                    .collect::<Vec<_>>()
+                    .join(" "),
+                "Initial residual breakdown (top 8 by magnitude)"
+            );
+        }
+
         tracing::debug!(iteration = 0, residual_norm = current_norm, "Initial state");
 
         // Check if already converged
@@ -458,6 +514,14 @@ impl Solver for NewtonConfig {
                     jacobian_matrix.update_from_builder(jacobian_builder.entries());
                 };
 
+                // Feed the LinearSolver backend with the freshly assembled
+                // matrix (assembles + factorizes; frozen iterations skip this).
+                linear_backend
+                    .set_matrix(jacobian_matrix.as_matrix())
+                    .map_err(|e| SolverError::InvalidSystem {
+                        message: format!("Failed to factorize Jacobian: {e}"),
+                    })?;
+
                 frozen_count = 0;
                 force_recompute = false;
 
@@ -487,27 +551,55 @@ impl Solver for NewtonConfig {
                 tracing::debug!(iteration, frozen_count, "Reusing frozen Jacobian");
             }
 
-            // Solve J·Δx = -r
-            let delta = match jacobian_matrix.solve(&residuals) {
-                Some(d) => d,
-                None => {
-                    let failure_diagnostics = self.finalize_failure_diagnostics(
-                        diagnostics.take(),
-                        iteration,
-                        current_norm,
-                        best_residual,
-                        start_time.elapsed().as_millis() as u64,
-                        cached_condition,
-                        Some(state.clone()),
-                    );
-                    return Err(SolverError::Divergence {
-                        reason: "Jacobian is singular".to_string(),
+            // Solve J·Δx = -r through the LinearSolver backend (FR3). The
+            // factorization is fresh when assembled above and reused when the
+            // Jacobian is frozen. Non-square systems keep the legacy
+            // least-squares path.
+            let delta = if n_equations != n_state {
+                match jacobian_matrix.solve(&residuals) {
+                    Some(d) => d,
+                    None => {
+                        let failure_diagnostics = self.finalize_failure_diagnostics(
+                            diagnostics.take(),
+                            iteration,
+                            current_norm,
+                            best_residual,
+                            start_time.elapsed().as_millis() as u64,
+                            cached_condition,
+                            Some(state.clone()),
+                        );
+                        return Err(SolverError::Divergence {
+                            reason: "Jacobian is singular".to_string(),
+                        }
+                        .with_optional_diagnostics(failure_diagnostics));
+                    }
+                }
+            } else {
+                for (d, r) in delta_buf.iter_mut().zip(residuals.iter()) {
+                    *d = -*r;
+                }
+                match linear_backend.solve_in_place(&mut delta_buf) {
+                    Ok(()) => delta_buf.clone(),
+                    Err(_) => {
+                        let failure_diagnostics = self.finalize_failure_diagnostics(
+                            diagnostics.take(),
+                            iteration,
+                            current_norm,
+                            best_residual,
+                            start_time.elapsed().as_millis() as u64,
+                            cached_condition,
+                            Some(state.clone()),
+                        );
+                        return Err(SolverError::Divergence {
+                            reason: "Jacobian is singular".to_string(),
+                        }
+                        .with_optional_diagnostics(failure_diagnostics));
                     }
-                    .with_optional_diagnostics(failure_diagnostics));
                 }
             };
 
             // Apply step with optional line search
+            let mut recoverable_events = 0u32;
             let alpha = if self.line_search {
                 match self.line_search(
                     system,
@@ -518,7 +610,8 @@ impl Solver for NewtonConfig {
                     &mut state_copy,
                     &mut new_residuals,
                     &clipping_mask,
-                ) {
+                    &mut recoverable_events,
+                )? {
                     Some(a) => a,
                     None => {
                         let failure_diagnostics = self.finalize_failure_diagnostics(
@@ -530,6 +623,17 @@ impl Solver for NewtonConfig {
                             cached_condition,
                             Some(state.clone()),
                         );
+                        // Backtracks exhausted on recoverable domain violations:
+                        // a typed outcome, not an untyped stall (AC #2).
+                        if recoverable_events > 0 {
+                            return Err(SolverError::DomainViolation(DomainViolation {
+                                component: None,
+                                detail:
+                                    "line search exhausted after recoverable domain violation(s)"
+                                        .into(),
+                            })
+                            .with_optional_diagnostics(failure_diagnostics));
+                        }
                         return Err(SolverError::Divergence {
                             reason: "Line search failed".to_string(),
                         }
@@ -543,8 +647,13 @@ impl Solver for NewtonConfig {
 
             system
                 .compute_residuals(&state, &mut residuals)
-                .map_err(|e| SolverError::InvalidSystem {
-                    message: format!("Failed to compute residuals: {:?}", e),
+                .map_err(|e| match e {
+                    ComponentError::DomainViolation(violation) => {
+                        SolverError::DomainViolation(violation)
+                    }
+                    fatal => SolverError::InvalidSystem {
+                        message: format!("Failed to compute residuals: {:?}", fatal),
+                    },
                 })?;
 
             previous_norm = current_norm;
@@ -606,6 +715,7 @@ impl Solver for NewtonConfig {
                     jacobian_condition: cached_condition,
                     max_residual_index,
                     max_residual,
+                    recoverable_events,
                 });
             }
 

crates/solver/src/strategies/pseudo_transient.rs

@@ -0,0 +1,334 @@
+//! Pseudo-transient continuation (PTC / Ψtc) globalization.
+//!
+//! When Newton-with-line-search stagnates at a local minimum of ‖F‖ — the
+//! classic failure on fully-coupled vapor-compression systems with two-phase
+//! pressure drops at high EXV opening — PTC walks the physical relaxation
+//! dynamics instead of minimizing ‖F‖ along a line. Each iteration is one
+//! linearized implicit Euler step of the fictitious dynamics `x′ = −F(x)`:
+//!
+//! ```text
+//!   (δ⁻¹·I + J(x))·s = −F(x),   x ← x + s
+//! ```
+//!
+//! with the timestep δ grown by the Switched Evolution Relaxation (SER) rule
+//! as the residual shrinks. Small δ far from the solution reconditions the
+//! near-singular momentum block (orifice ↔ two-phase ΔP) and bounds the step;
+//! δ → ∞ recovers exact Newton with quadratic terminal convergence.
+//!
+//! References: Kelley & Keyes, *Convergence Analysis of Pseudo-transient
+//! Continuation*, SIAM J. Numer. Anal. 35:508–523 (1998); Mulder & Van Leer
+//! (1985) for the SER rule; PETSc `TSPSEUDO` for the incremental SER variant.
+
+use crate::jacobian::JacobianMatrix;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{apply_newton_step, ConvergedState, ConvergenceStatus, Solver, SolverError};
+use crate::system::System;
+use entropyk_components::JacobianBuilder;
+use entropyk_solver_core::LinearSolver;
+use std::time::{Duration, Instant};
+
+/// Configuration for pseudo-transient continuation.
+///
+/// Used as a globalization stage of [`super::FallbackSolver`]: only invoked
+/// when the primary Newton/Picard stages fail, so it never perturbs problems
+/// that already converge.
+#[derive(Debug, Clone, PartialEq)]
+pub struct PtcConfig {
+    /// Maximum PTC iterations.
+    pub max_iterations: usize,
+    /// Convergence tolerance on the residual L2 norm.
+    pub tolerance: f64,
+    /// Initial timestep δ₀ (in scaled-residual units). Small δ₀ regularizes
+    /// the first steps; too small wastes iterations in the induction phase.
+    pub delta0: f64,
+    /// SER growth factor (> 1). Controls how fast δ grows as ‖F‖ shrinks.
+    pub growth: f64,
+    /// Maximum timestep before the diagonal shift is negligible (pure Newton).
+    pub delta_max: f64,
+    /// Timestep shrink factor on rejected steps (< 1).
+    pub shrink: f64,
+    /// Minimum timestep: if δ falls below this after repeated rejections the
+    /// method has failed (detectable failure mode, Kelley-Keyes Thm 2.3).
+    pub delta_min: f64,
+    /// Optional initial state (cold start shared with the primary solvers).
+    pub initial_state: Option<Vec<f64>>,
+    /// Optional timeout.
+    pub timeout: Option<Duration>,
+}
+
+impl Default for PtcConfig {
+    fn default() -> Self {
+        Self {
+            max_iterations: 200,
+            tolerance: 1e-6,
+            delta0: 1e-3,
+            growth: 1.5,
+            delta_max: 1e9,
+            shrink: 0.25,
+            delta_min: 1e-12,
+            initial_state: None,
+            timeout: None,
+        }
+    }
+}
+
+impl PtcConfig {
+    /// L2 residual norm, consistent with the other strategies.
+    fn residual_norm(residuals: &[f64]) -> f64 {
+        let norm = residuals.iter().map(|r| r * r).sum::<f64>().sqrt();
+        if norm.is_finite() {
+            norm
+        } else {
+            f64::MAX
+        }
+    }
+}
+
+impl Solver for PtcConfig {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+
+        let n_state = system.full_state_vector_len();
+        let n_equations: usize = system
+            .traverse_for_jacobian()
+            .map(|(_, c, _)| c.n_equations())
+            .sum::<usize>()
+            + system.constraints().count()
+            + system.coupling_residual_count()
+            + 2 * system.saturated_controller_count()
+            + system.mass_flow_closure_count();
+
+        if n_state == 0 || n_equations == 0 {
+            return Err(SolverError::InvalidSystem {
+                message: "Empty system has no state variables or equations".to_string(),
+            });
+        }
+
+        let mut state: Vec<f64> = match self.initial_state.as_ref() {
+            Some(s) if s.len() == n_state => s.clone(),
+            Some(s) => {
+                return Err(SolverError::InvalidSystem {
+                    message: format!(
+                        "initial_state length {} does not match system state length {}",
+                        s.len(),
+                        n_state
+                    ),
+                });
+            }
+            None => vec![0.0; n_state],
+        };
+        let mut residuals: Vec<f64> = vec![0.0; n_equations];
+        let mut saved_state: Vec<f64> = vec![0.0; n_state];
+        let mut jacobian_builder = JacobianBuilder::new();
+        let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state);
+
+        // LinearSolver backend (FR3): factorizes at each fresh assembly;
+        // bit-identical to JacobianMatrix::solve.
+        let mut linear_backend = crate::linear::NalgebraLuSolver::new();
+        linear_backend
+            .set_problem(n_equations, n_state)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to initialize linear backend: {e}"),
+            })?;
+        let mut step_buf: Vec<f64> = vec![0.0; n_state];
+
+        let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
+            .map(|i| system.get_solver_bounds_for_state_index(i))
+            .collect();
+
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute initial residuals: {:?}", e),
+            })?;
+        let mut current_norm = Self::residual_norm(&residuals);
+        let mut prev_norm = current_norm;
+
+        tracing::info!(
+            max_iterations = self.max_iterations,
+            tolerance = self.tolerance,
+            delta0 = self.delta0,
+            residual_norm = current_norm,
+            "Pseudo-transient continuation starting"
+        );
+
+        if current_norm < self.tolerance {
+            return Ok(ConvergedState::new(
+                state,
+                0,
+                current_norm,
+                ConvergenceStatus::Converged,
+                SimulationMetadata::new(system.input_hash()),
+            ));
+        }
+
+        let mut delta = self.delta0;
+
+        for iteration in 1..=self.max_iterations {
+            if let Some(timeout) = self.timeout {
+                if start_time.elapsed() > timeout {
+                    tracing::info!(iteration, "PTC timed out");
+                    return Err(SolverError::Timeout {
+                        timeout_ms: timeout.as_millis() as u64,
+                    });
+                }
+            }
+
+            // Fresh Jacobian each iteration (δ changes every step anyway).
+            jacobian_builder.clear();
+            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());
+
+            // Diagonal shift: (δ⁻¹·I + J)·s = −F.
+            let shift = 1.0 / delta;
+            {
+                let n = jacobian_matrix.as_matrix().nrows();
+                let m = jacobian_matrix.as_matrix_mut();
+                for i in 0..n {
+                    m[(i, i)] += shift;
+                }
+            }
+
+            // Solve through the LinearSolver backend (FR3); non-square keeps
+            // the legacy least-squares path.
+            let step = if n_equations == n_state {
+                linear_backend
+                    .set_matrix(jacobian_matrix.as_matrix())
+                    .map_err(|e| SolverError::InvalidSystem {
+                        message: format!("Failed to factorize Jacobian: {e}"),
+                    })?;
+                for (d, r) in step_buf.iter_mut().zip(residuals.iter()) {
+                    *d = -*r;
+                }
+                linear_backend
+                    .solve_in_place(&mut step_buf)
+                    .ok()
+                    .map(|_| step_buf.clone())
+            } else {
+                jacobian_matrix.solve(&residuals)
+            };
+            let Some(step) = step else {
+                // Singular even with the shift: shrink and retry without
+                // consuming an iteration of state progress.
+                delta *= self.shrink;
+                if delta < self.delta_min {
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+                continue;
+            };
+
+            saved_state.copy_from_slice(&state);
+            // No line search: the δ⁻¹ term already bounds the step.
+            apply_newton_step(&mut state, &step, &clipping_mask, 1.0);
+
+            if let Err(e) = system.compute_residuals(&state, &mut residuals) {
+                if !e.is_recoverable() {
+                    // Fatal evaluation error: abort immediately, do not burn
+                    // the shrink budget on it.
+                    return Err(SolverError::InvalidSystem {
+                        message: format!("Failed to compute residuals: {:?}", e),
+                    });
+                }
+                // Recoverable domain violation (KINSOL > 0): reject and shrink
+                // (residual undefined).
+                state.copy_from_slice(&saved_state);
+                delta *= self.shrink;
+                if delta < self.delta_min {
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+                continue;
+            }
+            current_norm = Self::residual_norm(&residuals);
+
+            tracing::debug!(
+                iteration,
+                residual_norm = current_norm,
+                delta,
+                "PTC iteration"
+            );
+
+            if current_norm < self.tolerance {
+                tracing::info!(
+                    iterations = iteration,
+                    final_residual = current_norm,
+                    "PTC converged"
+                );
+                return Ok(ConvergedState::new(
+                    state,
+                    iteration,
+                    current_norm,
+                    ConvergenceStatus::Converged,
+                    SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+
+            if current_norm > prev_norm {
+                // Reject: restore and cut the timestep.
+                state.copy_from_slice(&saved_state);
+                delta *= self.shrink;
+                tracing::debug!(iteration, delta, "PTC step rejected, shrinking timestep");
+                if delta < self.delta_min {
+                    tracing::warn!(
+                        iteration,
+                        final_residual = current_norm,
+                        "PTC timestep collapsed — globalization failed"
+                    );
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+                continue;
+            }
+
+            // SER incremental rule (PETSc form), bounded by delta_max.
+            delta = (self.growth * delta * prev_norm / current_norm).min(self.delta_max);
+            prev_norm = current_norm;
+        }
+
+        Err(SolverError::NonConvergence {
+            iterations: self.max_iterations,
+            final_residual: current_norm,
+        })
+    }
+
+    fn with_timeout(mut self, timeout: Duration) -> Self {
+        self.timeout = Some(timeout);
+        self
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_ptc_config_defaults() {
+        let cfg = PtcConfig::default();
+        assert_eq!(cfg.max_iterations, 200);
+        assert!(cfg.delta0 > 0.0 && cfg.delta0 < 1.0);
+        assert!(cfg.growth > 1.0);
+        assert!(cfg.delta_max > cfg.delta0);
+        assert!(cfg.shrink > 0.0 && cfg.shrink < 1.0);
+        assert!(cfg.delta_min > 0.0 && cfg.delta_min < cfg.delta0);
+    }
+
+    #[test]
+    fn test_ptc_rejects_empty_system() {
+        let mut system = System::new();
+        system.finalize().unwrap();
+        let mut solver = PtcConfig::default();
+        let result = solver.solve(&mut system);
+        assert!(matches!(result, Err(SolverError::InvalidSystem { .. })));
+    }
+}

crates/solver/src/strategies/sequential_substitution.rs

@@ -15,6 +15,7 @@ use crate::solver::{
     IterationDiagnostics, Solver, SolverError, SolverType, TimeoutConfig, VerboseConfig,
 };
 use crate::system::System;
+use entropyk_components::ComponentError;
 
 /// Configuration for the Sequential Substitution (Picard iteration) solver.
 ///
@@ -338,11 +339,17 @@ impl Solver for PicardConfig {
         let mut best_state: Vec<f64> = vec![0.0; n_state];
         let mut best_residual: f64;
 
-        // Initial residual computation
+        // Initial residual computation. A recoverable domain violation is a
+        // typed outcome (Picard has no step to shrink); fatal stays InvalidSystem.
         system
             .compute_residuals(&state, &mut residuals)
-            .map_err(|e| SolverError::InvalidSystem {
-                message: format!("Failed to compute initial residuals: {:?}", e),
+            .map_err(|e| match e {
+                ComponentError::DomainViolation(violation) => {
+                    SolverError::DomainViolation(violation)
+                }
+                fatal => SolverError::InvalidSystem {
+                    message: format!("Failed to compute initial residuals: {:?}", fatal),
+                },
             })?;
 
         let mut current_norm = Self::residual_norm(&residuals);
@@ -419,11 +426,17 @@ impl Solver for PicardConfig {
                 Self::apply_relaxation(&mut state, &residuals, self.relaxation_factor);
             }
 
-            // Compute new residuals
+            // Compute new residuals. Recoverable domain violation → typed
+            // outcome; fatal → InvalidSystem flattening.
             system
                 .compute_residuals(&state, &mut residuals)
-                .map_err(|e| SolverError::InvalidSystem {
-                    message: format!("Failed to compute residuals: {:?}", e),
+                .map_err(|e| match e {
+                    ComponentError::DomainViolation(violation) => {
+                        SolverError::DomainViolation(violation)
+                    }
+                    fatal => SolverError::InvalidSystem {
+                        message: format!("Failed to compute residuals: {:?}", fatal),
+                    },
                 })?;
 
             previous_norm = current_norm;
@@ -471,6 +484,7 @@ impl Solver for PicardConfig {
                     jacobian_condition: None, // No Jacobian in Picard
                     max_residual_index,
                     max_residual,
+                    recoverable_events: 0, // Picard has no step-reduction retry
                 });
             }
 

crates/solver/src/strategies/trust_region.rs

@@ -0,0 +1,717 @@
+//! Powell dogleg trust-region globalization.
+//!
+//! When Newton-with-line-search stagnates near a local minimum of ‖F‖ — the
+//! signature failure of fully-coupled vapor-compression systems with two-phase
+//! pressure drops at high EXV opening, where the Newton direction is almost
+//! orthogonal to the gradient of ‖F‖ and Armijo backtracking collapses to
+//! infinitesimal steps — trust-region globalization escapes by picking the
+//! step inside a ball of radius δ that best reduces the linear model
+//! m(p) = ½‖r + J·p‖². The classic Powell dogleg blends the Newton step n
+//! (J·n = −r) and the Cauchy (steepest-descent) step c inside the trust region;
+//! when the Jacobian is singular or ill-conditioned, the step falls back to a
+//! Levenberg–Marquardt solve (JᵀJ + λI)·p = −Jᵀr, clipped to the region.
+//! The radius δ is grown/shrunk from the gain ratio ρ of actual vs predicted
+//! ‖r‖ reduction. The method reuses Entropyk's analytical Jacobian — no finite
+//! differences — and the same Ruiz-equilibrated LU path as `NewtonConfig`.
+//!
+//! References:
+//! - Nocedal & Wright, *Numerical Optimization* (2nd ed.), ch. 4.6 (dogleg)
+//! - Madsen, Nielsen & Tingleff, *Methods for Non-Linear Least Squares
+//!   Problems* (2004) — LM with gain ratio
+//! - Dennis & Schnabel, *Numerical Methods for Unconstrained Optimization*
+//! - GSL `gsl_multifit_nlinear` hybrid scaled algorithm
+//! - gomez crate (MIT): <https://github.com/datamole-ai/gomez>
+
+use std::time::{Duration, Instant};
+
+use crate::jacobian::JacobianMatrix;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{apply_newton_step, ConvergedState, ConvergenceStatus, Solver, SolverError};
+use crate::system::System;
+use entropyk_components::JacobianBuilder;
+use nalgebra::DVector;
+
+/// Configuration for the Powell dogleg trust-region solver.
+///
+/// Used as a globalization stage of [`super::FallbackSolver`]: only invoked
+/// when the primary Newton/Picard stages fail, so it never perturbs problems
+/// that already converge. The defaults mirror the well-tested gomez/GSL
+/// parameterization.
+#[derive(Debug, Clone, PartialEq)]
+pub struct TrustRegionConfig {
+    /// Maximum trust-region iterations.
+    pub max_iterations: usize,
+    /// Convergence tolerance on the residual L2 norm.
+    pub tolerance: f64,
+    /// Initial trust-region radius δ₀ (in column-scaled variables).
+    pub delta_init: f64,
+    /// Minimum radius: if δ collapses below this the method has failed.
+    pub delta_min: f64,
+    /// Maximum radius (effectively unbounded above this).
+    pub delta_max: f64,
+    /// Gain ratio ρ below which a step is rejected and δ shrinks.
+    pub accept_threshold: f64,
+    /// ρ below which δ is shrunk even on accepted steps (poor model match).
+    pub shrink_threshold: f64,
+    /// ρ above which δ is expanded (only if ‖D·p‖ ≈ δ, i.e. the step hit the boundary).
+    pub expand_threshold: f64,
+    /// δ shrink multiplier (< 1) on rejected steps.
+    pub shrink_factor: f64,
+    /// δ expand multiplier (> 1) on the boundary with good gain.
+    pub expand_factor: f64,
+    /// Consecutive rejections before declaring failure at the current point.
+    pub max_rejections: usize,
+    /// Initial Levenberg–Marquardt λ for the singular-Newton fallback.
+    pub lm_lambda_init: f64,
+    /// LM λ upper bound (solves become pure gradient descent above this).
+    pub lm_lambda_max: f64,
+    /// LM λ growth on a rejected LM step.
+    pub lm_lambda_grow: f64,
+    /// LM λ shrink on an accepted LM step.
+    pub lm_lambda_shrink: f64,
+    /// Optional initial state (cold start shared with the primary solvers).
+    pub initial_state: Option<Vec<f64>>,
+    /// Optional timeout.
+    pub timeout: Option<Duration>,
+}
+
+impl Default for TrustRegionConfig {
+    fn default() -> Self {
+        Self {
+            max_iterations: 200,
+            tolerance: 1e-6,
+            delta_init: 1.0,
+            delta_min: f64::EPSILON.sqrt(), // ≈ 1.5e-8 (gomez)
+            delta_max: 1e9,
+            accept_threshold: 1e-4,
+            shrink_threshold: 0.25,
+            expand_threshold: 0.75,
+            shrink_factor: 0.5, // gomez / Madsen-Tingleff μ=0.5 on rejection
+            expand_factor: 2.0,
+            max_rejections: 20,
+            lm_lambda_init: 1e-3,
+            lm_lambda_max: 1e10,
+            lm_lambda_grow: 10.0,
+            lm_lambda_shrink: 0.1,
+            initial_state: None,
+            timeout: None,
+        }
+    }
+}
+
+impl TrustRegionConfig {
+    /// L2 residual norm, consistent with the other strategies. Falls back to
+    /// `f64::MAX` on non-finite values so the divergence guard fires cleanly
+    /// instead of poisoning the gain ratio with NaN.
+    fn residual_norm(residuals: &[f64]) -> f64 {
+        let norm = residuals.iter().map(|r| r * r).sum::<f64>().sqrt();
+        if norm.is_finite() {
+            norm
+        } else {
+            f64::MAX
+        }
+    }
+
+    /// Column-scaling diagonal `d_j = 1 / ‖col_j(J)‖` (GSL-style). Returns the
+    /// vector of `d_j` so the trust-region metric `‖D·p‖` treats every state
+    /// variable on the same scale — critical for thermodynamic systems that
+    /// mix P (~1e6 Pa), h (~1e5 J/kg) and dimensionless controls (~1) in the
+    /// same state vector. Zero columns fall back to `d_j = 1` to avoid div-by-0.
+    fn column_scaling(jacobian: &JacobianMatrix) -> Vec<f64> {
+        let m = jacobian.as_matrix();
+        let ncols = m.ncols();
+        let mut d = Vec::with_capacity(ncols);
+        for j in 0..ncols {
+            let col_norm = m.column(j).norm();
+            d.push(if col_norm > f64::EPSILON.sqrt() {
+                1.0 / col_norm
+            } else {
+                1.0
+            });
+        }
+        d
+    }
+
+    /// Weighted norm `‖D·v‖` used as the trust-region metric.
+    fn scaled_norm(v: &[f64], d: &[f64]) -> f64 {
+        let mut acc = 0.0_f64;
+        for (vi, di) in v.iter().zip(d.iter()) {
+            acc += (vi * di).powi(2);
+        }
+        acc.sqrt()
+    }
+
+    /// Computes the dogleg step blending the Cauchy point `c` and Newton step
+    /// `n` so that `‖D·p‖ ≤ δ`. Returns the step in the same coordinates as
+    /// the Newton step (i.e. in x-space, not D-scaled).
+    ///
+    /// Algorithm (Nocedal & Wright, Procedure 4.3, scaled form):
+    ///   * If `‖D·n‖ ≤ δ`, return `n` (full Newton step is inside the region).
+    ///   * Else if `‖D·c‖ ≥ δ`, clip the Cauchy step: `p = (δ/‖D·c‖)·c`.
+    ///   * Else find the unique τ ∈ (0,1] with `‖D·(c + τ(n − c))‖ = δ` by
+    ///     solving the quadratic `a·τ² + b·τ + c = 0` and return the linear
+    ///     blend on the segment from c to n.
+    ///
+    /// All inputs are in x-space; only the trust-region *radius* is measured
+    /// in D-scaled coordinates.
+    #[allow(clippy::too_many_arguments)]
+    fn dogleg_step(
+        newton: Option<&[f64]>,
+        cauchy: &[f64],
+        d: &[f64],
+        delta: f64,
+        buf_dp: &mut [f64],
+        buf_dn: &mut [f64],
+    ) -> Vec<f64> {
+        // Case 1: Newton step valid and inside the region.
+        if let Some(n) = newton {
+            let norm_dn = {
+                for (i, &ni) in n.iter().enumerate() {
+                    buf_dn[i] = ni * d[i];
+                }
+                buf_dn.iter().map(|x| x * x).sum::<f64>().sqrt()
+            };
+            if norm_dn <= delta {
+                return n.to_vec();
+            }
+        }
+
+        // Case 2: outside Newton — clip Cauchy if needed.
+        let norm_dc = {
+            for (i, &ci) in cauchy.iter().enumerate() {
+                buf_dp[i] = ci * d[i];
+            }
+            buf_dp.iter().map(|x| x * x).sum::<f64>().sqrt()
+        };
+        if norm_dc >= delta || newton.is_none() {
+            let alpha = if norm_dc > 0.0 { delta / norm_dc } else { 0.0 };
+            return cauchy.iter().map(|c| c * alpha).collect();
+        }
+
+        // Case 3: dogleg on the segment c → n. Solve the quadratic for τ ≥ 0
+        // such that ‖D·(c + τ(n − c))‖² = δ².
+        let n = newton.unwrap();
+        let mut a = 0.0_f64;
+        let mut b = 0.0_f64;
+        let mut c_coef = 0.0_f64;
+        for j in 0..n.len() {
+            let dj = d[j];
+            let cj = cauchy[j] * dj;
+            let nj = n[j] * dj;
+            let mj = nj - cj; // direction of the segment in scaled space
+            a += mj * mj;
+            b += 2.0 * cj * mj;
+            c_coef += cj * cj;
+        }
+        let rhs = c_coef - delta * delta;
+        // Solve a·τ² + b·τ + rhs = 0 for the smallest non-negative root. We
+        // want τ ∈ (0,1]: that is the intersection of the segment with the
+        // sphere surface. If a is ~0 (n ≈ c numerically) fall back to Cauchy.
+        let tau = if a.abs() < f64::EPSILON {
+            1.0
+        } else {
+            let disc = b * b - 4.0 * a * rhs;
+            if disc < 0.0 {
+                // Quadratic has no real root: the whole segment is inside or
+                // outside. Numerically degenerate — fall back to Cauchy clip.
+                return cauchy.to_vec();
+            }
+            let sqrt_disc = disc.sqrt();
+            let tau1 = (-b + sqrt_disc) / (2.0 * a);
+            let tau2 = (-b - sqrt_disc) / (2.0 * a);
+            // Pick the smallest non-negative root in (0,1].
+            let mut candidates = [tau1, tau2]
+                .iter()
+                .copied()
+                .filter(|&t| t.is_finite() && t > 0.0 && t <= 1.0 + 1e-12)
+                .collect::<Vec<_>>();
+            candidates.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
+            candidates.first().copied().unwrap_or(1.0).clamp(0.0, 1.0)
+        };
+        cauchy
+            .iter()
+            .zip(n.iter())
+            .map(|(c, n)| c + tau * (n - c))
+            .collect()
+    }
+
+    /// Solves the Levenberg–Marquardt normal equations `(JᵀJ + λI)·p = −Jᵀr`
+    /// via LU. Returns `None` if the regularized system is still singular.
+    fn lm_step(jacobian: &JacobianMatrix, residuals: &[f64], lambda: f64) -> Option<Vec<f64>> {
+        let m = jacobian.as_matrix();
+        let n = m.ncols();
+        let jt = m.transpose();
+        let mut regularized = &jt * m;
+        // Add λ on the diagonal of the normal-equations matrix. nalgebra
+        // exposes `diagonal_mut` only on owned `OMatrix`; we already have one
+        // (result of `&jt * m` is owned), so the call works through DMatrix.
+        for j in 0..n {
+            regularized[(j, j)] += lambda;
+        }
+        let jtr = DVector::from_iterator(
+            n,
+            (0..n).map(|j| {
+                let col = jt.column(j);
+                -(col.dot(&DVector::from_row_slice(residuals)))
+            }),
+        );
+        let lu = regularized.lu();
+        lu.solve(&jtr)
+            .map(|v| {
+                let p: Vec<f64> = v.iter().copied().collect();
+                if p.iter().all(|x| x.is_finite()) {
+                    Some(p)
+                } else {
+                    None
+                }
+            })
+            .flatten()
+    }
+}
+
+impl Solver for TrustRegionConfig {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+
+        let n_state = system.full_state_vector_len();
+        let n_equations: usize = system
+            .traverse_for_jacobian()
+            .map(|(_, c, _)| c.n_equations())
+            .sum::<usize>()
+            + system.constraints().count()
+            + system.coupling_residual_count()
+            + 2 * system.saturated_controller_count()
+            + system.mass_flow_closure_count();
+
+        if n_state == 0 || n_equations == 0 {
+            return Err(SolverError::InvalidSystem {
+                message: "Empty system has no state variables or equations".to_string(),
+            });
+        }
+
+        let mut state: Vec<f64> = match self.initial_state.as_ref() {
+            Some(s) if s.len() == n_state => s.clone(),
+            Some(s) => {
+                return Err(SolverError::InvalidSystem {
+                    message: format!(
+                        "initial_state length {} does not match system state length {}",
+                        s.len(),
+                        n_state
+                    ),
+                });
+            }
+            None => vec![0.0; n_state],
+        };
+        let mut residuals: Vec<f64> = vec![0.0; n_equations];
+        let mut saved_state: Vec<f64> = vec![0.0; n_state];
+        let mut trial_state: Vec<f64> = vec![0.0; n_state];
+        let mut trial_residuals: Vec<f64> = vec![0.0; n_equations];
+        let mut jacobian_builder = JacobianBuilder::new();
+        let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state);
+
+        // Scratch buffers for the dogleg computation — sized to n_state.
+        let mut buf_dp = vec![0.0; n_state];
+        let mut buf_dn = vec![0.0; n_state];
+
+        let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
+            .map(|i| system.get_solver_bounds_for_state_index(i))
+            .collect();
+
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute initial residuals: {:?}", e),
+            })?;
+        let mut current_norm = Self::residual_norm(&residuals);
+
+        tracing::info!(
+            max_iterations = self.max_iterations,
+            tolerance = self.tolerance,
+            delta_init = self.delta_init,
+            residual_norm = current_norm,
+            "Trust-region (Powell dogleg + LM) starting"
+        );
+
+        if current_norm < self.tolerance {
+            return Ok(ConvergedState::new(
+                state,
+                0,
+                current_norm,
+                ConvergenceStatus::Converged,
+                SimulationMetadata::new(system.input_hash()),
+            ));
+        }
+
+        let mut delta = self.delta_init;
+        let mut lm_lambda = self.lm_lambda_init;
+        let mut consecutive_rejections = 0usize;
+
+        for iteration in 1..=self.max_iterations {
+            if let Some(timeout) = self.timeout {
+                if start_time.elapsed() > timeout {
+                    tracing::info!(iteration, "Trust-region timed out");
+                    return Err(SolverError::Timeout {
+                        timeout_ms: timeout.as_millis() as u64,
+                    });
+                }
+            }
+
+            // Fresh analytical Jacobian every iteration — the trust-region
+            // accepts/rejects based on actual residuals, so a stale Jacobian
+            // would mis-predict the gain ratio.
+            jacobian_builder.clear();
+            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());
+
+            // Column scaling diagonal D (constant within an iteration).
+            let d = Self::column_scaling(&jacobian_matrix);
+
+            // Gradient of ½‖r‖²: g = Jᵀ·r.
+            let m = jacobian_matrix.as_matrix().clone();
+            let r_vec = DVector::from_row_slice(&residuals);
+            let grad = m.transpose() * &r_vec;
+
+            // ─── Newton step ──────────────────────────────────────────────
+            let newton_step = jacobian_matrix.solve(&residuals);
+
+            // ─── Cauchy (steepest-descent) step ───────────────────────────
+            // α = (g·g) / (g·JᵀJ·g) = ‖g‖² / ‖J·g‖². For F(x)=0 we use the
+            // convention that the descent direction is −g.
+            let grad_norm_sq = grad.dot(&grad);
+            let j_grad = &m * &grad; // J·g  (n_eq-vector)
+            let j_grad_norm_sq = j_grad.dot(&j_grad);
+            let cauchy_step: Vec<f64> = if j_grad_norm_sq > f64::EPSILON {
+                let alpha = grad_norm_sq / j_grad_norm_sq;
+                grad.iter().map(|gi| -alpha * gi).collect()
+            } else {
+                // Degenerate (gradient ≈ 0): no descent direction available.
+                grad.iter().map(|_| 0.0).collect()
+            };
+
+            // ─── Pick the trust-region step ───────────────────────────────
+            // Try the dogleg first; if the Newton step is missing (singular J),
+            // fall back to a Levenberg–Marquardt solve clipped to the region.
+            let step = if newton_step.is_some() || j_grad_norm_sq > f64::EPSILON {
+                Self::dogleg_step(
+                    newton_step.as_deref(),
+                    &cauchy_step,
+                    &d,
+                    delta,
+                    &mut buf_dp,
+                    &mut buf_dn,
+                )
+            } else {
+                // Both Newton and Cauchy are unusable — LM is the only hope.
+                match Self::lm_step(&jacobian_matrix, &residuals, lm_lambda) {
+                    Some(p) => {
+                        let norm_dp = Self::scaled_norm(&p, &d);
+                        if norm_dp > delta && norm_dp > 0.0 {
+                            p.iter().map(|pi| pi * (delta / norm_dp)).collect()
+                        } else {
+                            p
+                        }
+                    }
+                    None => {
+                        // Even LM failed: shrink δ and retry without consuming
+                        // an iteration of state progress.
+                        delta *= self.shrink_factor;
+                        if delta < self.delta_min {
+                            return Err(SolverError::NonConvergence {
+                                iterations: iteration - 1,
+                                final_residual: current_norm,
+                            });
+                        }
+                        continue;
+                    }
+                }
+            };
+
+            if step.iter().any(|s| !s.is_finite()) {
+                // Poisoned step: shrink and retry.
+                delta *= self.shrink_factor;
+                consecutive_rejections += 1;
+                if delta < self.delta_min || consecutive_rejections > self.max_rejections {
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+                continue;
+            }
+
+            // ─── Trial state ──────────────────────────────────────────────
+            saved_state.copy_from_slice(&state);
+            trial_state.copy_from_slice(&state);
+            apply_newton_step(&mut trial_state, &step, &clipping_mask, 1.0);
+
+            if let Err(e) = system.compute_residuals(&trial_state, &mut trial_residuals) {
+                if !e.is_recoverable() {
+                    // Fatal evaluation error: abort immediately, do not burn
+                    // the shrink budget on it.
+                    return Err(SolverError::InvalidSystem {
+                        message: format!("Failed to compute residuals: {:?}", e),
+                    });
+                }
+                // Recoverable domain violation (KINSOL > 0): reject and shrink.
+                delta *= self.shrink_factor;
+                consecutive_rejections += 1;
+                if delta < self.delta_min || consecutive_rejections > self.max_rejections {
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+                continue;
+            }
+            let trial_norm = Self::residual_norm(&trial_residuals);
+
+            // ─── Gain ratio ρ ─────────────────────────────────────────────
+            //   ρ = (½‖r‖² − ½‖r_trial‖²) / (½‖r‖² − ½‖r + J·p‖²)
+            // The denominator is the linear-model reduction; for square systems
+            // with an exact Newton step it equals ½‖r‖².
+            let actual_reduction = 0.5 * (current_norm * current_norm - trial_norm * trial_norm);
+            let predicted_reduction = {
+                let j_step = &m * DVector::from_row_slice(&step);
+                let predicted: DVector<f64> = &r_vec + &j_step;
+                0.5 * (current_norm * current_norm - predicted.norm().powi(2))
+            };
+            let rho = if predicted_reduction.abs() > f64::EPSILON {
+                actual_reduction / predicted_reduction
+            } else if actual_reduction > 0.0 {
+                // Model and reality both flat-ish; accept progress but do not
+                // expand the region (would lead to runaway δ growth).
+                1.0
+            } else {
+                // No predicted reduction and no actual reduction — reject.
+                -1.0
+            };
+
+            let step_norm_d = Self::scaled_norm(&step, &d);
+
+            if rho >= self.accept_threshold && trial_norm < current_norm {
+                // ─── Accept ───────────────────────────────────────────────
+                state.copy_from_slice(&trial_state);
+                residuals.copy_from_slice(&trial_residuals);
+                current_norm = trial_norm;
+                consecutive_rejections = 0;
+
+                // Grow δ on a boundary step with good model agreement. On a
+                // *very* good gain (ρ > 0.9) we also re-center δ at its init
+                // value when it has collapsed — a strong signal that we just
+                // escaped a tough spot and can afford bolder steps again.
+                if rho > self.expand_threshold
+                    && step_norm_d > 0.9 * delta
+                    && step_norm_d < 1.1 * delta
+                {
+                    delta = (delta * self.expand_factor).min(self.delta_max);
+                } else if rho > self.shrink_threshold {
+                    // Reasonable gain, no shrink. If δ had collapsed and we
+                    // just made real progress, lift it back up so subsequent
+                    // iterations don't get stuck at infinitesimal steps.
+                    if delta < self.delta_init {
+                        delta = self.delta_init.min(self.delta_max);
+                    }
+                } else {
+                    // Accepted but the model poorly predicted the reduction —
+                    // tighten the region for the next iteration.
+                    delta = (delta * self.shrink_factor).max(self.delta_min);
+                }
+
+                // LM λ decay on accepted LM steps.
+                if newton_step.is_none() {
+                    lm_lambda = (lm_lambda * self.lm_lambda_shrink).max(1e-12);
+                }
+
+                tracing::debug!(
+                    iteration,
+                    residual_norm = current_norm,
+                    delta,
+                    rho,
+                    step_norm_d,
+                    "Trust-region accepted step"
+                );
+
+                if current_norm < self.tolerance {
+                    tracing::info!(
+                        iterations = iteration,
+                        final_residual = current_norm,
+                        "Trust-region converged"
+                    );
+                    return Ok(ConvergedState::new(
+                        state,
+                        iteration,
+                        current_norm,
+                        ConvergenceStatus::Converged,
+                        SimulationMetadata::new(system.input_hash()),
+                    ));
+                }
+            } else {
+                // ─── Reject ───────────────────────────────────────────────
+                consecutive_rejections += 1;
+                delta = (delta * self.shrink_factor).max(self.delta_min);
+
+                // Grow LM λ on rejected LM steps so the next iteration is more
+                // gradient-descent-like.
+                if newton_step.is_none() {
+                    lm_lambda = (lm_lambda * self.lm_lambda_grow).min(self.lm_lambda_max);
+                }
+
+                tracing::debug!(
+                    iteration,
+                    residual_norm = current_norm,
+                    trial_norm,
+                    delta,
+                    rho,
+                    consecutive_rejections,
+                    "Trust-region rejected step"
+                );
+
+                if delta <= self.delta_min || consecutive_rejections > self.max_rejections {
+                    tracing::warn!(
+                        iteration,
+                        final_residual = current_norm,
+                        delta,
+                        consecutive_rejections,
+                        "Trust-region failed to make progress"
+                    );
+                    return Err(SolverError::NonConvergence {
+                        iterations: iteration - 1,
+                        final_residual: current_norm,
+                    });
+                }
+            }
+        }
+
+        tracing::warn!(
+            max_iterations = self.max_iterations,
+            final_residual = current_norm,
+            "Trust-region did not converge"
+        );
+        Err(SolverError::NonConvergence {
+            iterations: self.max_iterations,
+            final_residual: current_norm,
+        })
+    }
+
+    fn with_timeout(mut self, timeout: Duration) -> Self {
+        self.timeout = Some(timeout);
+        self
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::solver::Solver;
+    use crate::system::System;
+
+    #[test]
+    fn test_trust_region_defaults() {
+        let cfg = TrustRegionConfig::default();
+        assert_eq!(cfg.max_iterations, 200);
+        assert!(cfg.delta_init > 0.0);
+        assert!(cfg.delta_min > 0.0 && cfg.delta_min < cfg.delta_init);
+        assert!(cfg.delta_max > cfg.delta_init);
+        assert!(cfg.shrink_factor > 0.0 && cfg.shrink_factor < 1.0);
+        assert!(cfg.expand_factor > 1.0);
+        assert!(cfg.accept_threshold > 0.0 && cfg.accept_threshold < 1.0);
+    }
+
+    #[test]
+    fn test_trust_region_rejects_empty_system() {
+        let mut system = System::new();
+        system.finalize().unwrap();
+        let mut solver = TrustRegionConfig::default();
+        let result = solver.solve(&mut system);
+        assert!(matches!(result, Err(SolverError::InvalidSystem { .. })));
+    }
+
+    /// Smoke check on the dogleg step: when Newton is inside the region the
+    /// returned step must equal the Newton step exactly.
+    #[test]
+    fn test_dogleg_returns_newton_inside_region() {
+        // Newton step of norm 0.5 (in D=I), region radius 1.0 → inside.
+        let newton = vec![0.3, 0.4]; // ‖n‖ = 0.5
+        let cauchy = vec![0.1, 0.1];
+        let d = vec![1.0, 1.0];
+        let mut buf_dp = vec![0.0; 2];
+        let mut buf_dn = vec![0.0; 2];
+        let p = TrustRegionConfig::dogleg_step(
+            Some(&newton),
+            &cauchy,
+            &d,
+            1.0,
+            &mut buf_dp,
+            &mut buf_dn,
+        );
+        assert!((p[0] - 0.3).abs() < 1e-12);
+        assert!((p[1] - 0.4).abs() < 1e-12);
+    }
+
+    /// When the Newton step is outside the region but the Cauchy point is
+    /// inside, the dogleg must lie on the segment c → n at the boundary.
+    #[test]
+    fn test_dogleg_blends_to_boundary() {
+        // Newton ‖n‖ = 10, Cauchy ‖c‖ = 0.1, δ = 1.0 → τ must place the step
+        // on the segment c→n with scaled norm 1.
+        let newton = vec![10.0, 0.0];
+        let cauchy = vec![-0.1, 0.0];
+        let d = vec![1.0, 1.0];
+        let mut buf_dp = vec![0.0; 2];
+        let mut buf_dn = vec![0.0; 2];
+        let p = TrustRegionConfig::dogleg_step(
+            Some(&newton),
+            &cauchy,
+            &d,
+            1.0,
+            &mut buf_dp,
+            &mut buf_dn,
+        );
+        let norm = (p[0].powi(2) + p[1].powi(2)).sqrt();
+        // Step must hit the trust-region boundary (within numerical tolerance).
+        assert!(
+            (norm - 1.0).abs() < 1e-9,
+            "dogleg must hit the boundary, got norm = {}",
+            norm
+        );
+        // Step must lie on the segment c → n (x-component between −0.1 and 10).
+        assert!(p[0] >= -0.1 && p[0] <= 10.0);
+        assert!((p[1]).abs() < 1e-12);
+    }
+
+    /// When Newton is missing (singular Jacobian) and Cauchy is outside the
+    /// region, the step is the Cauchy direction clipped to the boundary.
+    #[test]
+    fn test_dogleg_clips_cauchy_when_newton_missing() {
+        let cauchy = vec![3.0, 4.0]; // ‖c‖ = 5
+        let d = vec![1.0, 1.0];
+        let mut buf_dp = vec![0.0; 2];
+        let mut buf_dn = vec![0.0; 2];
+        let p = TrustRegionConfig::dogleg_step(None, &cauchy, &d, 1.0, &mut buf_dp, &mut buf_dn);
+        let norm = (p[0].powi(2) + p[1].powi(2)).sqrt();
+        assert!((norm - 1.0).abs() < 1e-12);
+        // Same direction as Cauchy.
+        assert!(p[0] > 0.0 && p[1] > 0.0);
+    }
+
+    /// Column scaling must assign ~1 to near-zero columns and 1/‖col‖ otherwise.
+    #[test]
+    fn test_column_scaling_handles_zero_and_unit_columns() {
+        // J = [[10, 0], [0, 2]] → d = [1/10, 1/2].
+        let entries = vec![(0, 0, 10.0), (1, 1, 2.0)];
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+        let d = TrustRegionConfig::column_scaling(&j);
+        assert!((d[0] - 0.1).abs() < 1e-12);
+        assert!((d[1] - 0.5).abs() < 1e-12);
+
+        // Zero column → d = 1.
+        let entries = vec![(0, 0, 5.0), (1, 0, 5.0)]; // col 1 all zeros
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+        let d = TrustRegionConfig::column_scaling(&j);
+        assert!((d[1] - 1.0).abs() < 1e-12);
+    }
+}

crates/solver/src/system.rs

@@ -772,11 +772,7 @@ impl System {
                 tracing::debug!("{}", report.summary());
             }
             SystemDofBalance::UnderConstrained { free_dofs } => {
-                tracing::warn!(
-                    free_dofs,
-                    "{}",
-                    report.summary()
-                );
+                tracing::warn!(free_dofs, "{}", report.summary());
             }
             SystemDofBalance::OverConstrained { excess_equations } => {
                 tracing::error!(excess_equations, "{}", report.summary());
@@ -2339,9 +2335,7 @@ impl System {
             }
         }
         for id in self.inverse_control.linked_controls() {
-            unknowns.push(UnknownKind::InverseControl {
-                id: id.to_string(),
-            });
+            unknowns.push(UnknownKind::InverseControl { id: id.to_string() });
         }
         for index in 0..self.coupling_residual_count() {
             unknowns.push(UnknownKind::CouplingHeat { index });
@@ -2351,9 +2345,7 @@ impl System {
             unknowns.push(UnknownKind::SaturatedIntegrator { index });
         }
         for id in &self.free_actuators {
-            unknowns.push(UnknownKind::FreeActuator {
-                id: id.to_string(),
-            });
+            unknowns.push(UnknownKind::FreeActuator { id: id.to_string() });
         }
 
         if unknowns.len() != n_unknowns {