From aa096465a695cb0dbecdc1b4b4d882b7eff3d639 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 18:53:00 -0500 Subject: [PATCH] embedded3 S2-2a/S2-3: moving bodies on the device (phased tables, host rebuild per step, impose kernel; gate test), the per-span cut wall route, the flag-wake driver and its geometry pre-flight Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/src/kernels/cuda/e3_cut.cu | 20 ++ .../incompressible/embedded3/cutwall.rs | 37 ++- .../incompressible/embedded3/step/device.rs | 3 +- .../embedded3/step/device/cut.rs | 115 +++++++- .../incompressible/embedded3/step/mod.rs | 89 +++--- .../solvers/incompressible/embedded3/wall.rs | 6 + .../rtx-cfd/tests/embedded3_device_moving.rs | 118 ++++++++ .../rtx-cfd/tests/embedded3_flag_geometry.rs | 75 +++++ .../rtx-cfd/tests/embedded3_flag_wake.rs | 273 ++++++++++++++++++ 9 files changed, 687 insertions(+), 49 deletions(-) create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu index d367912..3796086 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -307,3 +307,23 @@ extern "C" __global__ void e3_cut_add_p(E3Params g, E3Ptrs f, E3Cut m) if (!m.active[t]) return; f.p[t] += f.pp[m.owner[t]]; } + +/* The prescribed interior faces of component c take the surface velocity + * (the host `impose` for a cut mask; `open` = the instantaneous kinds). */ +extern "C" __global__ void e3_cut_impose(E3Params g, E3Ptrs f, E3Cut m, int c) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + int ni = c == 0 ? g.nx + 1 : g.nx; + int nj = c == 1 ? g.ny + 1 : g.ny; + int nk = c == 2 ? g.nz + 1 : g.nz; + if (t >= ni * nj * nk) return; + int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj); + if (c == 0 && (i == 0 || i == g.nx)) return; + if (c == 1 && (j == 0 || j == g.ny)) return; + if (c == 2) { if (g.periodic_z) { if (k == g.nz) return; } else if (k == 0 || k == g.nz) return; } + const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w); + if (open[t]) return; + const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w); + double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w); + out[t] = ubt[t]; +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs index 33674e9..56e1182 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -472,6 +472,22 @@ impl Mask { Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]]) } + /// The cut-cell load route restricted to the cells (and faces) of the + /// planes `k0..k1`, divided by the slab's thickness: the load per unit + /// span on a body's mid-section. + pub fn cut_wall_force_per_span( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + (k0, k1): (usize, usize), + ) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_parts_in(body, f, mu, t, Some((k0, k1)))?; + let lz = (k1 - k0) as f64 * self.grid.dz; + Some([(p[0] + s[0]) / lz, (p[1] + s[1]) / lz, (p[2] + s[2]) / lz]) + } + /// The cut-cell load route split into its pressure and shear parts. pub fn cut_wall_force_parts( &self, @@ -479,14 +495,27 @@ impl Mask { f: &Field, mu: f64, t: f64, + ) -> Option<([f64; 3], [f64; 3])> { + self.cut_wall_force_parts_in(body, f, mu, t, None) + } + + fn cut_wall_force_parts_in( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + planes: Option<(usize, usize)>, ) -> Option<([f64; 3], [f64; 3])> { let cut = self.cut.as_ref()?; let g = self.grid; let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let (k0, k1) = planes.unwrap_or((0, nz)); let mut pressure = [0.0; 3]; let mut force = [0.0; 3]; for (idx, w) in cut.wall.iter().enumerate() { - if self.cell_fluid[idx] { + let k = g.kji(idx).0; + if self.cell_fluid[idx] && k >= k0 && k < k1 { for c in 0..3 { pressure[c] += f.p[idx] * w[c]; } @@ -497,9 +526,9 @@ impl Mask { let w_range = if self.periodic_z { 0..nz } else { 1..nz }; for c in 0..3 { let (ir, jr, kr) = match c { - 0 => (1..nx, 0..ny, 0..nz), - 1 => (0..nx, 1..ny, 0..nz), - _ => (0..nx, 0..ny, w_range.clone()), + 0 => (1..nx, 0..ny, k0..k1), + 1 => (0..nx, 1..ny, k0..k1), + _ => (0..nx, 0..ny, w_range.start.max(k0)..w_range.end.min(k1)), }; for k in kr { for j in jr.clone() { diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs index 1e3dfb2..14a813e 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs @@ -121,6 +121,7 @@ pub struct StepTimers { pub predictor_ns: u64, pub poisson_ns: u64, pub apply_ns: u64, + /// The moving body's host rebuild (mask, tables, transfers) per step. pub transfer_ns: u64, pub steps: u64, pub cg_iterations: u64, @@ -181,7 +182,7 @@ impl DeviceStep { let timers = std::env::var("RTX_PROFILE") .is_ok() .then(StepTimers::default); - let cut = cut::DeviceCut::build(&solver, grid); + let cut = cut::DeviceCut::build(&solver, grid, cut::Phase::Predictor, solver.time()); Self { solver, grid, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs index da82d3c..1b92ab3 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs @@ -6,9 +6,11 @@ use super::{DeviceStep, E3Params, E3Ptrs, StepResult}; use crate::solvers::incompressible::embedded3::Grid; +use crate::solvers::incompressible::embedded3::field::Field; use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime}; use crate::solvers::incompressible::embedded3::poisson::device_cg::DeviceCg; use crate::solvers::incompressible::embedded3::step::Solver; +use crate::solvers::incompressible::embedded3::wall::FaceKind; use crate::solvers::incompressible::poisson::MultigridParameters; use cudarc::driver::{ CudaFunction, CudaModule, CudaSlice, DevicePtr, DeviceRepr, PushKernelArg, ValidAsZeroBits, @@ -28,6 +30,7 @@ struct CutKernels { fold: CudaFunction, correct: CudaFunction, add_p: CudaFunction, + impose: CudaFunction, } static CUT_KERNELS_ONCE: OnceLock = OnceLock::new(); @@ -42,6 +45,7 @@ fn cut_kernels() -> &'static CutKernels { fold: f("e3_cut_fold"), correct: f("e3_cut_correct"), add_p: f("e3_cut_add_p"), + impose: f("e3_cut_impose"), _module: module, } }) @@ -71,12 +75,23 @@ pub(super) struct DeviceCut { pub(super) merged: usize, } +/// Which phase the tables serve: the predictor reads the instantaneous +/// apertures and kinds of the mask at its time; the projection reads the +/// step-averaged apertures and the space-time classification. +#[derive(Clone, Copy, PartialEq, Eq)] +pub(super) enum Phase { + Predictor, + Projection, +} + impl DeviceCut { - /// The tables of the solver's cut mask (`None` without one). - pub(super) fn build(solver: &Solver, g: Grid) -> Option { + /// The tables of the solver's cut mask (`None` without one) for + /// `phase`, the surface velocities at the mask's time `t`. + pub(super) fn build(solver: &Solver, g: Grid, phase: Phase, t: f64) -> Option { let mask = solver.mask()?; let cut = mask.cut()?; let body = solver.body()?; + let projection = phase == Phase::Projection; let rt = runtime(); let (nx, ny, nz) = (g.nx, g.ny, g.nz); let h = [g.dx, g.dy, g.dz]; @@ -112,11 +127,20 @@ impl DeviceCut { (j as f64 + if c == 1 { 0.0 } else { 0.5 }) * h[1], (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2], ]; - ubc[idx] = mask.surface_velocity_at(body, x, c, 0.0); - opc[idx] = i32::from(match c { - 0 => mask.u_open(idx), - 1 => mask.v_open(idx), - _ => mask.w_open(idx), + ubc[idx] = mask.surface_velocity_at(body, x, c, t); + opc[idx] = i32::from(if projection { + match c { + 0 => mask.u_open(idx), + 1 => mask.v_open(idx), + _ => mask.w_open(idx), + } + } else { + let kind = match c { + 0 => mask.u_kind(idx), + 1 => mask.v_kind(idx), + _ => mask.w_kind(idx), + }; + kind == FaceKind::Fluid }); } } @@ -124,9 +148,31 @@ impl DeviceCut { ub[c] = ubc; open[c] = opc; } - let (wall_flux, _) = mask.wall_flux_table(body, 0.0); + // The solver's current table (the GCL table on a moving body) when + // it has one, else the static porous table. + let wall_flux: Vec = if solver.wall_fluxes().len() == g.cells() { + solver.wall_fluxes().to_vec() + } else { + mask.wall_flux_table(body, t).0 + }; let nc = g.cells(); - let active: Vec = (0..nc).map(|i| i32::from(mask.cell_active(i))).collect(); + let active: Vec = (0..nc) + .map(|i| { + i32::from(if projection { + mask.cell_active(i) + } else { + mask.is_fluid_cell(i) + }) + }) + .collect(); + let step = if projection { + mask.step_apertures() + } else { + None + }; + let ap_u: &[f64] = step.map_or(&cut.a_u, |a| &a.0); + let ap_v: &[f64] = step.map_or(&cut.a_v, |a| &a.1); + let ap_w: &[f64] = step.map_or(&cut.a_w, |a| &a.2); let owner: Vec = (0..nc) .map(|i| mask.master(i).unwrap_or(i) as u32) .collect(); @@ -146,7 +192,7 @@ impl DeviceCut { fold_ptr.push(fold_idx.len() as u32); } Some(Self { - a: [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)], + a: [up_f(ap_u), up_f(ap_v), up_f(ap_w)], d: [up_f(&cut.d_u), up_f(&cut.d_v), up_f(&cut.d_w)], ub: [up_f(&ub[0]), up_f(&ub[1]), up_f(&ub[2])], wall_flux: up_f(&wall_flux), @@ -246,6 +292,30 @@ impl DeviceStep { .expect("w*"); rt.stream.synchronize().expect("sync"); let t_pred = t0.elapsed(); + // A moving body: the mask at the end-of-step geometry on the host + // (the predicted field down, the rebuilt one up), the projection + // tables and the operator rebuilt. + let t_rebuild = Instant::now(); + let mut fresh_cells = 0; + if self.solver.is_moving() { + let mut field = Field::new(g); + self.download(&mut field); + fresh_cells = self.solver.rebuild_moving_mask(&mut field, dt, t_new); + self.upload(&field); + self.cut = DeviceCut::build(&self.solver, g, Phase::Projection, t_new); + self.cg = None; + rt.stream + .memcpy_dtod(&self.u, &mut self.u_star) + .expect("u*"); + rt.stream + .memcpy_dtod(&self.v, &mut self.v_star) + .expect("v*"); + rt.stream + .memcpy_dtod(&self.w, &mut self.w_star) + .expect("w*"); + } + let cptrs = self.cut.as_ref().expect("cut").ptrs(); + let rebuild = t_rebuild.elapsed(); if self.cg.is_none() || self.cg_dt != dt { let problem = self.solver.poisson_operator(g, dt); let params = MultigridParameters { @@ -361,16 +431,39 @@ impl DeviceStep { .memcpy_dtod(&self.w, &mut self.w_star) .expect("w*"); } + // The prescribed faces take the surface velocity (the host's + // end-of-step impose), and the next predictor's tables. + if self.solver.is_moving() { + self.cut = DeviceCut::build(&self.solver, g, Phase::Predictor, t_new); + let cptrs = self.cut.as_ref().expect("cut").ptrs(); + for c in 0..3i32 { + unsafe { + rt.stream + .launch_builder(&k.impose) + .arg(&prm) + .arg(&ptrs) + .arg(&cptrs) + .arg(&c) + .launch(cfg(counts[c as usize])) + .expect("e3_cut_impose"); + } + } + if prm.periodic_z != 0 { + self.launch_sides(prm, &ptrs, 0); + } + rt.stream.synchronize().expect("sync"); + } self.solver.set_time(t_new); if let Some(tm) = self.timers.as_mut() { tm.predictor_ns += t_pred.as_nanos() as u64; tm.poisson_ns += t_poisson.as_nanos() as u64; tm.apply_ns += t_apply.as_nanos() as u64; + tm.transfer_ns += rebuild.as_nanos() as u64; tm.steps += 1; tm.cg_iterations += cg_iterations as u64; } StepResult { - fresh_cells: 0, + fresh_cells, converged: final_residual < self.solver.params.tolerance, corrector_steps_performed: total, final_residual, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index c40b88f..bfbe96c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -463,6 +463,56 @@ impl Solver { self.initialized = true; } + /// Whether the body moves (the mask is rebuilt every step). + #[must_use] + pub fn is_moving(&self) -> bool { + self.moving + } + + /// The current step's compatible wall-flux table (cut wall). + #[must_use] + pub fn wall_fluxes(&self) -> &[f64] { + &self.wall_fluxes + } + + /// The moving body's mask at the end-of-step geometry `t_new`: the + /// pressure of the cells that just became fluid refilled from their + /// neighbours (fluid in both masks), the new mask's prescribed and + /// ghost values imposed from the previous corrected field, the + /// step-averaged apertures and the GCL wall-flux table (cut wall). + /// Returns the fresh-cell count. `field` holds the predicted field. + pub fn rebuild_moving_mask(&mut self, field: &mut Field, dt: f64, t_new: f64) -> usize { + let Some(body) = &self.body else { + return 0; + }; + let mut fresh_cells = 0; + let mut new_mask = self.build_mask(body, field.grid, t_new); + if let Some(old_mask) = &self.mask { + fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); + new_mask.set_step_apertures(old_mask); + } + new_mask.impose_from( + body, + &field.u_old, + &field.v_old, + &field.w_old, + &mut field.u, + &mut field.v, + &mut field.w, + t_new, + ); + if new_mask.cut().is_some() { + let (table, correction) = match &self.mask { + Some(old_mask) => new_mask.gcl_flux_table(old_mask, dt), + None => new_mask.wall_flux_table(body, t_new), + }; + self.wall_fluxes = table; + self.last_ghost_correction = correction; + } + self.mask = Some(new_mask); + fresh_cells + } + /// One step of `dt`: predictor, correctors, clock. pub fn advance(&mut self, field: &mut Field, dt: f64) -> StepResult { assert!( @@ -477,39 +527,12 @@ impl Solver { field.update_old_values(); self.momentum_predictor(field, dt, t_old); self.apply_boundary_normals(field, t_new); - // A moving body: the mask at the end-of-step geometry, the pressure - // of the cells that just became fluid refilled from their - // neighbours (fluid in both masks), the new mask's prescribed and - // ghost values imposed from the previous corrected field. - let mut fresh_cells = 0; - if self.moving { - if let Some(body) = &self.body { - let mut new_mask = self.build_mask(body, field.grid, t_new); - if let Some(old_mask) = &self.mask { - fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); - new_mask.set_step_apertures(old_mask); - } - new_mask.impose_from( - body, - &field.u_old, - &field.v_old, - &field.w_old, - &mut field.u, - &mut field.v, - &mut field.w, - t_new, - ); - if new_mask.cut().is_some() { - let (table, correction) = match &self.mask { - Some(old_mask) => new_mask.gcl_flux_table(old_mask, dt), - None => new_mask.wall_flux_table(body, t_new), - }; - self.wall_fluxes = table; - self.last_ghost_correction = correction; - } - self.mask = Some(new_mask); - } - } + // A moving body: the mask at the end-of-step geometry. + let fresh_cells = if self.moving { + self.rebuild_moving_mask(field, dt, t_new) + } else { + 0 + }; field.copy_to_starred(); let mut cut_correction = None; if let (Some(body), Some(mask)) = (&self.body, &self.mask) { diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs index 3d5eff6..a3adfe5 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -505,6 +505,12 @@ impl Mask { }) } + /// The step-averaged apertures of a moving cut wall (`None` at rest). + #[must_use] + pub fn step_apertures(&self) -> Option<&(Vec, Vec, Vec)> { + self.step_apertures.as_ref() + } + /// The master of a virtually merged small cell. #[inline] #[must_use] diff --git a/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs b/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs new file mode 100644 index 0000000..a0807cb --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs @@ -0,0 +1,118 @@ +//! embedded3 S2-2a: a moving body on the device — the circle of the +//! falsifier (R 0.05 across a periodic 4-cell slab on h = 1/152, at 1 m/s +//! peak) on the cut wall, host vs device over 100 steps under tight +//! tolerances to `1e-9·scale`; the host rebuild's share of the step time +//! recorded. +//! +//! `RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_device_moving -- --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, +}; + +const N: usize = 152; +const AMP: f64 = 0.08; +const U: f64 = 1.0; + +fn circle(moving: bool) -> Body { + let yc = move |t: f64| { + if moving { + 0.5 + AMP * (U / AMP * t).sin() + } else { + 0.5 + } + }; + let vc = move |t: f64| if moving { U * (U / AMP * t).cos() } else { 0.0 }; + Body::from_sdf(move |x, y, _z, t| ((x - 0.5_f64).powi(2) + (y - yc(t)).powi(2)).sqrt() - 0.05) + .with_surface_velocity(move |_, _, _, t| (0.0, vc(t), 0.0)) +} + +fn make(tight: bool) -> (Solver, Grid) { + let (tolerance, inner_stop_factor) = if tight { (1e-12, 1e-6) } else { (1e-8, 1e-2) }; + let mut solver = Solver::new( + Fluid { + density: 1000.0, + viscosity: 1.0, + reference_velocity: 1.0, + reference_length: 0.1, + }, + Parameters { + corrector_steps: 2, + tolerance, + inner_stop_factor, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|_, _, _, _| (0.0, 0.0, 0.0)); + solver.set_moving_body(circle(true)); + (solver, Grid::cubic(N, N, 4, 1.0 / N as f64)) +} + +fn max_diff(a: &[f64], b: &[f64]) -> f64 { + a.iter() + .zip(b) + .fold(0.0_f64, |m, (&x, &y)| m.max((x - y).abs())) +} +fn scale(a: &[f64]) -> f64 { + a.iter().fold(0.0_f64, |m, &x| m.max(x.abs())) +} + +fn march(tight: bool, steps: usize, bound: Option) { + let dt = 3.24e-4; + let (mut host, g) = make(tight); + let mut fh = Field::new(g); + host.initialize(&mut fh); + let (mut ds, _) = make(tight); + let mut fd = Field::new(g); + ds.initialize(&mut fd); + unsafe { std::env::set_var("RTX_PROFILE", "1") }; + let mut device = DeviceStep::new(ds, g); + device.upload(&fd); + let start = std::time::Instant::now(); + let (mut differ, mut fresh_h, mut fresh_d) = (0, 0, 0); + for _ in 0..steps { + let rh = host.advance(&mut fh, dt); + let rd = device.advance(dt); + if rh.corrector_steps_performed != rd.corrector_steps_performed { + differ += 1; + } + fresh_h += rh.fresh_cells; + fresh_d += rd.fresh_cells; + } + let seconds = start.elapsed().as_secs_f64(); + device.download(&mut fd); + let du = max_diff(&fh.u, &fd.u) + .max(max_diff(&fh.v, &fd.v)) + .max(max_diff(&fh.w, &fd.w)); + let su = scale(&fh.u).max(scale(&fh.v)).max(scale(&fh.w)); + let dp = max_diff(&fh.p, &fd.p); + let sp = scale(&fh.p).max(1000.0); + let t = device.timers().expect("timers"); + println!( + " moving circle 152²×4 CutCell (tight {tight}): {steps} steps; host vs device max |Δu| {du:.3e} on {su:.3e}, max |Δp| {dp:.3e} on {sp:.3e}; fresh cells host {fresh_h} device {fresh_d}; corrector counts differ on {differ} steps; {:.1} ms per step (host + device), device step {:.1} ms of which rebuild {:.1} ms", + 1e3 * seconds / steps as f64, + 1e-6 * (t.predictor_ns + t.poisson_ns + t.apply_ns + t.transfer_ns) as f64 / steps as f64, + 1e-6 * t.transfer_ns as f64 / steps as f64 + ); + assert_eq!(fresh_h, fresh_d, "fresh-cell counts differ"); + if let Some(b) = bound { + assert!(du < b * su, "velocity differs: {du:.3e} on {su:.3e}"); + assert!(dp < b * sp, "pressure differs: {dp:.3e} on {sp:.3e}"); + } +} + +#[test] +fn moving_circle_host_equals_device() { + march(false, 100, None); + march(true, 100, Some(1e-9)); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs new file mode 100644 index 0000000..bd1e701 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs @@ -0,0 +1,75 @@ +//! S2-3 pre-flight (host): the flag + cylinder body at ny 62 — the mask +//! builds, its cut geometry closes, the counts and the build time (the +//! moving body's per-step host rebuild cost) are recorded. +use rtx_cfd::solvers::incompressible::embedded3::{Body, Boundaries, Grid, Mask, Side}; + +const H: f64 = 0.41; +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 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 flag_2d(x: f64, y: f64, phase: f64) -> f64 { + let n = 40; + let mut best = f64::INFINITY; + let point = |m: usize| { + let s = m as f64 / n as f64; + (FLAG_X0 + s * FLAG_LEN, 0.2 + AMP * mode(s) * phase) + }; + for m in 0..n { + let (ax, ay) = point(m); + let (bx, by) = point(m + 1); + 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(); + best = best.min(d); + } + best - FLAG_HALF +} + +fn flag_3d(x: f64, y: f64, z: f64, phase: f64, r: f64) -> f64 { + let d2 = flag_2d(x, y, phase); + let q1 = d2 + r; + let q2 = (z - 0.5 * H).abs() - 0.5 * FLAG_SPAN + r; + (q1.max(0.0).powi(2) + q2.max(0.0).powi(2)).sqrt() + q1.max(q2).min(0.0) - r +} + +#[test] +fn flag_body_builds_at_ny_62() { + let ny = 62; + let h = H / ny as f64; + let nx = (2.5 / h).round() as usize; + let g = Grid::cubic(nx, ny, ny, h); + let cyl = |x: f64, y: f64| ((x - 0.2_f64).powi(2) + (y - 0.2_f64).powi(2)).sqrt() - 0.05; + for phase in [0.0, 1.0] { + let body = Body::from_sdf(move |x, y, z, _t| cyl(x, y).min(flag_3d(x, y, z, phase, h))); + let b = Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + }; + let start = std::time::Instant::now(); + let mask = Mask::build_cut(&body, g, 0.0, b).expect("mask"); + let build = start.elapsed().as_secs_f64(); + let cut = mask.cut().unwrap(); + let (area, closure) = cut.wall_area_and_closure(); + let solid = g.cells() - mask.fluid_cells(); + println!( + " phase {phase}: {} cells, {} fluid, {solid} solid, {} merged; wall area {area:.4} m² (cylinder 0.129 + flag ~0.156), closure {:.2e}; build {build:.2} s", + g.cells(), + mask.fluid_cells(), + mask.merged_cells(), + (closure[0].powi(2) + closure[1].powi(2) + closure[2].powi(2)).sqrt() + ); + assert!(solid > 1000 && mask.fluid_cells() > g.cells() / 2); + assert!((closure[0].powi(2) + closure[1].powi(2) + closure[2].powi(2)).sqrt() < 1e-9); + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs new file mode 100644 index 0000000..970566d --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -0,0 +1,273 @@ +//! 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= RTX_E3_FLAG_CSV= \ +//! 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:?}"); + } +}