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

//! 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(())
}