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

//! embedded3 S2-3: the free-ended flag with prescribed motion — the first
//! honest 3D wake. The Turek–Hron channel (2.5 × 0.41) extruded to depth
//! 0.41 with the cylinder (D 0.1 at (0.2, 0.2)) across the width, the flag
//! 0.35 × 0.02 × 0.2 centred in z (z 0.105–0.305), its centreline deflected
//! as the first clamped-free beam mode with the 2D FSI2 flat-tip record's
//! tip amplitude 84 mm at 1.930 Hz (motion prescribed, no structure), the
//! flag's span edges rounded to one cell and its tip a semicircle (the
//! linear cut geometry needs smooth edges; disclosed). Inflow parabolic
//! in y and z with U_m 2.25 (Ū 1.0 = FSI2's mean, Re 100 on D), ρ 1000,
//! ν 1e-3. CutCell wall with merging, TVD, moving body on the device
//! (S2-2a: the host rebuild per step).
//!
//! Gate (`docs/embedded3_campaign.md` S2-3): ny 62, two full periods, no
//! death, mass residual ≤ 1e-8 every step; the mid-plane per-span loads
//! within 30 % of the 2D FSI2 record (drag mean 224.6 N/m, lift swing
//! ±215 flat tip / ±256 semicircle); 32 VTK phases of the last period.
//!
//! `RTX_E3_FLAG_NY=62 RTX_E3_FLAG_PERIODS=2 RTX_E3_FLAG_VTK=<dir> RTX_E3_FLAG_CSV=<path> \
//!  RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_flag_wake -- --ignored --nocapture`
#![cfg(feature = "cuda")]

use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, write_vtk,
};
use std::io::Write as _;

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;
const FLAG_X0: f64 = 0.6;
const FLAG_LEN: f64 = 0.35;
const FLAG_HALF: f64 = 0.01;
const FLAG_SPAN: f64 = 0.2;
const AMP: f64 = 0.084;
const FREQ: f64 = 1.930;
const U_M: f64 = 2.25;
const RHO: f64 = 1000.0;
const NU: f64 = 1e-3;
/// The first clamped-free beam mode's `β L`.
const BETA_L: f64 = 1.875_104_069;

fn env_f(name: &str, default: f64) -> f64 {
    std::env::var(name)
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(default)
}

/// The first mode shape normalised to 1 at the tip, `s ∈ [0, 1]`.
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)
}

/// Centreline deflection and its velocity at arc parameter `s`, time `t`.
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(),
    )
}

/// Signed distance to the deflected flag's cross-section (a capsule
/// around the centreline polyline of `n` segments) and the centreline's
/// transverse velocity at the closest point.
fn flag_2d(x: f64, y: f64, t: f64) -> (f64, f64) {
    let n = 40;
    let mut best = f64::INFINITY;
    let mut v_best = 0.0;
    let point = |m: usize| {
        let s = m as f64 / n as f64;
        let (d, v) = deflection(s, t);
        (FLAG_X0 + s * FLAG_LEN, CY + d, v)
    };
    for m in 0..n {
        let (ax, ay, av) = point(m);
        let (bx, by, bv) = point(m + 1);
        let (ex, ey) = (bx - ax, by - ay);
        let l2 = ex * ex + ey * ey;
        let u = (((x - ax) * ex + (y - ay) * ey) / l2).clamp(0.0, 1.0);
        let (px, py) = (ax + u * ex, ay + u * ey);
        let d = ((x - px).powi(2) + (y - py).powi(2)).sqrt();
        if d < best {
            best = d;
            v_best = av + u * (bv - av);
        }
    }
    (best - FLAG_HALF, v_best)
}

/// The flag in 3D: the extruded capsule cut to the span with edges
/// rounded to radius `r`.
fn flag_3d(x: f64, y: f64, z: f64, t: f64, r: f64) -> (f64, f64) {
    let (d2, v) = flag_2d(x, y, t);
    let zc = 0.5 * H;
    let q1 = d2 + r;
    let q2 = (z - zc).abs() - 0.5 * FLAG_SPAN + r;
    let outside = (q1.max(0.0).powi(2) + q2.max(0.0).powi(2)).sqrt();
    (outside + q1.max(q2).min(0.0) - r, v)
}

fn inflow(y: f64, z: f64) -> f64 {
    16.0 * U_M * y * z * (H - y) * (H - z) / (H * H * H * H)
}

#[test]
#[ignore = "S2-3: the flag wake on the device (about an hour at ny 62)"]
fn flag_wake_on_the_device() {
    let ny = env_f("RTX_E3_FLAG_NY", 62.0) as usize;
    let periods = env_f("RTX_E3_FLAG_PERIODS", 2.0);
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let nz = ny;
    let r_edge = h;
    let dt_cfl = 0.3 * h / (U_M.max(2.0 * std::f64::consts::PI * FREQ * AMP));
    let dt = dt_cfl.min(0.5 * h * h / (6.0 * NU));
    let period = 1.0 / FREQ;
    let t_end = periods * period;
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: RHO * NU,
            reference_velocity: 1.0,
            reference_length: 2.0 * R_CYL,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            convection_scheme: ConvectionScheme::TvdVanAlbada,
            wall_scheme: WallScheme::CutCell,
            boundaries: Boundaries {
                x1: Side::PressureOutlet,
                ..Boundaries::default()
            },
            ..Parameters::default()
        },
    );
    solver.set_boundary_velocity(|x, y, z, _t| {
        if x <= 0.0 {
            (inflow(y, z), 0.0, 0.0)
        } else {
            (0.0, 0.0, 0.0)
        }
    });
    let cyl = move |x: f64, y: f64| ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL;
    let body = Body::from_sdf(move |x, y, z, t| cyl(x, y).min(flag_3d(x, y, z, t, r_edge).0))
        .with_surface_velocity(move |x, y, z, t| {
            let (df, v) = flag_3d(x, y, z, t, r_edge);
            if df <= cyl(x, y) {
                (0.0, v, 0.0)
            } else {
                (0.0, 0.0, 0.0)
            }
        });
    solver.set_moving_body(body);
    let g = Grid::cubic(nx, ny, nz, h);
    let mut field = Field::new(g);
    for k in 0..nz {
        for j in 0..ny {
            let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
            for i in 0..=nx {
                field.u[g.uface(k, j, i)] = u0;
            }
        }
    }
    solver.initialize(&mut field);
    println!(
        "  flag wake ny {ny}: {nx}×{ny}×{nz} = {} cells, h {h:.4e}, dt {dt:.3e}, {periods} periods = {t_end:.3} s, {} steps",
        g.cells(),
        (t_end / dt).ceil() as usize
    );
    unsafe { std::env::set_var("RTX_PROFILE", "1") };
    let mut device = DeviceStep::new(solver, g);
    device.upload(&field);
    let steps = (t_end / dt).ceil() as usize;
    let mut csv = std::env::var("RTX_E3_FLAG_CSV").ok().map(|p| {
        let mut f = std::fs::File::create(p).expect("csv");
        writeln!(
            f,
            "t,tip,drag_span,lift_span,drag_total,lift_total,residual,cg,fresh"
        )
        .unwrap();
        f
    });
    let vtk_dir = std::env::var("RTX_E3_FLAG_VTK").ok();
    let phases = 32;
    let last_period_start = t_end - period;
    let mut next_phase = 0;
    let mid = nz / 2;
    let slab = (mid - 2, mid + 2);
    let start = std::time::Instant::now();
    let (mut drag_sum, mut lift_min, mut lift_max, mut samples) =
        (0.0, f64::INFINITY, f64::NEG_INFINITY, 0usize);
    let mut worst_residual = 0.0_f64;
    for step in 0..steps {
        let r = device.advance(dt);
        worst_residual = worst_residual.max(r.final_residual);
        assert!(r.final_residual.is_finite(), "death at step {step}");
        let t = device.solver.time();
        let sample = (step + 1) % 10 == 0 || step + 1 == steps;
        let phase_due = vtk_dir.is_some()
            && t >= last_period_start + next_phase as f64 * period / phases as f64
            && next_phase < phases;
        if sample || phase_due {
            device.download(&mut field);
            let solver = &device.solver;
            let mask = solver.mask().expect("mask");
            let body = solver.body().expect("body");
            let fs = mask
                .cut_wall_force_per_span(body, &field, RHO * NU, t, slab)
                .expect("wall");
            let ft = mask
                .cut_wall_force(body, &field, RHO * NU, t)
                .expect("wall");
            let tip = deflection(1.0, t).0;
            if sample {
                println!(
                    "  t {t:7.4} (tip {tip:+.4}): drag/span {:.1} lift/span {:+.1} N/m; total {:.3} {:+.3} N; residual {:.1e} CG {} fresh {}; [{:.0} s]",
                    fs[0],
                    fs[1],
                    ft[0],
                    ft[1],
                    r.final_residual,
                    r.poisson_iterations,
                    r.fresh_cells,
                    start.elapsed().as_secs_f64()
                );
                if let Some(f) = csv.as_mut() {
                    writeln!(
                        f,
                        "{t:.5},{tip:.5},{:.4},{:.4},{:.5},{:.5},{:.3e},{},{}",
                        fs[0],
                        fs[1],
                        ft[0],
                        ft[1],
                        r.final_residual,
                        r.poisson_iterations,
                        r.fresh_cells
                    )
                    .unwrap();
                }
                if t >= last_period_start {
                    drag_sum += fs[0];
                    lift_min = lift_min.min(fs[1]);
                    lift_max = lift_max.max(fs[1]);
                    samples += 1;
                }
            }
            if phase_due {
                let path = std::path::Path::new(vtk_dir.as_ref().unwrap())
                    .join(format!("flag_ny{ny}_phase{next_phase:02}.vtk"));
                write_vtk(&path, &field, Some(mask)).expect("vtk");
                next_phase += 1;
            }
        }
    }
    let drag_mean = drag_sum / samples.max(1) as f64;
    println!(
        "  FINAL ny {ny}: last period drag/span mean {drag_mean:.1} N/m (2D FSI2 224.6), lift/span {lift_min:+.1} … {lift_max:+.1} (2D ±215 flat tip, ±256 semicircle); worst residual {worst_residual:.1e}; {} phases written; {:.0} s",
        next_phase,
        start.elapsed().as_secs_f64()
    );
    if let Some(t) = device.timers() {
        println!("  timers: {t:?}");
    }
}