Entropyk/crates/components/src/pump.rs

//! Pump Component Implementation
//!
//! This module provides a pump component for hydraulic systems using
//! polynomial performance curves and affinity laws for variable speed operation.
//!
//! ## Performance Curves
//!
//! **Head Curve:** H = a₀ + a₁Q + a₂Q² + a₃Q³
//!
//! **Efficiency Curve:** η = b₀ + b₁Q + b₂Q²
//!
//! **Hydraulic Power:** P_hydraulic = ρ × g × Q × H / η
//!
//! ## Affinity Laws (Variable Speed)
//!
//! When operating at reduced speed (VFD):
//! - Q₂/Q₁ = N₂/N₁
//! - H₂/H₁ = (N₂/N₁)²
//! - P₂/P₁ = (N₂/N₁)³

use crate::polynomials::{AffinityLaws, PerformanceCurves, Polynomial1D};
use crate::port::{Connected, Disconnected, FluidId, Port};
use crate::state_machine::StateManageable;
use crate::{
    CircuitId, Component, ComponentError, ConnectedPort, JacobianBuilder, OperationalState,
    ResidualVector, StateSlice,
};
use entropyk_core::{MassFlow, Power};
use serde::{Deserialize, Serialize};
use std::marker::PhantomData;

/// Pump performance curve coefficients.
///
/// Defines the polynomial coefficients for the pump's head-flow curve
/// and efficiency curve.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct PumpCurves {
    /// Performance curves (head, efficiency, optional power)
    curves: PerformanceCurves,
}

impl PumpCurves {
    /// Creates pump curves from performance curves.
    pub fn new(curves: PerformanceCurves) -> Result<Self, ComponentError> {
        curves.validate()?;
        Ok(Self { curves })
    }

    /// Creates pump curves from polynomial coefficients.
    ///
    /// # Arguments
    ///
    /// * `head_coeffs` - Head curve coefficients [a0, a1, a2, ...] for H = a0 + a1*Q + a2*Q²
    /// * `eff_coeffs` - Efficiency coefficients [b0, b1, b2, ...] for η = b0 + b1*Q + b2*Q²
    ///
    /// # Units
    ///
    /// * Q (flow) in m³/s
    /// * H (head) in meters
    /// * η (efficiency) as decimal (0.0 to 1.0)
    pub fn from_coefficients(
        head_coeffs: Vec<f64>,
        eff_coeffs: Vec<f64>,
    ) -> Result<Self, ComponentError> {
        let head_curve = Polynomial1D::new(head_coeffs);
        let eff_curve = Polynomial1D::new(eff_coeffs);
        let curves = PerformanceCurves::simple(head_curve, eff_curve);
        Self::new(curves)
    }

    /// Creates a quadratic pump curve.
    ///
    /// H = a0 + a1*Q + a2*Q²
    /// η = b0 + b1*Q + b2*Q²
    pub fn quadratic(
        h0: f64,
        h1: f64,
        h2: f64,
        e0: f64,
        e1: f64,
        e2: f64,
    ) -> Result<Self, ComponentError> {
        Self::from_coefficients(vec![h0, h1, h2], vec![e0, e1, e2])
    }

    /// Creates a cubic pump curve (3rd-order polynomial for head).
    ///
    /// H = a0 + a1*Q + a2*Q² + a3*Q³
    /// η = b0 + b1*Q + b2*Q²
    pub fn cubic(
        h0: f64,
        h1: f64,
        h2: f64,
        h3: f64,
        e0: f64,
        e1: f64,
        e2: f64,
    ) -> Result<Self, ComponentError> {
        Self::from_coefficients(vec![h0, h1, h2, h3], vec![e0, e1, e2])
    }

    /// Returns the head at the given flow rate (at full speed).
    ///
    /// # Arguments
    ///
    /// * `flow_m3_per_s` - Volumetric flow rate in m³/s
    ///
    /// # Returns
    ///
    /// Head in meters
    pub fn head_at_flow(&self, flow_m3_per_s: f64) -> f64 {
        self.curves.head_curve.evaluate(flow_m3_per_s)
    }

    /// Returns the efficiency at the given flow rate (at full speed).
    ///
    /// # Arguments
    ///
    /// * `flow_m3_per_s` - Volumetric flow rate in m³/s
    ///
    /// # Returns
    ///
    /// Efficiency as decimal (0.0 to 1.0)
    pub fn efficiency_at_flow(&self, flow_m3_per_s: f64) -> f64 {
        let eta = self.curves.efficiency_curve.evaluate(flow_m3_per_s);
        // Clamp efficiency to valid range
        eta.clamp(0.0, 1.0)
    }

    /// Returns reference to the performance curves.
    pub fn curves(&self) -> &PerformanceCurves {
        &self.curves
    }
}

impl Default for PumpCurves {
    fn default() -> Self {
        Self::quadratic(30.0, 0.0, 0.0, 0.7, 0.0, 0.0).unwrap()
    }
}

/// A pump component with polynomial performance curves.
///
/// The pump uses the Type-State pattern to ensure ports are connected
/// before use in simulations.
///
/// # Example
///
/// ```ignore
/// use entropyk_components::pump::{Pump, PumpCurves};
/// use entropyk_components::port::{FluidId, Port};
/// use entropyk_core::{Pressure, Enthalpy};
///
/// // Create pump curves: H = 30 - 10*Q - 50*Q² (in m and m³/s)
/// let curves = PumpCurves::quadratic(30.0, -10.0, -50.0, 0.5, 0.3, -0.5).unwrap();
///
/// let inlet = Port::new(
///     FluidId::new("Water"),
///     Pressure::from_bar(1.0),
///     Enthalpy::from_joules_per_kg(100000.0),
/// );
/// let outlet = Port::new(
///     FluidId::new("Water"),
///     Pressure::from_bar(1.0),
///     Enthalpy::from_joules_per_kg(100000.0),
/// );
///
/// let pump = Pump::new(curves, inlet, outlet, 1000.0).unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct Pump<State> {
    /// Performance curves
    curves: PumpCurves,
    /// Inlet port
    port_inlet: Port<State>,
    /// Outlet port
    port_outlet: Port<State>,
    /// Fluid density in kg/m³
    fluid_density_kg_per_m3: f64,
    /// Speed ratio (0.0 to 1.0), default 1.0 (full speed)
    speed_ratio: f64,
    /// Circuit identifier
    circuit_id: CircuitId,
    /// Operational state
    operational_state: OperationalState,
    inlet_m_idx: Option<usize>,
    inlet_p_idx: Option<usize>,
    inlet_h_idx: Option<usize>,
    outlet_m_idx: Option<usize>,
    outlet_p_idx: Option<usize>,
    outlet_h_idx: Option<usize>,
    /// Phantom data for type state
    _state: PhantomData<State>,
}

impl Pump<Disconnected> {
    /// Creates a new disconnected pump.
    ///
    /// # Arguments
    ///
    /// * `curves` - Pump performance curves
    /// * `port_inlet` - Inlet port (disconnected)
    /// * `port_outlet` - Outlet port (disconnected)
    /// * `fluid_density` - Fluid density in kg/m³
    ///
    /// # Errors
    ///
    /// Returns an error if:
    /// - Ports have different fluid types
    /// - Fluid density is not positive
    pub fn new(
        curves: PumpCurves,
        port_inlet: Port<Disconnected>,
        port_outlet: Port<Disconnected>,
        fluid_density: f64,
    ) -> Result<Self, ComponentError> {
        if port_inlet.fluid_id() != port_outlet.fluid_id() {
            return Err(ComponentError::InvalidState(
                "Inlet and outlet ports must have the same fluid type".to_string(),
            ));
        }

        if fluid_density <= 0.0 {
            return Err(ComponentError::InvalidState(
                "Fluid density must be positive".to_string(),
            ));
        }

        Ok(Self {
            curves,
            port_inlet,
            port_outlet,
            fluid_density_kg_per_m3: fluid_density,
            speed_ratio: 1.0,
            circuit_id: CircuitId::default(),
            operational_state: OperationalState::default(),
            inlet_m_idx: None,
            inlet_p_idx: None,
            inlet_h_idx: None,
            outlet_m_idx: None,
            outlet_p_idx: None,
            outlet_h_idx: None,
            _state: PhantomData,
        })
    }

    /// Returns the fluid identifier.
    pub fn fluid_id(&self) -> &FluidId {
        self.port_inlet.fluid_id()
    }

    /// Returns the inlet port.
    pub fn port_inlet(&self) -> &Port<Disconnected> {
        &self.port_inlet
    }

    /// Returns the outlet port.
    pub fn port_outlet(&self) -> &Port<Disconnected> {
        &self.port_outlet
    }

    /// Returns the fluid density.
    pub fn fluid_density(&self) -> f64 {
        self.fluid_density_kg_per_m3
    }

    /// Returns the performance curves.
    pub fn curves(&self) -> &PumpCurves {
        &self.curves
    }

    /// Returns the speed ratio.
    pub fn speed_ratio(&self) -> f64 {
        self.speed_ratio
    }

    /// Sets the speed ratio (0.0 to 1.0).
    pub fn set_speed_ratio(&mut self, ratio: f64) -> Result<(), ComponentError> {
        if !(0.0..=1.0).contains(&ratio) {
            return Err(ComponentError::InvalidState(
                "Speed ratio must be between 0.0 and 1.0".to_string(),
            ));
        }
        self.speed_ratio = ratio;
        Ok(())
    }

    /// Connects the pump to inlet and outlet ports.
    ///
    /// This consumes the disconnected pump and returns a connected one,
    /// transitioning the state at compile time.
    pub fn connect(
        self,
        inlet: Port<Disconnected>,
        outlet: Port<Disconnected>,
    ) -> Result<Pump<Connected>, ComponentError> {
        let (p_in, _) = self
            .port_inlet
            .connect(inlet)
            .map_err(|e| ComponentError::InvalidState(e.to_string()))?;
        let (p_out, _) = self
            .port_outlet
            .connect(outlet)
            .map_err(|e| ComponentError::InvalidState(e.to_string()))?;

        Ok(Pump {
            curves: self.curves,
            port_inlet: p_in,
            port_outlet: p_out,
            fluid_density_kg_per_m3: self.fluid_density_kg_per_m3,
            speed_ratio: self.speed_ratio,
            circuit_id: self.circuit_id,
            operational_state: self.operational_state,
            inlet_m_idx: None,
            inlet_p_idx: None,
            inlet_h_idx: None,
            outlet_m_idx: None,
            outlet_p_idx: None,
            outlet_h_idx: None,
            _state: PhantomData,
        })
    }
}

impl Pump<Connected> {
    /// Returns the inlet port.
    pub fn port_inlet(&self) -> &Port<Connected> {
        &self.port_inlet
    }

    /// Returns the outlet port.
    pub fn port_outlet(&self) -> &Port<Connected> {
        &self.port_outlet
    }

    /// Calculates the pressure rise across the pump.
    ///
    /// Uses the head curve and converts to pressure:
    /// ΔP = ρ × g × H
    ///
    /// Applies affinity laws for variable speed operation.
    ///
    /// # Arguments
    ///
    /// * `flow_m3_per_s` - Volumetric flow rate in m³/s
    ///
    /// # Returns
    ///
    /// Pressure rise in Pascals
    pub fn pressure_rise(&self, flow_m3_per_s: f64) -> f64 {
        // Handle zero speed - pump produces no pressure
        if self.speed_ratio <= 0.0 {
            return 0.0;
        }

        // Handle negative flow gracefully by using a linear extrapolation from Q=0
        // to prevent polynomial extrapolation issues with quadratic/cubic terms
        if flow_m3_per_s < 0.0 {
            let h0 = self.curves.head_at_flow(0.0);
            let h_eps = self.curves.head_at_flow(1e-6);
            let dh_dq = (h_eps - h0) / 1e-6;

            let head_m = h0 + dh_dq * flow_m3_per_s;
            let actual_head = AffinityLaws::scale_head(head_m, self.speed_ratio);
            const G: f64 = 9.80665; // m/s²
            return self.fluid_density_kg_per_m3 * G * actual_head;
        }

        // Handle exactly zero flow
        if flow_m3_per_s == 0.0 {
            // At zero flow, use the shut-off head scaled by speed
            let head_m = self.curves.head_at_flow(0.0);
            let actual_head = AffinityLaws::scale_head(head_m, self.speed_ratio);
            const G: f64 = 9.80665; // m/s²
            return self.fluid_density_kg_per_m3 * G * actual_head;
        }

        // Apply affinity law to get equivalent flow at full speed
        let equivalent_flow = AffinityLaws::unscale_flow(flow_m3_per_s, self.speed_ratio);

        // Get head at equivalent flow
        let head_m = self.curves.head_at_flow(equivalent_flow);

        // Apply affinity law to scale head back to actual speed
        let actual_head = AffinityLaws::scale_head(head_m, self.speed_ratio);

        // Convert head to pressure: P = ρ × g × H
        const G: f64 = 9.80665; // m/s²
        self.fluid_density_kg_per_m3 * G * actual_head
    }

    /// Calculates the efficiency at the given flow rate.
    ///
    /// Applies affinity laws to find the equivalent operating point.
    pub fn efficiency(&self, flow_m3_per_s: f64) -> f64 {
        // Handle zero speed - pump is not running
        if self.speed_ratio <= 0.0 {
            return 0.0;
        }

        // Handle zero flow
        if flow_m3_per_s <= 0.0 {
            return self.curves.efficiency_at_flow(0.0);
        }

        let equivalent_flow = AffinityLaws::unscale_flow(flow_m3_per_s, self.speed_ratio);
        self.curves.efficiency_at_flow(equivalent_flow)
    }

    /// Calculates the hydraulic power consumption.
    ///
    /// P_hydraulic = Q × ΔP / η
    ///
    /// # Arguments
    ///
    /// * `flow_m3_per_s` - Volumetric flow rate in m³/s
    ///
    /// # Returns
    ///
    /// Power in Watts
    pub fn hydraulic_power(&self, flow_m3_per_s: f64) -> Power {
        if flow_m3_per_s <= 0.0 || self.speed_ratio <= 0.0 {
            return Power::from_watts(0.0);
        }

        let delta_p = self.pressure_rise(flow_m3_per_s);
        let eta = self.efficiency(flow_m3_per_s);

        if eta <= 0.0 {
            return Power::from_watts(0.0);
        }

        // P = Q × ΔP / η
        let power_w = flow_m3_per_s * delta_p / eta;
        Power::from_watts(power_w)
    }

    /// Calculates mass flow rate from volumetric flow.
    pub fn mass_flow_from_volumetric(&self, flow_m3_per_s: f64) -> MassFlow {
        MassFlow::from_kg_per_s(flow_m3_per_s * self.fluid_density_kg_per_m3)
    }

    /// Calculates volumetric flow rate from mass flow.
    pub fn volumetric_from_mass_flow(&self, mass_flow: MassFlow) -> f64 {
        mass_flow.to_kg_per_s() / self.fluid_density_kg_per_m3
    }

    /// Returns the fluid density.
    pub fn fluid_density(&self) -> f64 {
        self.fluid_density_kg_per_m3
    }

    /// Returns the performance curves.
    pub fn curves(&self) -> &PumpCurves {
        &self.curves
    }

    /// Returns the speed ratio.
    pub fn speed_ratio(&self) -> f64 {
        self.speed_ratio
    }

    /// Sets the speed ratio (0.0 to 1.0).
    pub fn set_speed_ratio(&mut self, ratio: f64) -> Result<(), ComponentError> {
        if !(0.0..=1.0).contains(&ratio) {
            return Err(ComponentError::InvalidState(
                "Speed ratio must be between 0.0 and 1.0".to_string(),
            ));
        }
        self.speed_ratio = ratio;
        Ok(())
    }

    /// Returns both ports as a slice for solver topology.
    pub fn get_ports_slice(&self) -> [&Port<Connected>; 2] {
        [&self.port_inlet, &self.port_outlet]
    }
}

impl Component for Pump<Connected> {
    fn set_system_context(
        &mut self,
        _state_offset: usize,
        external_edge_state_indices: &[(usize, usize, usize)],
    ) {
        if let Some(&(m, p, h)) = external_edge_state_indices.first() {
            self.inlet_m_idx = Some(m);
            self.inlet_p_idx = Some(p);
            self.inlet_h_idx = Some(h);
        }
        if let Some(&(m, p, h)) = external_edge_state_indices.get(1) {
            self.outlet_m_idx = Some(m);
            self.outlet_p_idx = Some(p);
            self.outlet_h_idx = Some(h);
        }
    }

    fn compute_residuals(
        &self,
        state: &StateSlice,
        residuals: &mut ResidualVector,
    ) -> Result<(), ComponentError> {
        if residuals.len() != self.n_equations() {
            return Err(ComponentError::InvalidResidualDimensions {
                expected: self.n_equations(),
                actual: residuals.len(),
            });
        }

        let (in_m, in_p, in_h, out_p, out_h) = match (
            self.inlet_m_idx,
            self.inlet_p_idx,
            self.inlet_h_idx,
            self.outlet_p_idx,
            self.outlet_h_idx,
        ) {
            (Some(in_m), Some(in_p), Some(in_h), Some(out_p), Some(out_h)) => {
                (in_m, in_p, in_h, out_p, out_h)
            }
            _ => {
                return Err(ComponentError::InvalidState(
                    "Pump requires live inlet and outlet edge state indices".to_string(),
                ));
            }
        };
        let max_idx = in_m.max(in_p).max(in_h).max(out_p).max(out_h);
        if max_idx >= state.len() {
            return Err(ComponentError::InvalidStateDimensions {
                expected: max_idx + 1,
                actual: state.len(),
            });
        }

        match self.operational_state {
            OperationalState::Off => {
                residuals[0] = state[out_p] - state[in_p];
                residuals[1] = state[out_h] - state[in_h];
                return Ok(());
            }
            OperationalState::Bypass => {
                residuals[0] = state[out_p] - state[in_p];
                residuals[1] = state[out_h] - state[in_h];
                return Ok(());
            }
            OperationalState::On => {}
        }

        let mass_flow_kg_s = state[in_m];
        let flow_m3_s = mass_flow_kg_s / self.fluid_density_kg_per_m3;
        let delta_p_calc = self.pressure_rise(flow_m3_s);
        residuals[0] = state[out_p] - (state[in_p] + delta_p_calc);
        let power_calc = self.hydraulic_power(flow_m3_s).to_watts();
        let enthalpy_rise_j_kg = if mass_flow_kg_s.abs() > 1e-12 {
            power_calc / mass_flow_kg_s
        } else {
            0.0
        };
        residuals[1] = state[out_h] - (state[in_h] + enthalpy_rise_j_kg);

        Ok(())
    }

    fn jacobian_entries(
        &self,
        state: &StateSlice,
        jacobian: &mut JacobianBuilder,
    ) -> Result<(), ComponentError> {
        let (in_m, in_p, in_h, out_p, out_h) = match (
            self.inlet_m_idx,
            self.inlet_p_idx,
            self.inlet_h_idx,
            self.outlet_p_idx,
            self.outlet_h_idx,
        ) {
            (Some(in_m), Some(in_p), Some(in_h), Some(out_p), Some(out_h)) => {
                (in_m, in_p, in_h, out_p, out_h)
            }
            _ => {
                return Err(ComponentError::InvalidState(
                    "Pump Jacobian requires live inlet and outlet edge state indices".to_string(),
                ));
            }
        };
        if in_m >= state.len() {
            return Err(ComponentError::InvalidStateDimensions {
                expected: in_m + 1,
                actual: state.len(),
            });
        }
        let mass_flow_kg_s = state[in_m];
        let flow_m3_s = mass_flow_kg_s / self.fluid_density_kg_per_m3;

        // Numerical derivative of pressure with respect to mass flow
        let h = 0.001;
        let p_plus = self.pressure_rise(flow_m3_s + h / self.fluid_density_kg_per_m3);
        let p_minus = self.pressure_rise(flow_m3_s - h / self.fluid_density_kg_per_m3);
        let dp_dm = (p_plus - p_minus) / (2.0 * h);

        // ∂r₀/∂ṁ = dΔP/dṁ
        jacobian.add_entry(0, in_m, -dp_dm);
        jacobian.add_entry(0, out_p, 1.0);
        jacobian.add_entry(0, in_p, -1.0);

        // Numerical derivative of power with respect to mass flow
        let pow_plus = self
            .hydraulic_power(flow_m3_s + h / self.fluid_density_kg_per_m3)
            .to_watts();
        let pow_minus = self
            .hydraulic_power(flow_m3_s - h / self.fluid_density_kg_per_m3)
            .to_watts();
        let dpow_dm = (pow_plus - pow_minus) / (2.0 * h);

        let dh_dm = if mass_flow_kg_s.abs() > 1e-12 {
            (dpow_dm * mass_flow_kg_s - self.hydraulic_power(flow_m3_s).to_watts())
                / (mass_flow_kg_s * mass_flow_kg_s)
        } else {
            0.0
        };
        jacobian.add_entry(1, in_m, -dh_dm);
        jacobian.add_entry(1, out_h, 1.0);
        jacobian.add_entry(1, in_h, -1.0);

        Ok(())
    }

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

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

    fn port_mass_flows(
        &self,
        state: &StateSlice,
    ) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
        let (Some(in_m), Some(out_m)) = (self.inlet_m_idx, self.outlet_m_idx) else {
            return Err(ComponentError::InvalidState(
                "Pump mass-flow reporting requires live inlet and outlet mass-flow indices"
                    .to_string(),
            ));
        };
        let max_idx = in_m.max(out_m);
        if max_idx >= state.len() {
            return Err(ComponentError::InvalidStateDimensions {
                expected: max_idx + 1,
                actual: state.len(),
            });
        }
        Ok(vec![
            entropyk_core::MassFlow::from_kg_per_s(state[in_m]),
            entropyk_core::MassFlow::from_kg_per_s(-state[out_m]),
        ])
    }

    fn port_enthalpies(
        &self,
        _state: &StateSlice,
    ) -> Result<Vec<entropyk_core::Enthalpy>, ComponentError> {
        let (Some(in_h), Some(out_h)) = (self.inlet_h_idx, self.outlet_h_idx) else {
            return Err(ComponentError::InvalidState(
                "Pump enthalpy reporting requires live inlet and outlet enthalpy indices"
                    .to_string(),
            ));
        };
        let max_idx = in_h.max(out_h);
        if max_idx >= _state.len() {
            return Err(ComponentError::InvalidStateDimensions {
                expected: max_idx + 1,
                actual: _state.len(),
            });
        }
        Ok(vec![
            entropyk_core::Enthalpy::from_joules_per_kg(_state[in_h]),
            entropyk_core::Enthalpy::from_joules_per_kg(_state[out_h]),
        ])
    }

    fn energy_transfers(
        &self,
        state: &StateSlice,
    ) -> Option<(entropyk_core::Power, entropyk_core::Power)> {
        match self.operational_state {
            OperationalState::Off | OperationalState::Bypass => Some((
                entropyk_core::Power::from_watts(0.0),
                entropyk_core::Power::from_watts(0.0),
            )),
            OperationalState::On => {
                let Some(in_m) = self.inlet_m_idx else {
                    return None;
                };
                let mass_flow_kg_s = *state.get(in_m)?;
                let flow_m3_s = mass_flow_kg_s / self.fluid_density_kg_per_m3;
                let power_calc = self.hydraulic_power(flow_m3_s).to_watts();
                Some((
                    entropyk_core::Power::from_watts(0.0),
                    entropyk_core::Power::from_watts(-power_calc),
                ))
            }
        }
    }

    fn signature(&self) -> String {
        format!("Pump(circuit={})", self.circuit_id.0)
    }

    fn to_params(&self) -> crate::ComponentParams {
        crate::ComponentParams::new("Pump")
            .with_param("circuitId", self.circuit_id.0)
            .with_param("fluidDensityKgPerM3", self.fluid_density_kg_per_m3)
            .with_param("speedRatio", self.speed_ratio)
    }
}

impl StateManageable for Pump<Connected> {
    fn state(&self) -> OperationalState {
        self.operational_state
    }

    fn set_state(&mut self, state: OperationalState) -> Result<(), ComponentError> {
        if self.operational_state.can_transition_to(state) {
            let from = self.operational_state;
            self.operational_state = state;
            self.on_state_change(from, state);
            Ok(())
        } else {
            Err(ComponentError::InvalidStateTransition {
                from: self.operational_state,
                to: state,
                reason: "Transition not allowed".to_string(),
            })
        }
    }

    fn can_transition_to(&self, target: OperationalState) -> bool {
        self.operational_state.can_transition_to(target)
    }

    fn circuit_id(&self) -> &CircuitId {
        &self.circuit_id
    }

    fn set_circuit_id(&mut self, circuit_id: CircuitId) {
        self.circuit_id = circuit_id;
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::port::FluidId;
    use approx::assert_relative_eq;
    use entropyk_core::{Enthalpy, Pressure};

    fn create_test_curves() -> PumpCurves {
        // Typical small pump:
        // H = 30 - 10*Q - 50*Q² (m, Q in m³/s)
        // η = 0.6 + 1.0*Q - 2.0*Q²
        PumpCurves::quadratic(30.0, -10.0, -50.0, 0.6, 1.0, -2.0).unwrap()
    }

    fn create_test_pump_disconnected() -> Pump<Disconnected> {
        let curves = create_test_curves();
        let inlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        let outlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        Pump::new(curves, inlet, outlet, 1000.0).unwrap()
    }

    fn create_test_pump_connected() -> Pump<Connected> {
        let curves = create_test_curves();
        let inlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        let outlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        let (inlet_conn, outlet_conn) = inlet.connect(outlet).unwrap();

        Pump {
            curves,
            port_inlet: inlet_conn,
            port_outlet: outlet_conn,
            fluid_density_kg_per_m3: 1000.0,
            speed_ratio: 1.0,
            circuit_id: CircuitId::default(),
            operational_state: OperationalState::default(),
            inlet_m_idx: None,
            inlet_p_idx: None,
            inlet_h_idx: None,
            outlet_m_idx: None,
            outlet_p_idx: None,
            outlet_h_idx: None,
            _state: PhantomData,
        }
    }

    #[test]
    fn test_pump_curves_creation() {
        let curves = create_test_curves();
        assert_eq!(curves.head_at_flow(0.0), 30.0);
        assert_relative_eq!(curves.efficiency_at_flow(0.0), 0.6);
    }

    #[test]
    fn test_pump_curves_head() {
        let curves = create_test_curves();
        // H = 30 - 10*0.5 - 50*0.25 = 30 - 5 - 12.5 = 12.5 m
        let head = curves.head_at_flow(0.5);
        assert_relative_eq!(head, 12.5, epsilon = 1e-10);
    }

    #[test]
    fn test_pump_curves_efficiency_clamped() {
        let curves = create_test_curves();
        // At very high flow, efficiency might go negative
        // Should be clamped to 0
        let eff = curves.efficiency_at_flow(10.0);
        assert!(eff >= 0.0);
    }

    #[test]
    fn test_pump_creation() {
        let pump = create_test_pump_disconnected();
        assert_eq!(pump.fluid_density(), 1000.0);
        assert_eq!(pump.speed_ratio(), 1.0);
    }

    #[test]
    fn test_pump_invalid_density() {
        let curves = create_test_curves();
        let inlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        let outlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );

        let result = Pump::new(curves, inlet, outlet, -1.0);
        assert!(result.is_err());
    }

    #[test]
    fn test_pump_different_fluids() {
        let curves = create_test_curves();
        let inlet = Port::new(
            FluidId::new("Water"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );
        let outlet = Port::new(
            FluidId::new("Glycol"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(100000.0),
        );

        let result = Pump::new(curves, inlet, outlet, 1000.0);
        assert!(result.is_err());
    }

    #[test]
    fn test_pump_set_speed_ratio() {
        let mut pump = create_test_pump_connected();
        assert!(pump.set_speed_ratio(0.8).is_ok());
        assert_eq!(pump.speed_ratio(), 0.8);
    }

    #[test]
    fn test_pump_set_speed_ratio_invalid() {
        let mut pump = create_test_pump_connected();
        assert!(pump.set_speed_ratio(1.5).is_err());
        assert!(pump.set_speed_ratio(-0.1).is_err());
    }

    #[test]
    fn test_pump_pressure_rise_full_speed() {
        let pump = create_test_pump_connected();
        // At Q=0: H=30m, P = 1000 * 9.8 * 30 ≈ 294200 Pa
        let delta_p = pump.pressure_rise(0.0);
        let expected = 1000.0 * 9.80665 * 30.0;
        assert_relative_eq!(delta_p, expected, epsilon = 100.0);
    }

    #[test]
    fn test_pump_pressure_rise_reduced_speed() {
        let mut pump = create_test_pump_connected();
        pump.set_speed_ratio(0.5).unwrap();

        // At 50% speed, shut-off head is 25% of full speed
        // H = 0.25 * 30 = 7.5 m
        let delta_p = pump.pressure_rise(0.0);
        let expected = 1000.0 * 9.80665 * 7.5;
        assert_relative_eq!(delta_p, expected, epsilon = 100.0);
    }

    #[test]
    fn test_pump_hydraulic_power() {
        let pump = create_test_pump_connected();

        // At Q=0.1 m³/s: H ≈ 30 - 1 - 0.5 = 28.5 m
        // η ≈ 0.6 + 0.1 - 0.02 = 0.68
        // P = 1000 * 9.8 * 0.1 * 28.5 / 0.68 ≈ 4110 W
        let power = pump.hydraulic_power(0.1);
        assert!(power.to_watts() > 0.0);
        assert!(power.to_watts() < 50000.0);
    }

    #[test]
    fn test_pump_affinity_laws_power() {
        let pump_full = create_test_pump_connected();

        let mut pump_half = create_test_pump_connected();
        pump_half.set_speed_ratio(0.5).unwrap();

        // Power at half speed should be ~12.5% of full speed (cube law)
        // At the same equivalent flow point
        let power_full = pump_full.hydraulic_power(0.1);
        let power_half = pump_half.hydraulic_power(0.05); // Half the flow

        // P_half / P_full ≈ 0.5³ = 0.125
        let ratio = power_half.to_watts() / power_full.to_watts();
        assert_relative_eq!(ratio, 0.125, epsilon = 0.05);
    }

    #[test]
    fn test_pump_component_n_equations() {
        let pump = create_test_pump_connected();
        assert_eq!(pump.n_equations(), 2);
    }

    #[test]
    fn test_pump_component_compute_residuals() {
        let mut pump = create_test_pump_connected();
        pump.set_system_context(0, &[(0, 1, 2), (3, 4, 5)]);
        let state = vec![50.0, 1.0e5, 100_000.0, 50.0, 4.0e5, 105_000.0];
        let mut residuals = vec![0.0; 2];

        let result = pump.compute_residuals(&state, &mut residuals);
        assert!(result.is_ok());
    }

    #[test]
    fn test_pump_state_manageable() {
        let pump = create_test_pump_connected();
        assert_eq!(pump.state(), OperationalState::On);
        assert!(pump.can_transition_to(OperationalState::Off));
    }
}