rtx-cfd: Poiseuille closed-form validation with exact-zero assertions

Plane channel driven by a uniform body force: u(y) = G/(2mu) y(1-y), v = 0,
p exactly constant. Convection vanishes identically, so this isolates
diffusion, the half-cell wall treatment and the pressure coupling — and two
of the three answers are exact zeros, which no benchmark comparison offers.

The ends are clamped to the profile the DISCRETISATION prefers — the 1-D
tridiagonal with half-cell wall closures, solved directly in the test —
rather than to the continuous parabola. That makes (u_hat, 0, const) an
exact fixed point of the 2-D discretisation, and the solver must sit on it:

    |u - u_hat| ~ 1e-10,  max |v| ~ 1e-10,  p spread ~ 8e-10   (16^2)

A first version clamped the ends to the continuous parabola instead; the
O(h^2) incompatibility between that profile and the discrete one drove a
weak secondary flow near the ends (max |v| = 1.3e-3) — a property of the
mismatched boundary data, not of the solver, recorded in the test docs so
nobody rediscovers it as a bug.

The wall treatment's own truncation is measured in isolation as
|u_hat - parabola|: 3.906e-3 at 16, 9.766e-4 at 32 — refinement ratio
exactly 4.00, second order, in closed form c h^2 / 4.

Co-Authored-By: Claude Fable 5 <[email protected]>

crates/specialized/rtx-cfd/tests/poiseuille.rs

@@ -0,0 +1,206 @@
+//! Plane Poiseuille flow: the closed-form channel profile.
+//!
+//! A channel of unit height with no-slip walls, driven by a uniform body
+//! force `G` in x (equivalent to a constant pressure gradient, but it needs
+//! no pressure boundary conditions, so it fits the closed staggered box this
+//! solver provides). The exact steady solution is
+//!
+//! ```text
+//!     u(y) = G / (2 mu) * y (1 - y),    v = 0,    p = constant
+//! ```
+//!
+//! Convection vanishes identically — `u` depends only on `y` and `v = 0` —
+//! so this isolates exactly the operators the cavity cannot: diffusion, the
+//! half-cell wall treatment, and the pressure coupling.
+//!
+//! # The discrete profile, and why the ends are clamped to it
+//!
+//! The interior three-point Laplacian is *exact* on a parabola, but the wall
+//! rows are not: the half-cell flux `mu (u_0 - u_wall) / (h/2)` evaluated on
+//! the sampled parabola leaves a residual of `G/4` in the wall row, so the
+//! discretisation prefers a profile `u_hat` that differs from the parabola by
+//! `O(h^2)`, concentrated at the walls. `u_hat` solves a tridiagonal system
+//! (`-3u_0 + u_1 = -G h^2 / mu` at the walls, the standard second difference
+//! inside) that this test solves directly.
+//!
+//! Clamping the inlet and outlet to `u_hat` makes `(u_hat, 0, const)` an
+//! exact fixed point of the whole 2-D discretisation, which buys the sharp
+//! assertions: the solved field must reproduce `u_hat` to solver tolerance,
+//! the transverse velocity must vanish, and the pressure must be *flat* —
+//! any structure in it is spurious coupling. A first version of this test
+//! clamped the ends to the continuous parabola instead; the `O(h^2)`
+//! incompatibility between that profile and the discrete one drove a weak
+//! secondary flow near the ends (max |v| = 1.3e-3 at 16^2), which is a
+//! property of the mismatched boundary data, not of the solver.
+//!
+//! The distance between `u_hat` and the true parabola is then measured
+//! separately: it is the wall treatment's truncation error in isolation, and
+//! it must fall at second order.
+
+use rtx_cfd::solvers::incompressible::{
+    BoundaryConditions, FlowField, SimpleParameters, SimpleSolver,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+
+const MU: f64 = 0.1;
+const G: f64 = 0.8; // gives u_max = G / (8 mu) = 1 at mid-channel
+
+fn u_exact(y: f64) -> f64 {
+    G / (2.0 * MU) * y * (1.0 - y)
+}
+
+/// The profile the discretisation converges to: the 1-D discrete channel
+/// equation with half-cell wall closures, solved by the Thomas algorithm.
+fn discrete_profile(n: usize) -> Vec<f64> {
+    let h = 1.0 / n as f64;
+    let rhs_value = -G * h * h / MU;
+
+    // Tridiagonal: diag[j] u_j + upper u_{j+1} + lower u_{j-1} = rhs.
+    let mut diag = vec![-2.0; n];
+    diag[0] = -3.0;
+    diag[n - 1] = -3.0;
+    let mut rhs = vec![rhs_value; n];
+
+    // Forward elimination (sub- and super-diagonals are all 1).
+    let mut upper = vec![1.0; n];
+    for j in 1..n {
+        let factor = 1.0 / diag[j - 1];
+        diag[j] -= factor * upper[j - 1];
+        rhs[j] -= factor * rhs[j - 1];
+    }
+    let mut u = vec![0.0; n];
+    u[n - 1] = rhs[n - 1] / diag[n - 1];
+    for j in (0..n - 1).rev() {
+        u[j] = (rhs[j] - upper[j] * u[j + 1]) / diag[j];
+    }
+    u
+}
+
+struct Measurement {
+    /// Solved field against the discrete profile — solver truncation only.
+    max_u_vs_discrete: f64,
+    /// Discrete profile against the closed form — wall truncation only.
+    max_discrete_vs_exact: f64,
+    max_v: f64,
+    p_spread: f64,
+}
+
+async fn measure(n: usize) -> CfdResult<Measurement> {
+    let dx = 1.0 / n as f64;
+    let dy = dx;
+    let u_hat = discrete_profile(n);
+
+    let config = CfdConfig::new()
+        .with_density(1.0)
+        .with_viscosity(MU)
+        .with_reference_velocity(1.0)
+        .with_reference_length(1.0);
+    let params = SimpleParameters::default()
+        .with_max_iterations(40000)
+        .with_tolerance(1e-10);
+    let mut solver = SimpleSolver::new(config, params)?;
+    solver.set_momentum_source(|_x, _y| (G, 0.0));
+    solver.set_wall_velocity(|_x, _y| (0.0, 0.0));
+
+    let mut field = FlowField::new(n, n, dx, dy)?;
+
+    // Inlet and outlet carry the discrete profile; the v faces on the walls
+    // stay at zero, which is what `FlowField::new` initialises.
+    for (j, &u_hat_j) in u_hat.iter().enumerate() {
+        field.u[(j, 0)] = u_hat_j;
+        field.u[(j, n)] = u_hat_j;
+    }
+
+    let empty = BoundaryConditions::new();
+    for _ in 0..40000 {
+        let (mass, momentum) = solver
+            .solve_simple_iteration(&mut field, &empty, 0.01)
+            .await?;
+        if (mass * mass + momentum * momentum).sqrt() < 1e-10 {
+            break;
+        }
+    }
+
+    let mut max_u_vs_discrete: f64 = 0.0;
+    for (j, &u_hat_j) in u_hat.iter().enumerate() {
+        for i in 1..n {
+            max_u_vs_discrete = max_u_vs_discrete.max((field.u[(j, i)] - u_hat_j).abs());
+        }
+    }
+
+    let mut max_discrete_vs_exact: f64 = 0.0;
+    for (j, &u_hat_j) in u_hat.iter().enumerate() {
+        let y = (j as f64 + 0.5) * dy;
+        max_discrete_vs_exact = max_discrete_vs_exact.max((u_hat_j - u_exact(y)).abs());
+    }
+
+    let mut max_v: f64 = 0.0;
+    for j in 1..n {
+        for i in 0..n {
+            max_v = max_v.max(field.v[(j, i)].abs());
+        }
+    }
+
+    let mut p_min = f64::INFINITY;
+    let mut p_max = f64::NEG_INFINITY;
+    for j in 0..n {
+        for i in 0..n {
+            p_min = p_min.min(field.p[(j, i)]);
+            p_max = p_max.max(field.p[(j, i)]);
+        }
+    }
+
+    Ok(Measurement {
+        max_u_vs_discrete,
+        max_discrete_vs_exact,
+        max_v,
+        p_spread: p_max - p_min,
+    })
+}
+
+#[tokio::test]
+async fn poiseuille_profile_matches_the_closed_form() -> CfdResult<()> {
+    let m16 = measure(16).await?;
+    let m32 = measure(32).await?;
+
+    for (n, m) in [(16, &m16), (32, &m32)] {
+        println!(
+            "  n = {n:2}   |u - u_hat| = {:.3e}   |u_hat - exact| = {:.3e}   \
+             max |v| = {:.3e}   p spread = {:.3e}",
+            m.max_u_vs_discrete, m.max_discrete_vs_exact, m.max_v, m.p_spread
+        );
+    }
+    println!(
+        "  wall-truncation refinement ratio = {:.2}",
+        m16.max_discrete_vs_exact / m32.max_discrete_vs_exact
+    );
+
+    // (u_hat, 0, const) is an exact fixed point of the discretisation, so
+    // the solved field must sit on it to solver tolerance — these three
+    // assertions have exact-zero answers, and any structure is a defect.
+    for m in [&m16, &m32] {
+        assert!(
+            m.max_u_vs_discrete < 1e-7,
+            "u departs from the discrete profile by {:.3e}",
+            m.max_u_vs_discrete
+        );
+        assert!(m.max_v < 1e-7, "spurious transverse flow {:.3e}", m.max_v);
+        assert!(m.p_spread < 1e-6, "spurious pressure {:.3e}", m.p_spread);
+    }
+
+    // The wall treatment's own truncation: second order, so the 16 -> 32
+    // ratio must sit near 4.
+    let ratio = m16.max_discrete_vs_exact / m32.max_discrete_vs_exact;
+    assert!(
+        m16.max_discrete_vs_exact < 5e-3,
+        "wall truncation {:.3e} at n = 16 on a profile of peak 1.0",
+        m16.max_discrete_vs_exact
+    );
+    assert!(
+        (3.3..4.7).contains(&ratio),
+        "wall truncation refines at ratio {ratio:.2}, not the ~4 of a \
+         second-order treatment"
+    );
+
+    Ok(())
+}