From f2cdda46910a38e73116036662e5ade8fe66ddbc Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 23:14:42 -0500 Subject: [PATCH] embedded3 S2-3c: the 2D solver on the flag wake's kinematics (the reference for the 3D full-span run) Co-Authored-By: Claude Fable 5.1 --- .../tests/embedded3_flag_reference_2d.rs | 240 ++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs new file mode 100644 index 0000000..24c3ef2 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs @@ -0,0 +1,240 @@ +//! 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; +const FLAG_X0: f64 = 0.6; +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 = (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::() / 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(()) +}