embedded3 S2-2b (host lever): narrow-band φ re-evaluation for moving bodies (max_surface_speed) — bit-identical, the flag body's ny 62 rebuild 3.54 → 0.30 s; the flag driver uses it

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs

@@ -28,6 +28,11 @@ pub struct CutGeometry {
     pub d_u: Vec<f64>,
     pub d_v: Vec<f64>,
     pub d_w: Vec<f64>,
+    /// 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<f64>,
 }
 
 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,
         }
     }
 

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<Self, String> {
-        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<Self, String> {
+        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);

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<f64>,
 }
 
 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);

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);
+}

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()
         },
     );