diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/mg_vcycle.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/mg_vcycle.cu index 837c9dd..f908ff2 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/mg_vcycle.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/mg_vcycle.cu @@ -110,3 +110,112 @@ extern "C" __global__ void mg_coarsest( } } } + +/* --------------------------------------------------------------------------- + * PERF-2 P3-iii (`MarchPlane`): the same maps over K LANES that each own + * their operator (their own mask, coefficients, index lists and work + * buffers). A LaneLevel is one lane's level: device pointers plus counts; + * the kernels take a [K] table of them per level and blockIdx.y = the lane. + * The per-cell arithmetic and its order are the K = 1 kernels' above, so a + * lane's correction does not depend on K or on its neighbours. The layout + * (14 pointers, then 6 ints) is mirrored by `LaneLevel` in device.rs. + * ------------------------------------------------------------------------- */ +struct LaneLevel { + const float* ae; const float* aw; const float* an; const float* as_; const float* ap; + float* b; float* x; float* r; + const unsigned int* cells; const unsigned int* red; const unsigned int* black; + const unsigned int* coarse_of; const unsigned int* children_ptr; const unsigned int* children_idx; + int n; int nx; int n_cells; int n_red; int n_black; int pad; +}; + +extern "C" __global__ void ml_rb_half(const LaneLevel* __restrict__ L, int colour) +{ + const LaneLevel& l = L[blockIdx.y]; + int t = blockIdx.x * blockDim.x + threadIdx.x; + const unsigned int* list = colour ? l.black : l.red; + int n_list = colour ? l.n_black : l.n_red; + if (t >= n_list) return; + int g = list[t]; + int nx = l.nx; + const float* x = l.x; + float s = 0.0f; + float e = l.ae[g]; if (e != 0.0f) s += e * x[g + 1]; + float w = l.aw[g]; if (w != 0.0f) s += w * x[g - 1]; + float nn = l.an[g]; if (nn != 0.0f) s += nn * x[g + nx]; + float ss = l.as_[g]; if (ss != 0.0f) s += ss * x[g - nx]; + l.x[g] = (l.b[g] + s) / l.ap[g]; +} + +extern "C" __global__ void ml_residual(const LaneLevel* __restrict__ L) +{ + const LaneLevel& l = L[blockIdx.y]; + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= l.n_cells) return; + int g = l.cells[t]; + int nx = l.nx; + const float* x = l.x; + float s = 0.0f; + float e = l.ae[g]; if (e != 0.0f) s += e * x[g + 1]; + float w = l.aw[g]; if (w != 0.0f) s += w * x[g - 1]; + float nn = l.an[g]; if (nn != 0.0f) s += nn * x[g + nx]; + float ss = l.as_[g]; if (ss != 0.0f) s += ss * x[g - nx]; + l.r[g] = l.b[g] - (l.ap[g] * x[g] - s); +} + +/* coarse.b[c] = sum of fine.r over the children of coarse cell c (fixed order). */ +extern "C" __global__ void ml_restrict(const LaneLevel* __restrict__ F, const LaneLevel* __restrict__ C) +{ + const LaneLevel& f = F[blockIdx.y]; + const LaneLevel& c = C[blockIdx.y]; + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= c.n_cells) return; + float s = 0.0f; + for (unsigned int p = f.children_ptr[t]; p < f.children_ptr[t + 1]; ++p) s += f.r[f.children_idx[p]]; + c.b[c.cells[t]] = s; +} + +/* fine.x += 2 coarse.x[coarse_of[idx]] */ +extern "C" __global__ void ml_prolong(const LaneLevel* __restrict__ F, const LaneLevel* __restrict__ C) +{ + const LaneLevel& f = F[blockIdx.y]; + const LaneLevel& c = C[blockIdx.y]; + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= f.n_cells) return; + int idx = f.cells[t]; + f.x[idx] += 2.0f * c.x[f.coarse_of[idx]]; +} + +extern "C" __global__ void ml_zero(const LaneLevel* __restrict__ L) +{ + const LaneLevel& l = L[blockIdx.y]; + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= l.n_cells) return; + l.x[l.cells[t]] = 0.0f; +} + +/* The coarsest level: one thread per lane, the CPU's red, black, black, red + * ordering from zero (see mg_coarsest). */ +extern "C" __global__ void ml_coarsest(const LaneLevel* __restrict__ L, int K, int sweeps) +{ + int k = blockIdx.x * blockDim.x + threadIdx.x; + if (k >= K) return; + const LaneLevel& l = L[k]; + int nx = l.nx; + float* x = l.x; + for (int t = 0; t < l.n_cells; ++t) x[l.cells[t]] = 0.0f; + for (int sw = 0; sw < sweeps; ++sw) { + for (int half = 0; half < 4; ++half) { + const unsigned int* list = (half == 0 || half == 3) ? l.red : l.black; + int n_list = (half == 0 || half == 3) ? l.n_red : l.n_black; + for (int t = 0; t < n_list; ++t) { + int g = list[t]; + float s = 0.0f; + float e = l.ae[g]; if (e != 0.0f) s += e * x[g + 1]; + float w = l.aw[g]; if (w != 0.0f) s += w * x[g - 1]; + float nn = l.an[g]; if (nn != 0.0f) s += nn * x[g + nx]; + float ss = l.as_[g]; if (ss != 0.0f) s += ss * x[g - nx]; + x[g] = (l.b[g] + s) / l.ap[g]; + } + } + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs index 5819242..02eb116 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs @@ -66,8 +66,9 @@ pub use piso::{PisoParameters, PisoResult, PisoSolver}; pub use piso_gpu::PisoGpuSolver; pub use poisson::{ LevelExport, MgPrecision, MgSmoother, MultigridParameters, PcgCache, PoissonProblem, - PoissonSolution, export_hierarchy, vcycle_f32_reference, vcycle_f32_work, - PoissonSolverKind, configure_threads, solve_multigrid_pcg, solve_multigrid_pcg_cached, + PoissonSolution, PoissonSolverKind, configure_threads, export_hierarchy, plane_counters, + set_plane_lane, solve_multigrid_pcg, solve_multigrid_pcg_cached, vcycle_f32_reference, + vcycle_f32_work, }; pub use polygon_sdf::PolygonSdf; pub use simple::{ConvectionScheme, SimpleParameters, SimpleResult, SimpleSolver}; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs index 5e29de2..73040b7 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs @@ -285,6 +285,18 @@ pub enum MgSmoother { #[cfg(feature = "cuda")] mod device; +/// PERF-2 P3-iii: mark the calling thread as plane lane `lane` (`None` = a +/// solo march). On a lane, every device V-cycle of a CG goes through the +/// plane's rendezvous. Without the `cuda` feature this is a no-op. +#[cfg(feature = "cuda")] +pub use device::{plane_counters, set_lane as set_plane_lane}; +#[cfg(not(feature = "cuda"))] +pub fn set_plane_lane(_lane: Option) {} +#[cfg(not(feature = "cuda"))] +pub fn plane_counters() -> (u64, u64) { + (0, 0) +} + /// Multigrid preconditioner parameters. #[derive(Debug, Clone)] pub struct MultigridParameters { @@ -1256,10 +1268,19 @@ fn run_pcg( let fine: &Level = fine; let cells: &[usize] = cells; let components: &Components = components; + // PERF-2 P3-iii: on a plane lane the device V-cycle goes through the + // rendezvous (one batched launch set for every lane inside a CG); the + // guard marks this CG's extent, whichever `return` below ends it. + #[cfg(feature = "cuda")] + let cg_guard = device.as_ref().and_then(|_| device::cg_enter()); let mut precond = |r: &[f64], z: &mut [f64]| { #[cfg(feature = "cuda")] if let Some(d) = device.as_mut() { - d.apply(r, z); + if cg_guard.is_some() { + device::apply_lane(d, r, z); + } else { + d.apply(r, z); + } return; } hier.apply_preconditioner(r, z); diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson/device.rs index 92050fb..b6e5421 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson/device.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson/device.rs @@ -6,10 +6,13 @@ use super::LevelExport; use cudarc::driver::{ - CudaContext, CudaFunction, CudaModule, CudaSlice, CudaStream, LaunchConfig, PushKernelArg, + CudaContext, CudaFunction, CudaModule, CudaSlice, CudaStream, DevicePtr, DeviceRepr, + LaunchConfig, PushKernelArg, ValidAsZeroBits, }; use cudarc::nvrtc::{CompileOptions, compile_ptx_with_opts}; -use std::sync::{Arc, OnceLock}; +use std::cell::Cell; +use std::collections::BTreeMap; +use std::sync::{Arc, Condvar, Mutex, OnceLock}; const KERNELS: &str = include_str!("../../../kernels/cuda/mg_vcycle.cu"); @@ -23,6 +26,13 @@ struct Runtime { f_prolong: CudaFunction, f_zero: CudaFunction, f_coarsest: CudaFunction, + // The lane-table kernels (P3-iii). + l_half: CudaFunction, + l_res: CudaFunction, + l_restrict: CudaFunction, + l_prolong: CudaFunction, + l_zero: CudaFunction, + l_coarsest: CudaFunction, } static RUNTIME: OnceLock = OnceLock::new(); @@ -49,6 +59,12 @@ fn runtime() -> &'static Runtime { f_prolong: f("mg_prolong"), f_zero: f("mg_zero"), f_coarsest: f("mg_coarsest"), + l_half: f("ml_rb_half"), + l_res: f("ml_residual"), + l_restrict: f("ml_restrict"), + l_prolong: f("ml_prolong"), + l_zero: f("ml_zero"), + l_coarsest: f("ml_coarsest"), _ctx: ctx, stream, _module: module, @@ -310,3 +326,553 @@ impl DeviceVcycle { } } } + +// --------------------------------------------------------------------------- +// PERF-2 P3-iii: the `MarchPlane` — K lanes' V-cycles in ONE batched launch +// set, and the rendezvous that gathers the lanes' preconditioner calls. +// --------------------------------------------------------------------------- + +/// One lane's level as the `ml_*` kernels read it: 14 device pointers then +/// 6 ints (`n, nx, n_cells, n_red, n_black, pad`), mirroring `struct +/// LaneLevel` in `mg_vcycle.cu`. +#[repr(C)] +#[derive(Clone, Copy, Default)] +struct LaneLevel { + ptrs: [u64; 14], + ints: [i32; 6], +} +unsafe impl DeviceRepr for LaneLevel {} +unsafe impl ValidAsZeroBits for LaneLevel {} + +impl DeviceVcycle { + pub(super) fn depth(&self) -> usize { + self.levels.len() + } + + fn lane_level(&self, l: usize, stream: &CudaStream) -> LaneLevel { + let lv = &self.levels[l]; + let pf = |s: &CudaSlice| s.device_ptr(stream).0; + let pu = |s: &CudaSlice| s.device_ptr(stream).0; + LaneLevel { + ptrs: [ + pf(&lv.ae), + pf(&lv.aw), + pf(&lv.an), + pf(&lv.as_), + pf(&lv.ap), + pf(&lv.b), + pf(&lv.x), + pf(&lv.r), + pu(&lv.cells), + pu(&lv.red), + pu(&lv.black), + pu(&lv.coarse_of), + pu(&lv.children_ptr), + pu(&lv.children_idx), + ], + ints: [ + lv.n as i32, + lv.nx, + lv.n_cells as i32, + lv.n_red as i32, + lv.n_black as i32, + 0, + ], + } + } +} + +/// One lane's preconditioner call: its operator on the device, the residual +/// in, the correction out. +pub(super) struct LaneJob<'a> { + pub dev: &'a mut DeviceVcycle, + pub r: &'a [f64], + pub z: &'a mut [f64], +} + +fn cfg_k(n_items: usize, k: usize) -> LaunchConfig { + LaunchConfig { + grid_dim: ((n_items as u32).div_ceil(256).max(1), k as u32, 1), + block_dim: (256, 1, 1), + shared_mem_bytes: 0, + } +} + +/// `z_k = M_k⁻¹ r_k` for every job, the lanes of equal hierarchy depth in +/// one batched launch set (grid.y = the lane). Each lane's result is the +/// K = 1 [`DeviceVcycle::apply`]'s, bit for bit (the same per-cell +/// arithmetic in the same order; gate G-P0). +pub(super) fn apply_batch(jobs: &mut [LaneJob<'_>]) { + let rt = runtime(); + let stream = &rt.stream; + let mut groups: BTreeMap> = BTreeMap::new(); + for (i, j) in jobs.iter().enumerate() { + groups.entry(j.dev.depth()).or_default().push(i); + } + for (depth, members) in groups { + let k = members.len(); + // Residuals up. + for &i in &members { + let job = &mut jobs[i]; + for (dst, &src) in job.dev.r_f32.iter_mut().zip(job.r) { + *dst = src as f32; + } + stream + .memcpy_htod(&job.dev.r_f32, &mut job.dev.levels[0].b) + .expect("upload r"); + } + // The lane tables, one per level. + let tables: Vec> = (0..depth) + .map(|l| { + let rows: Vec = members + .iter() + .map(|&i| jobs[i].dev.lane_level(l, stream)) + .collect(); + stream.memcpy_stod(&rows).expect("lane table") + }) + .collect(); + let max_of = |f: &dyn Fn(&DevLevel) -> usize, l: usize| -> usize { + members + .iter() + .map(|&i| f(&jobs[i].dev.levels[l])) + .max() + .unwrap_or(0) + }; + let sweeps = jobs[members[0]].dev.sweeps; + let smooth = |l: usize| { + let n_red = max_of(&|lv| lv.n_red, l); + let n_black = max_of(&|lv| lv.n_black, l); + let half = |colour: i32| { + let n = if colour == 0 { n_red } else { n_black }; + unsafe { + stream + .launch_builder(&rt.l_half) + .arg(&tables[l]) + .arg(&colour) + .launch(cfg_k(n, k)) + .expect("ml_rb_half"); + } + }; + for _ in 0..sweeps { + half(0); + half(1); + half(1); + half(0); + } + }; + // Down. + for l in 0..depth - 1 { + let n_cells = max_of(&|lv| lv.n_cells, l); + unsafe { + stream + .launch_builder(&rt.l_zero) + .arg(&tables[l]) + .launch(cfg_k(n_cells, k)) + .expect("ml_zero"); + } + smooth(l); + unsafe { + stream + .launch_builder(&rt.l_res) + .arg(&tables[l]) + .launch(cfg_k(n_cells, k)) + .expect("ml_residual"); + } + let n_coarse = max_of(&|lv| lv.n_cells, l + 1); + unsafe { + stream + .launch_builder(&rt.l_restrict) + .arg(&tables[l]) + .arg(&tables[l + 1]) + .launch(cfg_k(n_coarse, k)) + .expect("ml_restrict"); + } + } + // Coarsest. + { + let (k_i, sw_i) = (k as i32, 50i32); + unsafe { + stream + .launch_builder(&rt.l_coarsest) + .arg(&tables[depth - 1]) + .arg(&k_i) + .arg(&sw_i) + .launch(LaunchConfig { + grid_dim: ((k as u32).div_ceil(32), 1, 1), + block_dim: (32, 1, 1), + shared_mem_bytes: 0, + }) + .expect("ml_coarsest"); + } + } + // Up. + for l in (0..depth - 1).rev() { + let n_cells = max_of(&|lv| lv.n_cells, l); + unsafe { + stream + .launch_builder(&rt.l_prolong) + .arg(&tables[l]) + .arg(&tables[l + 1]) + .launch(cfg_k(n_cells, k)) + .expect("ml_prolong"); + } + smooth(l); + } + // Corrections down. + for &i in &members { + let job = &mut jobs[i]; + stream + .memcpy_dtoh(&job.dev.levels[0].x, &mut job.dev.z_f32) + .expect("download z"); + } + stream.synchronize().expect("sync"); + for &i in &members { + let job = &mut jobs[i]; + for &idx in &job.dev.fine_cells { + job.z[idx as usize] = job.dev.z_f32[idx as usize] as f64; + } + } + } +} + +// --- The rendezvous --------------------------------------------------------- +// +// K lane threads each run an unchanged march. A lane inside a CG (between +// `cg_enter` and the guard's drop) that calls the preconditioner submits its +// job and blocks; when every lane currently inside a CG has submitted, the +// last arrival serves the whole batch and wakes the rest. A lane LEAVING its +// CG can complete a batch too (the others were waiting for it). Lanes +// outside a CG are not counted, so they never hold a batch up. + +thread_local! { + static LANE: Cell> = const { Cell::new(None) }; +} + +/// Mark this thread as plane lane `lane` (`None` = a solo march: the K = 1 +/// path, unchanged). +pub fn set_lane(lane: Option) { + LANE.with(|l| l.set(lane)); +} + +pub fn lane() -> Option { + LANE.with(|l| l.get()) +} + +struct PendingJob { + dev: *mut DeviceVcycle, + r: *const f64, + r_len: usize, + z: *mut f64, + z_len: usize, +} +// The submitting thread blocks until its job is served, so the pointers +// outlive the batch. +unsafe impl Send for PendingJob {} + +#[derive(Default)] +struct PlaneState { + in_cg: usize, + pending: Vec, + generation: u64, + /// Batches served and lanes served (the plane's own bookkeeping). + batches: u64, + lanes_served: u64, +} + +static PLANE: Mutex = Mutex::new(PlaneState { + in_cg: 0, + pending: Vec::new(), + generation: 0, + batches: 0, + lanes_served: 0, +}); +static PLANE_CV: Condvar = Condvar::new(); + +fn plane_lock() -> std::sync::MutexGuard<'static, PlaneState> { + PLANE.lock().unwrap_or_else(|e| e.into_inner()) +} + +fn serve(st: &mut PlaneState) { + let pending = std::mem::take(&mut st.pending); + let mut jobs: Vec> = pending + .iter() + .map(|p| unsafe { + LaneJob { + dev: &mut *p.dev, + r: std::slice::from_raw_parts(p.r, p.r_len), + z: std::slice::from_raw_parts_mut(p.z, p.z_len), + } + }) + .collect(); + apply_batch(&mut jobs); + st.batches += 1; + st.lanes_served += pending.len() as u64; + st.generation += 1; + PLANE_CV.notify_all(); +} + +/// Held by a lane for the duration of one CG solve. +pub(super) struct CgGuard(()); + +impl Drop for CgGuard { + fn drop(&mut self) { + let mut st = plane_lock(); + st.in_cg -= 1; + if st.in_cg > 0 && st.pending.len() >= st.in_cg { + serve(&mut st); + } + } +} + +/// Enter a CG solve on this lane (`None` when this thread is not a lane). +pub(super) fn cg_enter() -> Option { + lane()?; + plane_lock().in_cg += 1; + Some(CgGuard(())) +} + +/// This lane's preconditioner call inside the plane: submit, and either +/// serve the batch (last arrival) or wait for it. +pub(super) fn apply_lane(dev: &mut DeviceVcycle, r: &[f64], z: &mut [f64]) { + let mut st = plane_lock(); + let my_generation = st.generation; + st.pending.push(PendingJob { + dev: std::ptr::from_mut(dev), + r: r.as_ptr(), + r_len: r.len(), + z: z.as_mut_ptr(), + z_len: z.len(), + }); + if st.pending.len() >= st.in_cg { + serve(&mut st); + return; + } + while st.generation == my_generation { + st = PLANE_CV.wait(st).unwrap_or_else(|e| e.into_inner()); + } +} + +/// `(batches served, lane calls served)` so far in this process. +pub fn plane_counters() -> (u64, u64) { + let st = plane_lock(); + (st.batches, st.lanes_served) +} + +#[cfg(test)] +mod tests { + //! G-P0 (`docs/perf2_campaign.md`, P3-iii): K distinct masked operators + //! through the batched path equal each one's K = 1 device apply bit for + //! bit, and the rendezvous serves uneven lanes exactly. GPU only: + //! `RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --lib device::tests`. + use super::*; + use crate::solvers::incompressible::poisson::{ + MgSmoother, MultigridParameters, PoissonProblem, export_hierarchy, + }; + + /// A masked channel operator with a hole at `(cx, cy)` and a seeded + /// right-hand side (the go/no-go benchmark's construction). + fn problem(nx: usize, ny: usize, cx: f64, cy: f64, seed: u64) -> PoissonProblem { + let mut p = PoissonProblem::new(nx, ny); + let (dx, dy, dt) = (2.5 / nx as f64, 0.41 / ny as f64, 3.24e-4); + let (ae, an) = (dt * dy / dx, dt * dx / dy); + let hole = |i: usize, j: usize| { + let (x, y) = ((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy); + (x - cx).powi(2) + (y - cy).powi(2) < 0.08 * 0.08 + }; + for j in 0..ny { + for i in 0..nx { + let idx = j * nx + i; + if hole(i, j) { + p.active[idx] = false; + continue; + } + if i + 1 < nx && !hole(i + 1, j) { + p.ae[idx] = ae; + } + if i > 0 && !hole(i - 1, j) { + p.aw[idx] = ae; + } + if j + 1 < ny && !hole(i, j + 1) { + p.an[idx] = an; + } + if j > 0 && !hole(i, j - 1) { + p.as_[idx] = an; + } + if i + 1 == nx { + p.extra_diag[idx] = 2.0 * ae; + } + } + } + let mut state = seed | 1; + for idx in 0..nx * ny { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + p.rhs[idx] = if p.active[idx] { + 1e-3 * ((state >> 11) as f64 / (1u64 << 53) as f64 - 0.5) + } else { + 0.0 + }; + } + p + } + + fn params() -> MultigridParameters { + MultigridParameters { + smoother: MgSmoother::RedBlack, + ..MultigridParameters::default() + } + } + + fn device_of(p: &PoissonProblem) -> DeviceVcycle { + DeviceVcycle::new( + &export_hierarchy(p, ¶ms()), + params().smoother_sweeps.max(1), + ) + } + + fn bits(v: &[f64]) -> Vec { + v.iter().map(|x| x.to_bits()).collect() + } + + /// K operators with different holes: solo applies, then one batch. + #[test] + fn batched_lanes_equal_the_solo_apply_bit_for_bit() { + let (nx, ny) = (186usize, 31usize); + let probs: Vec = [(0.2, 0.2, 3u64), (0.6, 0.15, 7), (1.1, 0.25, 11)] + .iter() + .map(|&(cx, cy, s)| problem(nx, ny, cx, cy, s)) + .collect(); + let n = nx * ny; + let mut solo: Vec> = Vec::new(); + for p in &probs { + let mut d = device_of(p); + let mut z = vec![0.0; n]; + d.apply(&p.rhs, &mut z); + solo.push(z); + } + assert!( + solo[0].iter().any(|v| *v != 0.0), + "the solo apply produced zeros" + ); + // K = 1 through the batched path (the lane kernels vs the K = 1 kernels). + { + let mut d = device_of(&probs[0]); + let mut z = vec![0.0; n]; + apply_batch(&mut [LaneJob { + dev: &mut d, + r: &probs[0].rhs, + z: &mut z, + }]); + assert_eq!( + bits(&z), + bits(&solo[0]), + "K = 1 through the lane kernels differs" + ); + } + // K = 3. + let mut devs: Vec = probs.iter().map(device_of).collect(); + let mut zs: Vec> = vec![vec![0.0; n]; 3]; + { + let mut jobs: Vec> = devs + .iter_mut() + .zip(probs.iter()) + .zip(zs.iter_mut()) + .map(|((dev, p), z)| LaneJob { dev, r: &p.rhs, z }) + .collect(); + apply_batch(&mut jobs); + } + for k in 0..3 { + assert_eq!( + bits(&zs[k]), + bits(&solo[k]), + "lane {k} differs from its solo apply" + ); + } + // A second batch on the same devices (persistent buffers) with the + // lanes in another order and one lane left out. + let mut z1 = vec![0.0; n]; + let mut z2 = vec![0.0; n]; + { + let (a, b) = devs.split_at_mut(2); + let mut jobs = [ + LaneJob { + dev: &mut b[0], + r: &probs[2].rhs, + z: &mut z2, + }, + LaneJob { + dev: &mut a[1], + r: &probs[1].rhs, + z: &mut z1, + }, + ]; + apply_batch(&mut jobs); + } + assert_eq!(bits(&z2), bits(&solo[2])); + assert_eq!(bits(&z1), bits(&solo[1])); + } + + /// Three lane threads with 2 / 4 / 5 preconditioner calls each inside + /// one CG guard: every call's correction equals the solo apply of the + /// same residual, the lanes leave at different times, and no thread + /// hangs. + #[test] + fn the_rendezvous_serves_uneven_lanes_exactly() { + let (nx, ny) = (124usize, 21usize); + let n = nx * ny; + let probs: Vec = [(0.2, 0.2, 5u64), (0.7, 0.2, 9), (1.3, 0.22, 13)] + .iter() + .map(|&(cx, cy, s)| problem(nx, ny, cx, cy, s)) + .collect(); + let calls = [2usize, 4, 5]; + // Expected: the solo apply of rhs × (1 + c) for call c. + let mut expected: Vec>> = Vec::new(); + for (k, p) in probs.iter().enumerate() { + let mut d = device_of(p); + let mut per_call = Vec::new(); + for c in 0..calls[k] { + let r: Vec = p.rhs.iter().map(|v| v * (1.0 + c as f64)).collect(); + let mut z = vec![0.0; n]; + d.apply(&r, &mut z); + per_call.push(bits(&z)); + } + expected.push(per_call); + } + let (b0, l0) = plane_counters(); + let handles: Vec<_> = probs + .into_iter() + .enumerate() + .map(|(k, p)| { + let calls_k = calls[k]; + std::thread::spawn(move || { + set_lane(Some(k)); + let mut d = device_of(&p); + let guard = cg_enter().expect("lane set"); + let mut out = Vec::new(); + for c in 0..calls_k { + let r: Vec = p.rhs.iter().map(|v| v * (1.0 + c as f64)).collect(); + let mut z = vec![0.0; n]; + apply_lane(&mut d, &r, &mut z); + out.push(bits(&z)); + } + drop(guard); + out + }) + }) + .collect(); + let results: Vec>> = handles + .into_iter() + .map(|h| h.join().expect("lane thread")) + .collect(); + for k in 0..3 { + assert_eq!(results[k].len(), calls[k]); + for c in 0..calls[k] { + assert_eq!(results[k][c], expected[k][c], "lane {k} call {c} differs"); + } + } + let (b1, l1) = plane_counters(); + assert_eq!(l1 - l0, 11, "every lane call served once"); + assert!(b1 - b0 >= 5 && b1 - b0 <= 11, "batches {}", b1 - b0); + } +} diff --git a/crates/specialized/rtx-cfd/tests/gpu_kernel_comprehensive_tests.rs b/crates/specialized/rtx-cfd/tests/gpu_kernel_comprehensive_tests.rs index c8f4325..0dfc988 100644 --- a/crates/specialized/rtx-cfd/tests/gpu_kernel_comprehensive_tests.rs +++ b/crates/specialized/rtx-cfd/tests/gpu_kernel_comprehensive_tests.rs @@ -420,7 +420,10 @@ mod gpu_kernel_tests { // Check that center has highest value let center_val = phi_result[center]; - assert!(center_val < 0.0, "Center value should be negative (∇²φ = f, a positive source)"); + assert!( + center_val < 0.0, + "Center value should be negative (∇²φ = f, a positive source)" + ); Ok(()) } diff --git a/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs b/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs index 3ec6a3d..c45393b 100644 --- a/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs +++ b/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs @@ -174,9 +174,21 @@ fn batched_device_vcycle_go_no_go() { cells: up_u(&l.cells), red: up_u(&l.red), black: up_u(&l.black), - coarse_of: up_u(if l.coarse_of.is_empty() { &[0u32][..] } else { &l.coarse_of }), - children_ptr: up_u(if l.children_ptr.is_empty() { &[0u32][..] } else { &l.children_ptr }), - children_idx: up_u(if l.children_idx.is_empty() { &[0u32][..] } else { &l.children_idx }), + coarse_of: up_u(if l.coarse_of.is_empty() { + &[0u32][..] + } else { + &l.coarse_of + }), + children_ptr: up_u(if l.children_ptr.is_empty() { + &[0u32][..] + } else { + &l.children_ptr + }), + children_idx: up_u(if l.children_idx.is_empty() { + &[0u32][..] + } else { + &l.children_idx + }), n_cells: l.cells.len(), n_red: l.red.len(), n_black: l.black.len(), @@ -206,7 +218,11 @@ fn batched_device_vcycle_go_no_go() { // Upload the residual into level 0's b. stream.memcpy_htod(r_f32, &mut dev[0].b).unwrap(); let half = |lv: &mut DeviceLevel, colour: u8| { - let (list, n_list) = if colour == 0 { (&lv.red, lv.n_red) } else { (&lv.black, lv.n_black) }; + let (list, n_list) = if colour == 0 { + (&lv.red, lv.n_red) + } else { + (&lv.black, lv.n_black) + }; let n_i = lv.n as i32; let n_list_i = n_list as i32; unsafe { @@ -371,7 +387,10 @@ fn batched_device_vcycle_go_no_go() { let (db, sb) = cmp(&b[..n], b_ref); let (dx, sx) = cmp(&x[..n], x_ref); let (dr, sr) = cmp(&r[..n], r_ref); - println!(" level {l} ({} cells): |Δb| {db:.3e} / {sb:.3e}, |Δx| {dx:.3e} / {sx:.3e}, |Δr| {dr:.3e} / {sr:.3e}", levels[l].cells.len()); + println!( + " level {l} ({} cells): |Δb| {db:.3e} / {sb:.3e}, |Δx| {dx:.3e} / {sx:.3e}, |Δr| {dr:.3e} / {sr:.3e}", + levels[l].cells.len() + ); } } let mut worst = 0.0_f64; @@ -381,8 +400,13 @@ fn batched_device_vcycle_go_no_go() { worst = worst.max(d); } } - println!(" K = {k}: device V-cycle vs CPU f32 reference: max |Δz| {worst:.3e} on a scale of {z_scale:.3e}"); - assert!(worst < 1e-4 * z_scale, "device V-cycle disagrees with the CPU one"); + println!( + " K = {k}: device V-cycle vs CPU f32 reference: max |Δz| {worst:.3e} on a scale of {z_scale:.3e}" + ); + assert!( + worst < 1e-4 * z_scale, + "device V-cycle disagrees with the CPU one" + ); // Timing: 50 V-cycles, per march. let reps = 50; diff --git a/crates/specialized/rtx-fsi/tests/fsi2_overset_plane.rs b/crates/specialized/rtx-fsi/tests/fsi2_overset_plane.rs new file mode 100644 index 0000000..ba3fac8 --- /dev/null +++ b/crates/specialized/rtx-fsi/tests/fsi2_overset_plane.rs @@ -0,0 +1,209 @@ +//! PERF-2 P3-iii (`docs/perf2_campaign.md` in omni-cortex): the +//! `MarchPlane` — K overset FSI marches in ONE process, one lane thread +//! each, the device V-cycles of every lane inside a CG gathered into one +//! batched launch set by the plane rendezvous in `rtx-cfd`'s Poisson +//! driver. Every lane runs the unchanged `run_march_overset`. +//! +//! Knobs (all the `RTX_FSI2O_*` knobs of `fsi2_on_the_overset` apply to +//! every lane; these add the plane): +//! - `RTX_FSI2O_PLANE_FILE`: the points, one per line, `u_mean e_s tag` +//! (the `p7_points.txt` format; `#` comments and blank lines skipped). +//! Without it the test returns. +//! - `RTX_FSI2O_PLANE_DIR`: each lane's CSV goes to `/.csv`. +//! - `RTX_FSI2O_PLANE_STACK_MB`: the lane threads' stack (default 64). +//! Requires `RTX_FSI2O_MG_RB=1 RTX_FSI2O_MG_DEVICE=1` on a +//! `--features rtx-cfd/cuda` build; the test fails loudly otherwise (a +//! CPU plane is K separate processes, not this driver). +//! +//! `RTX_FSI2O_PLANE_FILE=points.txt RTX_FSI2O_PLANE_DIR=$L RTX_FSI2O_MG_RB=1 RTX_FSI2O_MG_DEVICE=1 RTX_CUDA_ARCH=sm_120 fsi2_overset_plane --nocapture` + +mod fsi2_harness; + +use fsi2_harness::overset_march::{OversetMarchConfig, run_march_overset}; +use fsi2_harness::{BenchmarkCase, FSI2, FSI3}; +use rtx_cfd::solvers::incompressible::{plane_counters, set_plane_lane}; +use std::time::Instant; + +struct Point { + u_mean: f64, + e_s: f64, + tag: String, +} + +fn read_points(path: &str) -> Vec { + let text = std::fs::read_to_string(path).unwrap_or_else(|e| panic!("{path}: {e}")); + text.lines() + .map(str::trim) + .filter(|l| !l.is_empty() && !l.starts_with('#')) + .map(|l| { + let f: Vec<&str> = l.split_whitespace().collect(); + assert!( + f.len() >= 3, + "plane point line needs `u_mean e_s tag`: {l:?}" + ); + Point { + u_mean: f[0].parse().expect("u_mean"), + e_s: f[1].parse().expect("e_s"), + tag: f[2].to_string(), + } + }) + .collect() +} + +struct LaneSummary { + tag: String, + wall: f64, + steps: usize, + subit: f64, + retries: usize, + uy_amp_mm: f64, + uy_mid_mm: f64, + frequency: Option, + drag_mid: f64, + lift_amp: f64, + death: Option<(usize, f64, String)>, +} + +#[test] +fn fsi2_overset_plane() { + let Ok(points_path) = std::env::var("RTX_FSI2O_PLANE_FILE") else { + return; + }; + let dir = std::env::var("RTX_FSI2O_PLANE_DIR").unwrap_or_else(|_| ".".into()); + let stack_mb: usize = std::env::var("RTX_FSI2O_PLANE_STACK_MB") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(64); + let device_on = std::env::var("RTX_FSI2O_MG_DEVICE").is_ok_and(|v| v == "1") + && std::env::var("RTX_FSI2O_MG_RB").is_ok_and(|v| v == "1"); + assert!( + device_on, + "the plane needs RTX_FSI2O_MG_RB=1 RTX_FSI2O_MG_DEVICE=1 on a cuda build" + ); + let points = read_points(&points_path); + let k = points.len(); + assert!(k >= 1, "no plane points in {points_path}"); + + let config = OversetMarchConfig::from_env( + "FSI2O", + OversetMarchConfig { + ny: 41, + flag_nx: 35, + t_release: 6.0, + t_end: 7.0, + subcycle: 1, + tol_floor: 2e-4, + rtol: 1e-2, + stall_accept: 5.0, + max_subiterations: 12, + coupler: "aitken".into(), + reuse: 2, + initial_relaxation: 0.5, + c1_interface: false, + predictor: "structure".into(), + sweeps: 100, + max_rounds: 3, + csv_path: None, + trace_steps: 0, + }, + ); + let base = if std::env::var("RTX_FSI2O_CASE").as_deref() == Ok("FSI3") { + FSI3 + } else { + FSI2 + }; + println!( + " MARCH PLANE: {k} lanes in one process ({} base, ny {}, t_end {}, coupler {}, s = {}); points from {points_path}, CSVs under {dir}", + base.name, config.ny, config.t_end, config.coupler, config.subcycle + ); + for (i, p) in points.iter().enumerate() { + println!( + " lane {i}: {} u_mean {} (Re {:.0}) e_s {:.4e} (E/E0 {:.4})", + p.tag, + p.u_mean, + 100.0 * p.u_mean, + p.e_s, + p.e_s / base.e_s + ); + } + + let t_plane = Instant::now(); + let handles: Vec<_> = points + .into_iter() + .enumerate() + .map(|(i, p)| { + let mut cfg = config.clone(); + cfg.csv_path = Some(format!("{dir}/{}.csv", p.tag)); + let case = BenchmarkCase { + u_mean: p.u_mean, + e_s: p.e_s, + ..base + }; + std::thread::Builder::new() + .name(format!("lane-{i}-{}", p.tag)) + .stack_size(stack_mb << 20) + .spawn(move || { + set_plane_lane(Some(i)); + let t0 = Instant::now(); + let r = run_march_overset(case, &cfg); + let m = &r.result; + let w = m.window(3.0); + LaneSummary { + tag: p.tag, + wall: t0.elapsed().as_secs_f64(), + steps: m.coupled_steps, + subit: m.mean_subiterations, + retries: m.retried_steps, + uy_amp_mm: w.uy_amp * 1e3, + uy_mid_mm: w.uy_mid * 1e3, + frequency: w.frequency, + drag_mid: w.drag_mid, + lift_amp: w.lift_amp, + death: r.death, + } + }) + .expect("lane thread") + }) + .collect(); + let mut summaries = Vec::with_capacity(k); + let mut failures = Vec::new(); + for (i, h) in handles.into_iter().enumerate() { + match h.join() { + Ok(s) => summaries.push(s), + Err(_) => failures.push(i), + } + } + let plane_wall = t_plane.elapsed().as_secs_f64(); + let (batches, served) = plane_counters(); + println!( + " MARCH PLANE done: {k} lanes in {plane_wall:.0} s wall; {batches} batched V-cycle launches served {served} lane calls ({:.2} lanes per batch)", + if batches > 0 { + served as f64 / batches as f64 + } else { + 0.0 + } + ); + for s in &summaries { + println!( + " lane {}: {:.0} s wall, {} coupled steps, {:.1} subit/step, {} retries; last 3 s: uy(A) = {:.3} ± {:.3} mm, f = {:?} Hz, drag mid {:.2}, lift amp {:.2}; death {:?}", + s.tag, + s.wall, + s.steps, + s.subit, + s.retries, + s.uy_mid_mm, + s.uy_amp_mm, + s.frequency, + s.drag_mid, + s.lift_amp, + s.death + ); + } + assert!(failures.is_empty(), "lane threads panicked: {failures:?}"); + let dead: Vec<&str> = summaries + .iter() + .filter(|s| s.death.is_some()) + .map(|s| s.tag.as_str()) + .collect(); + assert!(dead.is_empty(), "lanes died: {dead:?}"); +}