embedded3 PERF-3 P1-3 (c)(d)(e): the moving path's host mirror persists between steps and only the arrays the rebuild reads come down (u, v, w, p, the previous velocities) and only those it writes go up (u, v, w, p — the correction stays on the device; the first form re-uploaded a stale one and the byte gate caught it); the operator assembly in parallel by plane with the predicates lifted off the solver; the fine export's casts as parallel maps — slab flag ny 62 CSV byte-identical, device moving/cg green; ny 124 rebuild block 3,386 → 2,680 ms per step

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

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

@@ -130,13 +130,15 @@ pub fn export_fine(problem: &Problem) -> LevelExport {
 /// is the identical export without a second clone of the problem and a
 /// second level build.
 pub(crate) fn export_fine_from(lv: &Level<f64>) -> LevelExport {
+    use rayon::prelude::*;
     let (_, coarse_of) = lv.coarsen();
+    // Per-entry maps (P1-3 (e)): the same values in the same order.
     let to_u32 = |v: &[usize]| {
-        v.iter()
+        v.par_iter()
             .map(|&i| if i == usize::MAX { u32::MAX } else { i as u32 })
             .collect::<Vec<u32>>()
     };
-    let cast = |v: &[f64]| v.iter().map(|&x| x as f32).collect::<Vec<f32>>();
+    let cast = |v: &[f64]| v.par_iter().map(|&x| x as f32).collect::<Vec<f32>>();
     LevelExport {
         nx: lv.problem.nx,
         ny: lv.problem.ny,

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

@@ -151,6 +151,9 @@ pub struct DeviceStep {
     n_blocks: usize,
     cg: Option<DeviceCg>,
     cg_dt: f64,
+    /// The host mirror a moving body's rebuild reads and writes, kept
+    /// between steps (P1-3 (c)).
+    host_field: Option<Field>,
     /// PERF-3 P2: per corrector index, the basis of the projected initial guess.
     pub(super) guess: Vec<crate::solvers::incompressible::embedded3::poisson::GuessBasis>,
     timers: Option<StepTimers>,
@@ -231,6 +234,7 @@ impl DeviceStep {
             n_blocks,
             cg: None,
             cg_dt: 0.0,
+            host_field: None,
             guess: Vec::new(),
             timers,
             initialized: false,
@@ -266,6 +270,40 @@ impl DeviceStep {
         rt.stream.synchronize().expect("sync");
     }
 
+    /// The arrays a moving body's mask rebuild writes (P1-3 (c)): the
+    /// velocities (imposed and refilled) and the pressure (refilled) — the
+    /// correction stays on the device, untouched by the rebuild.
+    pub fn upload_after_rebuild(&mut self, field: &Field) {
+        let rt = runtime();
+        assert_eq!(field.grid, self.grid);
+        rt.stream.memcpy_htod(&field.u, &mut self.u).expect("u");
+        rt.stream.memcpy_htod(&field.v, &mut self.v).expect("v");
+        rt.stream.memcpy_htod(&field.w, &mut self.w).expect("w");
+        rt.stream.memcpy_htod(&field.p, &mut self.p).expect("p");
+        rt.stream.synchronize().expect("sync");
+    }
+
+    /// The arrays a moving body's mask rebuild reads (P1-3 (c)): the
+    /// corrected velocities, the pressure and the previous velocities.
+    pub fn download_for_rebuild(&self, field: &mut Field) {
+        let rt = runtime();
+        assert_eq!(field.grid, self.grid);
+        rt.stream.memcpy_dtoh(&self.u, &mut field.u).expect("u");
+        rt.stream.memcpy_dtoh(&self.v, &mut field.v).expect("v");
+        rt.stream.memcpy_dtoh(&self.w, &mut field.w).expect("w");
+        rt.stream.memcpy_dtoh(&self.p, &mut field.p).expect("p");
+        rt.stream
+            .memcpy_dtoh(&self.u_old, &mut field.u_old)
+            .expect("u_old");
+        rt.stream
+            .memcpy_dtoh(&self.v_old, &mut field.v_old)
+            .expect("v_old");
+        rt.stream
+            .memcpy_dtoh(&self.w_old, &mut field.w_old)
+            .expect("w_old");
+        rt.stream.synchronize().expect("sync");
+    }
+
     /// The device field into the host mirror.
     pub fn download(&self, field: &mut Field) {
         let rt = runtime();

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

@@ -489,12 +489,17 @@ impl DeviceStep {
             // `RTX_E3_MOVING_PROFILE=1`: the laps of the moving path (S2-2b).
             let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok();
             let lap = Instant::now();
-            let mut field = Field::new(g);
-            self.download(&mut field);
+            // P1-3 (c): the host mirror persists (no 1.1 GB allocation per
+            // step) and only the arrays the rebuild reads come down — u, v,
+            // w, p (fresh cells take their neighbours' pressure) and the
+            // previous velocities (the imposition's source).
+            let mut field = self.host_field.take().unwrap_or_else(|| Field::new(g));
+            self.download_for_rebuild(&mut field);
             let l_down = lap.elapsed();
             fresh_cells = self.solver.rebuild_moving_mask(&mut field, dt, t_new);
             let l_mask = lap.elapsed();
-            self.upload(&field);
+            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);
             let l_tables = lap.elapsed();

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

@@ -25,12 +25,48 @@ impl Solver {
         let at_outlet = dt * (dx * dy) / (0.5 * dz);
         let mut problem = Problem::new(nx, ny, nz);
         problem.periodic_z = periodic;
-        for k in 0..nz {
+        // P1-3 (d): every entry depends on its own cell alone — assembled in
+        // parallel by plane with the serial loop's arithmetic.
+        use rayon::prelude::*;
+        let nxy = nx * ny;
+        // The predicates over the mask alone (the solver holds `Cell`s and
+        // cannot cross threads; these are its `u_is_unknown` etc. verbatim).
+        let mask = self.mask.as_ref();
+        let cell_is_active = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.cell_active(g.cell(k, j, i)));
+        let u_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.u_open(g.uface(k, j, i)));
+        let v_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.v_open(g.vface(k, j, i)));
+        let w_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.w_open(g.wface(k, j, i)));
+        let au = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.au_step(g.uface(k, j, i)));
+        let av = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.av_step(g.vface(k, j, i)));
+        let aw = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.aw_step(g.wface(k, j, i)));
+        let gu = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.grad_weight(0, g.uface(k, j, i)));
+        let gv = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.grad_weight(1, g.vface(k, j, i)));
+        let (mut p_active, mut p_ae, mut p_aw, mut p_an, mut p_as, mut p_at, mut p_ab, mut p_extra) = (
+            std::mem::take(&mut problem.active),
+            std::mem::take(&mut problem.ae),
+            std::mem::take(&mut problem.aw),
+            std::mem::take(&mut problem.an),
+            std::mem::take(&mut problem.as_),
+            std::mem::take(&mut problem.at),
+            std::mem::take(&mut problem.ab),
+            std::mem::take(&mut problem.extra_diag),
+        );
+        p_active
+            .par_chunks_mut(nxy)
+            .zip(p_ae.par_chunks_mut(nxy))
+            .zip(p_aw.par_chunks_mut(nxy))
+            .zip(p_an.par_chunks_mut(nxy))
+            .zip(p_as.par_chunks_mut(nxy))
+            .zip(p_at.par_chunks_mut(nxy))
+            .zip(p_ab.par_chunks_mut(nxy))
+            .zip(p_extra.par_chunks_mut(nxy))
+            .enumerate()
+            .for_each(|(k, (((((((c_active, c_ae), c_aw), c_an), c_as), c_at), c_ab), c_extra))| {
             for j in 0..ny {
                 for i in 0..nx {
-                    let idx = g.cell(k, j, i);
-                    if !self.cell_is_active(k, j, i) {
-                        problem.active[idx] = false;
+                    let idx = j * nx + i;
+                    if !cell_is_active(k, j, i) {
+                        c_active[idx] = false;
                         continue;
                     }
                     let mut extra = 0.0;
@@ -38,48 +74,56 @@ impl Solver {
                         if b.x1 == outlet {
                             extra += ae_outlet;
                         }
-                    } else if self.u_is_unknown(k, j, i + 1) {
-                        problem.ae[idx] = ae_interior * self.au(k, j, i + 1) * self.gu(k, j, i + 1);
+                    } else if u_is_unknown(k, j, i + 1) {
+                        c_ae[idx] = ae_interior * au(k, j, i + 1) * gu(k, j, i + 1);
                     }
                     if i == 0 {
                         if b.x0 == outlet {
                             extra += ae_outlet;
                         }
-                    } else if self.u_is_unknown(k, j, i) {
-                        problem.aw[idx] = ae_interior * self.au(k, j, i) * self.gu(k, j, i);
+                    } else if u_is_unknown(k, j, i) {
+                        c_aw[idx] = ae_interior * au(k, j, i) * gu(k, j, i);
                     }
                     if j + 1 == ny {
                         if b.y1 == outlet {
                             extra += an_outlet;
                         }
-                    } else if self.v_is_unknown(k, j + 1, i) {
-                        problem.an[idx] = an_interior * self.av(k, j + 1, i) * self.gv(k, j + 1, i);
+                    } else if v_is_unknown(k, j + 1, i) {
+                        c_an[idx] = an_interior * av(k, j + 1, i) * gv(k, j + 1, i);
                     }
                     if j == 0 {
                         if b.y0 == outlet {
                             extra += an_outlet;
                         }
-                    } else if self.v_is_unknown(k, j, i) {
-                        problem.as_[idx] = an_interior * self.av(k, j, i) * self.gv(k, j, i);
+                    } else if v_is_unknown(k, j, i) {
+                        c_as[idx] = an_interior * av(k, j, i) * gv(k, j, i);
                     }
                     if k + 1 == nz && !periodic {
                         if b.z1 == outlet {
                             extra += at_outlet;
                         }
-                    } else if self.w_is_unknown((k + 1) % nz, j, i) {
-                        problem.at[idx] = at_interior * self.aw((k + 1) % nz, j, i);
+                    } else if w_is_unknown((k + 1) % nz, j, i) {
+                        c_at[idx] = at_interior * aw((k + 1) % nz, j, i);
                     }
                     if k == 0 && !periodic {
                         if b.z0 == outlet {
                             extra += at_outlet;
                         }
-                    } else if self.w_is_unknown(k, j, i) {
-                        problem.ab[idx] = at_interior * self.aw(k, j, i);
+                    } else if w_is_unknown(k, j, i) {
+                        c_ab[idx] = at_interior * aw(k, j, i);
                     }
-                    problem.extra_diag[idx] = extra;
+                    c_extra[idx] = extra;
                 }
             }
-        }
+            });
+        problem.active = p_active;
+        problem.ae = p_ae;
+        problem.aw = p_aw;
+        problem.an = p_an;
+        problem.as_ = p_as;
+        problem.at = p_at;
+        problem.ab = p_ab;
+        problem.extra_diag = p_extra;
         self.merge_small_cells(&mut problem);
         problem
     }