diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs index 3c12b2a..ef67563 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs @@ -71,6 +71,9 @@ pub struct EmbeddedParameters { /// and the stop stay f64; only the preconditioner runs in single /// precision. No effect with [`PoissonSolverKind::Sor`]. pub poisson_precision: MgPrecision, + /// The multigrid smoother ordering (PERF-2; default lexicographic, the + /// recorded regime). + pub poisson_smoother: super::poisson::MgSmoother, /// Convective face values in the explicit predictor (default /// [`ConvectionScheme::Upwind`], which is bit-identical to the fixed-grid /// PISO). The TVD schemes add SIMPLE's limited correction to each @@ -93,6 +96,7 @@ impl Default for EmbeddedParameters { boundaries: AleBoundaries::default(), poisson_solver: PoissonSolverKind::Sor, poisson_precision: MgPrecision::F64, + poisson_smoother: super::poisson::MgSmoother::Lexicographic, convection_scheme: ConvectionScheme::Upwind, } } @@ -215,6 +219,12 @@ impl EmbeddedPisoSolver { self.parameters.poisson_precision = precision; } + /// The multigrid smoother ordering (PERF-2's red-black knob); the + /// cached operator is rebuilt on the next solve. + pub fn set_poisson_smoother(&mut self, smoother: super::poisson::MgSmoother) { + self.parameters.poisson_smoother = smoother; + } + /// Mask hysteresis for the moving-body rebuild, as a fraction of the /// min cell size (default 0, exactly the plain rebuild). With a band, /// a cell within `band * h_min` of the surface keeps the diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs index 9eb44c4..8ec7fd7 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs @@ -299,6 +299,7 @@ impl EmbeddedPisoSolver { &mut p_prime, &MultigridParameters { precision: self.parameters.poisson_precision, + smoother: self.parameters.poisson_smoother, ..MultigridParameters::default() }, inner_stop, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs index 4cccedb..d901b6f 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs @@ -65,8 +65,8 @@ pub use piso::{PisoParameters, PisoResult, PisoSolver}; #[cfg(feature = "cuda")] pub use piso_gpu::PisoGpuSolver; pub use poisson::{ - MgPrecision, MultigridParameters, PcgCache, PoissonProblem, PoissonSolution, PoissonSolverKind, - solve_multigrid_pcg, solve_multigrid_pcg_cached, + MgPrecision, MgSmoother, MultigridParameters, PcgCache, PoissonProblem, PoissonSolution, + PoissonSolverKind, solve_multigrid_pcg, solve_multigrid_pcg_cached, }; 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 cd2e78c..91d0b0c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs @@ -268,11 +268,27 @@ pub enum MgPrecision { F32, } +/// The V-cycle's smoother ordering (PERF-2, `docs/perf2_campaign.md`). +/// `Lexicographic` is the recorded regime (row-major symmetric +/// Gauss–Seidel, a dependency chain through the division per cell); +/// `RedBlack` updates the two colours of the five-point stencil in turn — +/// each colour a map of independent cells (threads, vectors, the GPU) — +/// and is a different preconditioner, gated by the noise probe and the +/// anchor's band, never bit-identical. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum MgSmoother { + #[default] + Lexicographic, + RedBlack, +} + /// Multigrid preconditioner parameters. #[derive(Debug, Clone)] pub struct MultigridParameters { /// Precision of the V-cycle (see [`MgPrecision`]). pub precision: MgPrecision, + /// Smoother ordering (see [`MgSmoother`]). + pub smoother: MgSmoother, /// Symmetric Gauss–Seidel sweeps before AND after the coarse correction /// (default 2; a value of 0 is treated as 1). One count for both on /// purpose: unequal pre/post counts make the V-cycle non-symmetric and @@ -290,6 +306,7 @@ impl Default for MultigridParameters { fn default() -> Self { Self { precision: MgPrecision::F64, + smoother: MgSmoother::Lexicographic, smoother_sweeps: 2, coarsest_cells: 32, max_iterations: 500, @@ -431,6 +448,10 @@ struct Level { ap: Vec, /// Row-major indices of the active cells. cells: Vec, + /// The active cells with `i + j` even / odd, each row-major (the + /// red-black smoother's two independent maps on the five-point stencil). + red: Vec, + black: Vec, /// Fine index → coarse index (empty on the coarsest level). coarse_of: Vec, } @@ -491,6 +512,17 @@ impl Level { } } let cells: Vec = (0..n).filter(|&idx| active[idx]).collect(); + let parity = |idx: usize| (idx % nx + idx / nx) % 2; + let red: Vec = cells + .iter() + .copied() + .filter(|&idx| parity(idx) == 0) + .collect(); + let black: Vec = cells + .iter() + .copied() + .filter(|&idx| parity(idx) == 1) + .collect(); let cast = |v: &[f64]| v.iter().map(|&x| T::from_f64(x)).collect::>(); Self { ae: cast(&problem.ae), @@ -501,6 +533,8 @@ impl Level { problem, active, cells, + red, + black, coarse_of: Vec::new(), } } @@ -559,6 +593,33 @@ impl Level { } } + /// One symmetric red-black Gauss–Seidel sweep: red, black, black, red — + /// within a colour every cell reads only the other colour, so each + /// half-sweep is a map (symmetric as a preconditioner, like the + /// lexicographic pair). + fn symmetric_gs_rb(&self, b: &[T], x: &mut [T]) { + for &idx in &self.red { + x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx]; + } + for &idx in &self.black { + x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx]; + } + for &idx in &self.black { + x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx]; + } + for &idx in &self.red { + x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx]; + } + } + + /// One symmetric sweep in the chosen ordering. + fn smooth(&self, b: &[T], x: &mut [T], smoother: MgSmoother) { + match smoother { + MgSmoother::Lexicographic => self.symmetric_gs(b, x), + MgSmoother::RedBlack => self.symmetric_gs_rb(b, x), + } + } + /// Galerkin coarsening: coarse face coefficient = sum of the fine /// coefficients across that coarse face, coarse `extra_diag` = sum of /// the children's. Returns the coarse problem and the parent map. @@ -601,6 +662,7 @@ pub(crate) struct Hierarchy { levels: Vec>, work: Vec>, sweeps: usize, + smoother: MgSmoother, } impl Hierarchy { @@ -630,6 +692,7 @@ impl Hierarchy { levels, work, sweeps: params.smoother_sweeps.max(1), + smoother: params.smoother, } } @@ -672,7 +735,7 @@ impl Hierarchy { wf.x[idx] = T::ZERO; } for _ in 0..self.sweeps { - fine.symmetric_gs(&wf.b, &mut wf.x); + fine.smooth(&wf.b, &mut wf.x, self.smoother); } fine.residual(&wf.b, &wf.x, &mut wf.r); for &idx in &coarse.cells { @@ -695,7 +758,7 @@ impl Hierarchy { wb.x[idx] = T::ZERO; } for _ in 0..COARSEST_SWEEPS { - bottom.symmetric_gs(&wb.b, &mut wb.x); + bottom.smooth(&wb.b, &mut wb.x, self.smoother); } } // Up: prolongate, smooth. @@ -707,7 +770,7 @@ impl Hierarchy { wf.x[idx] += T::from_f64(COARSE_CORRECTION) * wc.x[fine.coarse_of[idx]]; } for _ in 0..self.sweeps { - fine.symmetric_gs(&wf.b, &mut wf.x); + fine.smooth(&wf.b, &mut wf.x, self.smoother); } } for &idx in &levels[0].cells { @@ -762,6 +825,7 @@ struct OperatorKey { coefficients: Vec, smoother_sweeps: usize, coarsest_cells: usize, + smoother: MgSmoother, } impl OperatorKey { @@ -782,6 +846,7 @@ impl OperatorKey { coefficients, smoother_sweeps: params.smoother_sweeps, coarsest_cells: params.coarsest_cells, + smoother: params.smoother, } } @@ -790,6 +855,7 @@ impl OperatorKey { && self.ny == problem.ny && self.smoother_sweeps == params.smoother_sweeps && self.coarsest_cells == params.coarsest_cells + && self.smoother == params.smoother && self.active == problem.active && self.coefficients.iter().copied().eq(problem .ae diff --git a/crates/specialized/rtx-cfd/tests/poisson_redblack.rs b/crates/specialized/rtx-cfd/tests/poisson_redblack.rs new file mode 100644 index 0000000..482f4f8 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/poisson_redblack.rs @@ -0,0 +1,108 @@ +//! PERF-2 (`docs/perf2_campaign.md`): the red-black symmetric Gauss–Seidel +//! smoother is a different preconditioner, not a bit-identical one: the pin +//! is that it solves the same masked problem to the same stop, that its +//! solution agrees with the lexicographic one to the solver's tolerance, +//! and that the cached red-black solve is bit-identical to the uncached +//! red-black solve (the cache keys on the smoother). + +use rtx_cfd::solvers::incompressible::{ + MgSmoother, MultigridParameters, PcgCache, PoissonProblem, solve_multigrid_pcg, + solve_multigrid_pcg_cached, +}; + +fn problem(nx: usize, ny: usize, seed: u64) -> PoissonProblem { + let mut p = PoissonProblem::new(nx, ny); + let (dx, dy, dt) = (1.0 / nx as f64, 0.41 / ny as f64, 1e-3); + 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 - 0.2).powi(2) + (y - 0.2).powi(2) < 0.05 * 0.05 + }; + 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-6 * ((state >> 11) as f64 / (1u64 << 53) as f64 - 0.5) + } else { + 0.0 + }; + } + p +} + +#[test] +fn red_black_solves_the_masked_problem_and_caches_exactly() { + let (nx, ny) = (96, 40); + let prob = problem(nx, ny, 5); + let tol = 1e-12; + let lex = MultigridParameters::default(); + let rb = MultigridParameters { + smoother: MgSmoother::RedBlack, + ..MultigridParameters::default() + }; + let (mut p_lex, mut p_rb) = (vec![0.0; nx * ny], vec![0.0; nx * ny]); + let s_lex = solve_multigrid_pcg(&prob, &mut p_lex, &lex, tol, None); + let s_rb = solve_multigrid_pcg(&prob, &mut p_rb, &rb, tol, None); + assert!(s_lex.converged && s_rb.converged); + assert!( + prob.residual_l1(&p_rb) < tol, + "red-black residual {:.3e}", + prob.residual_l1(&p_rb) + ); + let scale = p_lex.iter().fold(0.0_f64, |m, v| m.max(v.abs())); + let diff = p_lex + .iter() + .zip(&p_rb) + .fold(0.0_f64, |m, (a, b)| m.max((a - b).abs())); + println!( + " lexicographic {} iterations vs red-black {} iterations; solutions differ by {:.3e} on a scale of {:.3e}", + s_lex.iterations, s_rb.iterations, diff, scale + ); + assert!( + diff < 1e-6 * scale, + "red-black and lexicographic disagree: {diff:.3e} of {scale:.3e}" + ); + // The cached red-black solve is the uncached one, bit for bit. + let mut cache = PcgCache::default(); + for k in 0..3u64 { + let mut q = prob.clone(); + q.rhs = problem(nx, ny, 20 + k).rhs; + let (mut a, mut b) = (vec![0.0; nx * ny], vec![0.0; nx * ny]); + let sa = solve_multigrid_pcg(&q, &mut a, &rb, tol, None); + let sb = solve_multigrid_pcg_cached(&q, &mut b, &rb, tol, None, &mut cache); + assert!(a.iter().zip(&b).all(|(x, y)| x.to_bits() == y.to_bits())); + assert_eq!(sa.iterations, sb.iterations); + } + // Switching the smoother is a cache miss (the key carries it), still exact. + let (mut a, mut b) = (vec![0.0; nx * ny], vec![0.0; nx * ny]); + let sa = solve_multigrid_pcg(&prob, &mut a, &lex, tol, None); + let sb = solve_multigrid_pcg_cached(&prob, &mut b, &lex, tol, None, &mut cache); + assert!(a.iter().zip(&b).all(|(x, y)| x.to_bits() == y.to_bits())); + assert_eq!(sa.iterations, sb.iterations); +} diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/march.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/march.rs index 097e59b..f6ae5b2 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/march.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/march.rs @@ -384,6 +384,14 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult { solver.set_poisson_precision(rtx_cfd::solvers::incompressible::MgPrecision::F32); println!(" poisson V-cycle precision: F32 (M1 probe)"); } + // PERF-2: `RTX_FSI2_MG_RB=1` — the red-black smoother on the embedded + // solver (the noise probe's regime gate). + if std::env::var("RTX_FSI2_MG_RB").is_ok_and(|v| v == "1") { + solver.set_poisson_smoother(rtx_cfd::solvers::incompressible::MgSmoother::RedBlack); + println!( + " multigrid smoother: RED-BLACK symmetric Gauss–Seidel (PERF-2 regime, RTX_FSI2_MG_RB)" + ); + } let dt_fluid = harness.dt_fluid; let dt = dt_fluid * subcycle as f64; let interface = &harness.interface; diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs index 21eeec5..e44fa3c 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs @@ -300,6 +300,16 @@ impl OversetFluid { }, poisson_solver: PoissonSolverKind::Multigrid, poisson_precision: MgPrecision::F64, + // PERF-2: `RTX_FSI2O_MG_RB=1` selects the red-black smoother + // (a regime change, band-gated; default = the recorded regime). + poisson_smoother: if std::env::var("RTX_FSI2O_MG_RB").is_ok_and(|v| v == "1") { + println!( + " multigrid smoother: RED-BLACK symmetric Gauss–Seidel (PERF-2 regime, RTX_FSI2O_MG_RB)" + ); + rtx_cfd::solvers::incompressible::MgSmoother::RedBlack + } else { + rtx_cfd::solvers::incompressible::MgSmoother::Lexicographic + }, convection_scheme: bg_convection(), }, )?;