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

//! S2-5 instrument: where does the cut wall sit? Poiseuille flow along z
//! (periodic, body force `f`) between an EMBEDDED flat wall at `y = y_w`
//! (cut at a chosen fraction θ of a cell) and the domain's top wall. The
//! flow rate per unit width is `f (H − y_w)³ / (12 μ)`, so the measured
//! rate gives the effective wall position `y_eff`; the offset
//! `(y_eff − y_w)/h` must vanish at second order and is the number the
//! DFG ladder reads as an effective radius ≈ 0.2 h short.
use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
};

const MU: f64 = 0.1;
const F: f64 = 1.0;
const HY: f64 = 1.0;

fn offset(ny: usize, theta: f64) -> (f64, Vec<(f64, f64, f64)>) {
    let h = HY / ny as f64;
    let (nx, nz) = (3 * ny, 2);
    let y_w = (ny as f64 / 4.0).floor() * h + theta * h;
    let exact = move |y: f64| {
        if y > y_w {
            F / (2.0 * MU) * (y - y_w) * (HY - y)
        } else {
            0.0
        }
    };
    let mut solver = Solver::new(
        Fluid {
            density: 1.0,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-10,
            convection_scheme: ConvectionScheme::Upwind,
            wall_scheme: WallScheme::CutCell,
            boundaries: Boundaries {
                z0: Side::Periodic,
                z1: Side::Periodic,
                ..Boundaries::default()
            },
            ..Parameters::default()
        },
    );
    solver.set_boundary_velocity(move |_x, y, _z, _t| (0.0, 0.0, exact(y)));
    solver.set_momentum_source(|_, _, _, _| (0.0, 0.0, F));
    solver.set_body(Body::from_sdf(move |_x, y, _z, _t| y - y_w));
    let g = Grid::cubic(nx, ny, nz, h);
    let mut field = Field::new(g);
    for k in 0..=nz {
        for j in 0..ny {
            for i in 0..nx {
                let idx = g.wface(k, j, i);
                field.w[idx] = exact((j as f64 + 0.5) * h);
            }
        }
    }
    solver.initialize(&mut field);
    let dt = 0.5 * h * h / (6.0 * MU);
    let steps = (3.0 / dt).ceil() as usize;
    for _ in 0..steps {
        solver.advance(&mut field, dt);
    }
    let mask = solver.mask().expect("mask");
    let i = nx / 2;
    let mut q = 0.0;
    let mut profile = Vec::new();
    for j in 0..ny {
        let f = g.wface(0, j, i);
        let a = mask.a_w(f);
        q += a * field.w[f] * h;
        let y = (j as f64 + 0.5) * h;
        if a > 0.0 && profile.len() < 4 {
            profile.push((a, field.w[f], exact(y)));
        }
    }
    let y_eff = HY - (12.0 * MU * q / F).cbrt();
    ((y_eff - y_w) / h, profile)
}

#[test]
#[ignore = "S2-5 instrument: the cut wall's effective position on a flat wall (a minute on the host)"]
fn flat_wall_effective_position() {
    for ny in [16usize, 32] {
        for theta in [0.05, 0.25, 0.5, 0.75, 0.95] {
            let (off, profile) = offset(ny, theta);
            let p: Vec<String> = profile
                .iter()
                .map(|(a, w, e)| format!("α {a:.2} w {w:.5} (exact at the face centre {e:.5})"))
                .collect();
            println!(
                "  ny {ny} θ {theta:.2}: effective wall offset {off:+.4} h (positive = the wall sits inside the fluid); first open faces: {}",
                p.join("; ")
            );
        }
    }
}

/// The gate: with the centroid diffusion (the default) the effective wall
/// position is second order — within 0.08 h at ny 16, halving in units of
/// h at ny 32, the same at every cut fraction (before S2-5: −½θ h at
/// every resolution).
#[test]
fn flat_wall_position_is_second_order() {
    if std::env::var("RTX_E3_DIFFUSION_CENTROID").is_ok_and(|v| v == "0") {
        return;
    }
    for theta in [0.25, 0.75] {
        let (coarse, _) = offset(16, theta);
        let (fine, _) = offset(32, theta);
        println!("  θ {theta}: offset {coarse:+.4} h at ny 16, {fine:+.4} h at ny 32");
        assert!(coarse.abs() < 0.08, "ny 16 offset {coarse}");
        assert!(
            fine.abs() < 0.6 * coarse.abs(),
            "the offset does not halve: {coarse} → {fine}"
        );
    }
}