hydra_engine_wds/hydraulics/
shared.rs

1/// Large conductance constant C∞ (§3.3, §3.5). Same as EPANET's CBIG.
2pub(super) const C_INF: f64 = 1.0e8;
3
4/// Minimum head-loss gradient g_min (§3.3). EPANET: CSMALL = 1e-6.
5pub(super) const G_MIN: f64 = 1.0e-6;
6
7/// Placeholder flow for closed links at initialisation (§3.10).
8pub(super) const Q_CLOSED: f64 = 1.0e-6;
9
10/// Specific weight of water, ρg (N/m³). Used for constant-HP pump.
11pub(super) const GAMMA_WATER: f64 = 9810.0;
12
13/// Power-function pump head-curve coefficients (§3.2.5).
14///
15/// For a pump at relative speed ω the head gain is:
16///   ΔH = ω² H₀ − r ω^(2−N) Q^N
17#[derive(Debug, Clone)]
18pub struct PumpCoeffs {
19    /// Shutoff head H₀ (m).
20    pub h0: f64,
21    /// Resistance coefficient r.
22    pub r: f64,
23    /// Flow exponent N.
24    pub n: f64,
25}
26
27impl PumpCoeffs {
28    /// Fits coefficients from three head-curve points `(0, h_shutoff)`,
29    /// `(q1, h1)`, `(q2, h2)` using the three-point formula (§3.2.5).
30    ///
31    /// Returns `None` when the constraints `h_shutoff > h1 > h2 > 0`,
32    /// `q2 > q1 > 0`, and `0 < N ≤ 20` are not all satisfied.
33    pub fn from_three_points(h_shutoff: f64, q1: f64, h1: f64, q2: f64, h2: f64) -> Option<Self> {
34        if h_shutoff <= h1 || h1 <= h2 || h2 < 0.0 {
35            return None;
36        }
37        if q1 <= 0.0 || q2 <= q1 {
38            return None;
39        }
40        let n = ((h_shutoff - h2) / (h_shutoff - h1)).ln() / (q2 / q1).ln();
41        if n <= 0.0 || n > 20.0 {
42            return None;
43        }
44        let r = (h_shutoff - h1) / q1.powf(n);
45        if r <= 0.0 {
46            return None;
47        }
48        Some(PumpCoeffs {
49            h0: h_shutoff,
50            r,
51            n,
52        })
53    }
54}
55
56/// Lower barrier: adds (Δh, Δg) to enforce q ≥ q₀ when called with δq = q − q₀.
57pub(super) fn lower_barrier(dq: f64) -> (f64, f64) {
58    let a = 1.0e9 * dq;
59    let b = (a * a + 1.0e-6_f64).sqrt();
60    ((a - b) / 2.0, 5.0e8 * (1.0 - a / b))
61}
62
63/// Upper barrier: adds (Δh, Δg) to enforce q ≤ q₁ when called with δq = q − q₁.
64pub(super) fn upper_barrier(dq: f64) -> (f64, f64) {
65    let a = 1.0e9 * dq;
66    let b = (a * a + 1.0e-6_f64).sqrt();
67    ((a + b) / 2.0, 5.0e8 * (1.0 + a / b))
68}
69
70/// Per-link P/Y coefficient pair (§3.3).
71pub(super) struct Py {
72    pub(super) p: f64,
73    pub(super) y: f64,
74}
75
76impl Py {
77    pub(super) fn closed(q: f64) -> Self {
78        Py {
79            p: 1.0 / C_INF,
80            y: q,
81        }
82    }
83
84    pub(super) fn from_hg(h: f64, g: f64, n_f: f64, q: f64) -> Self {
85        if g < G_MIN {
86            let g2 = G_MIN / n_f;
87            let h2 = g2 * q;
88            Py {
89                p: 1.0 / g2,
90                y: h2 / g2,
91            }
92        } else {
93            Py {
94                p: 1.0 / g,
95                y: h / g,
96            }
97        }
98    }
99}
100
101/// Outcome of a completed hydraulic solve.
102#[derive(Debug, Clone, Copy, PartialEq, Eq)]
103pub enum SolveResult {
104    /// All convergence criteria satisfied.
105    Converged,
106    /// `max_iter` reached; extra iterations run with frozen status. Results are
107    /// valid but potentially unbalanced.
108    Unbalanced,
109}
110
111/// Hydraulic solver error.
112#[derive(Debug, Clone)]
113pub enum HydraulicError {
114    /// Newton loop did not converge and `extra_iter = -1` (halt-on-failure mode).
115    NotConverged,
116    /// Sparse matrix numerically singular at the given permuted junction step.
117    SingularMatrix {
118        /// Zero-based permuted-junction index at which the singular pivot was detected.
119        junction_step: usize,
120    },
121    /// A pump that requires power-function coefficients has none fitted.
122    NoPumpCoeffs {
123        /// ID of the pump whose power-function coefficients could not be computed.
124        pump_id: String,
125    },
126}
127
128impl std::fmt::Display for HydraulicError {
129    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
130        match self {
131            Self::NotConverged => write!(f, "hydraulic solve did not converge"),
132            Self::SingularMatrix { junction_step } => {
133                write!(
134                    f,
135                    "sparse matrix singular at permuted junction step {junction_step}"
136                )
137            }
138            Self::NoPumpCoeffs { pump_id } => {
139                write!(
140                    f,
141                    "pump '{pump_id}' requires power-function coefficients but none were fitted"
142                )
143            }
144        }
145    }
146}
147
148#[cfg(test)]
149mod tests {
150    use super::{lower_barrier, upper_barrier, PumpCoeffs};
151
152    #[test]
153    fn pump_coeffs_three_point_fit() {
154        let c = PumpCoeffs::from_three_points(100.0, 1.0, 90.0, 2.0, 60.0).unwrap();
155        assert!((c.n - 2.0).abs() < 1.0e-10, "N={}", c.n);
156        assert!((c.r - 10.0).abs() < 1.0e-10, "r={}", c.r);
157        assert!((c.h0 - 100.0).abs() < 1.0e-10);
158    }
159
160    #[test]
161    fn pump_coeffs_rejects_bad_points() {
162        assert!(PumpCoeffs::from_three_points(90.0, 1.0, 90.0, 2.0, 60.0).is_none());
163        assert!(PumpCoeffs::from_three_points(100.0, 0.0, 90.0, 2.0, 60.0).is_none());
164        assert!(PumpCoeffs::from_three_points(100.0, 1.0, 99.99999, 2.0, 0.001).is_none());
165    }
166
167    #[test]
168    fn lower_barrier_feasible_near_zero() {
169        let (dh, _) = lower_barrier(1.0);
170        assert!(dh.abs() < 1.0e-3, "Δh={dh}");
171    }
172
173    #[test]
174    fn lower_barrier_violation_strongly_negative() {
175        let (dh, _) = lower_barrier(-1.0);
176        assert!(dh < -1.0e3, "Δh={dh}");
177    }
178
179    #[test]
180    fn upper_barrier_feasible_near_zero() {
181        let (dh, _) = upper_barrier(-1.0);
182        assert!(dh.abs() < 1.0e-3, "Δh={dh}");
183    }
184}
hydra_engine_wds/hydraulics/shared.rs

hydra_engine_wds/hydraulics/
shared.rs
