embedded3 PERF-3 P1-4 step e: Components::summary_planes — the component count and singular flags alone on the device path (no per-cell id / member arrays); slab CSV byte-identical, device moving/cg green; ny 124 rebuild block 2,119 → 1,995 ms per step (components 262 → 156)

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

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

@@ -254,14 +254,12 @@ impl DeviceCg {
         }
         let problem = &fine.problem;
         let l_level = lap.elapsed();
-        let components = Components::find_planes(problem, &fine.cells);
+        let (n_components, singular_flags) = Components::summary_planes(problem, &fine.cells);
         let l_components = lap.elapsed();
-        let singular_count = components.singular.iter().filter(|&&s| s).count();
+        let singular_count = singular_flags.iter().filter(|&&s| s).count();
         assert!(
-            singular_count == 0 || components.members.len() == 1,
-            "DeviceCg::refresh: {} components with {} singular",
-            components.members.len(),
-            singular_count
+            singular_count == 0 || n_components == 1,
+            "DeviceCg::refresh: {n_components} components with {singular_count} singular"
         );
         let to_u32 = |v: &[usize]| {
             v.iter()

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

@@ -718,6 +718,132 @@ impl Components {
         }
     }
 
+    /// P1-4 step e: the component count and the singular flags alone (what
+    /// the device path reads), by the plane search without the per-cell id
+    /// and member arrays. The flags are indexed by the same numbering as
+    /// `find_planes` (roots by their smallest cell).
+    pub(crate) fn summary_planes(problem: &Problem, cells: &[usize]) -> (usize, Vec<bool>) {
+        use rayon::prelude::*;
+        let (nx, ny, nz) = (problem.nx, problem.ny, problem.nz);
+        let nxy = nx * ny;
+        let links = problem.link_map();
+        let no_links: Vec<(usize, f64)> = Vec::new();
+        let links_of = |idx: usize| -> &Vec<(usize, f64)> { links.get(&idx).unwrap_or(&no_links) };
+        let active = |idx: usize| problem.active[idx];
+        // Per plane: local piece ids, the count, and per piece (its smallest
+        // cell, whether it has a Dirichlet cell).
+        let pieces: Vec<(Vec<usize>, Vec<(usize, bool)>)> = (0..nz)
+            .into_par_iter()
+            .map(|k| {
+                let base = k * nxy;
+                let mut local = vec![usize::MAX; nxy];
+                let mut info: Vec<(usize, bool)> = Vec::new();
+                let mut stack = Vec::new();
+                for seed in 0..nxy {
+                    let idx = base + seed;
+                    if local[seed] != usize::MAX || !active(idx) {
+                        continue;
+                    }
+                    let count = info.len();
+                    let mut dirichlet = false;
+                    local[seed] = count;
+                    stack.push(seed);
+                    while let Some(s) = stack.pop() {
+                        let idx = base + s;
+                        if problem.extra_diag[idx] > 0.0 {
+                            dirichlet = true;
+                        }
+                        let (j, i) = (s / nx, s % nx);
+                        let mut visit = |t: usize, coefficient: f64| {
+                            if coefficient > 0.0 && local[t] == usize::MAX && active(base + t) {
+                                local[t] = count;
+                                stack.push(t);
+                            }
+                        };
+                        if i + 1 < nx {
+                            visit(s + 1, problem.ae[idx]);
+                        }
+                        if i > 0 {
+                            visit(s - 1, problem.aw[idx]);
+                        }
+                        if j + 1 < ny {
+                            visit(s + nx, problem.an[idx]);
+                        }
+                        if j > 0 {
+                            visit(s - nx, problem.as_[idx]);
+                        }
+                        if !links.is_empty() {
+                            for &(other, c) in links_of(idx) {
+                                if other / nxy == k {
+                                    visit(other % nxy, c);
+                                }
+                            }
+                        }
+                    }
+                    // The seed is the piece's smallest cell (ascending scan).
+                    info.push((idx, dirichlet));
+                }
+                (local, info)
+            })
+            .collect();
+        let mut offset = vec![0usize; nz + 1];
+        for k in 0..nz {
+            offset[k + 1] = offset[k] + pieces[k].1.len();
+        }
+        let piece_of = |idx: usize| -> usize {
+            let k = idx / nxy;
+            offset[k] + pieces[k].0[idx % nxy]
+        };
+        let total = offset[nz];
+        let mut parent: Vec<usize> = (0..total).collect();
+        fn root(parent: &mut [usize], mut a: usize) -> usize {
+            while parent[a] != a {
+                parent[a] = parent[parent[a]];
+                a = parent[a];
+            }
+            a
+        }
+        let mut union = |parent: &mut [usize], a: usize, b: usize| {
+            let (ra, rb) = (root(parent, a), root(parent, b));
+            if ra != rb {
+                let (lo, hi) = if ra < rb { (ra, rb) } else { (rb, ra) };
+                parent[hi] = lo;
+            }
+        };
+        let z_pairs: Vec<(usize, usize)> = cells
+            .par_iter()
+            .filter_map(|&idx| {
+                let k = idx / nxy;
+                let t = problem.top(idx, k)?;
+                (problem.at[idx] > 0.0 && active(t)).then(|| (piece_of(idx), piece_of(t)))
+            })
+            .filter(|(a, b)| a != b)
+            .collect();
+        for (a, b) in z_pairs {
+            union(&mut parent, a, b);
+        }
+        for &(a, b, c) in &problem.links {
+            if c > 0.0 && active(a) && active(b) && a / nxy != b / nxy {
+                union(&mut parent, piece_of(a), piece_of(b));
+            }
+        }
+        // Per root: the smallest cell and the Dirichlet flag; roots ordered
+        // by their smallest cell.
+        let mut min_cell = vec![usize::MAX; total];
+        let mut has_dirichlet = vec![false; total];
+        for pc in 0..total {
+            let r = root(&mut parent, pc);
+            let (k, local) = (0..nz).find_map(|k| (pc >= offset[k] && pc < offset[k + 1]).then_some((k, pc - offset[k]))).expect("piece plane");
+            let (cell, d) = pieces[k].1[local];
+            min_cell[r] = min_cell[r].min(cell);
+            has_dirichlet[r] |= d;
+        }
+        let mut roots: Vec<usize> = (0..total).filter(|&r| min_cell[r] != usize::MAX && parent[r] == r).collect();
+        roots.sort_unstable_by_key(|&r| min_cell[r]);
+        let singular: Vec<bool> = roots.iter().map(|&r| !has_dirichlet[r]).collect();
+        (roots.len(), singular)
+    }
+
     /// The same components by a parallel hook-and-shortcut union-find
     /// (PERF-3 P1-3): every component's root is its smallest cell index, so
     /// the numbering (by root, ascending) equals the serial search's (whose