AI_Assistant/cad-editor
Watch
1
0
Fork
You've already forked cad-editor
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add Storm Sewer module: storm-drain network analysis (Rational + Manning + HGL)

Browse source 
Adds a "Storm Sewer" ribbon tab backed by a new in-repo `stormsewer` engine
crate (vendored under crates/). The engine implements standard public-domain
methods — Rational-method peak-flow accumulation, Manning circular-pipe
hydraulics, time-of-concentration / IDF intensity, and an HGL backwater pass
with junction losses — with 24 unit tests and no external dependencies.

Module (src/modules/storm_sewer/):
- Drafting: SS_INLET / SS_JUNCTION / SS_OUTFALL place structure markers and
  SS_PIPE chains pipe runs (interactive CadCommands).
- Analysis: SS_ANALYZE / SS_REPORT / SS_PROFILE open a .ssn network file, run
  the engine, draw the plan and HGL long-section, and print a report.
- SVG icons; SS_* commands registered for command-line autocomplete.

Integration is additive: build.rs auto-discovers the module and one dispatch
block is added to src/app/commands.rs. 7 module integration tests.

See src/modules/storm_sewer/INTEGRATION.md for the design and the optional
follow-up (drive analysis from canvas geometry via XDATA).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
Michael Flynn 2026-06-08 17:16:30 -04:00
commit 9b44eb6400
26 changed files with 2007 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

15
crates/stormsewer/Cargo.toml Normal file
View file

@ -0,0 +1,15 @@
[package]
name = "stormsewer"
version = "0.1.0"
edition = "2021"
description = "Native-Rust storm sewer network hydrology & hydraulics (Rational method, Manning, HGL). Engine core for the Open CAD Studio storm-sewer module and HydroComplete."
license = "GPL-3.0-or-later"
repository = "https://github.com/mf4633/stormsewer"
keywords = ["hydrology", "hydraulics", "stormwater", "civil-engineering", "cad"]
categories = ["science", "simulation"]
# Engine is std-only on purpose: dependency-light and WASM-friendly.
[dependencies]
[lib]
name = "stormsewer"

18
crates/stormsewer/examples/sample.ssn Normal file
View file

@ -0,0 +1,18 @@
# Sample storm-sewer network (stormsewer-cli examples/sample.ssn)
# A short, properly-sized trunk line: two catch basins + a junction to outfall.
IDF 60 10 0.8 # i = 60 / (t+10)^0.8 (in/hr, t in minutes) ~ 5 in/hr
TAILWATER 100.5 # outfall tailwater elevation (ft)
MINTC 10 # minimum time of concentration (min)
JUNCTIONK 0.5 # junction loss coefficient
# id kind x y invert rim area C tc_inlet
NODE N1 inlet 0 0 104.0 110.0 1.0 0.70 12
NODE N2 inlet 300 0 102.5 108.5 1.0 0.70 10
NODE N3 junction 550 0 101.2 107.0 0.5 0.80 8
NODE OUT outfall 730 0 100.0 106.0
# id from to length dia n (dia: 15in=1.25, 18in=1.5, 21in=1.75)
PIPE P1 N1 N2 300 1.25 0.013
PIPE P2 N2 N3 250 1.50 0.013
PIPE P3 N3 OUT 180 1.75 0.013

29
crates/stormsewer/src/bin/stormsewer-cli.rs Normal file
View file

@ -0,0 +1,29 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! `stormsewer-cli` — run a storm-sewer analysis from a `.ssn` network file.
//!
//! Usage: stormsewer-cli <network-file>
//!
//! See `stormsewer::parse` for the file format.
use std::process::exit;
use stormsewer::parse::parse_ssn;
use stormsewer::report::format_analysis;
fn die(msg: &str) -> ! {
eprintln!("error: {msg}");
exit(1);
}
fn main() {
let path = match std::env::args().nth(1) {
Some(p) => p,
None => die("usage: stormsewer-cli <network-file>"),
};
let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("cannot read {path}: {e}")));
let parsed = parse_ssn(&text).unwrap_or_else(|e| die(&e));
match parsed.network.analyze(&parsed.idf, &parsed.options) {
Ok(a) => print!("{}", format_analysis(&a)),
Err(e) => die(&e.to_string()),
}
}

288
crates/stormsewer/src/drawing.rs Normal file
View file

@ -0,0 +1,288 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Convert an analyzed network into drawable primitives (CAD-agnostic).
//!
//! Produces a **plan** view (pipes as lines, structures as markers, flow/HGL
//! labels) and a **profile** (HGL long-section of the main stem: ground,
//! invert, and hydraulic-grade-line polylines). Coordinates are plain `f64`
//! drawing units; the host CAD turns these into its own entities.
use crate::network::{Analysis, Network, NodeKind};
use std::collections::HashMap;
/// A pipe drawn in plan.
#[derive(Clone, Debug)]
pub struct PlanPipe {
pub x1: f64,
pub y1: f64,
pub x2: f64,
pub y2: f64,
pub surcharged: bool,
}
/// A structure marker drawn in plan.
#[derive(Clone, Debug)]
pub struct PlanNode {
pub x: f64,
pub y: f64,
pub radius: f64,
pub kind: NodeKind,
}
/// A text label placed at a point.
#[derive(Clone, Debug)]
pub struct Label {
pub x: f64,
pub y: f64,
pub text: String,
pub height: f64,
}
/// Which line of the profile a polyline represents.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum ProfileRole {
Ground,
Invert,
Hgl,
}
/// A polyline (sequence of points) with a role.
#[derive(Clone, Debug)]
pub struct Polyline {
pub pts: Vec<(f64, f64)>,
pub role: ProfileRole,
}
/// The full set of primitives for a network drawing.
#[derive(Clone, Debug, Default)]
pub struct NetworkDrawing {
pub plan_pipes: Vec<PlanPipe>,
pub plan_nodes: Vec<PlanNode>,
pub plan_labels: Vec<Label>,
pub profile_lines: Vec<Polyline>,
pub profile_labels: Vec<Label>,
}
/// Layout / styling knobs for [`draw_network`].
#[derive(Clone, Debug)]
pub struct DrawConfig {
pub text_height: f64,
pub node_radius: f64,
/// Plan X mapped to profile X=0 (station origin in drawing units).
pub profile_origin_x: f64,
/// Drawing Y at which the profile datum elevation is plotted.
pub profile_origin_y: f64,
/// Horizontal scale: station feet → drawing units.
pub h_scale: f64,
/// Vertical exaggeration applied to elevations in the profile.
pub v_exag: f64,
}
impl Default for DrawConfig {
fn default() -> Self {
Self {
text_height: 5.0,
node_radius: 3.0,
profile_origin_x: 0.0,
profile_origin_y: -200.0,
h_scale: 1.0,
v_exag: 10.0,
}
}
}
/// Build plan + profile primitives for an analyzed network.
pub fn draw_network(net: &Network, a: &Analysis, cfg: &DrawConfig) -> NetworkDrawing {
let mut d = NetworkDrawing::default();
let pos: HashMap<&str, (f64, f64)> =
net.nodes.iter().map(|n| (n.id.as_str(), (n.x, n.y))).collect();
let hgl: HashMap<&str, f64> = a.nodes.iter().map(|n| (n.id.as_str(), n.hgl)).collect();
// ── Plan: pipes + labels ────────────────────────────────────────────────
for pr in &a.pipes {
let (x1, y1) = pos[pr.from.as_str()];
let (x2, y2) = pos[pr.to.as_str()];
d.plan_pipes.push(PlanPipe { x1, y1, x2, y2, surcharged: pr.surcharged });
d.plan_labels.push(Label {
x: (x1 + x2) / 2.0,
y: (y1 + y2) / 2.0 + cfg.text_height,
text: format!("{}: {:.1} cfs {:.0}%", pr.id, pr.design_q, pr.pct_full * 100.0),
height: cfg.text_height,
});
}
// ── Plan: structure markers + labels ────────────────────────────────────
for n in &net.nodes {
d.plan_nodes.push(PlanNode { x: n.x, y: n.y, radius: cfg.node_radius, kind: n.kind });
let h = hgl.get(n.id.as_str()).copied().unwrap_or(f64::NAN);
let label = if h.is_finite() {
format!("{} HGL {:.1}", n.id, h)
} else {
n.id.clone()
};
d.plan_labels.push(Label {
x: n.x + cfg.node_radius,
y: n.y + cfg.node_radius,
text: label,
height: cfg.text_height,
});
}
// ── Profile of the main stem ────────────────────────────────────────────
let stem = main_stem(net);
if stem.len() >= 2 {
let nidx: HashMap<&str, usize> =
net.nodes.iter().enumerate().map(|(i, n)| (n.id.as_str(), i)).collect();
// Stations (ft) from the upstream end.
let mut stations = vec![0.0f64; stem.len()];
for k in 1..stem.len() {
let len = pipe_between(net, stem[k - 1], stem[k]).map(|p| p.length).unwrap_or(0.0);
stations[k] = stations[k - 1] + len;
}
let datum = stem.iter().map(|&i| net.nodes[i].invert).fold(f64::INFINITY, f64::min);
let px = |st: f64| cfg.profile_origin_x + st * cfg.h_scale;
let py = |elev: f64| cfg.profile_origin_y + (elev - datum) * cfg.v_exag;
let mut ground = Vec::new();
let mut invert = Vec::new();
let mut hgl_line = Vec::new();
for (k, &i) in stem.iter().enumerate() {
let n = &net.nodes[i];
let st = stations[k];
ground.push((px(st), py(n.rim)));
invert.push((px(st), py(n.invert)));
if let Some(&h) = hgl.get(n.id.as_str()) {
if h.is_finite() {
hgl_line.push((px(st), py(h)));
}
}
d.profile_labels.push(Label {
x: px(st),
y: py(n.rim) + cfg.text_height,
text: n.id.clone(),
height: cfg.text_height,
});
let _ = nidx; // (reserved for future cross-refs)
}
d.profile_lines.push(Polyline { pts: ground, role: ProfileRole::Ground });
d.profile_lines.push(Polyline { pts: invert, role: ProfileRole::Invert });
if hgl_line.len() >= 2 {
d.profile_lines.push(Polyline { pts: hgl_line, role: ProfileRole::Hgl });
}
}
d
}
/// The main trunk, upstream-first: walk from the outfall up the incoming pipe
/// whose upstream node carries the most accumulated drainage area.
fn main_stem(net: &Network) -> Vec<usize> {
let n = net.nodes.len();
if n == 0 {
return Vec::new();
}
let nidx: HashMap<&str, usize> =
net.nodes.iter().enumerate().map(|(i, nd)| (nd.id.as_str(), i)).collect();
let mut incoming: Vec<Vec<usize>> = vec![Vec::new(); n]; // upstream node indices
let mut has_out = vec![false; n];
for p in &net.pipes {
if let (Some(&u), Some(&v)) = (nidx.get(p.from.as_str()), nidx.get(p.to.as_str())) {
incoming[v].push(u);
has_out[u] = true;
}
}
let ca = net.accumulate_ca().unwrap_or_default();
let size = |i: usize| ca.get(net.nodes[i].id.as_str()).copied().unwrap_or(0.0);
// Start at the outfall (kind, else the first node with no outgoing pipe).
let start = net
.nodes
.iter()
.position(|nd| nd.kind == NodeKind::Outfall)
.or_else(|| (0..n).find(|&i| !has_out[i]))
.unwrap_or(0);
let mut stem = vec![start];
let mut cur = start;
let mut guard = 0;
while guard < n {
guard += 1;
match incoming[cur].iter().copied().max_by(|&a, &b| {
size(a).partial_cmp(&size(b)).unwrap_or(std::cmp::Ordering::Equal)
}) {
Some(up) => {
stem.push(up);
cur = up;
}
None => break,
}
}
stem.reverse(); // upstream-first
stem
}
fn pipe_between<'a>(net: &'a Network, up: usize, dn: usize) -> Option<&'a crate::network::Pipe> {
let up_id = net.nodes[up].id.as_str();
let dn_id = net.nodes[dn].id.as_str();
net.pipes.iter().find(|p| p.from == up_id && p.to == dn_id)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::idf::IdfCurve;
use crate::network::{AnalysisOptions, Node, Pipe};
fn sample() -> Network {
Network {
nodes: vec![
Node::inlet("N1", 104.0, 110.0, 1.0, 0.70).at(0.0, 0.0),
Node::inlet("N2", 102.5, 108.5, 1.0, 0.70).at(300.0, 0.0),
Node::junction("N3", 101.2, 107.0, 0.5, 0.80).at(550.0, 0.0),
Node::outfall("OUT", 100.0, 106.0).at(730.0, 0.0),
],
pipes: vec![
Pipe::new("P1", "N1", "N2", 300.0, 1.25, 0.013),
Pipe::new("P2", "N2", "N3", 250.0, 1.50, 0.013),
Pipe::new("P3", "N3", "OUT", 180.0, 1.75, 0.013),
],
}
}
fn analyzed() -> (Network, Analysis) {
let net = sample();
let a = net.analyze(&IdfCurve::new(60.0, 10.0, 0.8), &AnalysisOptions { tailwater: Some(100.5), ..Default::default() }).unwrap();
(net, a)
}
#[test]
fn plan_has_one_line_per_pipe_and_marker_per_node() {
let (net, a) = analyzed();
let d = draw_network(&net, &a, &DrawConfig::default());
assert_eq!(d.plan_pipes.len(), net.pipes.len());
assert_eq!(d.plan_nodes.len(), net.nodes.len());
assert!(d.plan_labels.len() >= net.pipes.len() + net.nodes.len());
}
#[test]
fn main_stem_is_full_trunk_upstream_first() {
let net = sample();
let stem = main_stem(&net);
let ids: Vec<&str> = stem.iter().map(|&i| net.nodes[i].id.as_str()).collect();
assert_eq!(ids, vec!["N1", "N2", "N3", "OUT"]);
}
#[test]
fn profile_has_ground_invert_and_hgl() {
let (net, a) = analyzed();
let d = draw_network(&net, &a, &DrawConfig::default());
let roles: Vec<ProfileRole> = d.profile_lines.iter().map(|p| p.role).collect();
assert!(roles.contains(&ProfileRole::Ground));
assert!(roles.contains(&ProfileRole::Invert));
assert!(roles.contains(&ProfileRole::Hgl));
for pl in &d.profile_lines {
assert!(pl.pts.len() >= 2, "{:?} too short", pl.role);
}
}
}

208
crates/stormsewer/src/hydraulics.rs Normal file
View file

@ -0,0 +1,208 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Open-channel & partial-flow hydraulics for circular closed conduits.
//!
//! US customary units throughout (feet, seconds, cfs). Pass `K_MANNING_SI`
//! and `G_SI` for metric. All geometry is exact (no table lookups).
use std::f64::consts::PI;
/// Manning conversion factor, US customary (1.486 rounded to 1.49 by convention).
pub const K_MANNING_US: f64 = 1.49;
/// Manning conversion factor, SI.
pub const K_MANNING_SI: f64 = 1.0;
/// Gravitational acceleration, US customary (ft/s^2).
pub const G_US: f64 = 32.2;
/// Gravitational acceleration, SI (m/s^2).
pub const G_SI: f64 = 9.81;
/// Central angle (radians) subtended at the pipe centre by the water surface,
/// for flow depth `y` in a circular pipe of diameter `d`.
///
/// `0` when empty, `2*PI` when full.
pub fn circular_angle(y: f64, d: f64) -> f64 {
let y = y.clamp(0.0, d);
2.0 * (1.0 - 2.0 * y / d).acos()
}
/// Flow geometry at depth `y` in a circular pipe of diameter `d`.
///
/// Returns `(area, wetted_perimeter, hydraulic_radius, top_width)`.
pub fn circular_geometry(y: f64, d: f64) -> (f64, f64, f64, f64) {
let theta = circular_angle(y, d);
let area = d * d / 8.0 * (theta - theta.sin());
let perim = d * theta / 2.0;
let radius = if perim > 0.0 { area / perim } else { 0.0 };
let top = d * (theta / 2.0).sin();
(area, perim, radius, top)
}
/// Full cross-sectional area of a circular pipe.
pub fn full_area(d: f64) -> f64 {
PI * d * d / 4.0
}
/// Manning discharge from flow area and hydraulic radius.
///
/// `Q = (k/n) * A * R^(2/3) * S^(1/2)`
pub fn manning_q(n: f64, s: f64, area: f64, radius: f64, k: f64) -> f64 {
if n <= 0.0 {
return 0.0;
}
k / n * area * radius.powf(2.0 / 3.0) * s.max(0.0).sqrt()
}
/// Just-full (geometrically full, free-surface) capacity of a circular pipe.
///
/// Note: the *maximum* open-channel discharge actually occurs near `0.94 d`
/// and slightly exceeds this value — see [`max_capacity`].
pub fn full_flow_capacity(n: f64, s: f64, d: f64, k: f64) -> f64 {
manning_q(n, s, full_area(d), d / 4.0, k)
}
/// Discharge at partial depth `y` in a circular pipe.
pub fn circular_q(n: f64, s: f64, d: f64, y: f64, k: f64) -> f64 {
let (area, _p, radius, _t) = circular_geometry(y, d);
manning_q(n, s, area, radius, k)
}
/// Maximum open-channel discharge of a circular pipe and the depth at which it
/// occurs (~`0.938 d`). Returns `(q_max, y_at_max)`.
pub fn max_capacity(n: f64, s: f64, d: f64, k: f64) -> (f64, f64) {
// Q(y) is smooth with a single interior maximum; a fine scan is exact
// enough and cheap for the depths involved.
let steps = 2000;
let mut best_q = 0.0;
let mut best_y = 0.0;
for i in 1..=steps {
let y = d * i as f64 / steps as f64;
let q = circular_q(n, s, d, y, k);
if q > best_q {
best_q = q;
best_y = y;
}
}
(best_q, best_y)
}
/// Normal (uniform-flow) depth for a target discharge in a circular pipe.
///
/// Returns `None` when `q_target` exceeds the pipe's maximum open-channel
/// capacity — i.e. the pipe must surcharge / flow under pressure.
pub fn normal_depth(q_target: f64, n: f64, s: f64, d: f64, k: f64) -> Option<f64> {
if q_target <= 0.0 {
return Some(0.0);
}
let (q_max, y_max) = max_capacity(n, s, d, k);
if q_target > q_max {
return None;
}
// Q is monotonic increasing on (0, y_max]; bisect there.
let (mut lo, mut hi) = (0.0_f64, y_max);
for _ in 0..200 {
let mid = 0.5 * (lo + hi);
if circular_q(n, s, d, mid, k) > q_target {
hi = mid;
} else {
lo = mid;
}
if hi - lo < 1e-7 {
break;
}
}
Some(0.5 * (lo + hi))
}
/// Critical depth for a target discharge in a circular pipe.
///
/// Solves `Q^2 * T = g * A^3` (Froude = 1). Valid for `0 < y < d`.
pub fn critical_depth(q: f64, d: f64, g: f64) -> f64 {
if q <= 0.0 {
return 0.0;
}
let resid = |y: f64| {
let (area, _p, _r, top) = circular_geometry(y, d);
if area <= 0.0 || top <= 0.0 {
return -q * q;
}
q * q * top - g * area * area * area
};
// resid > 0 for small y, < 0 as y grows: bracket the sign change.
let (mut lo, mut hi) = (1e-6 * d, 0.999 * d);
for _ in 0..200 {
let mid = 0.5 * (lo + hi);
if resid(mid) > 0.0 {
lo = mid;
} else {
hi = mid;
}
if hi - lo < 1e-7 {
break;
}
}
0.5 * (lo + hi)
}
#[cfg(test)]
mod tests {
use super::*;
const K: f64 = K_MANNING_US;
#[test]
fn full_area_matches_geometry() {
// Geometry at y = d must equal the closed-form full area.
let d = 2.0;
let (area, perim, r, _t) = circular_geometry(d, d);
assert!((area - full_area(d)).abs() < 1e-9);
assert!((perim - PI * d).abs() < 1e-9); // full wetted perimeter = pi*d
assert!((r - d / 4.0).abs() < 1e-9); // hydraulic radius full = d/4
}
#[test]
fn half_depth_is_half_area() {
// A circle is exactly half-full at half-diameter depth.
let d = 3.0;
let (area, _p, _r, _t) = circular_geometry(d / 2.0, d);
assert!((area - full_area(d) / 2.0).abs() < 1e-9);
}
#[test]
fn manning_capacity_known_value() {
// 2-ft pipe, n=0.013, s=0.005 -> ~16.0 cfs full-flow (hand calc).
let q = full_flow_capacity(0.013, 0.005, 2.0, K);
assert!((q - 16.0).abs() < 0.2, "got {q}");
}
#[test]
fn normal_depth_round_trip() {
// Build a discharge from a depth, then recover the depth.
let (n, s, d) = (0.013, 0.01, 2.0);
let y0 = 1.2;
let q = circular_q(n, s, d, y0, K);
let y = normal_depth(q, n, s, d, K).expect("below capacity");
assert!((y - y0).abs() < 1e-3, "got {y}");
}
#[test]
fn surcharge_returns_none() {
// A tiny pipe cannot pass a huge flow as open channel.
let y = normal_depth(500.0, 0.013, 0.005, 1.0, K);
assert!(y.is_none());
}
#[test]
fn max_capacity_exceeds_full() {
// Peak open-channel flow occurs near 0.94 d and beats just-full.
let (qmax, ymax) = max_capacity(0.013, 0.005, 2.0, K);
let qfull = full_flow_capacity(0.013, 0.005, 2.0, K);
assert!(qmax > qfull);
assert!(ymax > 1.8 && ymax < 2.0, "y_at_max = {ymax}");
}
#[test]
fn critical_depth_in_range() {
let yc = critical_depth(10.0, 2.0, G_US);
assert!(yc > 0.0 && yc < 2.0, "yc = {yc}");
}
}

53
crates/stormsewer/src/idf.rs Normal file
View file

@ -0,0 +1,53 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Rainfall intensity–duration–frequency (IDF) curves.
//!
//! Uses the common three-parameter form
//!
//! ```text
//! i = a / (t + b)^c (i in in/hr, t in minutes)
//! ```
//!
//! which fits NOAA Atlas-14 partial-duration data well over the 5–180 min
//! range used in storm-sewer design.
/// A three-parameter IDF curve for a single return period.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct IdfCurve {
pub a: f64,
pub b: f64,
pub c: f64,
}
impl IdfCurve {
pub fn new(a: f64, b: f64, c: f64) -> Self {
Self { a, b, c }
}
/// Rainfall intensity (in/hr) for a storm duration `t_min` (minutes).
pub fn intensity(&self, t_min: f64) -> f64 {
let t = t_min.max(0.0);
self.a / (t + self.b).powf(self.c)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn intensity_decreases_with_duration() {
let idf = IdfCurve::new(120.0, 10.0, 0.8);
let short = idf.intensity(5.0);
let long = idf.intensity(60.0);
assert!(short > long, "short {short} long {long}");
assert!(long > 0.0);
}
#[test]
fn intensity_known_value() {
// a=120,b=10,c=0.8 ; t=15 -> 120 / 25^0.8 = 9.14 in/hr (hand calc).
let idf = IdfCurve::new(120.0, 10.0, 0.8);
assert!((idf.intensity(15.0) - 9.14).abs() < 0.05, "{}", idf.intensity(15.0));
}
}

45
crates/stormsewer/src/lib.rs Normal file
View file

@ -0,0 +1,45 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! # stormsewer
//!
//! Native-Rust storm-sewer network **hydrology & hydraulics** engine.
//!
//! It implements the standard, public-domain methods used by tools such as
//! Autodesk Hydraflow Storm Sewers:
//!
//! * **Rational method** peak-flow accumulation down a pipe network,
//! * **Manning** open-channel / partial-flow hydraulics for circular conduits,
//! * normal-depth, critical-depth and full-flow capacity,
//! * *(forthcoming)* **HEC-22** hydraulic-grade-line backwater with junction
//! and structure losses.
//!
//! This is an **engine only**: no GUI and no CAD dependencies, so it compiles
//! to a native library, to WASM (for hydrocomplete.com), and is consumable as
//! a module by an Open CAD Studio fork.
//!
//! ```
//! use stormsewer::{Network, Node, NodeKind, Pipe};
//! let net = Network {
//! nodes: vec![
//! Node::inlet("N1", 100.0, 105.0, 2.0, 0.7),
//! Node::outfall("OUT", 99.0, 104.0),
//! ],
//! pipes: vec![Pipe::new("P1", "N1", "OUT", 100.0, 1.5, 0.013)],
//! };
//! let results = net.analyze_rational(4.0).unwrap(); // i = 4 in/hr
//! assert_eq!(results.len(), 1);
//! assert!((results[0].design_q - 5.6).abs() < 1e-6); // 4 * (0.7*2.0)
//! ```
pub mod drawing;
pub mod hydraulics;
pub mod idf;
pub mod network;
pub mod parse;
pub mod report;
pub use drawing::*;
pub use hydraulics::*;
pub use idf::*;
pub use network::*;
pub use parse::*;

527
crates/stormsewer/src/network.rs Normal file
View file

@ -0,0 +1,527 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Storm-sewer network model: nodes (inlets / junctions / outfalls), pipes,
//! topology, Rational-method peak-flow accumulation, time-of-concentration /
//! IDF intensity, and a hydraulic-grade-line (HGL) backwater pass.
//!
//! The network is assumed **dendritic** (tree-like, draining to outfalls) as is
//! standard for gravity storm sewers. Looped networks are rejected by the
//! topological sort.
use crate::hydraulics::*;
use crate::idf::IdfCurve;
use std::collections::HashMap;
/// Kind of network node.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum NodeKind {
Inlet,
Junction,
Outfall,
}
/// A network node (structure).
#[derive(Clone, Debug)]
pub struct Node {
pub id: String,
pub kind: NodeKind,
/// Flowline / invert elevation (ft).
pub invert: f64,
/// Rim / ground elevation (ft) — used for surcharge-to-surface checks.
pub rim: f64,
/// Local contributing drainage area (acres).
pub area_ac: f64,
/// Rational runoff coefficient C for the local area (0–1).
pub c: f64,
/// Inlet time of concentration for the local catchment (minutes).
pub tc_inlet: f64,
/// Plan (map) X coordinate of the structure (drawing units).
pub x: f64,
/// Plan (map) Y coordinate of the structure (drawing units).
pub y: f64,
}
impl Node {
pub fn inlet(id: &str, invert: f64, rim: f64, area_ac: f64, c: f64) -> Self {
Self { id: id.into(), kind: NodeKind::Inlet, invert, rim, area_ac, c, tc_inlet: 10.0, x: 0.0, y: 0.0 }
}
pub fn junction(id: &str, invert: f64, rim: f64, area_ac: f64, c: f64) -> Self {
Self { id: id.into(), kind: NodeKind::Junction, invert, rim, area_ac, c, tc_inlet: 10.0, x: 0.0, y: 0.0 }
}
pub fn outfall(id: &str, invert: f64, rim: f64) -> Self {
Self { id: id.into(), kind: NodeKind::Outfall, invert, rim, area_ac: 0.0, c: 0.0, tc_inlet: 0.0, x: 0.0, y: 0.0 }
}
/// Builder: set the local inlet time (minutes).
pub fn with_tc_inlet(mut self, t_min: f64) -> Self {
self.tc_inlet = t_min;
self
}
/// Builder: set the plan (map) coordinates.
pub fn at(mut self, x: f64, y: f64) -> Self {
self.x = x;
self.y = y;
self
}
/// Local C*A product (acres).
pub fn ca(&self) -> f64 {
self.c * self.area_ac
}
}
/// A pipe (network link) carrying flow from `from` (upstream) to `to`.
#[derive(Clone, Debug)]
pub struct Pipe {
pub id: String,
pub from: String,
pub to: String,
pub length: f64, // ft
pub diameter: f64, // ft
pub n: f64, // Manning roughness
}
impl Pipe {
pub fn new(id: &str, from: &str, to: &str, length: f64, diameter: f64, n: f64) -> Self {
Self { id: id.into(), from: from.into(), to: to.into(), length, diameter, n }
}
}
/// A storm-sewer network.
#[derive(Clone, Debug, Default)]
pub struct Network {
pub nodes: Vec<Node>,
pub pipes: Vec<Pipe>,
}
/// Per-pipe analysis result.
#[derive(Clone, Debug)]
pub struct PipeResult {
pub id: String,
pub from: String,
pub to: String,
/// Pipe slope from upstream/downstream inverts (ft/ft).
pub slope: f64,
/// Accumulated C*A draining into this pipe (acres).
pub total_ca: f64,
/// Design rainfall intensity used (in/hr).
pub intensity: f64,
/// Time of concentration at the upstream end (minutes).
pub tc: f64,
/// Flow travel time through this pipe at design flow (minutes).
pub travel_time: f64,
/// Rational peak design discharge (cfs).
pub design_q: f64,
/// Just-full capacity (cfs).
pub capacity: f64,
/// Maximum open-channel capacity (~0.94 d) (cfs).
pub max_capacity: f64,
/// True if the design flow exceeds open-channel capacity (pipe surcharges).
pub surcharged: bool,
/// Normal depth (ft), or `None` if surcharged.
pub normal_depth: Option<f64>,
/// Critical depth (ft).
pub critical_depth: f64,
/// Actual velocity at design flow (ft/s).
pub velocity: f64,
/// Full-flow velocity (ft/s).
pub velocity_full: f64,
/// Design flow as a fraction of just-full capacity.
pub pct_full: f64,
/// HGL elevation at the upstream structure (ft), if an HGL pass was run.
pub hgl_up: Option<f64>,
/// HGL elevation at the downstream structure (ft), if an HGL pass was run.
pub hgl_dn: Option<f64>,
}
/// Per-node analysis result (populated by the HGL pass).
#[derive(Clone, Debug)]
pub struct NodeResult {
pub id: String,
pub tc: f64,
pub hgl: f64,
pub rim: f64,
/// True if the HGL rises above the rim (flooding).
pub surcharge_to_surface: bool,
}
/// Full network analysis (pipes + nodes).
#[derive(Clone, Debug)]
pub struct Analysis {
pub pipes: Vec<PipeResult>,
pub nodes: Vec<NodeResult>,
}
/// Options for [`Network::analyze`].
#[derive(Clone, Debug)]
pub struct AnalysisOptions {
/// Minimum time of concentration (minutes) — floors the IDF duration.
pub min_tc: f64,
/// Tailwater elevation at outfalls (ft). `None` → free outfall at invert.
pub tailwater: Option<f64>,
/// Junction/structure loss coefficient K in `H = K * V^2 / 2g`.
pub junction_k: f64,
/// If set, use this constant intensity (in/hr) instead of the IDF curve.
pub intensity_override: Option<f64>,
}
impl Default for AnalysisOptions {
fn default() -> Self {
Self { min_tc: 10.0, tailwater: None, junction_k: 0.5, intensity_override: None }
}
}
/// Errors raised during analysis.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum NetworkError {
UnknownNode(String),
CyclicNetwork,
}
impl std::fmt::Display for NetworkError {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
NetworkError::UnknownNode(id) => write!(f, "pipe references unknown node '{id}'"),
NetworkError::CyclicNetwork => write!(f, "network contains a cycle"),
}
}
}
impl std::error::Error for NetworkError {}
/// Conveyance K = (k/n) * A * R^(2/3), so that Q = K * sqrt(S).
fn conveyance(n: f64, area: f64, radius: f64, k: f64) -> f64 {
if n <= 0.0 {
return 0.0;
}
k / n * area * radius.powf(2.0 / 3.0)
}
impl Network {
fn index(&self) -> HashMap<&str, usize> {
self.nodes.iter().enumerate().map(|(i, nd)| (nd.id.as_str(), i)).collect()
}
/// Topological node order, most-upstream first (Kahn's algorithm).
fn topo_order(&self, idx: &HashMap<&str, usize>) -> Result<Vec<usize>, NetworkError> {
let n = self.nodes.len();
let mut indeg = vec![0usize; n];
let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
for p in &self.pipes {
let u = *idx.get(p.from.as_str()).ok_or_else(|| NetworkError::UnknownNode(p.from.clone()))?;
let v = *idx.get(p.to.as_str()).ok_or_else(|| NetworkError::UnknownNode(p.to.clone()))?;
adj[u].push(v);
indeg[v] += 1;
}
let mut queue: Vec<usize> = (0..n).filter(|&i| indeg[i] == 0).collect();
let mut order = Vec::with_capacity(n);
while let Some(u) = queue.pop() {
order.push(u);
for &v in &adj[u] {
indeg[v] -= 1;
if indeg[v] == 0 {
queue.push(v);
}
}
}
if order.len() != n {
return Err(NetworkError::CyclicNetwork);
}
Ok(order)
}
/// Accumulated C*A (acres) reaching each node, keyed by node id.
pub fn accumulate_ca(&self) -> Result<HashMap<String, f64>, NetworkError> {
let idx = self.index();
let order = self.topo_order(&idx)?;
let mut feeders: Vec<Vec<usize>> = vec![Vec::new(); self.nodes.len()];
for p in &self.pipes {
feeders[idx[p.to.as_str()]].push(idx[p.from.as_str()]);
}
let mut total = vec![0.0f64; self.nodes.len()];
for &i in &order {
let mut acc = self.nodes[i].ca();
for &u in &feeders[i] {
acc += total[u];
}
total[i] = acc;
}
Ok(self.nodes.iter().enumerate().map(|(i, nd)| (nd.id.clone(), total[i])).collect())
}
/// Simple Rational analysis with one constant intensity (no Tc/HGL).
/// Retained for quick checks; [`analyze`](Self::analyze) is the full method.
pub fn analyze_rational(&self, intensity: f64) -> Result<Vec<PipeResult>, NetworkError> {
let opts = AnalysisOptions { intensity_override: Some(intensity), ..Default::default() };
// A flat IDF is unused when intensity_override is set.
Ok(self.analyze(&IdfCurve::new(0.0, 1.0, 1.0), &opts)?.pipes)
}
/// Full analysis: Tc accumulation, per-pipe IDF intensity, Rational design
/// flows, pipe hydraulics, and an HGL backwater pass with junction losses.
pub fn analyze(&self, idf: &IdfCurve, opts: &AnalysisOptions) -> Result<Analysis, NetworkError> {
let k = K_MANNING_US;
let idxm = self.index();
let order = self.topo_order(&idxm)?;
let n_nodes = self.nodes.len();
let n_pipes = self.pipes.len();
// Resolve pipe endpoints to node indices once.
let pe: Vec<(usize, usize)> = self
.pipes
.iter()
.map(|p| (idxm[p.from.as_str()], idxm[p.to.as_str()]))
.collect();
// incoming[node] = (pipe_idx, upstream_node_idx)
let mut incoming: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n_nodes];
// outgoing[node] = (pipe_idx, downstream_node_idx)
let mut outgoing: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n_nodes];
for (pi, &(u, v)) in pe.iter().enumerate() {
outgoing[u].push((pi, v));
incoming[v].push((pi, u));
}
let ca = self.accumulate_ca()?;
let total_ca: Vec<f64> = self.nodes.iter().map(|nd| ca[&nd.id]).collect();
// Per-pipe scratch.
let mut p_slope = vec![0.0f64; n_pipes];
let mut p_intensity = vec![0.0f64; n_pipes];
let mut p_q = vec![0.0f64; n_pipes];
let mut p_vel = vec![0.0f64; n_pipes];
let mut p_surch = vec![false; n_pipes];
let mut p_yn: Vec<Option<f64>> = vec![None; n_pipes];
let mut p_travel = vec![0.0f64; n_pipes];
let mut tc_node = vec![0.0f64; n_nodes];
// ── Forward pass (upstream → downstream): Tc, intensity, design Q, vel.
for &i in &order {
let nd = &self.nodes[i];
let mut tc = if nd.kind == NodeKind::Outfall { 0.0 } else { nd.tc_inlet };
for &(pi, _u) in &incoming[i] {
tc = tc.max(tc_node[/*upstream*/ pe[pi].0] + p_travel[pi]);
}
tc = tc.max(opts.min_tc);
tc_node[i] = tc;
for &(pi, v) in &outgoing[i] {
let p = &self.pipes[pi];
let slope = if p.length > 0.0 {
(self.nodes[i].invert - self.nodes[v].invert) / p.length
} else {
0.0
};
let intensity = opts.intensity_override.unwrap_or_else(|| idf.intensity(tc));
let q = intensity * total_ca[i];
let (q_max, _) = max_capacity(p.n, slope, p.diameter, k);
let surcharged = q > q_max;
let yn = normal_depth(q, p.n, slope, p.diameter, k);
let area = if surcharged {
full_area(p.diameter)
} else {
circular_geometry(yn.unwrap_or(p.diameter), p.diameter).0
};
let vel = if area > 0.0 { q / area } else { 0.0 };
let travel = if vel > 0.0 { p.length / vel / 60.0 } else { 0.0 };
p_slope[pi] = slope;
p_intensity[pi] = intensity;
p_q[pi] = q;
p_vel[pi] = vel;
p_surch[pi] = surcharged;
p_yn[pi] = yn;
p_travel[pi] = travel;
}
}
// ── Backward pass (downstream → upstream): HGL with junction losses.
let mut hgl = vec![f64::NAN; n_nodes];
let mut pipe_hgl_up = vec![None; n_pipes];
let mut pipe_hgl_dn = vec![None; n_pipes];
for &d in order.iter().rev() {
// Seed outfall HGL from tailwater (or free outfall at invert).
if self.nodes[d].kind == NodeKind::Outfall && hgl[d].is_nan() {
hgl[d] = opts.tailwater.unwrap_or(self.nodes[d].invert);
}
if hgl[d].is_nan() {
hgl[d] = self.nodes[d].invert; // disconnected fallback
}
for &(pi, u) in &incoming[d] {
let p = &self.pipes[pi];
let dia = p.diameter;
let inv_d = self.nodes[d].invert;
let inv_u = self.nodes[u].invert;
let q = p_q[pi];
let (ws_d, hf) = if p_surch[pi] {
let conv_full = conveyance(p.n, full_area(dia), dia / 4.0, k);
let sf = if conv_full > 0.0 { (q / conv_full).powi(2) } else { 0.0 };
let crown_d = inv_d + dia;
(hgl[d].max(crown_d), sf * p.length)
} else {
let yn = p_yn[pi].unwrap_or(0.0);
let (a, _pp, r, _t) = circular_geometry(yn, dia);
let conv = conveyance(p.n, a, r, k);
let sf = if conv > 0.0 { (q / conv).powi(2) } else { p_slope[pi].max(0.0) };
(hgl[d].max(inv_d + yn), sf * p.length)
};
let yn_u = p_yn[pi].unwrap_or(0.0);
let hgl_us_pipe = (ws_d + hf).max(inv_u + yn_u);
let hj = opts.junction_k * p_vel[pi].powi(2) / (2.0 * G_US);
let hgl_u = hgl_us_pipe + hj;
hgl[u] = if hgl[u].is_nan() { hgl_u } else { hgl[u].max(hgl_u) };
pipe_hgl_dn[pi] = Some(hgl[d]);
pipe_hgl_up[pi] = Some(hgl_u);
}
}
// ── Assemble results.
let pipes = self
.pipes
.iter()
.enumerate()
.map(|(pi, p)| {
let capacity = full_flow_capacity(p.n, p_slope[pi], p.diameter, k);
let velocity_full = if full_area(p.diameter) > 0.0 {
capacity / full_area(p.diameter)
} else {
0.0
};
PipeResult {
id: p.id.clone(),
from: p.from.clone(),
to: p.to.clone(),
slope: p_slope[pi],
total_ca: total_ca[pe[pi].0],
intensity: p_intensity[pi],
tc: tc_node[pe[pi].0],
travel_time: p_travel[pi],
design_q: p_q[pi],
capacity,
max_capacity: max_capacity(p.n, p_slope[pi], p.diameter, k).0,
surcharged: p_surch[pi],
normal_depth: p_yn[pi],
critical_depth: critical_depth(p_q[pi], p.diameter, G_US),
velocity: p_vel[pi],
velocity_full,
pct_full: if capacity > 0.0 { p_q[pi] / capacity } else { 0.0 },
hgl_up: pipe_hgl_up[pi],
hgl_dn: pipe_hgl_dn[pi],
}
})
.collect();
let nodes = self
.nodes
.iter()
.enumerate()
.map(|(i, nd)| NodeResult {
id: nd.id.clone(),
tc: tc_node[i],
hgl: hgl[i],
rim: nd.rim,
surcharge_to_surface: nd.kind != NodeKind::Outfall && hgl[i] > nd.rim,
})
.collect();
Ok(Analysis { pipes, nodes })
}
}
#[cfg(test)]
mod tests {
use super::*;
/// N1 (inlet) -> N2 (inlet) -> OUT, two pipes.
fn sample() -> Network {
Network {
nodes: vec![
Node::inlet("N1", 100.0, 105.0, 2.0, 0.7), // CA = 1.4
Node::inlet("N2", 99.0, 104.0, 3.0, 0.8), // CA = 2.4
Node::outfall("OUT", 98.0, 103.0),
],
pipes: vec![
Pipe::new("P1", "N1", "N2", 100.0, 1.5, 0.013),
Pipe::new("P2", "N2", "OUT", 100.0, 1.5, 0.013),
],
}
}
#[test]
fn ca_accumulates_downstream() {
let ca = sample().accumulate_ca().unwrap();
assert!((ca["N1"] - 1.4).abs() < 1e-9);
assert!((ca["N2"] - 3.8).abs() < 1e-9);
assert!((ca["OUT"] - 3.8).abs() < 1e-9);
}
#[test]
fn rational_design_flows() {
let r = sample().analyze_rational(4.0).unwrap();
let p1 = r.iter().find(|x| x.id == "P1").unwrap();
let p2 = r.iter().find(|x| x.id == "P2").unwrap();
assert!((p1.design_q - 5.6).abs() < 1e-6, "P1 {}", p1.design_q);
assert!((p2.design_q - 15.2).abs() < 1e-6, "P2 {}", p2.design_q);
}
#[test]
fn slope_from_inverts() {
let r = sample().analyze_rational(4.0).unwrap();
assert!((r[0].slope - 0.01).abs() < 1e-9);
}
#[test]
fn small_pipe_surcharges_under_heavy_flow() {
let r = sample().analyze_rational(4.0).unwrap();
let p1 = r.iter().find(|x| x.id == "P1").unwrap();
let p2 = r.iter().find(|x| x.id == "P2").unwrap();
assert!(!p1.surcharged, "P1 should fit");
assert!(p1.normal_depth.is_some());
assert!(p2.surcharged, "P2 should surcharge");
assert!(p2.normal_depth.is_none());
}
#[test]
fn cycle_is_rejected() {
let net = Network {
nodes: vec![Node::junction("A", 10.0, 20.0, 0.0, 0.0), Node::junction("B", 9.0, 19.0, 0.0, 0.0)],
pipes: vec![Pipe::new("AB", "A", "B", 50.0, 1.0, 0.013), Pipe::new("BA", "B", "A", 50.0, 1.0, 0.013)],
};
assert_eq!(net.accumulate_ca().unwrap_err(), NetworkError::CyclicNetwork);
}
#[test]
fn unknown_node_is_rejected() {
let net = Network {
nodes: vec![Node::outfall("OUT", 98.0, 103.0)],
pipes: vec![Pipe::new("P", "GHOST", "OUT", 100.0, 1.5, 0.013)],
};
assert!(matches!(net.accumulate_ca().unwrap_err(), NetworkError::UnknownNode(_)));
}
#[test]
fn tc_accumulates_downstream() {
// Downstream pipe's Tc must include upstream inlet time + travel time.
let idf = IdfCurve::new(120.0, 10.0, 0.8);
let a = sample().analyze(&idf, &AnalysisOptions::default()).unwrap();
let p1 = a.pipes.iter().find(|x| x.id == "P1").unwrap();
let p2 = a.pipes.iter().find(|x| x.id == "P2").unwrap();
assert!(p2.tc >= p1.tc, "p2.tc {} p1.tc {}", p2.tc, p1.tc);
// Intensity falls as Tc grows downstream.
assert!(p2.intensity <= p1.intensity, "i2 {} i1 {}", p2.intensity, p1.intensity);
}
#[test]
fn hgl_rises_upstream() {
// With a tailwater, HGL must be monotonically higher going upstream.
let idf = IdfCurve::new(120.0, 10.0, 0.8);
let opts = AnalysisOptions { tailwater: Some(100.0), ..Default::default() };
let a = sample().analyze(&idf, &opts).unwrap();
let h = |id: &str| a.nodes.iter().find(|n| n.id == id).unwrap().hgl;
assert!(h("OUT") <= h("N2"), "OUT {} N2 {}", h("OUT"), h("N2"));
assert!(h("N2") <= h("N1"), "N2 {} N1 {}", h("N2"), h("N1"));
assert!((h("OUT") - 100.0).abs() < 1e-9, "tailwater seeded");
}
}

122
crates/stormsewer/src/parse.rs Normal file
View file

@ -0,0 +1,122 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Parser for the `.ssn` storm-sewer network text format.
//!
//! Whitespace-delimited; `#` starts a comment. Keywords:
//!
//! ```text
//! IDF <a> <b> <c> # IDF curve i = a/(t+b)^c
//! INTENSITY <in/hr> # OR a constant intensity (overrides IDF)
//! TAILWATER <elev_ft> # outfall tailwater elevation
//! MINTC <min> # minimum time of concentration
//! JUNCTIONK <k> # junction loss coefficient
//! NODE <id> <inlet|junction|outfall> <x> <y> <invert> <rim> [area] [C] [tc_inlet]
//! PIPE <id> <from> <to> <length> <dia> <n>
//! ```
//!
//! For `outfall` nodes the trailing `area`, `C`, and `tc_inlet` are omitted.
use crate::idf::IdfCurve;
use crate::network::{AnalysisOptions, Network, Node, Pipe};
/// Result of parsing a `.ssn` document.
#[derive(Clone, Debug)]
pub struct ParsedNetwork {
pub network: Network,
pub idf: IdfCurve,
pub options: AnalysisOptions,
}
/// Parse a `.ssn` document. Returns a human-readable error with the line number.
pub fn parse_ssn(text: &str) -> Result<ParsedNetwork, String> {
let mut network = Network::default();
let mut idf = IdfCurve::new(120.0, 10.0, 0.8);
let mut options = AnalysisOptions::default();
for (lineno, raw) in text.lines().enumerate() {
let line = raw.split('#').next().unwrap_or("").trim();
if line.is_empty() {
continue;
}
let t: Vec<&str> = line.split_whitespace().collect();
let ln = lineno + 1;
let num = |i: usize| -> Result<f64, String> {
t.get(i)
.ok_or_else(|| format!("line {ln}: missing value #{} in `{line}`", i + 1))?
.parse::<f64>()
.map_err(|_| format!("line {ln}: `{}` is not a number", t[i]))
};
match t[0].to_ascii_uppercase().as_str() {
"IDF" => idf = IdfCurve::new(num(1)?, num(2)?, num(3)?),
"INTENSITY" => options.intensity_override = Some(num(1)?),
"TAILWATER" => options.tailwater = Some(num(1)?),
"MINTC" => options.min_tc = num(1)?,
"JUNCTIONK" => options.junction_k = num(1)?,
"NODE" => {
let id = *t.get(1).ok_or_else(|| format!("line {ln}: NODE needs an id"))?;
let kind = t.get(2).ok_or_else(|| format!("line {ln}: NODE needs a kind"))?.to_ascii_lowercase();
let (x, y, invert, rim) = (num(3)?, num(4)?, num(5)?, num(6)?);
let node = match kind.as_str() {
"outfall" => Node::outfall(id, invert, rim),
"junction" => Node::junction(id, invert, rim, num(7)?, num(8)?),
"inlet" => Node::inlet(id, invert, rim, num(7)?, num(8)?),
other => return Err(format!("line {ln}: unknown node kind `{other}`")),
}
.at(x, y);
let node = if t.len() > 9 { node.with_tc_inlet(num(9)?) } else { node };
network.nodes.push(node);
}
"PIPE" => {
let id = *t.get(1).ok_or_else(|| format!("line {ln}: PIPE needs an id"))?;
let from = *t.get(2).ok_or_else(|| format!("line {ln}: PIPE needs a from-node"))?;
let to = *t.get(3).ok_or_else(|| format!("line {ln}: PIPE needs a to-node"))?;
network.pipes.push(Pipe::new(id, from, to, num(4)?, num(5)?, num(6)?));
}
other => return Err(format!("line {ln}: unknown keyword `{other}`")),
}
}
Ok(ParsedNetwork { network, idf, options })
}
#[cfg(test)]
mod tests {
use super::*;
const DOC: &str = "\
# demo
IDF 60 10 0.8
TAILWATER 100.5
NODE N1 inlet 0 0 104 110 1.0 0.70 12
NODE N2 inlet 300 0 102 108 1.0 0.70
NODE OUT outfall 600 0 100 106
PIPE P1 N1 N2 300 1.25 0.013
PIPE P2 N2 OUT 300 1.50 0.013
";
#[test]
fn parses_nodes_pipes_and_params() {
let p = parse_ssn(DOC).unwrap();
assert_eq!(p.network.nodes.len(), 3);
assert_eq!(p.network.pipes.len(), 2);
assert_eq!(p.options.tailwater, Some(100.5));
let n2 = p.network.nodes.iter().find(|n| n.id == "N2").unwrap();
assert_eq!((n2.x, n2.y), (300.0, 0.0));
assert_eq!(n2.tc_inlet, 10.0); // defaulted (no trailing value)
let n1 = p.network.nodes.iter().find(|n| n.id == "N1").unwrap();
assert_eq!(n1.tc_inlet, 12.0);
}
#[test]
fn analyzes_after_parse() {
let p = parse_ssn(DOC).unwrap();
let a = p.network.analyze(&p.idf, &p.options).unwrap();
assert_eq!(a.pipes.len(), 2);
}
#[test]
fn reports_bad_line() {
let err = parse_ssn("NODE N1 inlet 0 0 oops 110 1 0.7").unwrap_err();
assert!(err.contains("line 1"), "{err}");
}
}

93
crates/stormsewer/src/report.rs Normal file
View file

@ -0,0 +1,93 @@
// SPDX-License-Identifier: GPL-3.0-or-later
//! Plain-text report tables for an [`Analysis`], in the spirit of Hydraflow
//! Storm Sewers' pipe and HGL summaries.
use crate::network::Analysis;
fn f(x: f64, w: usize, p: usize) -> String {
format!("{:>w$.p$}", x, w = w, p = p)
}
/// Format the pipe hydraulics table.
pub fn pipe_table(a: &Analysis) -> String {
let mut s = String::new();
s.push_str(
"Pipe From To Slope Tc i Q Cap %Full V yn HGLup HGLdn Status\n",
);
s.push_str(
" ft/ft min in/hr cfs cfs % ft/s ft ft ft\n",
);
s.push_str(&"-".repeat(110));
s.push('\n');
for p in &a.pipes {
let yn = p.normal_depth.map(|y| f(y, 6, 2)).unwrap_or_else(|| " full".into());
let hup = p.hgl_up.map(|h| f(h, 8, 2)).unwrap_or_else(|| " --".into());
let hdn = p.hgl_dn.map(|h| f(h, 8, 2)).unwrap_or_else(|| " --".into());
let status = if p.surcharged { "SURCHARGED" } else { "ok" };
s.push_str(&format!(
"{:<6} {:<6} {:<6} {} {} {} {} {} {} {} {} {} {} {}\n",
p.id,
p.from,
p.to,
f(p.slope, 7, 4),
f(p.tc, 5, 1),
f(p.intensity, 6, 2),
f(p.design_q, 6, 2),
f(p.capacity, 6, 2),
f(p.pct_full * 100.0, 6, 1),
f(p.velocity, 6, 2),
yn,
hup,
hdn,
status,
));
}
s
}
/// Format the node / HGL table.
pub fn node_table(a: &Analysis) -> String {
let mut s = String::new();
s.push_str("Node Tc(min) Rim(ft) HGL(ft) Freeboard Status\n");
s.push_str(&"-".repeat(60));
s.push('\n');
for n in &a.nodes {
let fb = n.rim - n.hgl;
let status = if n.surcharge_to_surface { "FLOODING" } else { "ok" };
s.push_str(&format!(
"{:<6} {} {} {} {} {}\n",
n.id,
f(n.tc, 7, 1),
f(n.rim, 8, 2),
f(n.hgl, 8, 2),
f(fb, 9, 2),
status,
));
}
s
}
/// Full report: pipe table followed by node/HGL table.
pub fn format_analysis(a: &Analysis) -> String {
let mut s = String::new();
s.push_str("=== STORM SEWER ANALYSIS ===\n\n");
s.push_str(&pipe_table(a));
s.push('\n');
s.push_str(&node_table(a));
// Summary flags.
let surcharged: Vec<&str> = a.pipes.iter().filter(|p| p.surcharged).map(|p| p.id.as_str()).collect();
let flooding: Vec<&str> = a.nodes.iter().filter(|n| n.surcharge_to_surface).map(|n| n.id.as_str()).collect();
s.push('\n');
if surcharged.is_empty() && flooding.is_empty() {
s.push_str("All pipes flow open-channel; no surface flooding.\n");
} else {
if !surcharged.is_empty() {
s.push_str(&format!("Surcharged pipes: {}\n", surcharged.join(", ")));
}
if !flooding.is_empty() {
s.push_str(&format!("Structures flooding (HGL > rim): {}\n", flooding.join(", ")));
}
}
s
}