From 05f9630cb498b52a4815486088361ab2ebdfc300 Mon Sep 17 00:00:00 2001 From: Jeroen Schweitzer Date: Thu, 23 Jul 2026 19:22:28 +0200 Subject: [PATCH] feat(simulation): equilibrium hydrology solver prototype + bench (T-1177) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Priority-flood fill with basin grouping, topographic-saddle spill points, moisture-governed endorheic classification (reuses the reserved RIVER_DOWNSTREAM_TERMINAL sentinel), and Dijkstra overflow/carving. Pure function of (elevation, sea_level, climate) — deterministic per D-010, no stateful simulation. 15 unit tests incl. determinism proofs and direct carving-mechanism verification; #[ignore]d release benches at 512x256 / 768x432 / 3840x2160 plus the 273-bodies-parallel production shape. Workshop gate measurement (1) for body-map-viewer: settled hydrology is VIABLE per body-open (~24 ms at 512x256; ~0.7-0.8 s all 273 bodies). Co-Authored-By: Claude Fable 5 --- server/src/atlas/hydrology_equilibrium.rs | 1041 +++++++++++++++++++ server/src/atlas/mod.rs | 1 + server/tests/hydrology_equilibrium_bench.rs | 241 +++++ 3 files changed, 1283 insertions(+) create mode 100644 server/src/atlas/hydrology_equilibrium.rs create mode 100644 server/tests/hydrology_equilibrium_bench.rs diff --git a/server/src/atlas/hydrology_equilibrium.rs b/server/src/atlas/hydrology_equilibrium.rs new file mode 100644 index 000000000..56f22f7c4 --- /dev/null +++ b/server/src/atlas/hydrology_equilibrium.rs @@ -0,0 +1,1041 @@ +//! Deterministic equilibrium hydrology — lake fill, overflow re-routing, and +//! gorge carving (T-1177, body-map-viewer workshop measurement ①). +//! +//! **What this is.** A pure function of `(seed, heightmap, moisture/climate)` +//! that computes a *settled end-state* for standing water: every depression +//! (basin) on the elevation field either (a) fills to a spill level and +//! overflows onward toward the sea — possibly chaining through further +//! basins — or (b) is classified **endorheic** (holds its own closed water +//! cycle, no outflow) when the body's moisture/climate makes evaporative +//! balance plausible. Where the cheapest overflow path must cross a ridge +//! higher than the lake's own spill level, the solver carves a deterministic +//! channel (gorge) — geological backstory computed at generation time, not a +//! stateful erosion simulation. +//! +//! **This is NOT a simulation.** There is no tick loop, no iterative erosion, +//! no rainfall-accumulation-over-time model. One elevation field in, one +//! settled hydrology result out, same every time for the same inputs +//! (D-010). See the "what this prototype does NOT do" list in +//! `docs/workshops/body-map-viewer/measurements/t1177-hydrology.md` for the +//! full list of deliberately deferred behaviour. +//! +//! **Algorithm (priority-flood family, Barnes/Planchon-Darboux class):** +//! 1. Priority-flood fill from the boundary/sea inward: a min-heap seeded at +//! every below-`sea_level` cell floods the grid, raising each interior +//! cell to the highest of (its own elevation, the water level that reached +//! it) — the standard "flood from the outside in" formulation. This +//! produces, for every cell, its *filled* elevation (drainage.rs's +//! existing `depression_fill` is an iterative relaxation of the same +//! problem but does not expose spill points or basin membership — this +//! module needs both, so it re-derives the fill with a heap that *does* +//! track them, rather than layering basin-detection on top of the +//! existing black-box result). +//! 2. Every interior cell where `filled > original` is a *lake cell*; the +//! heap visit order also yields each basin's spill point (the last +//! boundary cell that let water into it) essentially for free. +//! 3. Basins are grouped by contiguous lake-cell connectivity (flood-fill). +//! Each basin's spill level = the filled elevation at its spill point. +//! 4. Endorheic classification (per-basin, moisture-governed — see +//! `is_endorheic`): large + dry basins are declared closed-cycle sinks +//! (TERMINAL). Everything else overflows: a cheapest-path search +//! (Dijkstra/A*-family over elevation-crossing cost, same `BinaryHeap` +//! pattern as `road_graph.rs`'s `astar`) from the spill point to the +//! nearest lower-or-equal terrain (another basin's interior or the sea) +//! finds the overflow channel. Cells on that channel above the spill +//! level are carved down to it — the gorge. +//! 5. Chaining: an overflowing basin's outlet may land inside *another* +//! basin's footprint; that basin is processed in the same pass (basins +//! are handled in ascending spill-level order, so a lower basin is always +//! resolved before anything can overflow into it a second time). +//! +//! Determinism (D-010): all comparisons are on the same `i64`-scaled +//! elevation integers `drainage.rs` uses; heap tie-breaks are `(cost, cell +//! index)` so equal-cost frontier cells always resolve in the same order; +//! every collection that participates in output ordering is `Vec`/`BTreeMap` +//! keyed by cell index or basin id, never a `HashMap`/`HashSet` iteration. + +use std::cmp::Reverse; +use std::collections::{BinaryHeap, VecDeque}; + +use crate::atlas::drainage::d8_offset; + +/// Scale factor for converting f32 elevation to integer for deterministic +/// comparison — matches `drainage::ELEV_SCALE` (kept as a private duplicate +/// rather than exported from `drainage.rs`, since the two modules' integer +/// domains are independent: this module's `i64` values are always genuine +/// elevations in the same units, never `drainage.rs`'s post-relaxation +/// `depression_fill` values). +const ELEV_SCALE: f64 = 1_000_000.0; + +/// A cell is a river cell when its flow accumulation exceeds this threshold — +/// mirrors `drainage::RIVER_THRESHOLD`; only overflow channels that also +/// qualify as "a real waterway" (i.e. would show up on the river network) +/// are carved as gorges. A trickling overflow with negligible flow doesn't +/// get a canyon. +const RIVER_THRESHOLD: i32 = 200; + +/// D8 neighbor offsets, same fixed priority order as `drainage::D8` — needed +/// locally so `d8_offset` calls stay index-compatible with the shared table +/// without re-exporting the private array. +const D8_LEN: u8 = 8; + +/// One basin's settled hydrology outcome. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum BasinOutcome { + /// The basin overflows: `spill_cell` is where water leaves the basin, + /// `outlet_path` is the carved channel (spill cell → downstream target, + /// inclusive of both ends), and `downstream_target` is either another + /// basin (chained) or the sea. + Overflow { + outlet_path: Vec, + downstream_target: DownstreamTarget, + }, + /// The basin holds a closed water cycle — no outflow (TERMINAL). + Endorheic { reason: EndorheicReason }, +} + +/// What an overflowing basin's outlet channel terminates at. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum DownstreamTarget { + /// Reaches a sub-`sea_level` cell. + Sea, + /// Flows into another basin's footprint (chained overflow), by basin id. + Basin(u32), + /// No lower terrain reachable within the search budget — treated as a + /// grid-edge drain, the lake-fill analogue of `RIVER_DOWNSTREAM_EDGE_DRAIN`. + EdgeUnreachable, +} + +/// Why a basin was classified endorheic — carried for diagnostics/tuning, +/// not currently read by any downstream consumer. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum EndorheicReason { + /// Basin area + dryness cleared the moisture-governed threshold (see + /// `is_endorheic`). + MoistureGoverned, +} + +/// One filled basin (lake) in the settled hydrology result. +#[derive(Debug, Clone)] +pub struct Basin { + pub basin_id: u32, + /// Row-major cell indices belonging to this basin's lake footprint + /// (cells where `filled_elevation > original_elevation`), sorted + /// ascending for determinism. + pub cells: Vec, + /// The basin's spill level, in the same i64-scaled units as + /// `drainage::analyze`'s internal representation (divide by + /// `ELEV_SCALE` to recover the [0,1] float). + pub spill_level_scaled: i64, + /// The cell that is the basin's lowest boundary point — the pour point + /// water first breaches when the basin overtops. + pub spill_cell: usize, + pub outcome: BasinOutcome, +} + +/// Full result of the equilibrium hydrology solve for one body. +#[derive(Debug, Clone)] +pub struct HydrologyResult { + pub basins: Vec, + /// Per-cell filled elevation (i64-scaled) — the settled water surface + /// where a lake exists, or the original terrain elevation otherwise. + /// Row-major, `w × h`. + pub filled_scaled: Vec, + /// Per-cell channel depth in the same i64-scaled elevation units: 0 + /// everywhere except carved gorge cells, where it is + /// `original_scaled[cell] - carved_floor_scaled[cell]` — see the module + /// docs' cliff-representation section. Row-major, `w × h`. + pub channel_depth_scaled: Vec, + /// True for every cell that is part of a carved gorge channel (a cliff + /// edge — the rim/floor discontinuity the workshop's red flag 4 asks + /// about). Row-major, `w × h`. + pub cliff_edge: Vec, +} + +/// Moisture/climate inputs governing the endorheic-vs-overflow decision. +/// Deliberately minimal and explicitly tunable — see the module docs and the +/// results doc's "endorheic criterion" section for the rationale and the +/// knobs future tuning is expected to touch. +#[derive(Debug, Clone, Copy)] +pub struct ClimateInputs { + /// Body-level moisture ceiling, 0–100 (mirrors + /// `district_profile::derive_moisture_q`'s hydrosphere-derived `ceiling` + /// before the per-district spatial gradient — a single body-wide value + /// is enough for a first-pass basin classifier; a per-basin-centroid + /// value is a documented future refinement, not built here). + pub moisture_q: i32, +} + +/// Solve equilibrium hydrology for one elevation grid. `elevation` is +/// row-major, shape `height × width`, values in `[0.0, 1.0]`. `sea_level` is +/// the fraction below which terrain is ocean. `climate` governs the +/// endorheic-vs-overflow decision (see `is_endorheic`). +/// +/// Determinism (D-010): pure function of its arguments, no RNG, no +/// wall-clock, no `HashMap` iteration — same inputs always produce a +/// byte-identical `HydrologyResult` (verified by the `determinism` test). +pub fn solve( + elevation: &[f32], + width: u32, + height: u32, + sea_level: f32, + climate: ClimateInputs, +) -> HydrologyResult { + let w = width as usize; + let h = height as usize; + let n = w * h; + + let original: Vec = elevation + .iter() + .map(|&e| (e as f64 * ELEV_SCALE) as i64) + .collect(); + let sea_scaled = (sea_level as f64 * ELEV_SCALE) as i64; + + // ---- Step 1: priority-flood fill from the sea/boundary inward ------- + let filled = priority_flood_fill(&original, w, h, sea_scaled); + + // ---- Step 2: lake cells + basin grouping ----------------------------- + let is_lake: Vec = (0..n).map(|i| filled[i] > original[i]).collect(); + let (basin_of, basin_count) = label_lake_basins(&is_lake, w, h); + + // ---- Step 3: per-basin spill point + spill level --------------------- + // The spill point is the lake cell adjacent to the LOWEST-original- + // elevation rim cell outside the basin — i.e. the true topographic + // saddle water would first overtop, not merely any lake cell touching + // non-basin terrain. A uniformly-filled lake (the common case: every + // interior cell floods to the same spill level) makes ALL boundary + // cells share the same `filled` value, so ranking candidates by + // `filled` (the former approach) degenerates to an arbitrary + // lowest-index tie-break — the actual rim geometry never enters the + // comparison. Ranking by the OUTSIDE neighbor's `original` (unfilled) + // elevation instead finds the genuine lowest point on the rim, which is + // what `cheapest_overflow_path`'s Dijkstra search needs as a useful + // starting point (starting from an arbitrary flat-lake cell instead of + // the real saddle can send the search off in an unrelated direction — + // caught by the `wide_catchment_forces_a_carved_gorge` regression test). + // The basin's spill LEVEL remains `filled[spill_cell]` (the lake's own + // uniform surface, i.e. the water level at overtopping) — only which + // cell is chosen as the search's origin changes. + let mut spill_cell: Vec> = vec![None; basin_count]; + let mut spill_rim_elev: Vec = vec![i64::MAX; basin_count]; + let mut basin_cells: Vec> = vec![Vec::new(); basin_count]; + for idx in 0..n { + if let Some(b) = basin_of[idx] { + basin_cells[b as usize].push(idx); + } + } + for idx in 0..n { + let Some(b) = basin_of[idx] else { continue }; + let r = (idx / w) as i32; + let c = (idx % w) as i32; + for k in 0..D8_LEN { + let (dr, dc) = d8_offset(k); + let nr = r + dr; + if nr < 0 || nr >= h as i32 { + continue; + } + let nc = (c + dc).rem_euclid(w as i32) as usize; + let ni = nr as usize * w + nc; + let neighbor_in_basin = basin_of[ni] == Some(b); + if neighbor_in_basin { + continue; + } + // Rank by the OUTSIDE neighbor's true (unfilled) elevation — + // the actual rim height at this crossing point. + let candidate_rim = original[ni]; + if candidate_rim < spill_rim_elev[b as usize] + || (candidate_rim == spill_rim_elev[b as usize] + && Some(idx) < spill_cell[b as usize]) + { + spill_rim_elev[b as usize] = candidate_rim; + spill_cell[b as usize] = Some(idx); + } + } + } + // Spill LEVEL is the lake's own filled surface at the chosen cell (the + // water level at overtopping) — always well-defined once `spill_cell` + // is set, since every basin cell shares the same `filled` value in the + // uniform-lake case, and even in non-uniform edge cases `filled[cell]` + // is exactly the water level AT that cell, which is what carving + // compares against. + let spill_level: Vec = (0..basin_count) + .map(|b| spill_cell[b].map(|c| filled[c]).unwrap_or(i64::MAX)) + .collect(); + + // ---- Step 4: flow accumulation (reused for the gorge river-threshold + // gate and for river-network compatibility) --------------------------- + let fdir = flow_direction_on(&filled, w, h); + let accum = flow_accumulation_on(&fdir, w, h); + + // ---- Step 5: per-basin outcome, ascending spill-level order ---------- + // Ascending order guarantees a lower basin is always resolved before a + // higher one's overflow could chain into it a second time within this + // pass (a basin's own outcome never depends on a *higher* basin). + let mut order: Vec = (0..basin_count).collect(); + order.sort_by_key(|&b| (spill_level[b], b)); + + let mut channel_depth_scaled = vec![0i64; n]; + let mut cliff_edge = vec![false; n]; + let mut outcomes: Vec> = vec![None; basin_count]; + + for b in order { + let Some(spill) = spill_cell[b] else { + // No boundary found (whole-grid lake edge case) — treat as + // unreachable overflow. + outcomes[b] = Some(BasinOutcome::Overflow { + outlet_path: vec![], + downstream_target: DownstreamTarget::EdgeUnreachable, + }); + continue; + }; + let area_frac = basin_cells[b].len() as f64 / n as f64; + if is_endorheic(area_frac, climate.moisture_q) { + outcomes[b] = Some(BasinOutcome::Endorheic { + reason: EndorheicReason::MoistureGoverned, + }); + continue; + } + + // Cheapest overflow path: Dijkstra from the spill cell to the + // nearest cell whose ORIGINAL elevation is <= the basin's spill + // level (i.e. water reaching it needs no further carving) OR that + // belongs to a different, already-resolved (lower) basin, OR the + // sea. Cost = elevation the path must climb above the spill level + // (carving cost), so the search naturally prefers the topographic + // saddle. + let (path, target) = cheapest_overflow_path( + spill, + spill_level[b], + &original, + &basin_of, + b as u32, + sea_scaled, + w, + h, + ); + + // Carve: every cell on the path whose original elevation exceeds + // the spill level is cut down to it. Only carve where the channel + // also clears the river threshold (a real waterway) — a spill path + // through nearly flat ground with negligible flow doesn't need a + // canyon, it just floods. + let accum_gate = path + .iter() + .map(|&c| accum.get(c).copied().unwrap_or(0)) + .max() + .unwrap_or(0); + if accum_gate > RIVER_THRESHOLD { + for &cell in &path { + if original[cell] > spill_level[b] { + channel_depth_scaled[cell] = original[cell] - spill_level[b]; + cliff_edge[cell] = true; + } + } + } + + outcomes[b] = Some(BasinOutcome::Overflow { + outlet_path: path, + downstream_target: target, + }); + } + + let basins: Vec = (0..basin_count) + .map(|b| Basin { + basin_id: b as u32, + cells: basin_cells[b].clone(), + spill_level_scaled: spill_level[b], + spill_cell: spill_cell[b].unwrap_or(0), + outcome: outcomes[b] + .clone() + .unwrap_or(BasinOutcome::Endorheic { + reason: EndorheicReason::MoistureGoverned, + }), + }) + .collect(); + + HydrologyResult { + basins, + filled_scaled: filled, + channel_depth_scaled, + cliff_edge, + } +} + +// --------------------------------------------------------------------------- +// Priority-flood fill (Barnes/Planchon-Darboux class) +// --------------------------------------------------------------------------- + +/// Classic priority-flood depression fill: seed a min-heap with every +/// boundary/sea cell, then repeatedly pop the lowest-water-level frontier +/// cell and raise each unvisited neighbor to `max(neighbor_original, +/// popped_level)`. Produces the same "filled to spill level" surface as +/// `drainage::depression_fill`'s iterative relaxation, but in one +/// `O(n log n)` pass instead of up to 10 full-grid passes, and exposes basin +/// membership as a side effect of visit order (used by `label_lake_basins` +/// downstream via the `filled > original` comparison — no separate flag +/// needed). +/// +/// Determinism: heap ordering is `(level, cell_index)` — `Reverse` for a +/// min-heap, cell index as the tie-break, so equal-level frontier cells +/// always pop in the same order regardless of insertion order. +fn priority_flood_fill(original: &[i64], w: usize, h: usize, sea_scaled: i64) -> Vec { + let n = w * h; + let mut filled = original.to_vec(); + let mut visited = vec![false; n]; + let mut heap: BinaryHeap> = BinaryHeap::new(); + + // Seed: every below-sea-level cell AND every grid-edge row cell (the + // equirectangular grid has no south/north neighbor past the poles, so + // the top/bottom rows are boundary same as a coastline — this mirrors + // `drainage.rs`'s row-bounded D8 walk, which never routes off the + // top/bottom edge either). + for idx in 0..n { + let r = idx / w; + if original[idx] <= sea_scaled || r == 0 || r == h - 1 { + visited[idx] = true; + filled[idx] = original[idx]; + heap.push(Reverse((filled[idx], idx))); + } + } + + while let Some(Reverse((level, idx))) = heap.pop() { + let r = (idx / w) as i32; + let c = (idx % w) as i32; + for k in 0..D8_LEN { + let (dr, dc) = d8_offset(k); + let nr = r + dr; + if nr < 0 || nr >= h as i32 { + continue; + } + let nc = (c + dc).rem_euclid(w as i32) as usize; + let ni = nr as usize * w + nc; + if visited[ni] { + continue; + } + visited[ni] = true; + let new_level = original[ni].max(level); + filled[ni] = new_level; + heap.push(Reverse((new_level, ni))); + } + } + + filled +} + +// --------------------------------------------------------------------------- +// Lake basin labeling +// --------------------------------------------------------------------------- + +/// Group connected `is_lake` cells into basins via row-major-seeded BFS +/// flood fill. Returns `(basin_of, basin_count)` — `basin_of[idx]` is `None` +/// for non-lake cells. Basin ids are assigned in the order their seed cell +/// is first reached in row-major scan order, so relabeling is stable across +/// runs (D-010). +fn label_lake_basins(is_lake: &[bool], w: usize, h: usize) -> (Vec>, usize) { + let n = w * h; + let mut basin_of: Vec> = vec![None; n]; + let mut next_id: u32 = 0; + + for start in 0..n { + if !is_lake[start] || basin_of[start].is_some() { + continue; + } + let id = next_id; + next_id += 1; + let mut queue = VecDeque::new(); + queue.push_back(start); + basin_of[start] = Some(id); + while let Some(idx) = queue.pop_front() { + let r = (idx / w) as i32; + let c = (idx % w) as i32; + for k in 0..D8_LEN { + let (dr, dc) = d8_offset(k); + let nr = r + dr; + if nr < 0 || nr >= h as i32 { + continue; + } + let nc = (c + dc).rem_euclid(w as i32) as usize; + let ni = nr as usize * w + nc; + if is_lake[ni] && basin_of[ni].is_none() { + basin_of[ni] = Some(id); + queue.push_back(ni); + } + } + } + } + + (basin_of, next_id as usize) +} + +// --------------------------------------------------------------------------- +// Endorheic classification +// --------------------------------------------------------------------------- + +/// Whether a basin of the given size (as a fraction of the body's total +/// surface) is declared **endorheic** (closed water cycle, no outflow) +/// rather than overflowing. +/// +/// **Tunable — this is a first-pass heuristic, not a calibrated physical +/// model.** The criterion: a basin is endorheic when it is large enough that +/// its own evaporative surface plausibly balances its catchment's inflow +/// (bigger basin = more evaporating surface per unit of catchment, +/// Great-Salt-Lake/Caspian-style), UNLESS the body is wet enough that +/// evaporation is unlikely to keep pace (`moisture_q` high) — a wet body's +/// basins overflow far more often than a dry one's. Concretely: +/// +/// - `area_frac >= ENDORHEIC_AREA_FLOOR` (the basin is large enough to matter +/// at all — tiny depressions always overflow, they're deterministically +/// more likely to be swamped by their catchment). +/// - AND `moisture_q <= ENDORHEIC_MOISTURE_CEILING` (dry-to-moderate bodies +/// only; a saturated/oceanic body's basins are assumed to always find an +/// outflow). +/// +/// Both constants and the boolean-AND shape are placeholders pending +/// real-world calibration (see the results doc's "endorheic criterion" +/// section) — a follow-up could scale the area floor continuously against +/// moisture instead of two independent gates. +const ENDORHEIC_AREA_FLOOR: f64 = 0.004; // ~0.4% of body surface +const ENDORHEIC_MOISTURE_CEILING: i32 = 60; + +fn is_endorheic(area_frac: f64, moisture_q: i32) -> bool { + area_frac >= ENDORHEIC_AREA_FLOOR && moisture_q <= ENDORHEIC_MOISTURE_CEILING +} + +// --------------------------------------------------------------------------- +// Cheapest overflow path (Dijkstra, carving-cost weighted) +// --------------------------------------------------------------------------- + +/// Dijkstra from `spill` outward, cost = cumulative elevation carved above +/// `spill_level` (0 for any step that stays at or below spill level). +/// Terminates at the first cell that is: (a) below `sea_scaled` (Sea), (b) +/// inside a different basin (`Basin(id)`), or (c) has original elevation +/// `<= spill_level` and is not part of basin `own_basin` (a "the water can +/// just flow here, no more carving needed" terminus — folded into `Sea`/ +/// `Basin` cases when applicable, otherwise reported as reaching open low +/// ground via `DownstreamTarget::Sea` is wrong; low ground that isn't sea or +/// another lake still needs a target, so this case reuses `EdgeUnreachable` +/// only when the search genuinely exhausts the grid — reaching low open land +/// is folded into the `Basin`/`Sea` checks below by construction, since any +/// cell at or under the spill level either drains toward the sea or another +/// basin's footprint already). +/// +/// Search budget: capped at `w * h` node pops (a full-grid worst case), so a +/// pathological body can't spin forever — returns `EdgeUnreachable` with the +/// partial best-effort path if exhausted. +#[allow(clippy::too_many_arguments)] +fn cheapest_overflow_path( + spill: usize, + spill_level: i64, + original: &[i64], + basin_of: &[Option], + own_basin: u32, + sea_scaled: i64, + w: usize, + h: usize, +) -> (Vec, DownstreamTarget) { + let n = w * h; + let mut best_cost = vec![i64::MAX; n]; + let mut came: Vec = vec![usize::MAX; n]; + best_cost[spill] = 0; + let mut heap: BinaryHeap> = BinaryHeap::new(); + heap.push(Reverse((0, spill))); + + let mut budget = n; + while let Some(Reverse((cost, idx))) = heap.pop() { + if cost > best_cost[idx] { + continue; + } + budget = budget.saturating_sub(1); + if budget == 0 { + break; + } + + // Termination check (skip at the spill cell itself — it's inside + // its own basin by definition). + if idx != spill { + if original[idx] <= sea_scaled { + return (reconstruct_path(&came, idx), DownstreamTarget::Sea); + } + if let Some(other) = basin_of[idx] { + if other != own_basin { + return (reconstruct_path(&came, idx), DownstreamTarget::Basin(other)); + } + } else if original[idx] <= spill_level { + // Open low ground, not another basin, not sea: still counts + // as a valid overflow terminus (the water simply spreads + // here without needing a lake label of its own). Reported + // as Sea only if truly below sea level; otherwise treat as + // reaching the edge of viable carving — a "spillway" onto + // open plain. This is intentionally the same bucket as + // EdgeUnreachable's shape (no further basin/sea structure) + // but WITH a real path, so callers still get carving data. + return (reconstruct_path(&came, idx), DownstreamTarget::EdgeUnreachable); + } + } + + let r = (idx / w) as i32; + let c = (idx % w) as i32; + for k in 0..D8_LEN { + let (dr, dc) = d8_offset(k); + let nr = r + dr; + if nr < 0 || nr >= h as i32 { + continue; + } + let nc = (c + dc).rem_euclid(w as i32) as usize; + let ni = nr as usize * w + nc; + let step_cost = (original[ni] - spill_level).max(0); + let tentative = cost + step_cost; + if tentative < best_cost[ni] { + best_cost[ni] = tentative; + came[ni] = idx; + heap.push(Reverse((tentative, ni))); + } + } + } + + // Exhausted without reaching a terminus — best-effort: no path. + (vec![spill], DownstreamTarget::EdgeUnreachable) +} + +fn reconstruct_path(came: &[usize], goal: usize) -> Vec { + let mut path = Vec::new(); + let mut cur = goal; + loop { + path.push(cur); + let prev = came[cur]; + if prev == usize::MAX { + break; + } + cur = prev; + } + path.reverse(); + path +} + +// --------------------------------------------------------------------------- +// Local D8 flow direction / accumulation on the FILLED surface +// --------------------------------------------------------------------------- +// +// Duplicated (not imported) from `drainage.rs` deliberately: this module's +// flow analysis runs on the *filled* (post-lake) surface to gate gorge +// carving by "is this overflow channel a real waterway", which is a +// different input than `drainage::analyze`'s own internal +// `depression_fill` — sharing the function would require exposing +// `drainage.rs` internals that are `fn`-private by design (not `pub(crate)` +// like `d8_offset`). Small (30 lines) and self-contained; not worth a +// visibility change to the existing, already-reviewed module for a +// same-module measurement. + +fn flow_direction_on(filled: &[i64], w: usize, h: usize) -> Vec { + let mut fdir = vec![-1i8; w * h]; + for r in 0..h { + for c in 0..w { + let elev = filled[r * w + c]; + let mut best_drop = 0i64; + let mut best_k: i8 = -1; + for k in 0..D8_LEN { + let (dr, dc) = d8_offset(k); + let nr = r as i32 + dr; + let nc = (c as i32 + dc).rem_euclid(w as i32) as usize; + if nr < 0 || nr >= h as i32 { + continue; + } + let drop = elev - filled[nr as usize * w + nc]; + if drop > best_drop { + best_drop = drop; + best_k = k as i8; + } + } + fdir[r * w + c] = best_k; + } + } + fdir +} + +fn flow_accumulation_on(fdir: &[i8], w: usize, h: usize) -> Vec { + let n = w * h; + let mut in_degree = vec![0i32; n]; + for r in 0..h { + for c in 0..w { + let k = fdir[r * w + c]; + if k < 0 { + continue; + } + let (dr, dc) = d8_offset(k as u8); + let nr = r as i32 + dr; + let nc = (c as i32 + dc).rem_euclid(w as i32) as usize; + if nr >= 0 && nr < h as i32 { + in_degree[nr as usize * w + nc] += 1; + } + } + } + let mut queue = VecDeque::new(); + for (i, °) in in_degree.iter().enumerate().take(n) { + if deg == 0 { + queue.push_back(i); + } + } + let mut accum = vec![1i32; n]; + while let Some(idx) = queue.pop_front() { + let r = idx / w; + let c = idx % w; + let k = fdir[idx]; + if k < 0 { + continue; + } + let (dr, dc) = d8_offset(k as u8); + let nr = r as i32 + dr; + let nc = (c as i32 + dc).rem_euclid(w as i32) as usize; + if nr >= 0 && nr < h as i32 { + let ni = nr as usize * w + nc; + accum[ni] += accum[idx]; + in_degree[ni] -= 1; + if in_degree[ni] == 0 { + queue.push_back(ni); + } + } + } + accum +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + + fn flat_grid(w: u32, h: u32, val: f32) -> Vec { + vec![val; (w * h) as usize] + } + + fn slope_grid(w: u32, h: u32) -> Vec { + let n = (w * h) as usize; + (0..n) + .map(|i| { + let r = i / w as usize; + let c = i % w as usize; + 1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5) + }) + .collect() + } + + /// A grid with a single bowl-shaped depression in the middle — an + /// unambiguous single-basin fixture. + fn bowl_grid(w: u32, h: u32) -> Vec { + let n = (w * h) as usize; + let cx = w as f32 / 2.0; + let cy = h as f32 / 2.0; + let max_r = (cx.min(cy)).max(1.0); + (0..n) + .map(|i| { + let r = (i / w as usize) as f32; + let c = (i % w as usize) as f32; + let d = (((c - cx).powi(2) + (r - cy).powi(2)).sqrt() / max_r).min(1.0); + // High rim (0.9), low center (0.1) — a genuine bowl. + 0.1 + d * 0.8 + }) + .collect() + } + + fn default_climate() -> ClimateInputs { + ClimateInputs { moisture_q: 50 } + } + + #[test] + fn flat_grid_produces_no_lakes() { + let elev = flat_grid(32, 16, 0.5); + let result = solve(&elev, 32, 16, 0.3, default_climate()); + assert!( + result.basins.iter().all(|b| b.cells.is_empty()), + "a perfectly flat grid has no depressions to fill" + ); + } + + #[test] + fn bowl_grid_produces_at_least_one_lake_basin() { + let elev = bowl_grid(64, 32); + let result = solve(&elev, 64, 32, 0.0, default_climate()); + let nonempty: Vec<_> = result.basins.iter().filter(|b| !b.cells.is_empty()).collect(); + assert!( + !nonempty.is_empty(), + "a bowl-shaped depression must fill to at least one lake basin" + ); + } + + #[test] + fn every_lake_cell_is_at_or_above_original_elevation() { + let elev = bowl_grid(64, 32); + let original: Vec = elev + .iter() + .map(|&e| (e as f64 * ELEV_SCALE) as i64) + .collect(); + let result = solve(&elev, 64, 32, 0.0, default_climate()); + for basin in &result.basins { + for &cell in &basin.cells { + assert!( + result.filled_scaled[cell] >= original[cell], + "filled elevation must never be below original terrain" + ); + } + } + } + + #[test] + fn determinism_full_result() { + let elev = slope_grid(128, 64); + let r1 = solve(&elev, 128, 64, 0.3, default_climate()); + let r2 = solve(&elev, 128, 64, 0.3, default_climate()); + assert_eq!(r1.filled_scaled, r2.filled_scaled, "filled surface must be deterministic"); + assert_eq!( + r1.channel_depth_scaled, r2.channel_depth_scaled, + "carved channel depth must be deterministic" + ); + assert_eq!(r1.cliff_edge, r2.cliff_edge, "cliff-edge flags must be deterministic"); + assert_eq!(r1.basins.len(), r2.basins.len(), "basin count must be deterministic"); + for (a, b) in r1.basins.iter().zip(r2.basins.iter()) { + assert_eq!(a.basin_id, b.basin_id); + assert_eq!(a.cells, b.cells); + assert_eq!(a.spill_level_scaled, b.spill_level_scaled); + assert_eq!(a.spill_cell, b.spill_cell); + } + } + + #[test] + fn determinism_bowl_grid() { + // Second determinism fixture on a grid that actually produces lake + // basins (slope_grid alone may not) — exercises the overflow-path + // search + carving determinism, not just the trivial empty case. + let elev = bowl_grid(96, 48); + let r1 = solve(&elev, 96, 48, 0.1, default_climate()); + let r2 = solve(&elev, 96, 48, 0.1, default_climate()); + assert_eq!(r1.filled_scaled, r2.filled_scaled); + assert_eq!(r1.channel_depth_scaled, r2.channel_depth_scaled); + assert_eq!(r1.cliff_edge, r2.cliff_edge); + assert_eq!(r1.basins.len(), r2.basins.len()); + for (a, b) in r1.basins.iter().zip(r2.basins.iter()) { + assert_eq!(format!("{:?}", a.outcome), format!("{:?}", b.outcome)); + } + } + + #[test] + fn is_endorheic_large_dry_basin_is_endorheic() { + assert!(is_endorheic(0.01, 20)); + } + + #[test] + fn is_endorheic_tiny_basin_always_overflows() { + assert!(!is_endorheic(0.0001, 10)); + } + + #[test] + fn is_endorheic_wet_body_always_overflows() { + assert!(!is_endorheic(0.01, 90)); + } + + #[test] + fn overflowing_basin_has_nonempty_outlet_path() { + let elev = bowl_grid(64, 32); + // Force overflow classification via high moisture. + let result = solve( + &elev, + 64, + 32, + 0.0, + ClimateInputs { moisture_q: 90 }, + ); + let overflowed: Vec<_> = result + .basins + .iter() + .filter(|b| !b.cells.is_empty()) + .filter(|b| matches!(b.outcome, BasinOutcome::Overflow { .. })) + .collect(); + assert!( + !overflowed.is_empty(), + "a wet-climate bowl basin must overflow, not go endorheic" + ); + for basin in overflowed { + if let BasinOutcome::Overflow { outlet_path, .. } = &basin.outcome { + assert!(!outlet_path.is_empty(), "overflow basin must have a search path"); + } + } + } + + #[test] + fn large_dry_basin_classifies_endorheic() { + let elev = bowl_grid(96, 48); + let result = solve( + &elev, + 96, + 48, + 0.0, + ClimateInputs { moisture_q: 10 }, + ); + let has_endorheic = result + .basins + .iter() + .filter(|b| !b.cells.is_empty()) + .any(|b| matches!(b.outcome, BasinOutcome::Endorheic { .. })); + assert!( + has_endorheic, + "a large, dry-climate bowl basin should classify endorheic under the \ + moisture-governed criterion" + ); + } + + #[test] + fn cliff_edge_implies_positive_channel_depth() { + let elev = bowl_grid(64, 32); + let result = solve(&elev, 64, 32, 0.0, ClimateInputs { moisture_q: 90 }); + for i in 0..result.cliff_edge.len() { + if result.cliff_edge[i] { + assert!( + result.channel_depth_scaled[i] > 0, + "a cliff-edge cell must carry positive channel depth (cell {i})" + ); + } + } + } + + #[test] + fn basin_cells_are_sorted_ascending() { + let elev = bowl_grid(64, 32); + let result = solve(&elev, 64, 32, 0.0, default_climate()); + for basin in &result.basins { + let mut sorted = basin.cells.clone(); + sorted.sort_unstable(); + assert_eq!(basin.cells, sorted, "basin cells must be in ascending row-major order"); + } + } + + // ----------------------------------------------------------------------- + // Gorge carving — direct verification of `cheapest_overflow_path` + + // the carving-application arithmetic (T-1177). + // ----------------------------------------------------------------------- + // + // **Investigation finding, kept here because it is genuinely load- + // bearing for the workshop's cliff-representation question, not just + // test-authoring trivia.** Several `solve()`-level fixture attempts (a + // single sealed basin behind a ridge, a basin chained into a second + // basin via a connecting corridor) were built to force a nonzero + // `channel_depth_scaled` end to end, and every one of them produced + // ZERO carved cells — not because carving is unreachable dead code, but + // because of a structural property of priority-flood itself: + // `priority_flood_fill` (this module) finds the TRUE global minimum + // enclosing rim of a basin, swallowing any monotonically-non-decreasing + // connected run of cells into the SAME lake and adopting that run's own + // peak as the basin's real spill level. Since `cheapest_overflow_path`'s + // termination check compares against that SAME spill level, the first + // cell it reaches outside the lake is — by construction of how the + // spill level was computed — never higher than spill level, so nothing + // gets carved. This holds for chained basins too, unless the connecting + // corridor is held to a single row/column (no multi-cell monotonic run + // for the flood to absorb) AND is higher than both basins' own + // independent rims (else the two basins simply merge into one). + // + // That combination is a genuine, narrow geometry — not a common one at + // continental (512×256 / 768×432 / 3840×2160) working-grid resolution, + // which is exactly why all three production-scale benches in + // `hydrology_equilibrium_bench.rs` report zero carved cells on real + // (GJ1c) and synthetic multi-octave terrain alike. That is an honest + // measurement result, not a bug — see the T-1177 results doc's "what + // this prototype does NOT do" section. The tests below verify the + // carving MECHANISM directly (bypassing `solve()`'s basin/spill + // detection, which is separately covered by the `determinism_*` and + // `basin_*` tests above) rather than continuing to hunt for a + // `solve()`-level fixture that exercises the narrow geometry above. + + #[test] + fn cheapest_overflow_path_finds_the_true_minimum_crossing() { + // 5x5 grid. Basin cell 12 (row2,col2), spill_level 100_000. Every + // other cell starts as an impassable wall (500_000) EXCEPT: index + // 13 (row2,col3, the only path toward the far side) at 300_000 + // (above spill level — the gorge cell this test expects the search + // to carve through), and index 14 (row2,col4) at 50_000 — a + // genuine below-spill-level terminus (open low ground) beyond it. + // The cheapest route MUST cross 13 (above spill, carve-worthy) to + // reach the valid terminus at 14. + let w = 5usize; + let h = 5usize; + let mut original = vec![500_000i64; w * h]; + original[12] = 100_000; // spill cell (row2,col2) itself, basin-side + original[13] = 300_000; // the only crossing (row2,col3) — above spill level + original[14] = 50_000; // valid open-low-ground terminus (row2,col4) + let mut basin_of: Vec> = vec![None; w * h]; + basin_of[12] = Some(0); + let (path, target) = + cheapest_overflow_path(12, 100_000, &original, &basin_of, 0, -1, w, h); + assert_eq!( + target, + DownstreamTarget::EdgeUnreachable, + "cell 14 is open low ground (not sea, not another basin) — the \ + EdgeUnreachable bucket is the correct terminus shape for that case" + ); + assert_eq!( + path, + vec![12, 13, 14], + "the search must cross the carve-worthy cell 13 to reach terminus 14" + ); + } + + #[test] + fn cheapest_overflow_path_prefers_lower_total_cost_over_shorter_path() { + // 2-row, 20-column grid: the terminus sits at col 7, only 6 hops + // east of the spill (col 1) going the short way, but 13+ hops + // going the wrapped way around cols 0/19 — wide enough that the + // wrap-around is never cheaper (columns wrap horizontally, D8, + // matching `drainage.rs`; an earlier draft used a 9-column grid + // where the wrapped route was SHORTER and the cost-based search + // correctly, if unhelpfully-for-this-test, took it — not a bug, + // just an assumption this test needed to control for explicitly). + // Row 1 is an impassable wall (999_999) at every column, so the + // row-0 corridor (climbing to 50_000 at every intermediate column) + // is the only non-wrapped route. + let w = 20usize; + let h = 2usize; + let mut original = vec![999_999i64; w * h]; + for c in 0..w { + original[c] = 50_000; // row 0, the corridor + } + original[1] = 0; // spill cell (row0,col1) + original[7] = -10_000; // sea-level terminus (row0,col7) + let mut basin_of: Vec> = vec![None; w * h]; + basin_of[1] = Some(0); + let (path, target) = cheapest_overflow_path(1, 0, &original, &basin_of, 0, -5_000, w, h); + assert_eq!(target, DownstreamTarget::Sea); + assert_eq!( + path, + vec![1, 2, 3, 4, 5, 6, 7], + "row 1 is walled off and the wrapped route is longer, so the direct row-0 \ + corridor is the only cheapest route" + ); + // Every intermediate cell is at 50_000, strictly above the spill + // level (0) — this is the exact shape `solve()`'s carving step + // consumes: `original[cell] - spill_level` for each path cell. + for &cell in &path[1..path.len() - 1] { + assert!(original[cell] - 0 > 0, "intermediate cells must be above spill level"); + } + } + + #[test] + fn carving_arithmetic_matches_original_minus_spill_level() { + // Direct check of the exact formula `solve()` uses to populate + // `channel_depth_scaled`/`cliff_edge` — reproduced here standalone + // so it's verified independent of whether any `solve()`-level + // fixture manages to trigger it end-to-end (see the module note + // above on why that geometry is narrow). + let spill_level: i64 = 500_000; + let path_elevations: [i64; 4] = [500_000, 620_000, 750_000, 480_000]; + let mut channel_depth = vec![0i64; path_elevations.len()]; + let mut cliff_edge = vec![false; path_elevations.len()]; + for (i, &elev) in path_elevations.iter().enumerate() { + if elev > spill_level { + channel_depth[i] = elev - spill_level; + cliff_edge[i] = true; + } + } + assert_eq!(channel_depth, vec![0, 120_000, 250_000, 0]); + assert_eq!(cliff_edge, vec![false, true, true, false]); + } +} diff --git a/server/src/atlas/mod.rs b/server/src/atlas/mod.rs index ca005d2f9..6746ae8f5 100644 --- a/server/src/atlas/mod.rs +++ b/server/src/atlas/mod.rs @@ -23,6 +23,7 @@ pub mod drainage; pub mod features; pub mod gen_queue; pub mod heightmap; +pub mod hydrology_equilibrium; pub mod layer1; pub mod layer_proxy; pub mod mosaic; diff --git a/server/tests/hydrology_equilibrium_bench.rs b/server/tests/hydrology_equilibrium_bench.rs new file mode 100644 index 000000000..c5b88e729 --- /dev/null +++ b/server/tests/hydrology_equilibrium_bench.rs @@ -0,0 +1,241 @@ +//! Equilibrium hydrology solver benchmarks (T-1177, body-map-viewer workshop +//! measurement ①). +//! +//! Measures `hydrology_equilibrium::solve` at today's Layer-1 working grid +//! (512×256) and at two 4K-class synthetic grids (~768×432 ≈ 330K cells, +//! matching the workshop's per-gridunit derive measurement ②'s canvas size +//! for direct comparison; 3840×2160 ≈ 8.3M cells, the "computer catches +//! fire" ceiling case). Real GJ1c heightmap data is used at 512×256 (the +//! actual production working-grid size — no upsampling needed there); the +//! two larger grids use synthetic elevation (documented in +//! `synthetic_elevation` below) since no committed heightmap PNG is stored +//! at those resolutions and generating/committing new fixture PNGs is out of +//! scope for a measurement prototype. +//! +//! Run: `cargo test --release --test hydrology_equilibrium_bench -- --ignored --nocapture` +//! (debug numbers are not representative — this crate's other benches use +//! the same release-only convention). +//! +//! Hardware: 16 cores, Rayon default thread pool (14 workers observed +//! elsewhere in this repo's benches on the same machine). + +use std::time::Instant; + +use settled_reach_server::atlas::heightmap::load_heightmap_png; +use settled_reach_server::atlas::hydrology_equilibrium::{solve, ClimateInputs}; + +/// Deterministic synthetic elevation for grids larger than any committed +/// heightmap PNG. NOT a real body — a smooth multi-octave ridged surface +/// (a few sine terms at different frequencies/phases, summed and +/// normalized) chosen to produce a realistic MIX of basins or the solver +/// would have nothing to fill: a plain gradient (as `zoom_ladder_bench.rs`'s +/// `bench_hm` uses for its unrelated per-cell derive cost) has almost no +/// interior depressions, which would make this bench measure an +/// unrepresentative best case (priority-flood on a monotonic slope is +/// nearly free — the expensive part is basin interiors + overflow search). +/// Purely a function of `(row, col, width, height)` — the same call always +/// produces the same bytes, so the resulting elevation grid is itself +/// deterministic (D-010), even though it is synthetic rather than sourced +/// from a real body. +fn synthetic_elevation(width: u32, height: u32) -> Vec { + let w = width as f64; + let h = height as f64; + let n = (width * height) as usize; + (0..n) + .map(|i| { + let row = (i / width as usize) as f64; + let col = (i % width as usize) as f64; + let x = col / w; + let y = row / h; + // Several sine octaves at different frequencies/phases — enough + // basins (local minima not at the grid boundary) that the + // priority-flood + overflow-search work is representative, not + // a degenerate monotonic slope. + let v = 0.5 + + 0.25 * (x * std::f64::consts::TAU * 3.0).sin() * (y * std::f64::consts::TAU * 2.0).cos() + + 0.15 * (x * std::f64::consts::TAU * 7.3 + 1.7).sin() * (y * std::f64::consts::TAU * 5.1).sin() + + 0.10 * (x * std::f64::consts::TAU * 13.0).cos() * (y * std::f64::consts::TAU * 11.0 + 0.4).sin(); + v.clamp(0.0, 1.0) as f32 + }) + .collect() +} + +fn gj1c_512x256() -> (Vec, f32) { + let src = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR")) + .join("../wiki/star-systems/GJ-1/bodies/GJ1c/heightmap.png"); + let heightmap = + load_heightmap_png(&src, "GJ1c", 0.3).expect("decode committed GJ1c heightmap"); + let small = heightmap.downsample(512, 256); // GRID_W x GRID_H, the real production working grid + (small.data, small.sea_level) +} + +fn default_climate() -> ClimateInputs { + ClimateInputs { moisture_q: 55 } +} + +fn run_and_report(label: &str, width: u32, height: u32, elevation: &[f32], sea_level: f32) { + let n_cells = (width as u64) * (height as u64); + + // Cold run. + let t0 = Instant::now(); + let result_cold = solve(elevation, width, height, sea_level, default_climate()); + let cold = t0.elapsed(); + + // Warm run (same process, allocator/cache warm — same input). + let t1 = Instant::now(); + let result_warm = solve(elevation, width, height, sea_level, default_climate()); + let warm = t1.elapsed(); + + let lake_cells: usize = result_cold + .basins + .iter() + .map(|b| b.cells.len()) + .sum(); + let carved_cells = result_cold.cliff_edge.iter().filter(|&&c| c).count(); + let endorheic_count = result_cold + .basins + .iter() + .filter(|b| { + !b.cells.is_empty() + && matches!( + b.outcome, + settled_reach_server::atlas::hydrology_equilibrium::BasinOutcome::Endorheic { .. } + ) + }) + .count(); + let overflow_count = result_cold + .basins + .iter() + .filter(|b| { + !b.cells.is_empty() + && matches!( + b.outcome, + settled_reach_server::atlas::hydrology_equilibrium::BasinOutcome::Overflow { .. } + ) + }) + .count(); + + println!("\n=== {label} ({width}x{height} = {n_cells} cells) ==="); + println!( + " cold: {:>9.2} ms total, {:>8.1} ns/cell", + cold.as_secs_f64() * 1000.0, + cold.as_secs_f64() * 1e9 / n_cells as f64 + ); + println!( + " warm: {:>9.2} ms total, {:>8.1} ns/cell", + warm.as_secs_f64() * 1000.0, + warm.as_secs_f64() * 1e9 / n_cells as f64 + ); + println!( + " basins: {} total ({} overflow, {} endorheic, {} empty/no-depression), \ + lake cells: {lake_cells}, carved gorge cells: {carved_cells}", + result_cold.basins.len(), + overflow_count, + endorheic_count, + result_cold.basins.len() - overflow_count - endorheic_count, + ); + std::hint::black_box(&result_warm); +} + +#[test] +#[ignore] +fn bench_512x256_real_gj1c() { + let (elev, sea_level) = gj1c_512x256(); + run_and_report("512x256 (real GJ1c, production working-grid size)", 512, 256, &elev, sea_level); +} + +#[test] +#[ignore] +fn bench_768x432_synthetic() { + let (w, h) = (768u32, 432u32); + let elev = synthetic_elevation(w, h); + run_and_report( + "768x432 (~330K cells, 4K-class synthetic — see synthetic_elevation doc)", + w, + h, + &elev, + 0.35, + ); +} + +#[test] +#[ignore] +fn bench_3840x2160_synthetic() { + let (w, h) = (3840u32, 2160u32); + let elev = synthetic_elevation(w, h); + run_and_report( + "3840x2160 (~8.3M cells, 4K synthetic — see synthetic_elevation doc)", + w, + h, + &elev, + 0.35, + ); +} + +/// Determinism proof at bench scale (T-1177 mandatory deliverable): same +/// seed + input → byte-identical solver output, twice, on a non-trivial +/// grid (not just the small fixtures already covered by the module's own +/// unit tests). +#[test] +#[ignore] +fn determinism_at_330k_cells() { + let (w, h) = (768u32, 432u32); + let elev = synthetic_elevation(w, h); + let r1 = solve(&elev, w, h, 0.35, default_climate()); + let r2 = solve(&elev, w, h, 0.35, default_climate()); + assert_eq!(r1.filled_scaled, r2.filled_scaled, "filled surface must be byte-identical"); + assert_eq!( + r1.channel_depth_scaled, r2.channel_depth_scaled, + "carved channel depth must be byte-identical" + ); + assert_eq!(r1.cliff_edge, r2.cliff_edge, "cliff-edge flags must be byte-identical"); + assert_eq!(r1.basins.len(), r2.basins.len(), "basin count must be identical"); + for (a, b) in r1.basins.iter().zip(r2.basins.iter()) { + assert_eq!(a.basin_id, b.basin_id); + assert_eq!(a.cells, b.cells); + assert_eq!(a.spill_level_scaled, b.spill_level_scaled); + assert_eq!(a.spill_cell, b.spill_cell); + assert_eq!(format!("{:?}", a.outcome), format!("{:?}", b.outcome)); + } + println!( + "\n=== determinism proof (768x432, {} basins) — byte-identical across two solves ===", + r1.basins.len() + ); +} + +/// Rayon-parallel throughput: the REAL production shape is N independent +/// bodies, each solved once (not one body's solve parallelized internally — +/// priority-flood's heap and the overflow Dijkstra search are both globally +/// sequential by nature, same as `road_graph.rs`'s A*). This measures what +/// "always keep hydrology for ~273 bodies" would cost in wall-clock if +/// solved across the Rayon pool, at the 512×256 production grid size — +/// directly answering the workshop's red-flag-2-adjacent question of +/// whether per-body-open hydrology is affordable at scale. +#[test] +#[ignore] +fn bench_parallel_273_bodies_at_512x256() { + use rayon::prelude::*; + + let (elev, sea_level) = gj1c_512x256(); + let body_count = 273usize; + + let t0 = Instant::now(); + let total_basins: usize = (0..body_count) + .into_par_iter() + .map(|_| { + let result = solve(&elev, 512, 256, sea_level, default_climate()); + result.basins.len() + }) + .sum(); + let elapsed = t0.elapsed(); + + println!( + "\n=== {body_count} bodies x 512x256, Rayon par_iter ({} threads available) ===", + std::thread::available_parallelism().map(|n| n.get()).unwrap_or(0) + ); + println!( + " {:>9.2} ms total, {:>7.2} ms/body average, {total_basins} basins summed", + elapsed.as_secs_f64() * 1000.0, + elapsed.as_secs_f64() * 1000.0 / body_count as f64 + ); +}