embedded3 PERF-3 P1-5 (c) part 2 + (d): the device tables persist across steps and are updated in place within the band — one DeviceCut with both phases' aperture / open / activity buffers (the predictor phase is a flag); the faces of the changed cells for apertures, distances, shifts and open flags, the imposition band's faces for the surface velocities, the changed cells for activity and owner, the changed and wall cells for the wall flux, the merged-cell CSR whole; RTX_E3_BAND_CHECK=1 compares every table against a full build (passed); slab CSV byte-identical, device moving/cg green; ny 124 rebuild block 1,910 → 1,578 ms per step (tables 442 → 132)

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

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

@@ -65,6 +65,14 @@ unsafe impl ValidAsZeroBits for E3CutPtrs {}
 pub(super) struct DeviceCut {
     /// The surface-velocity tables (kept to share between the phases).
     ub_host: [Vec<f64>; 3],
+    /// Which phase's apertures / open flags / activity `ptrs` hands out
+    /// (P1-5: both sets live in the struct and persist across steps).
+    phase: Phase,
+    /// The predictor phase's instantaneous apertures, open flags (fluid
+    /// kind) and activity (fluid cells).
+    a_pred: [CudaSlice<f64>; 3],
+    open_pred: [CudaSlice<i32>; 3],
+    active_pred: CudaSlice<i32>,
     a: [CudaSlice<f64>; 3],
     d: [CudaSlice<f64>; 3],
     ub: [CudaSlice<f64>; 3],
@@ -307,8 +315,43 @@ impl DeviceCut {
         let l_cells = lap.elapsed();
         let ub_host = ub;
         let ub_dev = [up_f(&ub_host[0]), up_f(&ub_host[1]), up_f(&ub_host[2])];
+        // The predictor phase's set (instantaneous apertures, fluid kinds,
+        // fluid cells) — built alongside so that both persist (P1-5).
+        let open_pred = |c: usize| -> Vec<i32> {
+            (0..counts[c])
+                .into_par_iter()
+                .map(|idx| {
+                    let kind = match c {
+                        0 => mask.u_kind(idx),
+                        1 => mask.v_kind(idx),
+                        _ => mask.w_kind(idx),
+                    };
+                    i32::from(kind == FaceKind::Fluid)
+                })
+                .collect()
+        };
+        let active_pred: Vec<i32> = (0..nc).into_par_iter().map(|i| i32::from(mask.is_fluid_cell(i))).collect();
+        let (a_pred, open_pred_dev, active_pred_dev) = if projection {
+            (
+                [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)],
+                [up_i(&open_pred(0)), up_i(&open_pred(1)), up_i(&open_pred(2))],
+                up_i(&active_pred),
+            )
+        } else {
+            // A predictor-phase build IS the predictor set; the projection
+            // set is filled with the same arrays (replaced on the next update).
+            (
+                [up_f(ap_u), up_f(ap_v), up_f(ap_w)],
+                [up_i(&open[0]), up_i(&open[1]), up_i(&open[2])],
+                up_i(&active),
+            )
+        };
         let built = Self {
             ub_host,
+            phase,
+            a_pred,
+            open_pred: open_pred_dev,
+            active_pred: active_pred_dev,
             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: ub_dev,
@@ -347,42 +390,230 @@ impl DeviceCut {
     /// activity flags and the apertures (end-of-step instead of the step's)
     /// are recomputed and uploaded. Identical to a fresh `build_with` in
     /// every table.
-    pub(super) fn predictor_from(prev: Self, solver: &Solver, g: Grid) -> Option<Self> {
+    pub(super) fn predictor_from(mut prev: Self, solver: &Solver, g: Grid) -> Option<Self> {
+        // P1-5: the predictor set was built or updated with the projection's;
+        // switching phases is a flag.
+        let _ = (solver, g);
+        prev.phase = Phase::Predictor;
+        Some(prev)
+    }
+
+    /// P1-5: update the persistent tables in place for the projection phase
+    /// of the next step of a MOVING body: the per-face tables for the faces
+    /// of the changed cells and the surface velocities for the imposition
+    /// band's faces, the per-cell tables for the changed cells (the wall
+    /// flux for the wall cells too: its compatibility correction is
+    /// global). `false` when a full build is needed (no changed set, the
+    /// advancing-wall closure, a size change).
+    pub(super) fn update(&mut self, solver: &Solver, g: Grid, t: f64) -> bool {
+        use crate::solvers::incompressible::embedded3::poisson::device_cg::{scatter_f64, scatter_i32, scatter_u32};
         use rayon::prelude::*;
-        let mask = solver.mask()?;
-        let cut = mask.cut()?;
-        let rt = runtime();
-        let up_f = |v: &[f64]| -> CudaSlice<f64> {
-            rt.stream
-                .memcpy_stod(if v.is_empty() { &[0.0f64][..] } else { v })
-                .expect("upload")
-        };
-        let up_i = |v: &[i32]| -> CudaSlice<i32> { rt.stream.memcpy_stod(v).expect("upload") };
+        let Some(mask) = solver.mask() else { return false };
+        let Some(cut) = mask.cut() else { return false };
+        let Some(body) = solver.body() else { return false };
+        let Some(changed) = mask.changed_cells() else { return false };
+        if mask.wall_advancing || mask.wall_exchange_foot {
+            return false;
+        }
         let (nx, ny, nz) = (g.nx, g.ny, g.nz);
+        let nxy = nx * ny;
+        let nc = g.cells();
         let counts = [(nx + 1) * ny * nz, nx * (ny + 1) * nz, nx * ny * (nz + 1)];
-        let open = |c: usize| -> Vec<i32> {
-            (0..counts[c])
-                .into_par_iter()
-                .map(|idx| {
+        if self.active.len() != nc || self.a[0].len() != counts[0] {
+            return false;
+        }
+        let lap = Instant::now();
+        let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok();
+        // The faces of the changed cells, per component, ascending.
+        let mut marks: [Vec<bool>; 3] = [vec![false; counts[0]], vec![false; counts[1]], vec![false; counts[2]]];
+        for &idx in changed {
+            let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
+            marks[0][g.uface(k, j, i)] = true;
+            marks[0][g.uface(k, j, i + 1)] = true;
+            marks[1][g.vface(k, j, i)] = true;
+            marks[1][g.vface(k, j + 1, i)] = true;
+            marks[2][g.wface(k, j, i)] = true;
+            marks[2][g.wface(k + 1, j, i)] = true;
+        }
+        let touched: [Vec<u32>; 3] = [0, 1, 2].map(|c| {
+            (0..counts[c]).into_par_iter().filter(|&f| marks[c][f]).map(|f| f as u32).collect()
+        });
+        let cells_u32: Vec<u32> = changed.iter().map(|&i| i as u32).collect();
+        // The imposition band's faces (their surface velocity moves with t).
+        let band = mask.impose_band().unwrap_or(f64::INFINITY);
+        let dists: [&[f64]; 3] = [&cut.d_u, &cut.d_v, &cut.d_w];
+        let band_faces: [Vec<u32>; 3] = [0, 1, 2].map(|c| {
+            (0..counts[c]).into_par_iter().filter(|&f| dists[c][f].abs() <= band).map(|f| f as u32).collect()
+        });
+        let h = [g.dx, g.dy, g.dz];
+        let shifts = mask.face_shift_tables();
+        let step = mask.step_apertures();
+        let step_open = mask.step_open_flags();
+        for c in 0..3 {
+            let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1));
+            let pos = |idx: usize| -> [f64; 3] {
+                let (k, j, i) = (idx / (nj * ni), (idx / ni) % nj, idx % ni);
+                [
+                    (i as f64 + if c == 0 { 0.0 } else { 0.5 }) * h[0],
+                    (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],
+                ]
+            };
+            let a_inst: &[f64] = match c { 0 => &cut.a_u, 1 => &cut.a_v, _ => &cut.a_w };
+            let a_step: &[f64] = match (step, c) {
+                (Some(a), 0) => &a.0,
+                (Some(a), 1) => &a.1,
+                (Some(a), _) => &a.2,
+                (None, _) => a_inst,
+            };
+            let tf = &touched[c];
+            let pick = |v: &[f64]| tf.iter().map(|&f| v[f as usize]).collect::<Vec<f64>>();
+            scatter_f64(tf, &pick(a_step), &mut self.a[c]);
+            scatter_f64(tf, &pick(a_inst), &mut self.a_pred[c]);
+            scatter_f64(tf, &pick(dists[c]), &mut self.d[c]);
+            let open_step: Vec<i32> = tf.iter().map(|&f| i32::from(step_open.map_or(a_inst[f as usize] > 0.0, |o| o[c][f as usize]))).collect();
+            scatter_i32(tf, &open_step, &mut self.open[c]);
+            let open_inst: Vec<i32> = tf
+                .iter()
+                .map(|&f| {
                     let kind = match c {
-                        0 => mask.u_kind(idx),
-                        1 => mask.v_kind(idx),
-                        _ => mask.w_kind(idx),
+                        0 => mask.u_kind(f as usize),
+                        1 => mask.v_kind(f as usize),
+                        _ => mask.w_kind(f as usize),
                     };
                     i32::from(kind == FaceKind::Fluid)
                 })
-                .collect()
-        };
-        let active: Vec<i32> = (0..g.cells())
-            .into_par_iter()
-            .map(|i| i32::from(mask.is_fluid_cell(i)))
                 .collect();
-        Some(Self {
-            a: [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)],
-            open: [up_i(&open(0)), up_i(&open(1)), up_i(&open(2))],
-            active: up_i(&active),
-            ..prev
+            scatter_i32(tf, &open_inst, &mut self.open_pred[c]);
+            if let Some(sh) = shifts {
+                let idx3: Vec<u32> = tf.iter().flat_map(|&f| [3 * f, 3 * f + 1, 3 * f + 2]).collect();
+                let val3: Vec<f64> = tf.iter().flat_map(|&f| [sh[c][3 * f as usize], sh[c][3 * f as usize + 1], sh[c][3 * f as usize + 2]]).collect();
+                scatter_f64(&idx3, &val3, &mut self.shift[c]);
+            }
+            // The surface velocity at the foot for the band's faces.
+            let bf = &band_faces[c];
+            let ubv: Vec<f64> = bf
+                .par_iter()
+                .map(|&f| {
+                    let x = pos(f as usize);
+                    let xf = match (shifts, mask.wall_foot_centroid) {
+                        (Some(sh), true) => [x[0] + sh[c][3 * f as usize], x[1] + sh[c][3 * f as usize + 1], x[2] + sh[c][3 * f as usize + 2]],
+                        _ => x,
+                    };
+                    mask.surface_velocity_at(body, xf, c, t)
                 })
+                .collect();
+            for (&f, &v) in bf.iter().zip(&ubv) {
+                self.ub_host[c][f as usize] = v;
+            }
+            scatter_f64(bf, &ubv, &mut self.ub[c]);
+        }
+        let l_faces = lap.elapsed();
+        // Per-cell tables for the changed cells.
+        let active_step: Vec<i32> = changed.iter().map(|&i| i32::from(mask.cell_active(i))).collect();
+        let active_inst: Vec<i32> = changed.iter().map(|&i| i32::from(mask.is_fluid_cell(i))).collect();
+        let owner: Vec<u32> = changed.iter().map(|&i| mask.master(i).unwrap_or(i) as u32).collect();
+        scatter_i32(&cells_u32, &active_step, &mut self.active);
+        scatter_i32(&cells_u32, &active_inst, &mut self.active_pred);
+        scatter_u32(&cells_u32, &owner, &mut self.owner);
+        // The wall flux: the changed cells and every wall cell (the correction).
+        let wall_flux: &[f64] = solver.wall_fluxes();
+        if wall_flux.len() != nc {
+            return false;
+        }
+        let mut wmark = vec![false; nc];
+        for &i in changed {
+            wmark[i] = true;
+        }
+        let wcells: Vec<u32> = (0..nc)
+            .into_par_iter()
+            .filter(|&i| wmark[i] || cut.wall[i] != [0.0; 3])
+            .map(|i| i as u32)
+            .collect();
+        let wvals: Vec<f64> = wcells.iter().map(|&i| wall_flux[i as usize]).collect();
+        scatter_f64(&wcells, &wvals, &mut self.wall_flux);
+        // The merged-cell CSR whole (its lengths change).
+        let mut fold_ptr = vec![0u32; nc + 1];
+        let mut merged = 0;
+        for i in 0..nc {
+            if let Some(m) = mask.master(i) {
+                fold_ptr[m + 1] += 1;
+                merged += 1;
+            }
+        }
+        for i in 0..nc {
+            fold_ptr[i + 1] += fold_ptr[i];
+        }
+        let mut fold_idx = vec![0u32; merged];
+        let mut cursor: Vec<u32> = fold_ptr[..nc].to_vec();
+        for i in 0..nc {
+            if let Some(m) = mask.master(i) {
+                fold_idx[cursor[m] as usize] = i as u32;
+                cursor[m] += 1;
+            }
+        }
+        let rt = runtime();
+        let up_u = |v: &[u32]| -> CudaSlice<u32> {
+            rt.stream
+                .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v })
+                .expect("upload")
+        };
+        self.fold_ptr = up_u(&fold_ptr);
+        self.fold_idx = up_u(&fold_idx);
+        self.merged = merged;
+        self.phase = Phase::Projection;
+        rt.stream.synchronize().expect("sync");
+        if profile {
+            let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3;
+            eprintln!(
+                "    table laps (update): faces {:.0} ms ({} / {} / {} touched, {} / {} / {} band), cells {:.0} ms ({} changed)",
+                ms(l_faces),
+                touched[0].len(),
+                touched[1].len(),
+                touched[2].len(),
+                band_faces[0].len(),
+                band_faces[1].len(),
+                band_faces[2].len(),
+                ms(lap.elapsed() - l_faces),
+                changed.len()
+            );
+        }
+        true
+    }
+
+    /// `RTX_E3_BAND_CHECK=1`: every device table against a fresh full build.
+    pub(super) fn check_against(&self, full: &Self) {
+        let rt = runtime();
+        let same_f = |name: &str, a: &CudaSlice<f64>, b: &CudaSlice<f64>| {
+            let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl"));
+            assert!(x.len() == y.len() && x.iter().zip(&y).all(|(p, q)| p.to_bits() == q.to_bits()), "band tables: {name} differs");
+        };
+        let same_i = |name: &str, a: &CudaSlice<i32>, b: &CudaSlice<i32>| {
+            let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl"));
+            assert!(x == y, "band tables: {name} differs");
+        };
+        let same_u = |name: &str, a: &CudaSlice<u32>, b: &CudaSlice<u32>| {
+            let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl"));
+            assert!(x == y, "band tables: {name} differs");
+        };
+        for c in 0..3 {
+            same_f("a", &self.a[c], &full.a[c]);
+            same_f("a_pred", &self.a_pred[c], &full.a_pred[c]);
+            same_f("d", &self.d[c], &full.d[c]);
+            same_i("open", &self.open[c], &full.open[c]);
+            same_i("open_pred", &self.open_pred[c], &full.open_pred[c]);
+            same_f("shift", &self.shift[c], &full.shift[c]);
+            // ub: the band's faces only (beyond it the full build writes 0, ours keeps stale — never read).
+            let (x, y) = (rt.stream.memcpy_dtov(&self.ub[c]).expect("dl"), rt.stream.memcpy_dtov(&full.ub[c]).expect("dl"));
+            let bad = x.iter().zip(&y).filter(|(p, q)| **q != 0.0 && p.to_bits() != q.to_bits()).count();
+            assert!(bad == 0, "band tables: ub[{c}] differs on {bad} band faces");
+        }
+        same_f("wall_flux", &self.wall_flux, &full.wall_flux);
+        same_i("active", &self.active, &full.active);
+        same_i("active_pred", &self.active_pred, &full.active_pred);
+        same_u("owner", &self.owner, &full.owner);
+        same_u("fold_ptr", &self.fold_ptr, &full.fold_ptr);
+        same_u("fold_idx", &self.fold_idx, &full.fold_idx);
     }
 
     fn ptrs(&self) -> E3CutPtrs {
@@ -391,11 +622,15 @@ impl DeviceCut {
         let pf = |x: &CudaSlice<f64>| x.device_ptr(s).0;
         let pi = |x: &CudaSlice<i32>| x.device_ptr(s).0;
         let pu = |x: &CudaSlice<u32>| x.device_ptr(s).0;
+        let (a, open, active) = match self.phase {
+            Phase::Projection => (&self.a, &self.open, &self.active),
+            Phase::Predictor => (&self.a_pred, &self.open_pred, &self.active_pred),
+        };
         E3CutPtrs {
             ptrs: [
-                pf(&self.a[0]),
-                pf(&self.a[1]),
-                pf(&self.a[2]),
+                pf(&a[0]),
+                pf(&a[1]),
+                pf(&a[2]),
                 pf(&self.d[0]),
                 pf(&self.d[1]),
                 pf(&self.d[2]),
@@ -403,10 +638,10 @@ impl DeviceCut {
                 pf(&self.ub[1]),
                 pf(&self.ub[2]),
                 pf(&self.wall_flux),
-                pi(&self.open[0]),
-                pi(&self.open[1]),
-                pi(&self.open[2]),
-                pi(&self.active),
+                pi(&open[0]),
+                pi(&open[1]),
+                pi(&open[2]),
+                pi(active),
                 pu(&self.owner),
                 pu(&self.fold_ptr),
                 pu(&self.fold_idx),
@@ -501,7 +736,25 @@ impl DeviceStep {
             self.upload_after_rebuild(&field);
             self.host_field = Some(field);
             let l_up = lap.elapsed();
-            self.cut = DeviceCut::build(&self.solver, g, Phase::Projection, t_new);
+            // P1-5: the persistent tables updated in place when the mask has a
+            // changed set; a full build otherwise (and under the identity
+            // reference `RTX_E3_TABLES_REBUILD=1`).
+            let rebuild_all = std::env::var("RTX_E3_TABLES_REBUILD").is_ok();
+            let check = std::env::var("RTX_E3_BAND_CHECK").is_ok();
+            self.cut = match self.cut.take() {
+                Some(mut prev) if !rebuild_all => {
+                    if prev.update(&self.solver, g, t_new) {
+                        if check {
+                            let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new).expect("full");
+                            prev.check_against(&full);
+                        }
+                        Some(prev)
+                    } else {
+                        DeviceCut::build(&self.solver, g, Phase::Projection, t_new)
+                    }
+                }
+                _ => DeviceCut::build(&self.solver, g, Phase::Projection, t_new),
+            };
             let l_tables = lap.elapsed();
             if profile {
                 eprintln!(

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

@@ -578,6 +578,12 @@ impl Mask {
             _ => None,
         }
     }
+    /// The space-time open flags per u / v / w face (a moving cut wall).
+    #[must_use]
+    pub fn step_open_flags(&self) -> Option<[&[bool]; 3]> {
+        self.step_open.as_ref().map(|o| [o.0.as_slice(), o.1.as_slice(), o.2.as_slice()])
+    }
+
     /// The cells whose operator rows may have changed since the previous
     /// step (P1-4), or `None` when every row must be rebuilt.
     #[must_use]