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

//! S2-3c: the 2D solver on the flag wake's exact kinematics — the same
//! cylinder + capsule flag, the same first-mode motion (tip 84 mm at
//! 1.930 Hz), the same inflow mean (Ū 1.0, the 2D parabola), TVD, on the
//! ny 62 spacing (378 × 62, h 6.6 mm). Its loads are the 2D answer for
//! this kinematics; the 3D full-span run must reproduce them. Loads by
//! the traction sampler over the flag's and the cylinder's surface
//! samples per unit span, every 10 steps; the last period's mean drag and
//! median-filtered lift swing.
//!
//! `cargo test --release -p rtx-cfd --test embedded3_flag_reference_2d -- --ignored --nocapture`
use rtx_cfd::solvers::incompressible::{
    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
    FlowField, MgPrecision, MgSmoother, PoissonSolverKind, SideBoundary,
};
use rtx_cfd::{CfdConfig, CfdResult};

const H: f64 = 0.41;
const L: f64 = 2.5;
const CX: f64 = 0.2;
const CY: f64 = 0.2;
const R_CYL: f64 = 0.05;
/// The flag's ROOT: the cylinder's rear (the Turek–Hron flag runs from the
/// cylinder to its tip A at x = 0.6). Until 2026-09-18 this was 0.6 — the
/// flag sat DETACHED, its root where the benchmark's tip is; every flag
/// record before that date is of that geometry.
const FLAG_X0: f64 = 0.25;
const FLAG_LEN: f64 = 0.35;
const FLAG_HALF: f64 = 0.01;
const AMP: f64 = 0.084;
const FREQ: f64 = 1.930;
const U_MEAN: f64 = 1.0;
const RHO: f64 = 1000.0;
const NU: f64 = 1e-3;
const BETA_L: f64 = 1.875_104_069;

fn mode(s: f64) -> f64 {
    let b = BETA_L;
    let sigma = (b.sinh() - b.sin()) / (b.cosh() + b.cos());
    let phi = |s: f64| (b * s).cosh() - (b * s).cos() - sigma * ((b * s).sinh() - (b * s).sin());
    phi(s) / phi(1.0)
}

fn deflection(s: f64, t: f64) -> (f64, f64) {
    let w = 2.0 * std::f64::consts::PI * FREQ;
    (
        AMP * mode(s) * (w * t).sin(),
        AMP * mode(s) * w * (w * t).cos(),
    )
}

fn centreline(m: usize, n: usize, t: f64) -> (f64, f64, f64) {
    let s = m as f64 / n as f64;
    let (d, v) = deflection(s, t);
    (FLAG_X0 + s * FLAG_LEN, CY + d, v)
}

/// Distance to the capsule flag and the centreline velocity at the foot.
fn flag_sdf(x: f64, y: f64, t: f64) -> (f64, f64) {
    let n = 40;
    let mut best = f64::INFINITY;
    let mut v_best = 0.0;
    for m in 0..n {
        let (ax, ay, av) = centreline(m, n, t);
        let (bx, by, bv) = centreline(m + 1, n, t);
        let (ex, ey) = (bx - ax, by - ay);
        let u = (((x - ax) * ex + (y - ay) * ey) / (ex * ex + ey * ey)).clamp(0.0, 1.0);
        let d = ((x - ax - u * ex).powi(2) + (y - ay - u * ey).powi(2)).sqrt();
        if d < best {
            best = d;
            v_best = av + u * (bv - av);
        }
    }
    (best - FLAG_HALF, v_best)
}

fn cyl_sdf(x: f64, y: f64) -> f64 {
    ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL
}

/// Surface samples `(x, y, nx, ny, ds)` at `t`: the flag's two sides and
/// tip along the deflected centreline, the cylinder's circle; samples
/// inside the other body are dropped.
fn samples(t: f64, ds: f64) -> Vec<(f64, f64, f64, f64, f64)> {
    let mut out = Vec::new();
    let n = ((FLAG_LEN / ds).ceil() as usize).max(8);
    for m in 0..n {
        let (ax, ay, _) = centreline(m, n, t);
        let (bx, by, _) = centreline(m + 1, n, t);
        let (ex, ey) = (bx - ax, by - ay);
        let len = (ex * ex + ey * ey).sqrt();
        let (tx, ty) = (ex / len, ey / len);
        let (nx, ny) = (-ty, tx);
        let (mx, my) = (0.5 * (ax + bx), 0.5 * (ay + by));
        for sign in [1.0, -1.0] {
            let (px, py) = (mx + sign * FLAG_HALF * nx, my + sign * FLAG_HALF * ny);
            if cyl_sdf(px, py) > 0.0 {
                out.push((px, py, sign * nx, sign * ny, len));
            }
        }
    }
    // The tip: a semicircle around the last centreline point.
    let (tx0, ty0, _) = centreline(n, n, t);
    let (px, py, _) = centreline(n - 1, n, t);
    let ang0 = (ty0 - py).atan2(tx0 - px);
    let n_arc = ((std::f64::consts::PI * FLAG_HALF / ds).ceil() as usize).max(4);
    for k in 0..n_arc {
        let a = ang0 - std::f64::consts::FRAC_PI_2
            + (k as f64 + 0.5) / n_arc as f64 * std::f64::consts::PI;
        out.push((
            tx0 + FLAG_HALF * a.cos(),
            ty0 + FLAG_HALF * a.sin(),
            a.cos(),
            a.sin(),
            std::f64::consts::PI * FLAG_HALF / n_arc as f64,
        ));
    }
    let n_c = ((2.0 * std::f64::consts::PI * R_CYL / ds).ceil() as usize).max(16);
    for k in 0..n_c {
        let a = (k as f64 + 0.5) / n_c as f64 * 2.0 * std::f64::consts::PI;
        let (px, py) = (CX + R_CYL * a.cos(), CY + R_CYL * a.sin());
        if flag_sdf(px, py, t).0 > 0.0 {
            out.push((
                px,
                py,
                a.cos(),
                a.sin(),
                2.0 * std::f64::consts::PI * R_CYL / n_c as f64,
            ));
        }
    }
    out
}

fn median(v: &[f64]) -> f64 {
    let mut s = v.to_vec();
    s.sort_by(|a, b| a.partial_cmp(b).unwrap());
    s[s.len() / 2]
}

#[tokio::test]
#[ignore = "S2-3c: the 2D reference for the flag wake's kinematics (minutes on the host)"]
async fn flag_wake_2d_reference() -> CfdResult<()> {
    let ny = 62;
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let dt = 8.817e-4;
    let period = 1.0 / FREQ;
    let t_end = 2.0 * period;
    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(RHO * NU)
        .with_reference_velocity(U_MEAN)
        .with_reference_length(2.0 * R_CYL);
    let mut solver = EmbeddedPisoSolver::new(
        config,
        EmbeddedParameters {
            corrector_steps: 3,
            tolerance: 1e-8,
            boundaries: AleBoundaries {
                right: SideBoundary::PressureOutlet,
                ..AleBoundaries::default()
            },
            poisson_solver: PoissonSolverKind::Multigrid,
            poisson_precision: MgPrecision::F64,
            poisson_smoother: MgSmoother::Lexicographic,
            convection_scheme: ConvectionScheme::TvdVanAlbada,
        },
    )?;
    solver.set_boundary_velocity(|x, y, _t| {
        if x <= 0.0 {
            (1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2), 0.0)
        } else {
            (0.0, 0.0)
        }
    });
    let flag = EmbeddedBody::from_sdf(|x, y, t| flag_sdf(x, y, t).0)
        .with_surface_velocity(|x, y, t| (0.0, flag_sdf(x, y, t).1));
    solver.set_moving_body(EmbeddedBody::union(
        EmbeddedBody::circle(CX, CY, R_CYL),
        flag,
    ));
    let mut field = FlowField::new(nx, ny, h, h)?;
    for j in 0..ny {
        let u0 =
            1.5 * U_MEAN * ((j as f64 + 0.5) * h) * (H - (j as f64 + 0.5) * h) / (0.5 * H).powi(2);
        for i in 0..=nx {
            field.u[(j, i)] = u0;
        }
    }
    solver.initialize(&mut field)?;
    let steps = (t_end / dt).ceil() as usize;
    println!("  2D reference: {nx}×{ny}, h {h:.4e}, dt {dt:.3e}, {steps} steps");
    let mu = RHO * NU;
    let start = std::time::Instant::now();
    let mut drag = Vec::new();
    let mut lift = Vec::new();
    let mut worst = 0.0_f64;
    for step in 0..steps {
        let r = solver.advance(&mut field, dt).await?;
        worst = worst.max(r.solver_result.final_residual);
        let t = (step + 1) as f64 * dt;
        if (step + 1) % 10 == 0 {
            let mask = solver.mask().expect("mask");
            let body = solver.body().expect("body");
            let (mut fx, mut fy, mut skipped) = (0.0, 0.0, 0usize);
            for (x, y, nx_, ny_, ds) in samples(t, 0.5 * h) {
                match mask.traction_at(body, &field.u, &field.v, &field.p, mu, t, x, y, nx_, ny_) {
                    Some((tx, ty)) => {
                        fx += tx * ds;
                        fy += ty * ds;
                    }
                    None => skipped += 1,
                }
            }
            if t >= t_end - period {
                drag.push(fx);
                lift.push(fy);
            }
            if (step + 1) % 100 == 0 {
                println!(
                    "  t {t:6.3} tip {:+.4}: drag {fx:7.1} lift {fy:+8.1} N/m (skipped {skipped}); residual {:.1e}; [{:.0} s]",
                    deflection(1.0, t).0,
                    r.solver_result.final_residual,
                    start.elapsed().as_secs_f64()
                );
            }
        }
    }
    let filt: Vec<f64> = (0..lift.len())
        .map(|i| median(&lift[i.saturating_sub(5)..(i + 6).min(lift.len())]))
        .collect();
    let lo = filt.iter().cloned().fold(f64::INFINITY, f64::min);
    let hi = filt.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    println!(
        "  FINAL 2D reference: last period drag mean {:.1} N/m, lift swing {:+.1} … {:+.1} (raw {:+.1} … {:+.1}); worst residual {worst:.1e}; {:.0} s",
        drag.iter().sum::<f64>() / drag.len() as f64,
        lo,
        hi,
        lift.iter().cloned().fold(f64::INFINITY, f64::min),
        lift.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
        start.elapsed().as_secs_f64()
    );
    Ok(())
}