From fe76413c669a112f1364a8a4b0b7a93e949f4af7 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sun, 20 Sep 2026 13:46:33 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20PERF-3=20P1-4=20(steps=201=E2=80=93?= =?UTF-8?q?3=20+=20d):=20the=20band-persistent=20operator=20=E2=80=94=20Cu?= =?UTF-8?q?tGeometry=20marks=20re-evaluated=20corners;=20Mask::changed=5Fc?= =?UTF-8?q?ells=20(touched=20by=20either=20build,=20dilated=20by=20one);?= =?UTF-8?q?=20Level::patch=20re-derives=20the=20changed=20rows=20with=20th?= =?UTF-8?q?e=20constructor's=20formulas=20on=20the=20kept=20level;=20the?= =?UTF-8?q?=20fine=20export=20patched=20likewise=20with=20the=20parent=20m?= =?UTF-8?q?ap=20as=20a=20parallel=20per-cell=20map;=20components=20plane?= =?UTF-8?q?=20by=20plane=20with=20a=20small=20union=20across=20planes;=20P?= =?UTF-8?q?roblem::link=5Fmap=20(sparse)=20replaces=20the=20per-cell=20lin?= =?UTF-8?q?k=20lists=20on=20the=20per-step=20path;=20RTX=5FE3=5FBAND=5FCHE?= =?UTF-8?q?CK=3D1=20compares=20against=20full=20rebuilds=20(passed=20on=20?= =?UTF-8?q?the=20slab)=20=E2=80=94=20slab=20CSV=20byte-identical=20at=20ev?= =?UTF-8?q?ery=20step,=20device=20moving/cg=20green;=20ny=20124=20rebuild?= =?UTF-8?q?=20block=202,569=20=E2=86=92=202,119=20ms=20per=20step?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../solvers/incompressible/embedded3/cut.rs | 34 ++- .../incompressible/embedded3/cutwall.rs | 31 ++ .../embedded3/poisson/device_cg.rs | 87 +++++- .../embedded3/poisson/export.rs | 30 +- .../embedded3/poisson/hierarchy.rs | 280 ++++++++++++++++++ .../embedded3/poisson/problem.rs | 12 + .../embedded3/step/device/cut.rs | 3 +- .../solvers/incompressible/embedded3/wall.rs | 12 + 8 files changed, 477 insertions(+), 12 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs index 414c119..ba192a2 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs @@ -33,6 +33,10 @@ pub struct CutGeometry { /// motion, comes within the band; far corners keep a stale value with /// the right sign, which is all their cells use). pub bound: Vec, + /// P1-4: the corners re-evaluated by this build (all of them without a + /// narrow band) — a cell whose corners are all untouched has unchanged + /// apertures, volume and activity. + pub touched: Vec, } impl CutGeometry { @@ -64,10 +68,12 @@ impl CutGeometry { use rayon::prelude::*; let mut phi = vec![0.0; n_nodes]; let mut bound = vec![0.0; n_nodes]; + let mut touched = vec![true; n_nodes]; phi.par_iter_mut() .zip(bound.par_iter_mut()) + .zip(touched.par_iter_mut()) .enumerate() - .for_each(|(n, (phi_n, bound_n))| { + .for_each(|(n, ((phi_n, bound_n), touched_n))| { let (k, j, i) = ( n / ((ny + 1) * (nx + 1)), (n / (nx + 1)) % (ny + 1), @@ -78,6 +84,7 @@ impl CutGeometry { if b > band { *phi_n = p.phi[n]; *bound_n = b; + *touched_n = false; return; } } @@ -196,9 +203,34 @@ impl CutGeometry { d_v, d_w, bound, + touched, } } + /// The cells with a touched corner (P1-4). + #[must_use] + pub fn touched_cells(&self) -> Vec { + use rayon::prelude::*; + let g = self.grid; + let (nx, ny) = (g.nx, g.ny); + (0..g.cells()) + .into_par_iter() + .map(|idx| { + let (k, j, i) = g.kji(idx); + let mut t = false; + for dk in 0..2 { + for dj in 0..2 { + for di in 0..2 { + t |= self.touched[Self::node(g, k + dk, j + dj, i + di)]; + } + } + } + let _ = (nx, ny); + t + }) + .collect() + } + /// φ at the corner `(k, j, i)` of the corner lattice. #[inline] #[must_use] diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs index f916c9c..b167cfa 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -211,6 +211,7 @@ impl Mask { wall_advancing: false, exchange_convection_off: false, wall_exchange_axis: false, + changed_cells: None, cv_sides_exact: false, wall_order2_centroid: false, wall_exchange_foot: false, @@ -318,7 +319,37 @@ impl Mask { .collect(); self.step_open = Some((open(&au), open(&av), open(&aw), active)); self.step_apertures = Some((au, av, aw)); + // P1-4: the changed set — cells touched by either build (this step's + // apertures average the two geometries), dilated by one. + let t_new = cut.touched_cells(); + let t_old = old_cut.touched_cells(); self.compute_merging(Some(old)); + { + use rayon::prelude::*; + let g = self.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let nxy = nx * ny; + let periodic = self.periodic_z; + let changed: Vec = (0..g.cells()) + .into_par_iter() + .filter(|&idx| { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + let t = |q: usize| t_new[q] || t_old[q]; + if t(idx) { + return true; + } + (i + 1 < nx && t(idx + 1)) + || (i > 0 && t(idx - 1)) + || (j + 1 < ny && t(idx + nx)) + || (j > 0 && t(idx - nx)) + || (k + 1 < nz && t(idx + nxy)) + || (k > 0 && t(idx - nxy)) + || (periodic && nz > 1 && k + 1 == nz && t(idx - (nz - 1) * nxy)) + || (periodic && nz > 1 && k == 0 && t(idx + (nz - 1) * nxy)) + }) + .collect(); + self.changed_cells = Some(changed); + } } pub(super) fn lattice(&self) -> Lattice { 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 537bf7f..522fb2b 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 @@ -89,6 +89,10 @@ pub struct DeviceCg { /// The singular component's members (all cells when singular; empty /// otherwise) and whether the system is singular. singular: bool, + /// The fine level kept between refreshes (P1-4: patched on the changed set). + fine: Option>, + /// Its export, kept likewise. + fine_export: Option, active_host: Vec, max_iterations: usize, scalar_host: Vec, @@ -104,7 +108,7 @@ impl DeviceCg { let lap = std::time::Instant::now(); let fine = Level::::new(problem.clone()); let l_level = lap.elapsed(); - let components = Components::find_parallel(problem, &fine.cells); + let components = Components::find_planes(problem, &fine.cells); let l_components = lap.elapsed(); let singular_count = components.singular.iter().filter(|&&s| s).count(); assert!( @@ -174,6 +178,8 @@ impl DeviceCg { scalar: rt.stream.alloc_zeros::(1).expect("alloc"), vcycle, singular: singular_count > 0, + fine: None, + fine_export: None, active_host: fine.active.clone(), max_iterations: params.max_iterations, scalar_host: vec![0.0], @@ -194,16 +200,61 @@ impl DeviceCg { /// 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) { + self.refresh_with(problem, params, None); + } + + /// As [`Self::refresh`]; with `changed` (P1-4) the previous fine level is + /// patched on those rows instead of rebuilt. `RTX_E3_BAND_CHECK=1` + /// compares the patched level against a full rebuild (the gate). + pub fn refresh_with(&mut self, problem: Problem, params: &MultigridParameters, changed: Option<&[usize]>) { let rt = runtime(); let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); let lap = std::time::Instant::now(); // 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 check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); + let full = if check { Some(Level::::new(problem.clone())) } else { None }; + let mut patched = false; + let fine = match (self.fine.take(), changed) { + (Some(mut prev), Some(rows)) => { + prev.patch(problem, rows); + patched = true; + prev + } + _ => Level::::new(problem), + }; + if let Some(full) = full { + let same = |a: &[f64], b: &[f64]| a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()); + assert!(fine.active == full.active, "band patch: active differs"); + assert!(fine.cells == full.cells && fine.red == full.red && fine.black == full.black, "band patch: lists differ"); + assert!(fine.top == full.top && fine.bot == full.bot, "band patch: z links differ"); + for (name, a, b) in [ + ("ae", &fine.ae, &full.ae), + ("aw", &fine.aw, &full.aw), + ("an", &fine.an, &full.an), + ("as", &fine.as_, &full.as_), + ("at", &fine.at, &full.at), + ("ab", &fine.ab, &full.ab), + ("ap", &fine.ap, &full.ap), + ] { + assert!(same(a, b), "band patch: {name} differs"); + } + for (name, a, b) in [ + ("p.ae", &fine.problem.ae, &full.problem.ae), + ("p.aw", &fine.problem.aw, &full.problem.aw), + ("p.an", &fine.problem.an, &full.problem.an), + ("p.as", &fine.problem.as_, &full.problem.as_), + ("p.at", &fine.problem.at, &full.problem.at), + ("p.ab", &fine.problem.ab, &full.problem.ab), + ] { + assert!(same(a, b), "band patch: {name} differs"); + } + assert!(fine.problem.active == full.problem.active, "band patch: problem.active differs"); + } let problem = &fine.problem; let l_level = lap.elapsed(); - let components = Components::find_parallel(problem, &fine.cells); + let components = Components::find_planes(problem, &fine.cells); let l_components = lap.elapsed(); let singular_count = components.singular.iter().filter(|&&s| s).count(); assert!( @@ -223,15 +274,17 @@ impl DeviceCg { .expect("upload") }; let up_f = |v: &[f64]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; - let lists = problem.link_lists(); + let map = problem.link_map(); let mut link_ptr = Vec::with_capacity(self.n + 1); let mut link_idx = Vec::new(); let mut link_coef = Vec::new(); link_ptr.push(0u32); - for list in &lists { - for &(other, c) in list { - link_idx.push(other as u32); - link_coef.push(c); + for idx in 0..self.n { + if let Some(list) = map.get(&idx) { + for &(other, c) in list { + link_idx.push(other as u32); + link_coef.push(c); + } } link_ptr.push(link_idx.len() as u32); } @@ -265,9 +318,25 @@ impl DeviceCg { let l_key = lap.elapsed(); // The V-cycle's finest level follows the operator (its coarser // levels stay): the fine level alone, no hierarchy build. - let fine_level = super::export::export_fine_from(&fine); + // P1-4 step 3: the export patched on the changed rows when the level was. + let fine_level = match (self.fine_export.take(), changed, patched) { + (Some(mut prev), Some(rows), true) => { + super::export::patch_fine_export(&mut prev, &fine, rows); + prev + } + _ => super::export::export_fine_from(&fine), + }; + if check { + let full = super::export::export_fine_from(&fine); + let same = |a: &[f32], b: &[f32]| a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()); + assert!(fine_level.cells == full.cells && fine_level.red == full.red && fine_level.black == full.black, "band export: lists differ"); + assert!(fine_level.top == full.top && fine_level.bot == full.bot && fine_level.coarse_of == full.coarse_of, "band export: maps differ"); + assert!(same(&fine_level.ae, &full.ae) && same(&fine_level.aw, &full.aw) && same(&fine_level.an, &full.an) && same(&fine_level.as_, &full.as_) && same(&fine_level.at, &full.at) && same(&fine_level.ab, &full.ab) && same(&fine_level.ap, &full.ap), "band export: coefficients differ"); + } let l_export = lap.elapsed(); self.vcycle.refresh_fine(&fine_level); + self.fine_export = Some(fine_level); + self.fine = Some(fine); if profile { let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; eprintln!( diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs index fb14cf0..cd22114 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs @@ -131,7 +131,7 @@ pub fn export_fine(problem: &Problem) -> LevelExport { /// second level build. pub(crate) fn export_fine_from(lv: &Level) -> LevelExport { use rayon::prelude::*; - let (_, coarse_of) = lv.coarsen(); + let coarse_of = lv.coarse_of_only(); // Per-entry maps (P1-3 (e)): the same values in the same order. let to_u32 = |v: &[usize]| { v.par_iter() @@ -171,3 +171,31 @@ pub fn vcycle_f32_reference( let mut hier = Hierarchy::::build(problem, params); hier.apply_preconditioner(r, z); } + +/// P1-4 step 3: the previous fine export with the rows of `changed` +/// re-cast from `lv` and the lists and parent map rebuilt (they follow the +/// activity, ascending) — equal to [`export_fine_from`] when every other +/// row of `lv` is unchanged. +pub(crate) fn patch_fine_export(prev: &mut LevelExport, lv: &Level, changed: &[usize]) { + use rayon::prelude::*; + let to_u32 = |v: &[usize]| { + v.par_iter() + .map(|&i| if i == usize::MAX { u32::MAX } else { i as u32 }) + .collect::>() + }; + for &idx in changed { + prev.ae[idx] = lv.ae[idx] as f32; + prev.aw[idx] = lv.aw[idx] as f32; + prev.an[idx] = lv.an[idx] as f32; + prev.as_[idx] = lv.as_[idx] as f32; + prev.at[idx] = lv.at[idx] as f32; + prev.ab[idx] = lv.ab[idx] as f32; + prev.ap[idx] = lv.ap[idx] as f32; + } + prev.cells = to_u32(&lv.cells); + prev.red = to_u32(&lv.red); + prev.black = to_u32(&lv.black); + prev.top = to_u32(&lv.top); + prev.bot = to_u32(&lv.bot); + prev.coarse_of = to_u32(&lv.coarse_of_only()); +} 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 62fc3f1..746a33d 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 @@ -139,6 +139,116 @@ impl Level { } } + /// P1-4: re-derive the rows of `changed` (ascending cell indices) from + /// `new` with [`Self::new`]'s formulas, every other row kept — exact + /// when every row outside `changed` is a function of unchanged inputs. + /// `new` replaces the stored problem (its arrays patched to the masked + /// values as the constructor leaves them). + pub(crate) fn patch(&mut self, mut new: Problem, changed: &[usize]) { + let _ = &mut new; + let (nx, ny, nz) = (new.nx, new.ny, new.nz); + let nxy = nx * ny; + let n = nxy * nz; + // The raw diagonal of a changed row, as `Problem::diagonals` sums it. + let mut link_sum = vec![0.0; n]; + for &(a, b, c) in &new.links { + link_sum[a] += c; + link_sum[b] += c; + } + let ap_of = |idx: usize| { + let stencil = new.ae[idx] + + new.aw[idx] + + new.an[idx] + + new.as_[idx] + + new.at[idx] + + new.ab[idx] + + new.extra_diag[idx]; + stencil + link_sum[idx] + }; + // Pass 1: the activity of the changed rows (their neighbours' new + // activity is either recomputed here or unchanged). + let new_active: Vec<(usize, bool, f64)> = changed + .par_iter() + .map(|&idx| { + let ap = ap_of(idx); + (idx, new.active[idx] && ap > 0.0, ap) + }) + .collect(); + for &(idx, a, ap) in &new_active { + self.active[idx] = a; + self.ap[idx] = T::from_f64(ap); + } + // Pass 2: the masked couplings and the z links of the changed rows. + let active = &self.active; + let link = |idx: usize, up: bool| -> Option { + let k = idx / nxy; + let o = if up { new.top(idx, k) } else { new.bottom(idx, k) }; + o.filter(|&t| active[t]) + }; + let rows: Vec<(usize, [f64; 6], usize, usize)> = changed + .par_iter() + .map(|&idx| { + let a = active[idx]; + let m = |v: f64, ok: bool| if a && ok { v } else { 0.0 }; + let (i, j) = (idx % nx, (idx % nxy) / nx); + let vals = [ + m(new.ae[idx], i + 1 < nx && active[idx + 1]), + m(new.aw[idx], i > 0 && active[idx - 1]), + m(new.an[idx], j + 1 < ny && active[idx + nx]), + m(new.as_[idx], j > 0 && active[idx - nx]), + m(new.at[idx], link(idx, true).is_some()), + m(new.ab[idx], link(idx, false).is_some()), + ]; + let top = if a { link(idx, true).unwrap_or(usize::MAX) } else { usize::MAX }; + let bot = if a { link(idx, false).unwrap_or(usize::MAX) } else { usize::MAX }; + (idx, vals, top, bot) + }) + .collect(); + // The stored problem keeps its masked arrays; the changed rows take + // the new masked values, their raw activity and Dirichlet part; the + // links are the new step's. + for &(idx, v, top, bot) in &rows { + self.problem.ae[idx] = v[0]; + self.problem.aw[idx] = v[1]; + self.problem.an[idx] = v[2]; + self.problem.as_[idx] = v[3]; + self.problem.at[idx] = v[4]; + self.problem.ab[idx] = v[5]; + self.problem.active[idx] = new.active[idx]; + self.problem.extra_diag[idx] = new.extra_diag[idx]; + self.ae[idx] = T::from_f64(v[0]); + self.aw[idx] = T::from_f64(v[1]); + self.an[idx] = T::from_f64(v[2]); + self.as_[idx] = T::from_f64(v[3]); + self.at[idx] = T::from_f64(v[4]); + self.ab[idx] = T::from_f64(v[5]); + self.top[idx] = top; + self.bot[idx] = bot; + } + self.problem.links = std::mem::take(&mut new.links); + self.problem.periodic_z = new.periodic_z; + drop(new); + // The lists from the activity flags (ascending, the constructor's). + let active = &self.active; + self.cells = (0..n).into_par_iter().filter(|&idx| active[idx]).collect(); + let parity = |idx: usize| (idx % nx + (idx % nxy) / nx + idx / nxy) % 2; + self.red = self.cells.par_iter().copied().filter(|&idx| parity(idx) == 0).collect(); + self.black = self.cells.par_iter().copied().filter(|&idx| parity(idx) == 1).collect(); + self.links = if self.problem.links.is_empty() { + Vec::new() + } else { + let map = self.problem.link_map(); + (0..n) + .into_par_iter() + .map(|idx| { + map.get(&idx).map_or_else(Vec::new, |l| { + l.iter().map(|&(o, c)| (o, T::from_f64(c))).collect() + }) + }) + .collect() + }; + } + /// `Σ a_nb x_nb`: the 2D order (e, w, n, s) with t, b appended. #[inline] fn neighbour_sum(&self, x: &[T], idx: usize) -> T { @@ -217,6 +327,29 @@ impl Level { } /// Galerkin coarsening by 2 in every direction (the 2D rule with k). + /// The parent map alone (P1-4 step 3): `coarsen`'s `coarse_of` as a + /// parallel per-cell map (the same values; the coarse operator it also + /// builds is what the export never used). + pub(crate) fn coarse_of_only(&self) -> Vec { + let p = &self.problem; + let (nx, ny, nz) = (p.nx, p.ny, p.nz); + let nxc = (nx / 2).max(1); + let nyc = (ny / 2).max(1); + let nzc = (nz / 2).max(1); + let active = &self.active; + (0..nx * ny * nz) + .into_par_iter() + .map(|idx| { + if !active[idx] { + return usize::MAX; + } + let (k, j, i) = (idx / (nx * ny), (idx % (nx * ny)) / nx, idx % nx); + let (ic, jc, kc) = ((i / 2).min(nxc - 1), (j / 2).min(nyc - 1), (k / 2).min(nzc - 1)); + (kc * nyc + jc) * nxc + ic + }) + .collect() + } + pub(crate) fn coarsen(&self) -> (Problem, Vec) { let p = &self.problem; let (nx, ny, nz) = (p.nx, p.ny, p.nz); @@ -438,6 +571,153 @@ impl Components { } } + /// The same components plane by plane (P1-4 step 1): a serial search + /// inside every z plane in parallel over the planes (x / y faces and + /// in-plane links), then one small union-find over the plane pieces + /// through the z faces and the cross-plane links. Exact: a component + /// is a union of plane pieces joined by those edges. Numbering by the + /// smallest cell index (the serial search's), members ascending. + pub(crate) fn find_planes(problem: &Problem, cells: &[usize]) -> Self { + use rayon::prelude::*; + let (nx, ny, nz) = (problem.nx, problem.ny, problem.nz); + let nxy = nx * ny; + let n = nxy * nz; + 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 (usize::MAX = not active), the count. + let pieces: Vec<(Vec, usize)> = (0..nz) + .into_par_iter() + .map(|k| { + let base = k * nxy; + let mut local = vec![usize::MAX; nxy]; + let mut count = 0; + let mut stack = Vec::new(); + for seed in 0..nxy { + let idx = base + seed; + if local[seed] != usize::MAX || !active(idx) { + continue; + } + local[seed] = count; + stack.push(seed); + while let Some(s) = stack.pop() { + let idx = base + s; + 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); + } + } + } + } + count += 1; + } + (local, count) + }) + .collect(); + // Global piece numbering and the union through z faces and links. + let mut offset = vec![0usize; nz + 1]; + for k in 0..nz { + offset[k + 1] = offset[k] + pieces[k].1; + } + 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; + } + }; + // The z edges (each once, from the cell to its top neighbour) — the + // candidate pairs gathered in parallel, joined serially. + 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)); + } + } + // Components numbered by their smallest cell: the first cell (in + // ascending order) of every root's set, i.e. the root's first piece + // in plane order — gather (root, cell) and sort. + let mut root_of_piece = vec![0usize; total]; + for pc in 0..total { + root_of_piece[pc] = root(&mut parent, pc); + } + // `cells` is ascending, so the pairs are ascending in idx already. + let pairs: Vec<(usize, usize)> = cells + .par_iter() + .map(|&idx| (root_of_piece[piece_of(idx)], idx)) + .collect(); + // Roots ordered by their smallest cell (pairs are ascending in idx). + let mut order: Vec = vec![usize::MAX; total]; + let mut next = 0; + for &(r, _) in &pairs { + if order[r] == usize::MAX { + order[r] = next; + next += 1; + } + } + let mut id = vec![usize::MAX; n]; + let mut members: Vec> = vec![Vec::new(); next]; + let mut singular = vec![true; next]; + for &(r, idx) in &pairs { + let c = order[r]; + id[idx] = c; + members[c].push(idx); + if problem.extra_diag[idx] > 0.0 { + singular[c] = false; + } + } + Self { + id, + members, + 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 diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs index a9e9276..4abef32 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs @@ -54,6 +54,18 @@ impl Problem { } /// The link coefficients per cell (`(other, coefficient)` lists). + /// The links per cell as a sparse map (P1-4: `link_lists` allocates a + /// vector per cell — 280 MB at ny 124 — for a few hundred links). + #[must_use] + pub fn link_map(&self) -> std::collections::HashMap> { + let mut out: std::collections::HashMap> = std::collections::HashMap::new(); + for &(a, b, c) in &self.links { + out.entry(a).or_default().push((b, c)); + out.entry(b).or_default().push((a, c)); + } + out + } + #[must_use] pub fn link_lists(&self) -> Vec> { let mut out = vec![Vec::new(); self.nx * self.ny * self.nz]; 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 4ba9126..56895ea 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 @@ -537,7 +537,8 @@ impl DeviceStep { smoother: self.solver.params.poisson_smoother, ..MultigridParameters::default() }; - cg.refresh(problem, ¶ms); + let changed: Option> = self.solver.mask().and_then(|m| m.changed_cells().map(|c| c.to_vec())); + cg.refresh_with(problem, ¶ms, changed.as_deref()); } rt.stream .memcpy_dtod(&self.u, &mut self.u_star) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs index f0e0aae..0ec01dd 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -110,6 +110,10 @@ pub struct Mask { pub(super) wall_exchange_axis: bool, /// S2-7: the momentum control volumes' side apertures from the /// interpolant on the sides' own corners (host prototype). + /// P1-4: the cells whose operator rows may differ from the previous + /// step's (touched by either mask's build, dilated by one), ascending; + /// `None` on a wall at rest or before the first step. + pub(super) changed_cells: Option>, pub(super) cv_sides_exact: bool, /// S2-7: the quadratic wall gradient's second point at the neighbour's /// own centroid distance (host prototype). @@ -546,6 +550,7 @@ impl Mask { wall_advancing: false, exchange_convection_off: false, wall_exchange_axis: false, + changed_cells: None, cv_sides_exact: false, wall_order2_centroid: false, wall_exchange_foot: false, @@ -573,6 +578,13 @@ impl Mask { _ => None, } } + /// The cells whose operator rows may have changed since the previous + /// step (P1-4), or `None` when every row must be rebuilt. + #[must_use] + pub fn changed_cells(&self) -> Option<&[usize]> { + self.changed_cells.as_deref() + } + #[must_use] pub fn merged_cells(&self) -> usize { self.merge_master