feat(simulation): equilibrium hydrology solver prototype + bench (T-1177)

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 <noreply@anthropic.com>

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;

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<f32> {
+    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>, 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
+    );
+}