overset: the energy pin (P5-3) — the rigid cylinder + flag outline pitching about the cylinder's centre in still fluid on the moving patch at the FSI2 fluid and grid; per-cycle net fluid work on the body (must be ≤ 0 at every h), gross +/− split, and the last-two-period moment fit → added inertia + rotational damping; knobs RTX_OVERSET_NY / _EP_PERIODS / _EP_TIP / _EP_F; ny 41 one-period smoke: net −5.62 J/m, damping 56.9 N·m·s/m (π c θ0² ω = 5.6 J, self-consistent), 21 s

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01LzcjQX7tvgn87CQCyg9Cfr

crates/specialized/rtx-cfd/tests/overset_energy_pin.rs

@@ -0,0 +1,222 @@
+//! P5-3's ENERGY pin: the rigid cylinder + flag outline PITCHING about
+//! the cylinder's centre in still fluid, on the moving patch, at the FSI2
+//! fluid (ρ 1000, ν 1e-3) and the FSI2 background grid — the fluid's net
+//! work on the body per cycle must be ≤ 0 (viscous) at every h. No wake
+//! ambiguity, no coupling: a positive net work anywhere on the ladder
+//! locates an energy source in the composite's moving-wall solution.
+//!
+//! Registered before the run (2026-09-13): |W| per cycle negative at
+//! ny = 41 / 62 / 82 and shrinking slowly with h (a better-resolved Stokes
+//! layer); the added-mass-like moment coefficient (in phase with θ̈)
+//! should converge, not drift like the cylinder pin's C_m (1.069 / 1.097 /
+//! 1.127 at n 32 / 64 / 128).
+//!
+//! `RTX_OVERSET_NY` (41), `RTX_OVERSET_EP_PERIODS` (4), `RTX_OVERSET_EP_TIP`
+//! (tip amplitude, m, 0.02), `RTX_OVERSET_EP_F` (Hz, 2.0).
+
+use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
+use rtx_cfd::mesh::{PatchMesh, PatchSide};
+use rtx_cfd::solvers::incompressible::{
+    ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
+    EmbeddedPisoSolver, FlowField, NormalDiffusion, OversetField, OversetParameters,
+    OversetPisoSolver, PatchConvection, PatchField, PoissonSolverKind,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+use std::f64::consts::PI;
+use std::sync::{Arc, RwLock};
+
+const RHO: f64 = 1000.0;
+const NU: f64 = 1e-3;
+const CX: f64 = 0.2;
+const CY: f64 = 0.2;
+const R: f64 = 0.05;
+const T: f64 = 0.01;
+const X_TIP: f64 = 0.6;
+const H_BOX: f64 = 0.41;
+const L_BOX: f64 = 1.2;
+const FILLET: f64 = 0.5 * 0.41 / 41.0;
+
+fn env_f(k: &str, d: f64) -> f64 {
+    std::env::var(k)
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(d)
+}
+
+/// The rigid patch rotated by `theta` about the cylinder's centre.
+fn rotated(base: &PatchMesh, theta: f64) -> CfdResult<PatchMesh> {
+    let (ns, nn) = (base.ns(), base.nn());
+    let (c, s) = (theta.cos(), theta.sin());
+    let mut xs = Vec::with_capacity((ns + 1) * (nn + 1));
+    let mut ys = Vec::with_capacity((ns + 1) * (nn + 1));
+    for k in 0..=nn {
+        for i in 0..=ns {
+            let p = base.node_xy(base.node(k, i));
+            let (dx, dy) = (p[0] - CX, p[1] - CY);
+            xs.push(CX + c * dx - s * dy);
+            ys.push(CY + s * dx + c * dy);
+        }
+    }
+    PatchMesh::from_nodes(ns, nn, xs, ys, Some([0.0, 0.0]))
+}
+
+#[tokio::test]
+async fn pitching_flag_in_still_fluid_extracts_energy_at_every_h() -> CfdResult<()> {
+    let ny = env_f("RTX_OVERSET_NY", 41.0) as usize;
+    let periods = env_f("RTX_OVERSET_EP_PERIODS", 4.0) as usize;
+    let tip_amp = env_f("RTX_OVERSET_EP_TIP", 0.02);
+    let f = env_f("RTX_OVERSET_EP_F", 2.0);
+    let h = H_BOX / ny as f64;
+    let nx = (L_BOX / h).round() as usize;
+    let omega = 2.0 * PI * f;
+    let theta0 = tip_amp / (X_TIP - CX);
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(RHO * NU)
+        .with_reference_velocity(tip_amp * omega)
+        .with_reference_length(2.0 * R);
+    let mut background = EmbeddedPisoSolver::new(
+        config.clone(),
+        EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            poisson_solver: PoissonSolverKind::Multigrid,
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+            ..EmbeddedParameters::default()
+        },
+    )?;
+    background.set_boundary_velocity(|_, _, _| (0.0, 0.0));
+    let (base, _) = cylinder_flag_patch([CX, CY], R, T, X_TIP, h, FILLET, 6.0 * h, 12, 4.0, 500)?;
+    let mut patch = CurvilinearPisoSolver::new(
+        config,
+        CurvilinearParameters {
+            convection: PatchConvection::TvdVanAlbada,
+            normal_diffusion: NormalDiffusion::LineImplicit,
+            ..CurvilinearParameters::default()
+        },
+        base.clone(),
+    )?;
+    // The wall's exact velocity: θ̇ × (r − c); the patch's clock is the step's end.
+    let theta_dot = Arc::new(RwLock::new(0.0_f64));
+    let td = theta_dot.clone();
+    patch.set_side_velocity(PatchSide::Inner, move |x, y, _| {
+        let w = *td.read().unwrap();
+        (-w * (y - CY), w * (x - CX))
+    });
+    let mut patch_field = PatchField::new(patch.mesh());
+    patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));
+    let params = OversetParameters {
+        stall_rounds: 0,
+        max_rounds: 3,
+        ..OversetParameters::default()
+    };
+    let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
+    let mut field = OversetField {
+        background: FlowField::new(nx, ny, h, h)?,
+        patch: patch_field,
+    };
+    solver.initialize(&mut field)?;
+
+    let mut hs = f64::INFINITY;
+    for c in 0..solver.patch().mesh().cell_count() {
+        for (fc, _) in solver.patch().mesh().cell_faces(c) {
+            if solver.patch().mesh().is_sface(fc) {
+                let d = solver.patch().mesh().faces()[fc].d;
+                hs = hs.min((d[0] * d[0] + d[1] * d[1]).sqrt());
+            }
+        }
+    }
+    let u_peak = tip_amp * omega;
+    let dt_bg = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+    let dt_patch = 0.4 * (hs * hs / (4.0 * NU)).min(hs / u_peak);
+    let period = 1.0 / f;
+    let steps_per_period = (period / dt_bg.min(dt_patch)).ceil() as usize;
+    let dt = period / steps_per_period as f64;
+    println!(
+        "  ENERGY PIN ny = {ny} (h {h:.4}, nx {nx}, patch ns {} nn {}, inner edge {hs:.4}): tip ± {tip_amp} m at {f} Hz (θ0 {theta0:.4} rad), dt {dt:.3e} ({steps_per_period} per period), {periods} periods",
+        base.ns(), base.nn()
+    );
+    let start = std::time::Instant::now();
+    let mut per_period: Vec<(f64, f64, f64)> = Vec::new(); // (work, +part, −part)
+    let (mut w_acc, mut wp, mut wn) = (0.0, 0.0, 0.0);
+    let mut fit: Vec<(f64, f64)> = Vec::new(); // (t, moment about the centre) for the last two periods
+    for step in 0..periods * steps_per_period {
+        let t_new = (step + 1) as f64 * dt;
+        let theta = theta0 * (omega * t_new).sin();
+        *theta_dot.write().unwrap() = theta0 * omega * (omega * t_new).cos();
+        solver.set_patch_mesh(rotated(&base, theta)?)?;
+        solver.advance(&mut field, dt).await?;
+        let w = *theta_dot.read().unwrap();
+        let (mut power, mut moment) = (0.0, 0.0);
+        for (centre, _normal, len, traction) in
+            solver
+                .patch()
+                .wall_tractions(&field.patch, PatchSide::Inner, solver.time())
+        {
+            let (vx, vy) = (-w * (centre[1] - CY), w * (centre[0] - CX));
+            power += (traction[0] * vx + traction[1] * vy) * len;
+            moment += ((centre[0] - CX) * traction[1] - (centre[1] - CY) * traction[0]) * len;
+        }
+        w_acc += power * dt;
+        if power > 0.0 {
+            wp += power * dt
+        } else {
+            wn += power * dt
+        }
+        if step >= periods.saturating_sub(2) * steps_per_period {
+            fit.push((t_new, moment));
+        }
+        if (step + 1) % steps_per_period == 0 {
+            per_period.push((w_acc, wp, wn));
+            println!(
+                "    period {}: net work {w_acc:+.4e} J/m (+{wp:.3e} / {wn:+.3e}), {:.0} s",
+                (step + 1) / steps_per_period,
+                start.elapsed().as_secs_f64()
+            );
+            w_acc = 0.0;
+            wp = 0.0;
+            wn = 0.0;
+        }
+    }
+    // Moment fit over the last two periods: M = a sin ωt + b cos ωt + c.
+    // θ̈ = −θ0 ω² sin ωt → the added-inertia reaction is +I_a θ0 ω² sin ωt: I_a = a / (θ0 ω²).
+    let n = fit.len() as f64;
+    let (mut ss, mut sc, mut cc, mut s1, mut c1, mut ys, mut yc, mut y1) =
+        (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+    for &(t, m) in &fit {
+        let (s, c) = ((omega * t).sin(), (omega * t).cos());
+        ss += s * s;
+        sc += s * c;
+        cc += c * c;
+        s1 += s;
+        c1 += c;
+        ys += s * m;
+        yc += c * m;
+        y1 += m;
+    }
+    let mat = [[ss, sc, s1], [sc, cc, c1], [s1, c1, n]];
+    let rhs = [ys, yc, y1];
+    let det = |a: [[f64; 3]; 3]| {
+        a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
+            - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
+            + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
+    };
+    let d = det(mat);
+    let mut sol = [0.0; 3];
+    for k in 0..3 {
+        let mut mk = mat;
+        for i in 0..3 {
+            mk[i][k] = rhs[i];
+        }
+        sol[k] = det(mk) / d;
+    }
+    let inertia = sol[0] / (theta0 * omega * omega);
+    let damping = -sol[1] / (theta0 * omega);
+    let last = per_period.last().copied().unwrap_or((f64::NAN, 0.0, 0.0));
+    println!(
+        "  ENERGY PIN ny = {ny}: last-period net work {:+.4e} J/m (gross +{:.3e} / {:+.3e}); added inertia {inertia:.4e} kg·m²/m, rotational damping {damping:.4e} N·m·s/m; {} steps in {:.0} s",
+        last.0, last.1, last.2, periods * steps_per_period, start.elapsed().as_secs_f64()
+    );
+    assert!(last.0.is_finite(), "non-finite work");
+    Ok(())
+}