From c6a35a77e21ee6251ca865a6ce5feecfd86d350d Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sun, 20 Sep 2026 13:54:37 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20PERF-3=20P1-4=20step=20e:=20Compone?= =?UTF-8?q?nts::summary=5Fplanes=20=E2=80=94=20the=20component=20count=20a?= =?UTF-8?q?nd=20singular=20flags=20alone=20on=20the=20device=20path=20(no?= =?UTF-8?q?=20per-cell=20id=20/=20member=20arrays);=20slab=20CSV=20byte-id?= =?UTF-8?q?entical,=20device=20moving/cg=20green;=20ny=20124=20rebuild=20b?= =?UTF-8?q?lock=202,119=20=E2=86=92=201,995=20ms=20per=20step=20(component?= =?UTF-8?q?s=20262=20=E2=86=92=20156)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../embedded3/poisson/device_cg.rs | 10 +- .../embedded3/poisson/hierarchy.rs | 126 ++++++++++++++++++ 2 files changed, 130 insertions(+), 6 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs index 522fb2b..6a4f4fb 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs +++ b/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() diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/hierarchy.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/hierarchy.rs index 746a33d..3dfaa6d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/hierarchy.rs +++ b/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) { + 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, 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 = (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 = (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 = 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