facett-cfd/src/pipeflow.rs
//! **pipeflow** — *real liquid* in a pipe graph: a **gravity-driven, capacity-limited
//! 1-D free-surface solver** constrained to the tubes.
//!
//! This is the additive sibling of [`crate::fluid`]. Where [`FluidField`](crate::FluidField)
//! advects a **dye** at a velocity the host sets (a pretty flow indicator, no physics),
//! `PipeFlow` simulates a **liquid with volume** that
//!
//! 1. **falls** — each pipe's axial acceleration is *gravity projected onto that pipe's
//! own direction*, `a_p = |g| · (ĝ · d̂_p)`, so tipping the map re-aims every pipe's
//! flow at once;
//! 2. **stays inside the pipe** — a cell's state is a **fill fraction** `f ∈ [0, 1]` of
//! the tube's own cross-section, and every transfer is limited by both the source's
//! liquid *and* the destination's free space, so `f` can never leave `[0, 1]`;
//! 3. **is conserved** — every movement of liquid in the whole step is a **paired
//! transfer of one computed scalar** (`x` off the source, the same `x` onto the
//! destination). Conservation is therefore *structural*, not something the clamps
//! have to be careful about: no clamp can create or destroy liquid, only refuse it;
//! 4. **pools in the tanks** — the graph's nodes are reservoirs, so liquid genuinely
//! pours out of a downhill pipe into a junction and on out of that junction's other
//! downhill pipes.
//!
//! ## Why 1-D along the pipe and not a 3-D solver
//! A pipe network is a **graph of tubes**. A free-surface SPH/FLIP solver would spend
//! all its work rediscovering the boundary we already know exactly — and would need
//! collision geometry to keep particles inside, i.e. containment becomes a *tolerance*
//! instead of an *invariant*. Constrained to the tube axis, containment is a type-level
//! fact (`f ≤ 1`) and the whole solve is `O(pipes · cells)` flops. See
//! `facett/.nornir/systemmap-fluid-reasoning.md`.
//!
//! ## Determinism
//! No wall clock, no RNG at runtime: [`step`](PipeFlow::step) takes a host-supplied
//! `dt`, gravity is set by the host from its camera, and the initial state comes from a
//! **seeded** splitmix64 ([`seed_fill`](PipeFlow::seed_fill)). The gatling fan-out
//! writes disjoint per-pipe rows, so it is bit-identical to
//! [`step_serial`](PipeFlow::step_serial) (a test asserts exactly that).
//!
//! ## Cost
//! The advection fans over pipes through the znippy **gatling** fork-join (rayon is
//! forbidden in facett — ROOT LAW #0). The junction pass is a scatter onto shared nodes,
//! so it stays **serial in pipe-index order** — `O(pipes)` scalar work, far below the
//! dispatch cost of parallelising it, and order-independent determinism comes free.
use znippy_zoomies::gatling_forkjoin::gatling_run;
/// Viscous damping (per second) applied to each pipe's axial velocity, so liquid
/// **settles** instead of sloshing forever.
const DRAG: f32 = 2.6;
/// Hard cap on the per-step Courant number `|v·dt|` in cells. Below 1 the upwind
/// transfer can never skip a cell, which is what makes the fill bound inductive.
const MAX_COURANT: f32 = 0.45;
/// Hard cap on |axial velocity| (cells/second) — a numerical blip can never run away.
const MAX_VEL: f32 = 40.0;
/// **Levelling** rate (per second): how fast neighbouring cells equalise their depth.
/// Liquid in a tube finds its own level — without this a pipe keeps whatever surface
/// profile it was charged with for ever, which is both wrong and looks like hatching
/// rather than like water.
const LEVEL_RELAX: f32 = 4.0;
/// Junction exchange gain (per second) between a pipe end cell and its node tank.
const JUNCTION_GAIN: f32 = 5.0;
/// Default node tank capacity, in cell-volumes.
const DEFAULT_TANK_CAP: f32 = 6.0;
/// **Total cell count at/above which [`step`](PipeFlow::step) fans out** across cores
/// via gatling instead of running the serial kernel. A *measured* threshold, not a
/// guess: `gatling_run` is `std::thread::scope` with one spawn per core per call, so it
/// carries a fixed dispatch cost of roughly a millisecond on a 32-core host. Below this
/// size that dispatch is hundreds of times the entire solve — at the System Map's 528
/// cells it is ~1.9 ms of thread churn against ~1.8 µs of work. See
/// `facett-cfd/examples/pipeflow_bench.rs` and `facett/.nornir/systemmap-fluid-design.md`
/// f6. Both lanes are bit-identical (a test asserts it), so this only ever changes speed.
const PARALLEL_MIN_CELLS: usize = 500_000;
/// A `Send`+`Sync` raw base pointer into the `back` buffer so gatling workers each
/// write their **own** disjoint pipe row without a lock. Soundness: worker `p` writes
/// only `back[p*cells .. (p+1)*cells]` and every `p` is claimed by exactly one worker,
/// so the writes never alias — the same disjoint-index argument [`crate::fluid`] and
/// `facett_core::barnes_hut` use.
#[derive(Clone, Copy)]
struct RowsMut(*mut f32);
// SAFETY: only ever used to write disjoint, per-worker-unique index ranges.
unsafe impl Send for RowsMut {}
unsafe impl Sync for RowsMut {}
/// **Liquid in a pipe graph.** `pipes` tubes of `cells` cells each, joined at `nodes`
/// junction tanks, under a gravity vector the host supplies in **model space**.
pub struct PipeFlow {
pipes: usize,
cells: usize,
nodes: usize,
/// `(from, to)` node index per pipe. Cell `0` sits at `from`, cell `cells-1` at `to`.
ends: Vec<(u32, u32)>,
/// Unit direction `from → to` per pipe, in model space.
axis: Vec<[f32; 3]>,
/// Node positions in model space (for the geometry + the fall moment).
pos: Vec<[f32; 3]>,
/// Fill fraction per cell, row-major `[p*cells + i]`, invariant `0 ≤ f ≤ 1`.
front: Vec<f32>,
/// Scratch back buffer; swapped with `front` (double buffer, pointer swap).
back: Vec<f32>,
/// Signed axial velocity per pipe, in cells/second (`> 0` flows `from → to`).
vel: Vec<f32>,
/// Liquid held at each junction, in cell-volumes.
tank: Vec<f32>,
/// Junction capacity, in cell-volumes.
tank_cap: Vec<f32>,
/// Gravity **direction** in model space (unit).
g: [f32; 3],
/// Gravity magnitude, in cells/second².
g_mag: f32,
/// Accumulated sim time (seconds) — host-supplied `dt` only, never a wall clock.
t: f32,
}
impl PipeFlow {
/// Allocate a dry network: `nodes` junctions, `edges` pipes of `cells` cells, at the
/// given model-space node positions. **The only allocation** — [`step`](Self::step)
/// never allocates.
pub fn new(pos: Vec<[f32; 3]>, edges: &[(usize, usize)], cells: usize) -> Self {
let cells = cells.max(2);
let nodes = pos.len();
let pipes = edges.len();
let mut f = Self {
pipes,
cells,
nodes,
ends: edges.iter().map(|&(a, b)| (a as u32, b as u32)).collect(),
axis: vec![[0.0, -1.0, 0.0]; pipes],
pos,
front: vec![0.0; pipes * cells],
back: vec![0.0; pipes * cells],
vel: vec![0.0; pipes],
tank: vec![0.0; nodes],
tank_cap: vec![DEFAULT_TANK_CAP; nodes],
g: [0.0, -1.0, 0.0],
g_mag: 9.0,
t: 0.0,
};
f.refresh_geometry();
f
}
/// Recompute every pipe's unit axis from the current node positions. Call after the
/// host re-lays-out the graph (the pipes moved, so gravity hits them differently).
pub fn refresh_geometry(&mut self) {
for p in 0..self.pipes {
let (a, b) = self.ends[p];
let (a, b) = (a as usize, b as usize);
self.axis[p] = if a < self.nodes && b < self.nodes {
unit(sub(self.pos[b], self.pos[a]))
} else {
[0.0, -1.0, 0.0]
};
}
}
/// Move a node (and re-derive the incident pipe axes on the next
/// [`refresh_geometry`](Self::refresh_geometry)).
pub fn set_node_pos(&mut self, n: usize, p: [f32; 3]) {
if let Some(s) = self.pos.get_mut(n) {
*s = p;
}
}
/// Set gravity's **direction in model space** (normalised; a zero vector is ignored).
/// This is the whole trick: the host rotates its camera, expresses "screen-down" in
/// the model's own frame, and hands it here — every pipe then re-derives its own
/// axial acceleration from it.
pub fn set_gravity(&mut self, dir: [f32; 3]) {
let n = norm(dir);
if n > 1e-6 {
self.g = [dir[0] / n, dir[1] / n, dir[2] / n];
}
}
/// Set gravity's magnitude, in cells/second².
pub fn set_gravity_magnitude(&mut self, g: f32) {
self.g_mag = g.max(0.0);
}
/// Gravity direction in model space (unit).
#[must_use]
pub fn gravity(&self) -> [f32; 3] {
self.g
}
/// **Deterministic seeded initial charge.** Fills every pipe to roughly `level`
/// (± a seeded jitter) and every tank to `level` of its capacity, from a splitmix64
/// stream — no `rand`, no clock, so two runs with the same seed are bit-identical.
pub fn seed_fill(&mut self, seed: u64, level: f32) {
let level = level.clamp(0.0, 1.0);
let mut s = seed;
let mut next = move || -> f32 {
// splitmix64 → [0,1)
s = s.wrapping_add(0x9E37_79B9_7F4A_7C15);
let mut z = s;
z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
z ^= z >> 31;
(z >> 40) as f32 / 16_777_216.0
};
for v in self.front.iter_mut() {
*v = (level + (next() - 0.5) * 0.5 * level).clamp(0.0, 1.0);
}
self.back.copy_from_slice(&self.front);
for n in 0..self.nodes {
self.tank[n] = (self.tank_cap[n] * level * (0.6 + 0.8 * next())).clamp(0.0, self.tank_cap[n]);
}
for v in self.vel.iter_mut() {
*v = 0.0;
}
}
/// Charge cell `i` of pipe `p` to fill fraction `f` (clamped to `0..=1`). The host's
/// way to place liquid explicitly; out-of-range is ignored.
pub fn set_fill(&mut self, p: usize, i: usize, f: f32) {
if p < self.pipes && i < self.cells {
self.front[p * self.cells + i] = f.clamp(0.0, 1.0);
}
}
/// Set junction `n`'s liquid volume in cell-volumes (clamped to its capacity).
pub fn set_tank(&mut self, n: usize, vol: f32) {
if n < self.nodes {
self.tank[n] = vol.clamp(0.0, self.tank_cap[n]);
}
}
/// Set junction `n`'s capacity in cell-volumes (its current volume is re-clamped).
pub fn set_tank_capacity(&mut self, n: usize, cap: f32) {
if n < self.nodes {
self.tank_cap[n] = cap.max(0.0);
self.tank[n] = self.tank[n].min(self.tank_cap[n]);
}
}
/// Empty the whole network (pipes, tanks and velocities) — the dry start a host
/// charges by hand with [`set_fill`](Self::set_fill) / [`set_tank`](Self::set_tank).
pub fn drain(&mut self) {
self.front.fill(0.0);
self.back.fill(0.0);
self.tank.fill(0.0);
self.vel.fill(0.0);
}
/// Number of pipes.
#[must_use]
pub fn pipes(&self) -> usize {
self.pipes
}
/// Cells per pipe.
#[must_use]
pub fn cells(&self) -> usize {
self.cells
}
/// Number of junction tanks.
#[must_use]
pub fn nodes(&self) -> usize {
self.nodes
}
/// Accumulated sim time (seconds).
#[must_use]
pub fn time(&self) -> f32 {
self.t
}
/// **Borrow** pipe `p`'s fill row (`cells` fractions in `0..=1`) — the zero-copy read
/// the renderer samples along the tube. Empty slice if `p` is out of range.
#[must_use]
pub fn fill(&self, p: usize) -> &[f32] {
if p < self.pipes {
&self.front[p * self.cells..(p + 1) * self.cells]
} else {
&[]
}
}
/// Junction `n`'s level as a fraction of its capacity (`0..=1`).
#[must_use]
pub fn tank_level(&self, n: usize) -> f32 {
match (self.tank.get(n), self.tank_cap.get(n)) {
(Some(&v), Some(&c)) if c > 0.0 => (v / c).clamp(0.0, 1.0),
_ => 0.0,
}
}
/// Pipe `p`'s signed axial velocity (cells/second), `0.0` out of range.
#[must_use]
pub fn velocity(&self, p: usize) -> f32 {
self.vel.get(p).copied().unwrap_or(0.0)
}
/// Pipe `p`'s **axial gravity component** `|g|·(ĝ·d̂_p)` — the acceleration the
/// liquid in that pipe actually feels. Positive drives `from → to`.
#[must_use]
pub fn axial_gravity(&self, p: usize) -> f32 {
self.axis.get(p).map(|d| self.g_mag * dot(self.g, *d)).unwrap_or(0.0)
}
/// **Total liquid volume** in the sealed system (pipes + tanks), in cell-volumes.
/// This is the conserved quantity: it must not drift.
#[must_use]
pub fn total_volume(&self) -> f32 {
self.front.iter().copied().sum::<f32>() + self.tank.iter().copied().sum::<f32>()
}
/// The largest fill fraction anywhere. **The containment invariant is `≤ 1.0`** —
/// a cell holding more than its tube's cross-section *is* liquid outside the pipe.
#[must_use]
pub fn max_fill(&self) -> f32 {
self.front.iter().copied().fold(0.0f32, f32::max)
}
/// The smallest fill fraction anywhere (must stay `≥ 0`: liquid is never negative).
#[must_use]
pub fn min_fill(&self) -> f32 {
self.front.iter().copied().fold(1.0f32, f32::min)
}
/// **How far the liquid has fallen** — the volume-weighted distance along gravity,
/// `Σ vol · (ĝ · x)`. Falling *increases* it; it is (minus) the gravitational
/// potential energy. A settling run under fixed gravity must end higher than it
/// started, which is the "it really falls" property.
#[must_use]
pub fn fall_moment(&self) -> f32 {
let mut m = 0.0f32;
for p in 0..self.pipes {
let (a, b) = self.ends[p];
let (a, b) = (a as usize, b as usize);
if a >= self.nodes || b >= self.nodes {
continue;
}
let (pa, pb) = (self.pos[a], self.pos[b]);
for i in 0..self.cells {
let f = self.front[p * self.cells + i];
if f == 0.0 {
continue;
}
let t = (i as f32 + 0.5) / self.cells as f32;
let x = [pa[0] + (pb[0] - pa[0]) * t, pa[1] + (pb[1] - pa[1]) * t, pa[2] + (pb[2] - pa[2]) * t];
m += f * dot(self.g, x);
}
}
for n in 0..self.nodes {
m += self.tank[n] * dot(self.g, self.pos[n]);
}
m
}
/// **The step (gatling · zero-alloc · zero-copy · zero-wait).** Advance the whole
/// network by `dt` host-supplied seconds:
///
/// 1. **gravity → velocity** (serial, `O(pipes)`): `v ← (v + |g|(ĝ·d̂)·dt) / (1 + drag·dt)`;
/// 2. **capacity-limited upwind advection** along each pipe (gatling fan-out over
/// pipes — disjoint rows), then the double-buffer pointer swap;
/// 3. **junction exchange** (serial, `O(pipes)`): the downhill end drains into its
/// node tank, the uphill tank feeds the pipe — each a paired transfer.
///
/// **The lane is chosen by measurement.** Below [`PARALLEL_MIN_CELLS`] this is
/// [`step_serial`](Self::step_serial); above it, [`step_gatling`](Self::step_gatling).
/// The two are bit-identical, so the choice only ever changes speed — and at any size
/// a UI pane holds, the fork-join's per-call thread spawn costs orders of magnitude
/// more than the whole solve.
pub fn step(&mut self, dt: f32) {
if self.pipes * self.cells >= PARALLEL_MIN_CELLS {
self.step_gatling(dt);
} else {
self.step_serial(dt);
}
}
/// The **fan-out lane** — the per-pipe advection spread across every core through the
/// znippy gatling fork-join, each worker writing its own disjoint back row (`rayon`
/// is forbidden in facett, ROOT LAW #0). Bit-identical to
/// [`step_serial`](Self::step_serial); worth its dispatch only on a very large
/// network, which is what [`step`](Self::step) decides.
pub fn step_gatling(&mut self, dt: f32) {
if self.pipes == 0 || dt <= 0.0 {
self.t += dt.max(0.0);
return;
}
self.t += dt;
self.integrate_velocity(dt);
let (cells, pipes) = (self.cells, self.pipes);
let front: &[f32] = &self.front;
let vel: &[f32] = &self.vel;
let back = RowsMut(self.back.as_mut_ptr());
gatling_run(pipes, 0, |p| {
// Force whole-struct capture of the `Sync` wrapper (not the raw `*mut`
// field) so the closure stays `Sync` under edition-2024 disjoint capture.
#[allow(clippy::redundant_locals)] // not redundant: it moves the WHOLE wrapper in (see above)
let back = back;
let src = &front[p * cells..p * cells + cells];
// SAFETY: pipe `p` writes only its own disjoint row `back[p*cells..]`; each
// `p` is claimed by exactly one gatling worker, so these writes never alias.
let dst = unsafe { std::slice::from_raw_parts_mut(back.0.add(p * cells), cells) };
advect_pipe(src, dst, vel[p], dt);
});
std::mem::swap(&mut self.front, &mut self.back);
self.exchange_junctions(dt);
}
/// The **serial kernel** — identical math, no fan-out, and the lane a pane actually
/// runs. The determinism test asserts [`step_gatling`](Self::step_gatling) equals this
/// bit-for-bit; since the advected pipes are independent that must hold, so a broken
/// fan-out (aliasing writes, a dropped pipe) diverges → RED. Also the honest 1-core /
/// wasm lane.
pub fn step_serial(&mut self, dt: f32) {
if self.pipes == 0 || dt <= 0.0 {
self.t += dt.max(0.0);
return;
}
self.t += dt;
self.integrate_velocity(dt);
let cells = self.cells;
for p in 0..self.pipes {
let (lo, hi) = (p * cells, p * cells + cells);
let src: &[f32] = &self.front[lo..hi];
let dst_ptr = self.back[lo..hi].as_mut_ptr();
// SAFETY: `dst` is a unique slice into `back`, disjoint from `src` (`front`).
let dst = unsafe { std::slice::from_raw_parts_mut(dst_ptr, cells) };
advect_pipe(src, dst, self.vel[p], dt);
}
std::mem::swap(&mut self.front, &mut self.back);
self.exchange_junctions(dt);
}
/// Gravity projected on every pipe's own axis, integrated with viscous drag.
fn integrate_velocity(&mut self, dt: f32) {
for p in 0..self.pipes {
let a = self.g_mag * dot(self.g, self.axis[p]);
let v = (self.vel[p] + a * dt) / (1.0 + DRAG * dt);
self.vel[p] = v.clamp(-MAX_VEL, MAX_VEL);
}
}
/// Junction pass: the **downhill** end of each pipe drains into its node tank, and
/// the **uphill** node tank feeds the pipe's uphill end. Both are paired transfers
/// of one scalar limited by source liquid AND destination free space, so volume is
/// conserved and no bound can be broken. Serial in pipe-index order (the scatter
/// touches shared nodes) — deterministic by construction.
fn exchange_junctions(&mut self, dt: f32) {
let cells = self.cells;
for p in 0..self.pipes {
let (na, nb) = self.ends[p];
let (na, nb) = (na as usize, nb as usize);
if na >= self.nodes || nb >= self.nodes {
continue;
}
let ga = self.g_mag * dot(self.g, self.axis[p]);
// |cos θ| between gravity and the pipe: a level pipe barely exchanges, a
// vertical one exchanges hard. Keeps the junction from pumping uphill.
let gate = if self.g_mag > 0.0 { (ga / self.g_mag).abs() } else { 0.0 };
let k = (JUNCTION_GAIN * gate * dt).clamp(0.0, 0.5);
if k <= 0.0 {
continue;
}
// Which end is downhill, in cell-index terms.
let (down_cell, down_node, up_cell, up_node) =
if ga >= 0.0 { (cells - 1, nb, 0usize, na) } else { (0usize, na, cells - 1, nb) };
// (a) the downhill end pours into its tank.
let di = p * cells + down_cell;
let room = (self.tank_cap[down_node] - self.tank[down_node]).max(0.0);
let x = (k * self.front[di]).min(self.front[di]).min(room);
if x > 0.0 {
self.front[di] -= x;
self.tank[down_node] += x;
}
// (b) the uphill tank feeds the pipe.
let ui = p * cells + up_cell;
let room = (1.0 - self.front[ui]).max(0.0);
let y = (k * self.tank[up_node]).min(self.tank[up_node]).min(room);
if y > 0.0 {
self.tank[up_node] -= y;
self.front[ui] += y;
}
}
}
/// The observable surface (for `state_json` / a robot proof): sizes, gravity, the
/// conserved volume, the containment bound, the per-pipe axial gravity + velocity,
/// and the tank levels.
#[must_use]
pub fn state_json(&self) -> serde_json::Value {
let axial: Vec<f32> = (0..self.pipes).map(|p| self.axial_gravity(p)).collect();
let levels: Vec<f32> = (0..self.nodes).map(|n| self.tank_level(n)).collect();
serde_json::json!({
"pipes": self.pipes,
"cells": self.cells,
"nodes": self.nodes,
"time": self.t,
"gravity": self.g,
"gravity_magnitude": self.g_mag,
"total_volume": self.total_volume(),
"max_fill": self.max_fill(),
"min_fill": self.min_fill(),
"fall_moment": self.fall_moment(),
"axial_gravity": axial,
"velocity": self.vel,
"tank_level": levels,
})
}
}
/// **One pipe's capacity-limited transport.** `src → dst` (same length), no allocation,
/// pure. Two effects, both expressed as a single signed face flux `Φ(j)` across face
/// `j|j+1` (positive = `j → j+1`), computed **entirely from `src`**:
///
/// - **advection** — bulk transport at the pipe's axial velocity (the Courant fraction
/// `c` of the upwind cell's liquid);
/// - **levelling** — liquid in a tube finds its own level, so a depth difference between
/// neighbours drives a flux from the deeper to the shallower at [`LEVEL_RELAX`].
/// Without it a pipe keeps its charged surface profile for ever.
///
/// Applying `dst[i] = src[i] − Φ(i) + Φ(i−1)` (with `Φ(−1) = Φ(n−1) = 0`) is then
///
/// - **conservative** — the sum telescopes; every `Φ` is added exactly where it was
/// subtracted, so no clamp below can create or destroy liquid, only refuse to move it;
/// - **bounded** — each face flux is limited to **half** the source cell's liquid and
/// **half** the destination cell's free space. A cell has two faces, so even when
/// levelling drains it through *both* at once (a local maximum — which advection alone
/// can never do, since a uniform velocity gives all fluxes one sign) the worst case is
/// `dst[i] ≥ src[i] − ½src[i] − ½src[i] = 0` and
/// `dst[i] ≤ src[i] + ½(1−src[i]) + ½(1−src[i]) = 1`. **That is why liquid can never
/// appear outside the tube**, and the halving is what keeps the bound inductive once
/// two-sided draining is possible.
#[inline]
fn advect_pipe(src: &[f32], dst: &mut [f32], vel: f32, dt: f32) {
let n = src.len();
if n == 0 {
return;
}
let c = (vel * dt).clamp(-MAX_COURANT, MAX_COURANT);
let relax = (LEVEL_RELAX * dt).clamp(0.0, 0.5);
// Φ(j), positive = from j to j+1.
let flux = |j: usize| -> f32 {
if j + 1 >= n {
return 0.0;
}
let (a, b) = (src[j], src[j + 1]);
// Bulk transport: the upwind cell donates the Courant fraction of its liquid.
let adv = if c >= 0.0 { c * a } else { -(-c * b) };
// Levelling: the deeper neighbour donates toward the shallower one.
let want = adv + relax * (a - b);
// Availability: at most half the source, at most half the destination's room.
if want >= 0.0 {
want.min(0.5 * a).min(0.5 * (1.0 - b).max(0.0))
} else {
-((-want).min(0.5 * b).min(0.5 * (1.0 - a).max(0.0)))
}
};
let mut inflow = 0.0f32; // Φ(i−1)
for i in 0..n {
let out = if i + 1 < n { flux(i) } else { 0.0 };
dst[i] = (src[i] - out + inflow).clamp(0.0, 1.0);
inflow = out;
}
}
#[inline]
fn dot(a: [f32; 3], b: [f32; 3]) -> f32 {
a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}
#[inline]
fn sub(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}
#[inline]
fn norm(a: [f32; 3]) -> f32 {
dot(a, a).sqrt()
}
#[inline]
fn unit(a: [f32; 3]) -> [f32; 3] {
let n = norm(a);
if n > 1e-6 { [a[0] / n, a[1] / n, a[2] / n] } else { [0.0, -1.0, 0.0] }
}
#[cfg(test)]
mod tests {
use super::*;
/// A little four-node ladder with a genuine height spread, so gravity has something
/// to do: two pipes down, one across, one back up.
fn ladder() -> PipeFlow {
let pos = vec![[0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]];
let edges = [(0usize, 1usize), (1, 2), (2, 3), (3, 0)];
let mut f = PipeFlow::new(pos, &edges, 32);
f.seed_fill(0xFACE_7700, 0.4);
f
}
// ── CONSERVATION — the sealed system never gains or loses liquid ────────────────
#[test]
fn volume_is_conserved_under_a_long_tumble() {
let mut f = ladder();
let v0 = f.total_volume();
assert!(v0 > 1.0, "the seeded charge is real liquid (v0={v0})");
// Tumble gravity right round twice while stepping — the hardest case for a
// conservative scheme (the flow direction reverses through every pipe).
for k in 0..2_400 {
let a = k as f32 * 0.005;
f.set_gravity([a.sin(), -a.cos(), 0.3 * (a * 0.7).sin()]);
f.step(1.0 / 60.0);
let v = f.total_volume();
assert!(
(v - v0).abs() / v0 < 1e-3,
"volume conserved at step {k}: {v} vs {v0} (drift {:.3e})",
(v - v0).abs() / v0
);
}
}
// ── CONTAINMENT — the RED-when-broken assertion the owner named ───────────────────
#[test]
fn liquid_never_leaves_the_pipe_interior() {
let mut f = ladder();
// Charge it hard: nearly full pipes AND full tanks, so any scheme that lets a
// cell over-fill will do it here.
f.seed_fill(0x51DE1, 0.95);
for k in 0..3_000 {
let a = k as f32 * 0.011;
f.set_gravity([a.sin(), -a.cos(), (a * 1.3).cos()]);
f.set_gravity_magnitude(9.0 + 6.0 * (a * 2.0).sin());
f.step(1.0 / 60.0);
let hi = f.max_fill();
let lo = f.min_fill();
assert!(hi <= 1.0 + 1e-6, "no cell exceeds its tube's cross-section at step {k} (max_fill={hi})");
assert!(lo >= -1e-6, "no cell holds negative liquid at step {k} (min_fill={lo})");
for n in 0..f.nodes() {
assert!(f.tank_level(n) <= 1.0 + 1e-6, "tank {n} never overflows at step {k}");
}
}
}
// ── GRAVITY — flow direction agrees with gravity projected on the pipe ──────────
#[test]
fn flow_direction_follows_gravity_projected_onto_each_pipe() {
// One pipe, from (0,1,0) down to (0,0,0): its axis is straight down. Charge a
// SLUG in the middle so the centre of mass has somewhere to go (a uniform column
// advects into itself and would prove nothing).
let mut f = PipeFlow::new(vec![[0.0, 1.0, 0.0], [0.0, 0.0, 0.0]], &[(0, 1)], 40);
f.drain();
f.set_tank_capacity(0, 0.0); // no tanks in this arm: pipe physics only
f.set_tank_capacity(1, 0.0);
for i in 16..24 {
f.set_fill(0, i, 1.0);
}
// Gravity down ⇒ axial gravity positive (from → to) ⇒ the slug moves toward the
// high cell index.
f.set_gravity([0.0, -1.0, 0.0]);
assert!(f.axial_gravity(0) > 0.0, "gravity down a downward pipe drives from→to");
let c0 = com(&f, 0);
for _ in 0..120 {
f.step(1.0 / 60.0);
}
let c1 = com(&f, 0);
assert!(c1 > c0 + 1.0, "liquid ran downhill (com {c0:.2} → {c1:.2})");
// Flip gravity: the SAME pipe must now drive the other way.
f.set_gravity([0.0, 1.0, 0.0]);
assert!(f.axial_gravity(0) < 0.0, "flipped gravity flips the axial component");
for _ in 0..300 {
f.step(1.0 / 60.0);
}
let c2 = com(&f, 0);
assert!(c2 < c1 - 1.0, "liquid ran back the other way (com {c1:.2} → {c2:.2})");
}
/// The **projection** is the physics: the same pipe, rotated, feels `|g|·cos θ`.
/// A 45° pipe must feel `1/√2` of what a vertical one feels, and a pipe tilted the
/// other way must feel the negation — this is the "as the camera rotates, the flow
/// re-aims" claim reduced to one arithmetic assertion.
#[test]
fn axial_gravity_is_the_cosine_projection_onto_the_pipe() {
let s = std::f32::consts::FRAC_1_SQRT_2;
let pos = vec![[0.0, 0.0, 0.0], [0.0, -1.0, 0.0], [s, -s, 0.0], [s, s, 0.0], [1.0, 0.0, 0.0]];
let f = {
let mut f = PipeFlow::new(pos, &[(0, 1), (0, 2), (0, 3), (0, 4)], 8);
f.set_gravity([0.0, -1.0, 0.0]);
f.set_gravity_magnitude(10.0);
f
};
assert!((f.axial_gravity(0) - 10.0).abs() < 1e-4, "straight down: the full |g|");
assert!((f.axial_gravity(1) - 10.0 * s).abs() < 1e-4, "45° down: |g|/√2");
assert!((f.axial_gravity(2) + 10.0 * s).abs() < 1e-4, "45° up: −|g|/√2 (drives to→from)");
assert!(f.axial_gravity(3).abs() < 1e-4, "level: nothing");
}
/// A pipe with **no** gravity component along it must not **pump**: its liquid may
/// level out, but it must not travel. Centre of mass and total are both pinned.
#[test]
fn a_level_pipe_does_not_pump() {
let mut f = PipeFlow::new(vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], &[(0, 1)], 24);
f.drain();
f.set_tank_capacity(0, 0.0); // no junction exchange in this arm
f.set_tank_capacity(1, 0.0);
for i in 0..24 {
f.set_fill(0, i, 0.5);
}
f.set_gravity([0.0, -1.0, 0.0]); // perpendicular to the pipe
assert!(f.axial_gravity(0).abs() < 1e-6, "gravity has no component along a level pipe");
let before: Vec<f32> = f.fill(0).to_vec();
for _ in 0..300 {
f.step(1.0 / 60.0);
}
assert_eq!(before.as_slice(), f.fill(0), "an already-level pipe's liquid just sits there");
}
/// **Liquid finds its own level.** Charge one pipe with a lumpy surface, hold it
/// level, and the surface flattens — without moving the liquid anywhere (the centre
/// of mass and the total are both unchanged). RED if the levelling flux is removed
/// (the profile freezes) or if it is not conservative (the total drifts).
#[test]
fn a_level_pipe_finds_its_own_level() {
let mut f = PipeFlow::new(vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], &[(0, 1)], 24);
f.drain();
f.set_tank_capacity(0, 0.0);
f.set_tank_capacity(1, 0.0);
f.set_gravity([0.0, -1.0, 0.0]); // level pipe: no axial forcing at all
// A mound in the middle, mirror-symmetric about cell 11.5 — so the level state
// it relaxes to has the SAME centre of mass, and any drift in the assertion
// below is levelling biasing the flow rather than the profile simply changing.
for i in 0..24 {
f.set_fill(0, i, if (i as f32 - 11.5).abs() < 4.0 { 0.8 } else { 0.3 });
}
let spread = |f: &PipeFlow| {
let d = f.fill(0);
let mean = d.iter().sum::<f32>() / d.len() as f32;
d.iter().map(|v| (v - mean) * (v - mean)).sum::<f32>()
};
let com = |f: &PipeFlow| {
let d = f.fill(0);
let (n, m): (f32, f32) = d.iter().enumerate().fold((0.0, 0.0), |(n, m), (i, &v)| (n + i as f32 * v, m + v));
n / m
};
let (s0, c0, v0) = (spread(&f), com(&f), f.total_volume());
for _ in 0..400 {
f.step(1.0 / 60.0);
}
let (s1, c1, v1) = (spread(&f), com(&f), f.total_volume());
assert!(s1 < s0 * 0.05, "the surface flattened (spread {s0:.4} → {s1:.4})");
assert!((c1 - c0).abs() < 1e-2, "levelling moves no liquid along the pipe (com {c0:.4} → {c1:.4})");
assert!((v1 - v0).abs() / v0 < 1e-4, "levelling is conservative ({v0:.5} → {v1:.5})");
}
// ── SETTLING — under fixed gravity the liquid ends up lower than it started ─────
#[test]
fn liquid_settles_downhill_and_comes_to_rest() {
let mut f = ladder();
f.set_gravity([0.0, -1.0, 0.0]);
let m0 = f.fall_moment();
for _ in 0..3_000 {
f.step(1.0 / 60.0);
}
let m1 = f.fall_moment();
// `fall_moment` is Σ vol·(ĝ·x): falling increases it.
assert!(m1 > m0 + 1e-3, "the liquid genuinely descended (fall moment {m0:.3} → {m1:.3})");
// And it stops: two consecutive frames are (nearly) the same field.
let a: Vec<f32> = f.fill(0).to_vec();
f.step(1.0 / 60.0);
let d: f32 = a.iter().zip(f.fill(0)).map(|(x, y)| (x - y).abs()).sum();
assert!(d < 1e-2, "the settled field is at rest (frame delta {d:.5})");
}
// ── the gatling fan-out equals the serial kernel, bit-for-bit ───────────────────
// Drives `step_gatling` EXPLICITLY: `step` dispatches on size, and this fixture is
// deliberately small, so calling `step` here would silently compare the serial lane
// with itself and prove nothing.
#[test]
fn gatling_step_matches_serial_step_bit_for_bit() {
let (mut g, mut s) = (ladder(), ladder());
for k in 0..600 {
let a = k as f32 * 0.01;
let dir = [a.sin(), -a.cos(), 0.2];
g.set_gravity(dir);
s.set_gravity(dir);
g.step_gatling(0.016);
s.step_serial(0.016);
}
assert_eq!(g.front, s.front, "gatling == serial (RED if the fan-out aliases or drops a pipe)");
assert_eq!(g.tank, s.tank, "the junction pass is order-deterministic");
}
/// `step` picks a lane by size and **must not change the answer doing it**: a small
/// network stepped through `step` matches the same network stepped through
/// `step_serial`, and a big one matches `step_gatling`.
#[test]
fn step_dispatches_on_size_without_changing_the_answer() {
let small = || {
let mut f = ladder();
for _ in 0..50 {
f.step(0.016);
}
f.front.clone()
};
let small_serial = || {
let mut f = ladder();
for _ in 0..50 {
f.step_serial(0.016);
}
f.front.clone()
};
const _: () = assert!(4 * 32 < PARALLEL_MIN_CELLS, "the ladder is below the threshold");
assert_eq!(small(), small_serial(), "below the threshold `step` IS the serial kernel");
// A network above the threshold: 12 000 pipes x 48 cells = 576 000 cells.
let big = |gatling: bool| {
let pos: Vec<[f32; 3]> = (0..4_000)
.map(|i| {
let a = i as f32 * 0.017;
[a.cos(), (a * 2.0).sin(), a.sin()]
})
.collect();
let edges: Vec<(usize, usize)> = (0..12_000).map(|i| (i % 4_000, (i * 7 + 3) % 4_000)).collect();
let mut f = PipeFlow::new(pos, &edges, 48);
f.seed_fill(5, 0.5);
for _ in 0..3 {
if gatling { f.step_gatling(0.016) } else { f.step(0.016) }
}
f.front.clone()
};
assert_eq!(big(false), big(true), "above the threshold `step` IS the gatling lane");
}
// ── determinism: same seed + same host dt ⇒ same field ──────────────────────────
#[test]
fn same_seed_and_same_dt_give_the_identical_field() {
let run = || {
let mut f = ladder();
for k in 0..400 {
let a = k as f32 * 0.013;
f.set_gravity([a.sin(), -a.cos(), 0.0]);
f.step(1.0 / 60.0);
}
(f.fill(0).to_vec(), f.total_volume(), f.fall_moment())
};
assert_eq!(run(), run(), "no wall clock, no RNG at runtime — the run is reproducible");
}
// ── ZERO-ALLOC guard: the double buffer ping-pongs in place ─────────────────────
#[test]
fn step_is_zero_alloc_double_buffer_ping_pong() {
let mut f = ladder();
let (pa, ca) = (f.front.as_ptr(), f.front.capacity());
let (pb, cb) = (f.back.as_ptr(), f.back.capacity());
for _ in 0..240 {
f.step(1.0 / 60.0);
}
assert_eq!(f.front.as_ptr(), pa, "front buffer reused in place (no realloc/clone)");
assert_eq!(f.back.as_ptr(), pb, "back buffer reused in place (no realloc/clone)");
assert_eq!(f.front.capacity(), ca);
assert_eq!(f.back.capacity(), cb);
let before = f.front.as_ptr();
f.step(1.0 / 60.0);
assert_eq!(f.front.as_ptr(), pb, "one step swaps the double buffer");
assert_ne!(f.front.as_ptr(), before);
}
#[test]
fn empty_and_degenerate_networks_are_safe() {
let mut f = PipeFlow::new(Vec::new(), &[], 8);
f.step(0.016); // no panic
assert_eq!(f.total_volume(), 0.0);
assert!(f.fill(0).is_empty());
assert_eq!(f.tank_level(3), 0.0);
// An edge pointing at a missing node must not index out of bounds.
let mut g = PipeFlow::new(vec![[0.0; 3]], &[(0, 9)], 8);
g.seed_fill(1, 0.5);
g.step(0.016);
assert!(g.max_fill() <= 1.0);
}
#[test]
fn state_json_publishes_the_invariants() {
let mut f = ladder();
for _ in 0..60 {
f.step(0.016);
}
let j = f.state_json();
assert_eq!(j["pipes"], 4);
assert_eq!(j["nodes"], 4);
assert!(j["total_volume"].as_f64().unwrap() > 0.0);
assert!(j["max_fill"].as_f64().unwrap() <= 1.0);
assert_eq!(j["axial_gravity"].as_array().unwrap().len(), 4);
assert_eq!(j["tank_level"].as_array().unwrap().len(), 4);
}
/// Fill-weighted centre of mass in cell-index units.
fn com(f: &PipeFlow, p: usize) -> f32 {
let d = f.fill(p);
let (mut num, mut den) = (0.0f32, 0.0f32);
for (i, &v) in d.iter().enumerate() {
num += i as f32 * v;
den += v;
}
if den > 0.0 { num / den } else { 0.0 }
}
}