From c30e3dba346632155c30e2710a21cb8c944b2ca1 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 19:30:54 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20S2-2b=20(host=20lever):=20narrow-ba?= =?UTF-8?q?nd=20=CF=86=20re-evaluation=20for=20moving=20bodies=20(max=5Fsu?= =?UTF-8?q?rface=5Fspeed)=20=E2=80=94=20bit-identical,=20the=20flag=20body?= =?UTF-8?q?'s=20ny=2062=20rebuild=203.54=20=E2=86=92=200.30=20s;=20the=20f?= =?UTF-8?q?lag=20driver=20uses=20it?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../solvers/incompressible/embedded3/cut.rs | 39 +++++- .../incompressible/embedded3/cutwall.rs | 15 ++- .../incompressible/embedded3/step/mod.rs | 21 ++- .../rtx-cfd/tests/embedded3_flag_geometry.rs | 124 ++++++++++++++++++ .../rtx-cfd/tests/embedded3_flag_wake.rs | 2 + 5 files changed, 192 insertions(+), 9 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs index 0a65ae0..a910718 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs @@ -28,6 +28,11 @@ pub struct CutGeometry { pub d_u: Vec, pub d_v: Vec, pub d_w: Vec, + /// A lower bound on |φ| per corner (the narrow band of a moving body: + /// a corner is re-evaluated only when its bound, decayed by the body's + /// motion, comes within the band; far corners keep a stale value with + /// the right sign, which is all their cells use). + pub bound: Vec, } impl CutGeometry { @@ -36,16 +41,41 @@ impl CutGeometry { (k * (g.ny + 1) + j) * (g.nx + 1) + i } - /// Build the cut data of `body` at time `t`. + /// Build the cut data of `body` at time `t` (every corner evaluated). pub fn build(body: &Body, grid: Grid, t: f64) -> Self { + Self::build_from(body, grid, t, None) + } + + /// Build the cut data of `body` at `t`, re-evaluating only the corners + /// of `prev` whose |φ| bound, decayed by `motion` (the body's largest + /// displacement since `prev`), falls within `band` of the surface. + /// Identical to [`Self::build`] in every cut cell. + pub fn build_from( + body: &Body, + grid: Grid, + t: f64, + prev: Option<(&CutGeometry, f64, f64)>, + ) -> Self { let g = grid; let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); - let mut phi = vec![0.0; (nx + 1) * (ny + 1) * (nz + 1)]; + let n_nodes = (nx + 1) * (ny + 1) * (nz + 1); + let mut phi = vec![0.0; n_nodes]; + let mut bound = vec![0.0; n_nodes]; for k in 0..=nz { for j in 0..=ny { for i in 0..=nx { - phi[Self::node(g, k, j, i)] = - body.phi(i as f64 * dx, j as f64 * dy, k as f64 * dz, t); + let n = Self::node(g, k, j, i); + if let Some((p, band, motion)) = prev { + let b = p.bound[n] - motion; + if b > band { + phi[n] = p.phi[n]; + bound[n] = b; + continue; + } + } + let v = body.phi(i as f64 * dx, j as f64 * dy, k as f64 * dz, t); + phi[n] = v; + bound[n] = v.abs(); } } } @@ -181,6 +211,7 @@ impl CutGeometry { d_u, d_v, d_w, + bound, } } 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 56e1182..a066067 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -108,7 +108,20 @@ pub(super) struct CvGeometry { impl Mask { /// Classify the grid against `body` at `t` by its cut geometry. pub fn build_cut(body: &Body, g: Grid, t: f64, b: Boundaries) -> Result { - let cut = CutGeometry::build(body, g, t); + Self::build_cut_from(body, g, t, b, None) + } + + /// As [`Self::build_cut`], re-evaluating φ only within `band` of the + /// previous geometry moved by at most `motion` (see + /// [`CutGeometry::build_from`]). + pub fn build_cut_from( + body: &Body, + g: Grid, + t: f64, + b: Boundaries, + prev: Option<(&CutGeometry, f64, f64)>, + ) -> Result { + let cut = CutGeometry::build_from(body, g, t, prev); let (nx, ny, nz) = (g.nx, g.ny, g.nz); let periodic = b.z0 == Side::Periodic; let allowed = |side: Side| matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall); 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 bfbe96c..6acf7c4 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 @@ -74,6 +74,11 @@ pub struct Parameters { pub inner_stop_factor: f64, /// The wall treatment of an embedded body. pub wall_scheme: WallScheme, + /// A moving body's largest surface speed: with it the cut geometry's + /// φ is re-evaluated only within three cells of the surface each step + /// (identical results; the far corners keep their sign). `None` = + /// every corner every step. + pub max_surface_speed: Option, } impl Default for Parameters { @@ -87,6 +92,7 @@ impl Default for Parameters { convection_scheme: ConvectionScheme::Upwind, inner_stop_factor: 1e-2, wall_scheme: WallScheme::GhostBinary, + max_surface_speed: None, } } } @@ -182,10 +188,17 @@ impl Solver { self.mask = None; } - fn build_mask(&self, body: &Body, g: Grid, t: f64) -> Mask { + fn build_mask(&self, body: &Body, g: Grid, t: f64, dt: f64) -> Mask { match self.params.wall_scheme { WallScheme::GhostBinary => Mask::build(body, g, t, self.params.boundaries), - WallScheme::CutCell => Mask::build_cut(body, g, t, self.params.boundaries), + WallScheme::CutCell => { + let h = g.dx.min(g.dy).min(g.dz); + let prev = match (self.params.max_surface_speed, &self.mask) { + (Some(speed), Some(m)) => m.cut().map(|c| (c, 3.0 * h, speed * dt)), + _ => None, + }; + Mask::build_cut_from(body, g, t, self.params.boundaries, prev) + } } .expect("embedded mask") } @@ -453,7 +466,7 @@ impl Solver { let t = self.time; if let Some(body) = &self.body { if self.mask.is_none() { - self.mask = Some(self.build_mask(body, field.grid, t)); + self.mask = Some(self.build_mask(body, field.grid, t, 0.0)); } } self.apply_boundary_normals(field, t); @@ -486,7 +499,7 @@ impl Solver { return 0; }; let mut fresh_cells = 0; - let mut new_mask = self.build_mask(body, field.grid, t_new); + let mut new_mask = self.build_mask(body, field.grid, t_new, dt); if let Some(old_mask) = &self.mask { fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); new_mask.set_step_apertures(old_mask); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs index bd1e701..b67ca35 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs @@ -73,3 +73,127 @@ fn flag_body_builds_at_ny_62() { assert!((closure[0].powi(2) + closure[1].powi(2) + closure[2].powi(2)).sqrt() < 1e-9); } } + +/// Where the per-step host rebuild's time goes at ny 62 (S2-2b's target). +#[test] +#[ignore = "profile: the rebuild's parts at ny 62 (seconds)"] +fn rebuild_profile_at_ny_62() { + use rtx_cfd::solvers::incompressible::embedded3::CutGeometry; + 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; + let body = Body::from_sdf(move |x, y, z, _t| cyl(x, y).min(flag_3d(x, y, z, 0.3, h))); + let b = Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + }; + let t0 = std::time::Instant::now(); + let _cut = CutGeometry::build(&body, g, 0.0); + let t_cut = t0.elapsed().as_secs_f64(); + let t1 = std::time::Instant::now(); + let nodes = (nx + 1) * (ny + 1) * (ny + 1); + let mut acc = 0.0; + for n in 0..nodes { + let (k, j, i) = ( + n / ((nx + 1) * (ny + 1)), + (n / (nx + 1)) % (ny + 1), + n % (nx + 1), + ); + acc += body.phi(i as f64 * h, j as f64 * h, k as f64 * h, 0.0); + } + let t_phi = t1.elapsed().as_secs_f64(); + let t2 = std::time::Instant::now(); + let mask = Mask::build_cut(&body, g, 0.0, b).expect("mask"); + let t_mask = t2.elapsed().as_secs_f64(); + let t3 = std::time::Instant::now(); + let (_table, _) = mask.wall_flux_table(&body, 0.0); + let t_flux = t3.elapsed().as_secs_f64(); + println!( + " ny 62 rebuild parts: cut geometry {t_cut:.2} s (of which φ at {nodes} nodes {t_phi:.2} s), whole mask build {t_mask:.2} s, wall-flux table {t_flux:.2} s (Σφ {acc:.1})" + ); +} + +/// The narrow band is invisible to the solution: the moving circle +/// marched with and without `max_surface_speed` gives bit-identical +/// fields; the flag body's rebuild time at ny 62 with the band recorded. +#[test] +fn narrow_band_is_bit_identical_and_fast() { + use rtx_cfd::solvers::incompressible::ConvectionScheme; + use rtx_cfd::solvers::incompressible::embedded3::{ + CutGeometry, Field, Fluid, Parameters, Solver, WallScheme, + }; + let n = 76; + let h = 1.0 / n as f64; + let dt = 3.24e-4; + let run = |band: bool| { + 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: 1e-8, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + max_surface_speed: band.then_some(1.0), + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|_, _, _, _| (0.0, 0.0, 0.0)); + let yc = |t: f64| 0.5 + 0.08 * (t / 0.08).sin(); + let vc = |t: f64| (t / 0.08).cos(); + solver.set_moving_body( + 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)), + ); + let g = Grid::cubic(n, n, 4, h); + let mut f = Field::new(g); + solver.initialize(&mut f); + let start = std::time::Instant::now(); + for _ in 0..40 { + solver.advance(&mut f, dt); + } + (f, start.elapsed().as_secs_f64()) + }; + let (full, t_full) = run(false); + let (band, t_band) = run(true); + let same = full.u == band.u && full.v == band.v && full.w == band.w && full.p == band.p; + println!( + " moving circle 76²×4, 40 steps: full {t_full:.2} s, band {t_band:.2} s; fields bit-identical: {same}" + ); + assert!(same, "the narrow band changed the solution"); + // The flag body at ny 62: a full build, then a banded rebuild after a small motion. + let ny = 62; + let hh = H / ny as f64; + let nx = (2.5 / hh).round() as usize; + let g = Grid::cubic(nx, ny, ny, hh); + let cyl = |x: f64, y: f64| ((x - 0.2_f64).powi(2) + (y - 0.2_f64).powi(2)).sqrt() - 0.05; + let body = Body::from_sdf(move |x, y, z, t| cyl(x, y).min(flag_3d(x, y, z, t, hh))); + let t0 = std::time::Instant::now(); + let prev = CutGeometry::build(&body, g, 0.3); + let t_prev = t0.elapsed().as_secs_f64(); + let t1 = std::time::Instant::now(); + let next = CutGeometry::build_from(&body, g, 0.31, Some((&prev, 3.0 * hh, 1.0e-3))); + let t_next = t1.elapsed().as_secs_f64(); + let exact = CutGeometry::build(&body, g, 0.31); + let vol_same = next.vol == exact.vol + && next.a_u == exact.a_u + && next.a_v == exact.a_v + && next.a_w == exact.a_w; + println!( + " flag ny 62: full build {t_prev:.2} s, banded rebuild {t_next:.2} s; volumes and apertures identical to a full build: {vol_same}" + ); + assert!(vol_same); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index 970566d..060942d 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -139,6 +139,8 @@ fn flag_wake_on_the_device() { x1: Side::PressureOutlet, ..Boundaries::default() }, + // The narrow band: the flag's tip speed bounds the surface motion. + max_surface_speed: Some(2.0 * std::f64::consts::PI * FREQ * AMP * 1.05), ..Parameters::default() }, );