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 f53deb6..537bf7f 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 @@ -104,7 +104,7 @@ impl DeviceCg { let lap = std::time::Instant::now(); let fine = Level::::new(problem.clone()); let l_level = lap.elapsed(); - let components = Components::find(problem, &fine.cells); + let components = Components::find_parallel(problem, &fine.cells); let l_components = lap.elapsed(); let singular_count = components.singular.iter().filter(|&&s| s).count(); assert!( @@ -193,13 +193,17 @@ impl DeviceCg { /// operator changes little per step; the hierarchy's rebuild is the /// cost). `z` is zeroed before every V-cycle scatter, so cells absent /// from the stale hierarchy get no correction rather than a stale one. - pub fn refresh(&mut self, problem: &Problem, params: &MultigridParameters) { + pub fn refresh(&mut self, problem: Problem, params: &MultigridParameters) { let rt = runtime(); let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); let lap = std::time::Instant::now(); - let fine = Level::::new(problem.clone()); + // PERF-3 P1-3: the operator moves into the level (no 8-array clone); + // the level's masked problem gives the same components (couplings + // toward inactive cells are zero either way) and the same links. + let fine = Level::::new(problem); + let problem = &fine.problem; let l_level = lap.elapsed(); - let components = Components::find(problem, &fine.cells); + let components = Components::find_parallel(problem, &fine.cells); let l_components = lap.elapsed(); let singular_count = components.singular.iter().filter(|&&s| s).count(); assert!( 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 0153475..62fc3f1 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 @@ -437,4 +437,106 @@ impl Components { 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 + /// seed is the first unlabelled cell in ascending order = the smallest + /// of its component); the member lists are ascending (the serial ones + /// are in visit order — only their sets are used by the device path). + pub(crate) fn find_parallel(problem: &Problem, cells: &[usize]) -> Self { + use rayon::prelude::*; + use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering}; + let (nx, ny, nz) = (problem.nx, problem.ny, problem.nz); + let n = nx * ny * nz; + let nxy = nx * ny; + let parent: Vec = (0..n).map(AtomicUsize::new).collect(); + let links = problem.link_lists(); + let is_active = |idx: usize| problem.active[idx]; + // The positive-direction edges of a cell (each edge once) plus its links. + let edges = |idx: usize, mut f: &mut dyn FnMut(usize)| { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + if i + 1 < nx && problem.ae[idx] > 0.0 && is_active(idx + 1) { + f(idx + 1); + } + if j + 1 < ny && problem.an[idx] > 0.0 && is_active(idx + nx) { + f(idx + nx); + } + if let Some(t) = problem.top(idx, k) { + if problem.at[idx] > 0.0 && is_active(t) { + f(t); + } + } + for &(other, c) in &links[idx] { + if c > 0.0 && is_active(other) { + f(other); + } + } + let _ = &mut f; + }; + let root_of = |idx: usize| { + let mut r = idx; + loop { + let p = parent[r].load(Ordering::Relaxed); + if p == r { + return r; + } + r = p; + } + }; + loop { + let changed = AtomicBool::new(false); + cells.par_iter().for_each(|&idx| { + edges(idx, &mut |nb: usize| { + let (ru, rv) = (root_of(idx), root_of(nb)); + if ru != rv { + let (lo, hi) = if ru < rv { (ru, rv) } else { (rv, ru) }; + // Hook the larger root under the smaller (atomic min). + let mut cur = parent[hi].load(Ordering::Relaxed); + while lo < cur { + match parent[hi].compare_exchange_weak(cur, lo, Ordering::Relaxed, Ordering::Relaxed) { + Ok(_) => { + changed.store(true, Ordering::Relaxed); + break; + } + Err(actual) => cur = actual, + } + } + } + }); + }); + // Shortcut every cell to its root. + cells.par_iter().for_each(|&idx| { + let r = root_of(idx); + parent[idx].store(r, Ordering::Relaxed); + }); + if !changed.load(Ordering::Relaxed) { + break; + } + } + let mut id = vec![usize::MAX; n]; + let mut pairs: Vec<(usize, usize)> = cells.par_iter().map(|&idx| (root_of(idx), idx)).collect(); + pairs.par_sort_unstable(); + let mut members: Vec> = Vec::new(); + let mut singular = Vec::new(); + let mut last_root = usize::MAX; + for &(root, idx) in &pairs { + if root != last_root { + members.push(Vec::new()); + singular.push(true); + last_root = root; + } + let c = members.len() - 1; + members[c].push(idx); + id[idx] = c; + if problem.extra_diag[idx] > 0.0 { + singular[c] = false; + } + } + Self { + id, + members, + singular, + } + } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs index b64c49f..cd85aeb 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs @@ -142,7 +142,8 @@ impl DeviceCut { // S2-7b: the wall velocity at the open part's centroid's foot (the // host predictor's `foot_of`), and the axis feet of the solid exchange. let centroid_foot = mask.wall_foot_centroid; - let shifts = mask.face_shift_tables().cloned(); + // Borrowed: a clone here was 830 MB per build at ny 124 (P1-3 found it). + let shifts = mask.face_shift_tables(); let axis_foot = mask.wall_exchange_foot; let mut foot: [Vec; 3] = [vec![0.0], vec![0.0], vec![0.0]]; for c in 0..3 { @@ -170,7 +171,7 @@ impl DeviceCut { *ub_out = match shared { Some(sh) => sh[c][idx], None if dists[c][idx].abs() <= band => { - let xf = match (&shifts, centroid_foot) { + let xf = match (shifts, centroid_foot) { (Some(sh), true) => [ x[0] + sh[c][3 * idx], x[1] + sh[c][3 * idx + 1], @@ -531,7 +532,7 @@ impl DeviceStep { smoother: self.solver.params.poisson_smoother, ..MultigridParameters::default() }; - cg.refresh(&problem, ¶ms); + cg.refresh(problem, ¶ms); } rt.stream .memcpy_dtod(&self.u, &mut self.u_star)