Entropyk/crates/solver/tests/recoverable_domain_violation.rs

//! Integration tests for Story 1.3: Recoverable Domain-Violation Errors.
//!
//! Covers:
//! - AC #2: a Newton trial step that triggers a recoverable domain violation
//!   is shrunk and retried by the line search; after exhaustion of backtracks
//!   the solve fails with `ConvergenceReason::DomainViolation` and diagnostics.
//! - AC #3: end-to-end mock proof — a component that throws a domain violation
//!   on oversized steps converges via step reduction, and the iteration
//!   history records the recoverable events.
//! - Fatal vs recoverable: a fatal `ComponentError` is never retried by the
//!   line search, and a recoverable violation outside the line search is a
//!   typed `SolverError::DomainViolation`, not an `InvalidSystem` flattening.
//!
//! Mock pattern follows `convergence_reason.rs` (`LinearSystem` /
//! `create_test_system` self-loop fixture): state layout `[ṁ, P, h]`, the
//! mass-flow pin row parks ṁ at `DEFAULT_MASS_FLOW_SEED_KG_S`, and the
//! enthalpy slot is parked at `MIN_SOLVER_PRESSURE_PA` (abstract unknown).

use std::sync::{
    atomic::{AtomicUsize, Ordering},
    Arc,
};

use entropyk_components::{
    Component, ComponentError, DomainViolation, JacobianBuilder, ResidualVector, StateSlice,
};
use entropyk_solver::solver::{NewtonConfig, Solver, SolverError, VerboseConfig};
use entropyk_solver::system::{System, DEFAULT_MASS_FLOW_SEED_KG_S, MIN_SOLVER_PRESSURE_PA};
use entropyk_solver::{ConvergenceReason, DomainViolation as SolverDomainViolation, SolveOutcome};

// ─────────────────────────────────────────────────────────────────────────────
// Mock component: DomainWallComponent
// ─────────────────────────────────────────────────────────────────────────────

/// Pressure-row residual shape.
///
/// A strictly linear residual can never overshoot with an exact Jacobian
/// (Newton solves an affine system exactly in one step), so AC#3's
/// "full step overshoots the wall" requirement forces a nonlinear row.
/// Both shapes below drive `P` monotonically toward the target with an
/// exactly computable Newton step `ΔP = -(P - T)/p`:
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum WallShape {
    /// `r = sign(P − T)·√|P − T|` (p = 1/2): the full Newton step is
    /// `ΔP = -2(P − T)`, i.e. it lands at `2T − P₀`, symmetrically past the
    /// target. With `T − P₀ = 256` (dyadic) the halved step lands on the
    /// target exactly — zero floating-point noise in the AC#3 assertion.
    Sqrt,
    /// `r = cbrt(P − T)` (p = 1/3): the full step is `ΔP = -3(P − T)`
    /// (trial at `3T − 2P₀`, Armijo-rejected), while the halved step shrinks
    /// `|P − T|` by exactly 1/2 per iteration — a slow, deterministic
    /// approach used to stage the AC#2 exhaustion scenario over several
    /// iterations before the wall activates.
    Cbrt,
}

/// What the wall returns when a trial pressure exceeds it.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum WallMode {
    /// Recoverable domain violation (KINSOL `> 0`): the line search must
    /// shrink the step and retry.
    Recoverable,
    /// Fatal calculation failure (KINSOL `< 0`): the line search must abort
    /// immediately without burning backtracks.
    Fatal,
}

/// A mock component whose pressure row drives `P` toward `target_p` and whose
/// "physical domain" ends at `wall_p`: any evaluation at `P > wall_p` fails.
///
/// Residuals (state layout `[ṁ, P, h]`, 3 equations):
/// - `r0 = shape(P − target_p)` — pressure row,
/// - `r1 = h − parked_h` — enthalpy pin,
/// - `r2 = ṁ − DEFAULT_MASS_FLOW_SEED_KG_S` — mass-flow pin (CM1.3).
///
/// The Jacobian is analytic and exact (project rule — no finite differences):
/// `∂r0/∂P = 1/(p·|P − T|^(1−p))`, `∂r1/∂h = 1`, `∂r2/∂ṁ = 1`.
///
/// `wall_after_evals` delays wall enforcement until the evaluation count
/// exceeds the given value, which stages the AC#2 scenario: the wall only
/// becomes active after a few iterations have been recorded in the
/// diagnostics history, so the exhausting iteration carries diagnostics.
struct DomainWallComponent {
    target_p: f64,
    wall_p: f64,
    parked_h: f64,
    shape: WallShape,
    mode: WallMode,
    wall_after_evals: usize,
    eval_count: Arc<AtomicUsize>,
}

impl DomainWallComponent {
    fn new(
        target_p: f64,
        wall_p: f64,
        shape: WallShape,
        mode: WallMode,
        wall_after_evals: usize,
        eval_count: Arc<AtomicUsize>,
    ) -> Self {
        Self {
            target_p,
            wall_p,
            parked_h: MIN_SOLVER_PRESSURE_PA,
            shape,
            mode,
            wall_after_evals,
            eval_count,
        }
    }
}

impl Component for DomainWallComponent {
    fn compute_residuals(
        &self,
        state: &StateSlice,
        residuals: &mut ResidualVector,
    ) -> Result<(), ComponentError> {
        let eval = self.eval_count.fetch_add(1, Ordering::SeqCst) + 1;
        let p = state[1];
        if eval > self.wall_after_evals && p > self.wall_p {
            let detail = format!(
                "trial pressure {p:.1} Pa exceeds domain wall at {:.1} Pa",
                self.wall_p
            );
            return Err(match self.mode {
                WallMode::Recoverable => ComponentError::DomainViolation(DomainViolation {
                    component: Some("domain-wall".into()),
                    detail,
                }),
                WallMode::Fatal => ComponentError::CalculationFailed(detail),
            });
        }
        let dp = p - self.target_p;
        residuals[0] = match self.shape {
            WallShape::Sqrt => dp.signum() * dp.abs().sqrt(),
            WallShape::Cbrt => dp.cbrt(),
        };
        residuals[1] = state[2] - self.parked_h;
        residuals[2] = state[0] - DEFAULT_MASS_FLOW_SEED_KG_S;
        Ok(())
    }

    fn jacobian_entries(
        &self,
        state: &StateSlice,
        jacobian: &mut JacobianBuilder,
    ) -> Result<(), ComponentError> {
        let dp_abs = (state[1] - self.target_p).abs();
        let dr_dp = match self.shape {
            // d/dP [sign(P−T)·√|P−T|] = 1/(2·√|P−T|)
            WallShape::Sqrt => 0.5 / dp_abs.sqrt(),
            // d/dP cbrt(P−T) = 1/(3·|P−T|^(2/3))
            WallShape::Cbrt => 1.0 / (3.0 * dp_abs.cbrt().powi(2)),
        };
        jacobian.add_entry(0, 1, dr_dp);
        jacobian.add_entry(1, 2, 1.0);
        jacobian.add_entry(2, 0, 1.0);
        Ok(())
    }

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

    fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
        &[]
    }
}

/// Creates a minimal self-loop system with a single mock component.
fn create_test_system(component: Box<dyn Component>) -> System {
    let mut system = System::new();
    let n0 = system.add_component(component);
    system.add_edge(n0, n0).unwrap();
    system.finalize().unwrap();
    system
}

/// Newton with line search and verbose diagnostics (iteration history is only
/// collected when verbose mode is active).
fn line_search_solver(initial_p: f64) -> NewtonConfig {
    NewtonConfig {
        line_search: true,
        initial_state: Some(vec![
            DEFAULT_MASS_FLOW_SEED_KG_S,
            initial_p,
            MIN_SOLVER_PRESSURE_PA,
        ]),
        verbose_config: VerboseConfig {
            enabled: true,
            log_residuals: true,
            ..VerboseConfig::default()
        },
        ..NewtonConfig::default()
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// AC #3: converges via step reduction, events recorded in iteration history
// ─────────────────────────────────────────────────────────────────────────────
//
// Geometry (sqrt shape, dyadic exact): target T = floor + 300, start
// P₀ = T − 256, wall = T + 100. The full Newton step ΔP = +512 proposes the
// trial 2T − P₀ = T + 256, which crosses the wall (T + 100) → recoverable
// violation. The first backtrack (alpha = 0.5) lands exactly on T, which sits
// below the wall, and is accepted: the solve converges in one iteration with
// `recoverable_events = 1` and `alpha = 0.5` recorded.
#[test]
fn line_search_converges_via_step_reduction_on_domain_violation() {
    let target_p = MIN_SOLVER_PRESSURE_PA + 300.0;
    let start_p = target_p - 256.0;
    let wall_p = target_p + 100.0;
    let eval_count = Arc::new(AtomicUsize::new(0));

    let mut system = create_test_system(Box::new(DomainWallComponent::new(
        target_p,
        wall_p,
        WallShape::Sqrt,
        WallMode::Recoverable,
        0,
        Arc::clone(&eval_count),
    )));
    let mut solver = line_search_solver(start_p);

    let converged = solver
        .solve(&mut system)
        .expect("line search must shrink the oversized step and converge");

    assert!(
        converged.is_converged(),
        "solve must report full convergence, got {:?}",
        converged.status
    );
    assert!(
        (converged.state[1] - target_p).abs() <= 1e-9 * target_p,
        "converged pressure {} must equal the target {}",
        converged.state[1],
        target_p
    );

    let diagnostics = converged
        .diagnostics
        .as_ref()
        .expect("verbose solve must attach iteration diagnostics");
    let event_iteration = diagnostics
        .iteration_history
        .iter()
        .find(|iteration| iteration.recoverable_events > 0)
        .expect("iteration history must record the recoverable domain violation");
    match event_iteration.alpha {
        Some(alpha) => assert!(
            alpha < 1.0,
            "the recovering iteration must use a reduced step, got alpha = {alpha}"
        ),
        None => panic!("the recovering iteration must record its line-search alpha"),
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// AC #2: backtrack exhaustion surfaces a typed DomainViolation outcome
// ─────────────────────────────────────────────────────────────────────────────
//
// Geometry (cbrt shape): target T = floor + 200, start P₀ = T − 100, wall at
// T − 50 (below the target: the solution is unreachable without crossing the
// wall). The wall activates only after 10 residual evaluations so three
// iterations are recorded first:
//   - initial eval (1); iter 1: alpha = 1 Armijo-rejected (2), alpha = 0.5
//     accepted (3), post-step eval (4) → P₁ = T + 50;
//   - iter 2: evals 5, 6, post-step 7 → P₂ = T − 25;
//   - iter 3: evals 8, 9, post-step 10 → P₃ = T + 12.5.
// Iteration 4 starts above the target with a downward Newton step of −37.5 Pa;
// every backtracked trial stays in [T − 25, T + 12.5], all above the wall, so
// all 20 trials fail recoverably and the line search exhausts.
#[test]
fn line_search_exhaustion_returns_typed_domain_violation_outcome() {
    let target_p = MIN_SOLVER_PRESSURE_PA + 200.0;
    let start_p = target_p - 100.0;
    let wall_p = target_p - 50.0;
    let build_system = |eval_count: Arc<AtomicUsize>| {
        create_test_system(Box::new(DomainWallComponent::new(
            target_p,
            wall_p,
            WallShape::Cbrt,
            WallMode::Recoverable,
            10,
            eval_count,
        )))
    };

    // Legacy `solve`: hard Err, but typed as DomainViolation with diagnostics.
    let mut system = build_system(Arc::new(AtomicUsize::new(0)));
    let mut solver = line_search_solver(start_p);
    let err = solver
        .solve(&mut system)
        .expect_err("an unreachable target must fail the solve");

    let violation = match err.base_error() {
        SolverError::DomainViolation(violation) => violation,
        other => panic!("expected SolverError::DomainViolation, got {other:?}"),
    };
    // The solver facade re-exports the shared solver-core type (Task 3).
    let _: &SolverDomainViolation = violation;
    assert!(
        violation.component.is_none(),
        "the exhaustion error is raised by the line search, not a component: {violation:?}"
    );
    assert!(
        violation.detail.contains("exhausted"),
        "exhaustion detail must say the line search exhausted, got: {}",
        violation.detail
    );
    assert_eq!(
        err.convergence_reason(),
        Some(ConvergenceReason::DomainViolation),
        "typed domain violation must classify as ConvergenceReason::DomainViolation"
    );
    let diagnostics = err
        .diagnostics()
        .expect("failure after recorded iterations must carry diagnostics");
    assert!(
        !diagnostics.iteration_history.is_empty(),
        "iterations completed before exhaustion must be in the history"
    );

    // `solve_outcome`: the same termination is outcome data, not a hard Err.
    let mut system = build_system(Arc::new(AtomicUsize::new(0)));
    let mut solver = line_search_solver(start_p);
    let outcome: SolveOutcome = solver
        .solve_outcome(&mut system)
        .expect("domain-violation exhaustion must be an outcome, not a hard error");
    assert_eq!(outcome.reason, ConvergenceReason::DomainViolation);
    assert!(!outcome.is_converged());
}

// ─────────────────────────────────────────────────────────────────────────────
// Fatal vs recoverable: fatal errors are never retried by the line search
// ─────────────────────────────────────────────────────────────────────────────
//
// Same geometry as AC#3 (the full step crosses the wall), but the wall returns
// `ComponentError::CalculationFailed` (fatal, KINSOL `< 0`). The line search
// must abort immediately: exactly 2 residual evaluations (initial + the fatal
// trial), no backtracking retries, and an `InvalidSystem`-shaped hard error
// whose `convergence_reason()` is `None`.
#[test]
fn fatal_error_aborts_line_search_without_retries() {
    let target_p = MIN_SOLVER_PRESSURE_PA + 300.0;
    let start_p = target_p - 256.0;
    let wall_p = target_p + 100.0;
    let eval_count = Arc::new(AtomicUsize::new(0));

    let mut system = create_test_system(Box::new(DomainWallComponent::new(
        target_p,
        wall_p,
        WallShape::Sqrt,
        WallMode::Fatal,
        0,
        Arc::clone(&eval_count),
    )));
    let mut solver = line_search_solver(start_p);

    let err = solver
        .solve(&mut system)
        .expect_err("a fatal component error must abort the solve");

    assert!(
        matches!(err.base_error(), SolverError::InvalidSystem { .. }),
        "fatal component errors keep the InvalidSystem flattening, got {:?}",
        err.base_error()
    );
    assert_eq!(
        err.convergence_reason(),
        None,
        "fatal errors are hard errors, not convergence outcomes"
    );
    assert_eq!(
        eval_count.load(Ordering::SeqCst),
        2,
        "fatal errors must not burn line-search backtracks (initial eval + fatal trial)"
    );
}

// ─────────────────────────────────────────────────────────────────────────────
// Recoverable violation without line search: typed error, no flattening
// ─────────────────────────────────────────────────────────────────────────────
//
// With `line_search: false` the full Newton step is applied unconditionally;
// the post-step residual evaluation then sits beyond the wall and fails
// recoverably. The uniform residual-eval rule maps this to a typed
// `SolverError::DomainViolation` (an outcome), not an `InvalidSystem`
// flattening (a hard error).
#[test]
fn recoverable_violation_without_line_search_is_typed_not_flattened() {
    let target_p = MIN_SOLVER_PRESSURE_PA + 300.0;
    let start_p = target_p - 256.0;
    let wall_p = target_p + 100.0;
    let eval_count = Arc::new(AtomicUsize::new(0));

    let mut system = create_test_system(Box::new(DomainWallComponent::new(
        target_p,
        wall_p,
        WallShape::Sqrt,
        WallMode::Recoverable,
        0,
        Arc::clone(&eval_count),
    )));
    let mut solver = NewtonConfig {
        line_search: false,
        initial_state: Some(vec![
            DEFAULT_MASS_FLOW_SEED_KG_S,
            start_p,
            MIN_SOLVER_PRESSURE_PA,
        ]),
        ..NewtonConfig::default()
    };

    let err = solver
        .solve(&mut system)
        .expect_err("the post-step recoverable violation must fail the solve");

    assert!(
        matches!(err.base_error(), SolverError::DomainViolation(_)),
        "post-step recoverable violation must be typed DomainViolation, got {:?}",
        err.base_error()
    );
    assert_eq!(
        err.convergence_reason(),
        Some(ConvergenceReason::DomainViolation),
        "typed domain violation must classify as ConvergenceReason::DomainViolation"
    );
    assert_eq!(
        eval_count.load(Ordering::SeqCst),
        2,
        "initial eval + post-step eval, no line-search trials"
    );
}