//! System-level regression test for the MSH tube-ΔP + fixed-opening EXV
//! solver-robustness fix (Epic-0 follow-up).
//!
//! Builds the exact user system that exhibited the Newton stall (4-component
//! emergent-pressure R134a chiller, `dp_model=msh` on both heat exchangers,
//! `fix_opening=true, opening=0.9`) at the exact CLI staged seed, and guards:
//!
//! 1. the cold-start residual signature (evaporator tube-ΔP row dominant),
//! 2. the **momentum-row Jacobians** (condenser + evaporator tube ΔP) against
//!    central finite differences — the NFR9 guard for the exact analytic
//!    `tube_dp` composition wired into both heat exchangers,
//! 3. the **totality / C¹** of the tube-ΔP residual when Newton iterates leave
//!    the saturation domain (no silent model switch, no error, smooth values).
//!
//! Run: cargo test -p entropyk-solver --features coolprop --test msh_tube_dp_robustness
#![cfg(feature = "coolprop")]

use std::sync::Arc;

use entropyk_components::heat_exchanger::two_phase_dp::{
    TubeChannelGeometry, TwoPhaseDpCorrelation,
};
use entropyk_components::{BrineSink, BrineSource, Condenser, Evaporator};
use entropyk_components::{ConnectedPort, FluidId as ComponentFluidId};
use entropyk_components::{IsenthalpicExpansionValve, IsentropicCompressor, JacobianBuilder, Port};
use entropyk_core::{Concentration, Enthalpy, Pressure, Temperature};
use entropyk_fluids::{CoolPropBackend, FluidBackend, FluidId, FluidState, Property};
use entropyk_solver::scaling::{equilibrate, unscale_dx};
use entropyk_solver::system::System;

fn water_h(backend: &Arc<dyn FluidBackend>, t_c: f64) -> f64 {
    backend
        .property(
            FluidId::new("Water"),
            Property::Enthalpy,
            FluidState::from_pt(Pressure::from_bar(2.0), Temperature::from_celsius(t_c)),
        )
        .expect("water h(P,T)")
}

fn water_port() -> ConnectedPort {
    let fluid = ComponentFluidId::new("Water");
    let a = Port::new(
        fluid.clone(),
        Pressure::from_bar(2.0),
        Enthalpy::from_joules_per_kg(100_000.0),
    );
    let b = Port::new(
        fluid,
        Pressure::from_bar(2.0),
        Enthalpy::from_joules_per_kg(100_000.0),
    );
    a.connect(b).expect("port connect").0
}

/// Builds the exact hang system (msh + fixed-opening EXV at `opening`) and
/// returns it with the exact CLI staged seed and per-index variable names.
fn build_hang_system_with_opening(opening: f64) -> (System, Vec<f64>, Vec<String>) {
    let backend: Arc<dyn FluidBackend> = Arc::new(CoolPropBackend::new());
    let fluid = "R134a";

    let geom = TubeChannelGeometry {
        length_m: 6.0,
        diameter_m: 0.0095,
        n_parallel: 2.0,
    };

    // comp: emergent metered-flow (energy-only) — EXV fixed orifice meters ṁ.
    let comp = Box::new(
        IsentropicCompressor::new(0.70, 318.15, 278.15, 5.0)
            .with_refrigerant(fluid)
            .with_fluid_backend(backend.clone())
            .with_emergent_metered_flow(),
    );

    // cond: UA=1500, msh tube ΔP, water 4-port, emergent with 5 K subcooling.
    let mut cond = Condenser::new(1500.0)
        .with_refrigerant(fluid)
        .with_fluid_backend(backend.clone())
        .with_emergent_pressure(5.0);
    cond.set_secondary_fluid("Water");
    cond.set_tube_pressure_drop(TwoPhaseDpCorrelation::MullerSteinhagenHeck1986, geom);
    cond.set_secondary_pressure_drop_coeff(5000.0);
    let cond = Box::new(cond);

    // exv: emergent + fixed orifice kv=2e-6.
    let exv = Box::new(
        IsenthalpicExpansionValve::new(278.15)
            .with_refrigerant(fluid)
            .with_fluid_backend(backend.clone())
            .with_emergent_pressure()
            .with_orifice_fixed(2e-6, opening),
    );

    // evap: UA=2500, msh tube ΔP, water 4-port, emergent (5 K superheat).
    let mut evap = Evaporator::new(2500.0)
        .with_refrigerant(fluid)
        .with_fluid_backend(backend.clone())
        .with_emergent_pressure();
    evap.set_secondary_fluid("Water");
    evap.set_tube_pressure_drop(TwoPhaseDpCorrelation::MullerSteinhagenHeck1986, geom);
    evap.set_secondary_pressure_drop_coeff(5000.0);
    let evap = Box::new(evap);

    // Water boundaries.
    let cwin = Box::new(
        BrineSource::new(
            "Water",
            Pressure::from_bar(2.0),
            Temperature::from_celsius(30.0),
            Concentration::from_percent(0.0),
            backend.clone(),
            water_port(),
        )
        .expect("BrineSource cond")
        .with_imposed_mass_flow(0.3583)
        .expect("imposed m"),
    );
    let cwout = Box::new(
        BrineSink::new(
            "Water",
            Pressure::from_bar(2.0),
            None,
            None,
            backend.clone(),
            water_port(),
        )
        .expect("BrineSink cond"),
    );
    let ewin = Box::new(
        BrineSource::new(
            "Water",
            Pressure::from_bar(2.0),
            Temperature::from_celsius(12.0),
            Concentration::from_percent(0.0),
            backend.clone(),
            water_port(),
        )
        .expect("BrineSource evap")
        .with_imposed_mass_flow(0.4778)
        .expect("imposed m"),
    );
    let ewout = Box::new(
        BrineSink::new(
            "Water",
            Pressure::from_bar(2.0),
            None,
            None,
            backend.clone(),
            water_port(),
        )
        .expect("BrineSink evap"),
    );

    let mut system = System::new();
    let n_comp = system.add_component(comp);
    let n_cond = system.add_component(cond);
    let n_exv = system.add_component(exv);
    let n_evap = system.add_component(evap);
    let n_cwin = system.add_component(cwin);
    let n_cwout = system.add_component(cwout);
    let n_ewin = system.add_component(ewin);
    let n_ewout = system.add_component(ewout);

    // Refrigerant loop (ports: inlet=0, outlet=1; HX secondary 2/3).
    system.add_edge_with_ports(n_comp, 1, n_cond, 0).unwrap(); // E0
    system.add_edge_with_ports(n_cond, 1, n_exv, 0).unwrap(); // E1
    system.add_edge_with_ports(n_exv, 1, n_evap, 0).unwrap(); // E2
    system.add_edge_with_ports(n_evap, 1, n_comp, 0).unwrap(); // E3
                                                               // Water loops.
    system.add_edge_with_ports(n_cwin, 1, n_cond, 2).unwrap(); // W0
    system.add_edge_with_ports(n_cond, 3, n_cwout, 0).unwrap(); // W1
    system.add_edge_with_ports(n_ewin, 1, n_evap, 2).unwrap(); // W2
    system.add_edge_with_ports(n_evap, 3, n_ewout, 0).unwrap(); // W3

    system.finalize().unwrap();

    let n_state = system.full_state_vector_len();
    assert_eq!(n_state, 19, "hang system must have 19 unknowns");

    // ── Exact CLI staged seed (validated against the production run: the
    // cold-start residual breakdown matches row-by-row) ─────────────────────
    let h_w30 = water_h(&backend, 30.0);
    let h_w12 = water_h(&backend, 12.0);
    let h_w20 = water_h(&backend, 20.0);
    let seed_refrig = [
        (1.159_924e6, 4.465_191e5), // E0 comp→cond
        (1.159_924e6, 2.589_429e5), // E1 cond→exv
        (3.496_586e5, 2.639_429e5), // E2 exv→evap
        (3.496_586e5, 4.060_707e5), // E3 evap→comp
    ];
    // Water loop ṁ slots are seeded at the generic 0.05 kg/s default (the
    // sink boundary seed overwrites the source's imposed flow).
    let seed_water = [
        (0.05, h_w30), // W0 source→cond
        (0.05, h_w20), // W1 cond→sink
        (0.05, h_w12), // W2 source→evap
        (0.05, h_w20), // W3 evap→sink
    ];

    let mut state = vec![0.0; n_state];
    let mut names: Vec<Option<String>> = vec![None; n_state];
    let edge_names = ["E0", "E1", "E2", "E3", "W0", "W1", "W2", "W3"];
    for (i, e) in system.edge_indices().enumerate() {
        let (mi, pi, hi) = system.edge_state_indices_full(e);
        if i < 4 {
            state[mi] = 0.05;
            state[pi] = seed_refrig[i].0;
            state[hi] = seed_refrig[i].1;
        } else {
            let (mw, hw) = seed_water[i - 4];
            state[mi] = mw;
            state[pi] = 2.0e5;
            state[hi] = hw;
        }
        let en = edge_names[i];
        for (idx, tag) in [(mi, "m"), (pi, "P"), (hi, "h")] {
            if names[idx].is_none() {
                names[idx] = Some(format!("{tag}({en})"));
            }
        }
    }
    let names: Vec<String> = names
        .into_iter()
        .enumerate()
        .map(|(i, n)| n.unwrap_or_else(|| format!("x[{i}]")))
        .collect();
    (system, state, names)
}

/// Row indices: node order comp(0) | cond(1..5) | exv(6..7) | evap(8..12) |
/// boundaries(13..18). Row 1 = condenser refrigerant momentum,
/// row 8 = evaporator refrigerant momentum.
const COND_MOMENTUM_ROW: usize = 1;
const EVAP_MOMENTUM_ROW: usize = 8;
const N_EQ: usize = 19;

/// The exact cold-start signature of the production stall (residual breakdown
/// matches the CLI run row-by-row): the evaporator tube-ΔP momentum row
/// dominates the cold residual.
#[test]
fn cold_start_residual_signature_matches_production() {
    let (system, state, _names) = build_hang_system_with_opening(0.9);
    let mut r = vec![0.0; N_EQ];
    system.compute_residuals(&state, &mut r).unwrap();
    let norm: f64 = r.iter().map(|v| v * v).sum::<f64>().sqrt();
    assert!(
        (norm - 44457.554).abs() < 0.01,
        "cold-start residual norm must match the production signature, got {norm}"
    );
    // Evaporator momentum row: the full tube ΔP is unbalanced at the seed
    // (uniform low-side pressure), ≈ 42 kPa.
    assert!(
        (r[EVAP_MOMENTUM_ROW] - 41_991.24).abs() < 1.0,
        "evap momentum residual: {}",
        r[EVAP_MOMENTUM_ROW]
    );
    // Condenser momentum row ≈ 7.7 kPa.
    assert!(
        (r[COND_MOMENTUM_ROW] - 7_673.32).abs() < 1.0,
        "cond momentum residual: {}",
        r[COND_MOMENTUM_ROW]
    );
}

/// NFR9 guard: the momentum-row Jacobian of the tube ΔP (now fully analytic
/// through `heat_exchanger::tube_dp`) must agree with central finite
/// differences of the residual at the cold seed, on every column.
#[test]
fn tube_dp_momentum_jacobian_matches_fd_at_cold_seed() {
    let (system, state, names) = build_hang_system_with_opening(0.9);
    let n_state = state.len();

    let mut jb = JacobianBuilder::new();
    system.assemble_jacobian(&state, &mut jb).unwrap();
    let mut analytic = vec![vec![0.0_f64; n_state]; N_EQ];
    for &(row, col, v) in jb.entries() {
        analytic[row][col] += v;
    }

    for row in [COND_MOMENTUM_ROW, EVAP_MOMENTUM_ROW] {
        for col in 0..n_state {
            let eps = (state[col].abs() * 1e-6).max(1e-7);
            let (mut sp, mut sm) = (state.clone(), state.clone());
            sp[col] += eps;
            sm[col] -= eps;
            let (mut rp, mut rm) = (vec![0.0; N_EQ], vec![0.0; N_EQ]);
            system.compute_residuals(&sp, &mut rp).unwrap();
            system.compute_residuals(&sm, &mut rm).unwrap();
            let fd = (rp[row] - rm[row]) / (2.0 * eps);
            let a = analytic[row][col];
            if a == 0.0 && fd == 0.0 {
                continue;
            }
            let tol = (1e-4 * fd.abs().max(a.abs())).max(1e-9);
            assert!(
                (a - fd).abs() <= tol,
                "momentum J[{row}][{}]: analytic={a} vs fd={fd}",
                names[col]
            );
        }
    }
}

/// Totality + C¹ guard: pushing the refrigerant pressures far outside the
/// saturation domain (as Newton iterates do when a step overshoots) must keep
/// the residual defined — no error, no panic, no silent ΔP-model switch — and
/// the evaporator momentum row must vary smoothly across the domain bound.
#[test]
fn tube_dp_residual_is_total_and_smooth_outside_sat_domain() {
    let (system, state, _names) = build_hang_system_with_opening(0.9);
    // Find the E2 (exv→evap) pressure index.
    let e2 = system.edge_indices().nth(2).unwrap();
    let (_, p_e2, _) = system.edge_state_indices_full(e2);

    let r8_at = |p: f64| -> f64 {
        let mut s = state.clone();
        s[p_e2] = p;
        let mut r = vec![0.0; N_EQ];
        system
            .compute_residuals(&s, &mut r)
            .expect("residual must stay defined outside the saturation domain");
        r[EVAP_MOMENTUM_ROW]
    };

    // Deep outside the R134a saturation domain in both directions: the tube
    // ΔP saturates to the bound's value (constant continuation), and the
    // residual stays finite and equal across the far exterior.
    let p_nominal = state[p_e2];
    let r_nominal = r8_at(p_nominal);
    assert!(r_nominal.is_finite());
    for p_extreme in [1.0, 10.0, 1.0e9, 1.0e12] {
        let r = r8_at(p_extreme);
        assert!(r.is_finite(), "residual not finite at P={p_extreme}");
    }
    // Constant continuation far outside: two far-exterior points give the
    // same clamped ΔP (isolate it from the row: P_out − P_in + ΔP_sat).
    let d1 = r8_at(1.0e9) + 1.0e9;
    let d2 = r8_at(1.0e10) + 1.0e10;
    assert!(
        (d1 - d2).abs() < 1e-6 * d1.abs().max(1.0),
        "clamped ΔP must be constant far outside the domain: {d1} vs {d2}"
    );
    // C¹ across the upper domain bound: FD slope just inside vs just outside
    // the saturation-domain top must not jump (smooth clamp, Story 0.2).
    let (p_min, p_max) = {
        let backend: Arc<dyn FluidBackend> = Arc::new(CoolPropBackend::new());
        entropyk_components::heat_exchanger::sat_domain::saturation_pressure_domain(
            &backend, "R134a",
        )
        .expect("R134a domain")
    };
    let _ = p_min;
    let h = p_max * 1e-5;
    let slope_in = (r8_at(p_max - h) - r8_at(p_max - 2.0 * h)) / h;
    let slope_out = (r8_at(p_max + 2.0 * h) - r8_at(p_max + h)) / h;
    let denom = slope_in.abs().max(1.0); // the explicit −P_in term dominates
    assert!(
        (slope_in - slope_out).abs() / denom < 0.2,
        "C¹ slope across domain bound: in={slope_in} out={slope_out}"
    );
}

/// Newton-step structure at the cold seed (documentation of the stall
/// mechanism): the first full Newton step must be finite and the scaled
/// Jacobian must be non-singular — the production stall is a *damping/seed
/// distance* problem, not a singular or misassembled Jacobian.
#[test]
fn cold_seed_jacobian_is_nonsingular_and_step_finite() {
    let (system, state, _names) = build_hang_system_with_opening(0.9);
    let n_state = state.len();
    let mut r = vec![0.0; N_EQ];
    system.compute_residuals(&state, &mut r).unwrap();

    let mut jb = JacobianBuilder::new();
    system.assemble_jacobian(&state, &mut jb).unwrap();
    let mut jm = nalgebra::DMatrix::<f64>::zeros(N_EQ, n_state);
    for &(row, col, v) in jb.entries() {
        jm[(row, col)] += v;
    }
    let (d_r, d_c) = equilibrate(&jm);
    let mut js = jm.clone();
    for i in 0..N_EQ {
        for j in 0..n_state {
            js[(i, j)] *= d_r[i] * d_c[j];
        }
    }
    let b: nalgebra::DVector<f64> =
        nalgebra::DVector::from_iterator(N_EQ, (0..N_EQ).map(|i| -d_r[i] * r[i]));
    let y = js
        .clone()
        .lu()
        .solve(&b)
        .expect("scaled Jacobian must solve");
    let delta = unscale_dx(y.as_slice(), &d_c);
    assert!(
        delta.iter().all(|v| v.is_finite()),
        "Newton step must be finite"
    );
    // The step is huge (stiff emergent-pressure mode) but the scaled Jacobian
    // is invertible — σ_min > 0.
    let sigma_min = js
        .svd(false, false)
        .singular_values
        .iter()
        .copied()
        .fold(f64::INFINITY, f64::min);
    assert!(sigma_min > 0.0, "scaled Jacobian must be non-singular");
}