embedded3 PERF-3 P1-3 (f): the mask build's whole-grid scans as parallel per-entry maps (cell classification with the anchor as the smallest fluid index, the three face-kind tables, the merging scan's small flags, the centroid-shift table per face) — slab flag ny 62 CSV byte-identical, device moving/cg green; ny 124 rebuild block 2,680 → 2,569 ms per step (build_mask 409 → 309)

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

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

@@ -48,42 +48,38 @@ impl Mask {
             vec![0.0; 3 * sizes[1]],
             vec![0.0; 3 * sizes[2]],
         ];
+        // P1-3 (f): per face in parallel (each entry from its own geometry).
+        use rayon::prelude::*;
         for c in 0..3 {
-            let (ni, nj, nk) = (
+            let (ni, nj) = (g.nx + usize::from(c == 0), g.ny + usize::from(c == 1));
-                g.nx + usize::from(c == 0),
+            tables[c]
-                g.ny + usize::from(c == 1),
+                .par_chunks_mut(3)
-                g.nz + usize::from(c == 2),
+                .enumerate()
-            );
+                .for_each(|(f, out)| {
-            for k in 0..nk {
+                    let (k, j, i) = (f / (nj * ni), (f / ni) % nj, f % ni);
-                for j in 0..nj {
+                    let p = [i as i64, j as i64, k as i64];
-                    for i in 0..ni {
+                    debug_assert_eq!(lat.face(c, p), Some(f));
-                        let p = [i as i64, j as i64, k as i64];
+                    let Some(alpha) = self.aperture(c, p) else { return };
-                        let Some(f) = lat.face(c, p) else { continue };
+                    if alpha <= 0.0 || alpha >= 1.0 {
-                        let Some(alpha) = self.aperture(c, p) else {
+                        return;
-                            continue;
-                        };
-                        if alpha <= 0.0 || alpha >= 1.0 {
-                            continue;
-                        }
-                        // Interior faces only (a control volume needs both cells).
-                        let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c];
-                        if on_side && !(c == 2 && self.periodic_z) {
-                            continue;
-                        }
-                        let cv = self.cv_geometry(c, p);
-                        let mut n = cv.wall;
-                        n[c] = 0.0;
-                        let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
-                        if a == 0.0 {
-                            continue;
-                        }
-                        for d in 0..3 {
-                            // `wall` points into the body: the open part lies the other way.
-                            tables[c][3 * f + d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a;
-                        }
                     }
-                }
+                    // Interior faces only (a control volume needs both cells).
-            }
+                    let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c];
+                    if on_side && !(c == 2 && self.periodic_z) {
+                        return;
+                    }
+                    let cv = self.cv_geometry(c, p);
+                    let mut n = cv.wall;
+                    n[c] = 0.0;
+                    let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
+                    if a == 0.0 {
+                        return;
+                    }
+                    for d in 0..3 {
+                        // `wall` points into the body: the open part lies the other way.
+                        out[d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a;
+                    }
+                });
         }
         self.face_shifts = Some(tables);
     }

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

@@ -128,35 +128,31 @@ impl Mask {
         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);
-        let mut cell_fluid = vec![true; g.cells()];
+        // P1-3 (f): the whole-grid classifications as parallel per-entry
-        let mut fluid_cells = 0;
+        // maps (the same values; the anchor is the smallest fluid index, the
-        let mut anchor = None;
+        // serial loop's first).
-        for k in 0..nz {
+        use rayon::prelude::*;
-            for j in 0..ny {
+        let nxy = nx * ny;
-                for i in 0..nx {
+        let cell_fluid: Vec<bool> = cut.vol.par_iter().map(|&v| v > 0.0).collect();
-                    let idx = g.cell(k, j, i);
+        let fluid_cells = cell_fluid.par_iter().filter(|&&f| f).count();
-                    let fluid = cut.vol[idx] > 0.0;
+        let anchor = (0..g.cells()).into_par_iter().find_first(|&idx| cell_fluid[idx]);
-                    cell_fluid[idx] = fluid;
+        let touching = (0..g.cells()).into_par_iter().find_first(|&idx| {
-                    if fluid {
+            if cell_fluid[idx] {
-                        fluid_cells += 1;
+                return false;
-                        if anchor.is_none() {
-                            anchor = Some(idx);
-                        }
-                    } else {
-                        let touches = (i == 0 && !allowed(b.x0))
-                            || (i + 1 == nx && !allowed(b.x1))
-                            || (j == 0 && !allowed(b.y0))
-                            || (j + 1 == ny && !allowed(b.y1))
-                            || (k == 0 && !allowed(b.z0))
-                            || (k + 1 == nz && !allowed(b.z1));
-                        if touches {
-                            return Err(format!(
-                                "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})"
-                            ));
-                        }
-                    }
-                }
             }
+            let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
+            (i == 0 && !allowed(b.x0))
+                || (i + 1 == nx && !allowed(b.x1))
+                || (j == 0 && !allowed(b.y0))
+                || (j + 1 == ny && !allowed(b.y1))
+                || (k == 0 && !allowed(b.z0))
+                || (k + 1 == nz && !allowed(b.z1))
+        });
+        if let Some(idx) = touching {
+            let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
+            return Err(format!(
+                "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})"
+            ));
         }
         let Some(anchor) = anchor else {
             return Err("embedded body covers the whole domain".into());
@@ -168,32 +164,28 @@ impl Mask {
                 FaceKind::Solid
             }
         };
-        let mut u_kind = vec![FaceKind::Fluid; g.n_ufaces()];
+        // Domain-side faces keep `Fluid` (the serial loops skipped them).
-        let mut v_kind = vec![FaceKind::Fluid; g.n_vfaces()];
+        let u_kind: Vec<FaceKind> = (0..g.n_ufaces())
-        let mut w_kind = vec![FaceKind::Fluid; g.n_wfaces()];
+            .into_par_iter()
-        for k in 0..nz {
+            .map(|f| {
-            for j in 0..ny {
+                let i = f % (nx + 1);
-                for i in 1..nx {
+                if i == 0 || i == nx { FaceKind::Fluid } else { kind(cut.a_u[f]) }
-                    let f = g.uface(k, j, i);
+            })
-                    u_kind[f] = kind(cut.a_u[f]);
+            .collect();
-                }
+        let v_kind: Vec<FaceKind> = (0..g.n_vfaces())
-            }
+            .into_par_iter()
-            for j in 1..ny {
+            .map(|f| {
-                for i in 0..nx {
+                let j = (f / nx) % (ny + 1);
-                    let f = g.vface(k, j, i);
+                if j == 0 || j == ny { FaceKind::Fluid } else { kind(cut.a_v[f]) }
-                    v_kind[f] = kind(cut.a_v[f]);
+            })
-                }
+            .collect();
-            }
+        let w_kind: Vec<FaceKind> = (0..g.n_wfaces())
-        }
+            .into_par_iter()
-        let w_range = if periodic { 0..nz + 1 } else { 1..nz };
+            .map(|f| {
-        for k in w_range {
+                let k = f / nxy;
-            for j in 0..ny {
+                if !periodic && (k == 0 || k == nz) { FaceKind::Fluid } else { kind(cut.a_w[f]) }
-                for i in 0..nx {
+            })
-                    let f = g.wface(k, j, i);
+            .collect();
-                    w_kind[f] = kind(cut.a_w[f]);
-                }
-            }
-        }
         let mut mask = Self {
             grid: g,
             periodic_z: periodic,
@@ -254,9 +246,13 @@ impl Mask {
             .ok()
             .and_then(|v| v.parse().ok())
             .unwrap_or(MERGE_FRACTION);
-        let small: Vec<bool> = (0..n)
+        let small: Vec<bool> = {
-            .map(|idx| self.cell_active(idx) && frac(idx) < threshold)
+            use rayon::prelude::*;
-            .collect();
+            (0..n)
+                .into_par_iter()
+                .map(|idx| self.cell_active(idx) && frac(idx) < threshold)
+                .collect()
+        };
         let mut master = vec![usize::MAX; n];
         for idx in (0..n).filter(|&i| small[i]) {
             let (k, j, i) = g.kji(idx);