From da2945fcf5be82851facb286ffbffb56d6f27e94 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 25 Sep 2026 04:09:51 -0500 Subject: [PATCH] embedded3 R6-3: the hierarchy rebuild's level 0 on the device (RTX_E3_POISSON_DEVICE=1) DeviceCg::new_device: level 0 by refresh_device, level 1's operator coarsened on the device (Level::coarsen per coarse cell, children in ascending index order) and brought down, levels 1.. by Hierarchy's loop (Hierarchy::chain_below, export_levels), level 0's children lists on the device. Falls back to DeviceCg::new when the hierarchy would stop at level 0. RTX_E3_BAND_CHECK=1: every level against export_hierarchy, the fine structures as refresh_device's check, the raw components' singular flag. DeviceCg::new gains profile laps; DeviceVcycle::new split into upload_level + from_levels, export_hierarchy into export_levels (same code). Co-Authored-By: Claude Opus 5.5 (1M context) --- .../rtx-cfd/src/kernels/cuda/e3_pset.cu | 105 +++++ .../embedded3/poisson/device.rs | 84 ++-- .../embedded3/poisson/device_cg.rs | 19 +- .../embedded3/poisson/device_cg/setup.rs | 377 +++++++++++++++++- .../embedded3/poisson/export.rs | 13 +- .../embedded3/poisson/hierarchy.rs | 23 ++ .../embedded3/step/device/cut.rs | 21 +- 7 files changed, 581 insertions(+), 61 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_pset.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_pset.cu index 6a5636e..829f8e4 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_pset.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_pset.cu @@ -461,3 +461,108 @@ extern "C" __global__ void e3_ps_scan_values(long long n, const u32* __restrict_ u32 ex = ps_block_excl_scan(i < n ? vals[i] : 0u, &total); if (i < n) out[i] = offsets[blockIdx.x] + ex; } + +/* ---- the hierarchy's level 1 from the device level 0 (Level::coarsen, per coarse cell) ---- */ + +struct PsCoarse { + int nxc, nyc, nzc, pad; +}; + +/* The fine children box of coarse cell (ic, jc, kc): i in [2 ic, ic == nxc - 1 ? nx - 1 : 2 ic + 1]. */ +__device__ __forceinline__ void ps_box(const PsGrid& g, const PsCoarse& q, int ic, int jc, int kc, + int* i0, int* i1, int* j0, int* j1, int* k0, int* k1) +{ + *i0 = 2 * ic; *i1 = ic == q.nxc - 1 ? g.nx - 1 : 2 * ic + 1; + *j0 = 2 * jc; *j1 = jc == q.nyc - 1 ? g.ny - 1 : 2 * jc + 1; + *k0 = 2 * kc; *k1 = kc == q.nzc - 1 ? g.nz - 1 : 2 * kc + 1; +} + +/* The coarse operator: per coarse cell, its level-active fine children in ascending index + * order (the serial loop's order over `cells`), each coefficient summed when its face leaves + * the coarse cell; the Dirichlet parts summed; active when any child is. */ +extern "C" __global__ void e3_ps_coarsen( + PsGrid g, PsCoarse q, PsOp op, const u8* __restrict__ lact, const u32* __restrict__ top, + const u32* __restrict__ bot, PsOp co, u8* __restrict__ cact) +{ + long long c = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long nc = (long long) q.nxc * q.nyc * q.nzc; + if (c >= nc) return; + int ic = (int) (c % q.nxc); + int jc = (int) ((c / q.nxc) % q.nyc); + int kc = (int) (c / ((long long) q.nxc * q.nyc)); + int i0, i1, j0, j1, k0, k1; + ps_box(g, q, ic, jc, kc, &i0, &i1, &j0, &j1, &k0, &k1); + long long nxy = (long long) g.nx * g.ny; + double ae = 0.0, aw = 0.0, an = 0.0, as_ = 0.0, at = 0.0, ab = 0.0, extra = 0.0; + u8 act = 0; +#define CX(v) ((v) / 2 < q.nxc - 1 ? (v) / 2 : q.nxc - 1) +#define CY(v) ((v) / 2 < q.nyc - 1 ? (v) / 2 : q.nyc - 1) +#define CZ(v) ((v) / 2 < q.nzc - 1 ? (v) / 2 : q.nzc - 1) + for (int k = k0; k <= k1; ++k) { + for (int j = j0; j <= j1; ++j) { + for (int i = i0; i <= i1; ++i) { + long long idx = k * nxy + (long long) j * g.nx + i; + if (!lact[idx]) continue; + act = 1; + extra += op.extra[idx]; + double v; + v = op.ae[idx]; if (v != 0.0 && CX(i + 1) != ic) ae += v; + v = op.aw[idx]; if (v != 0.0 && CX(i - 1) != ic) aw += v; + v = op.an[idx]; if (v != 0.0 && CY(j + 1) != jc) an += v; + v = op.as_[idx]; if (v != 0.0 && CY(j - 1) != jc) as_ += v; + v = op.at[idx]; if (v != 0.0 && CZ((int) (top[idx] / nxy)) != kc) at += v; + v = op.ab[idx]; if (v != 0.0 && CZ((int) (bot[idx] / nxy)) != kc) ab += v; + } + } + } +#undef CX +#undef CY +#undef CZ + co.ae[c] = ae; co.aw[c] = aw; co.an[c] = an; co.as_[c] = as_; + co.at[c] = at; co.ab[c] = ab; co.extra[c] = extra; + cact[c] = act; +} + +/* Level 0's children lists for level 1's cells (in its `cells` order): the count per position. */ +extern "C" __global__ void e3_ps_children_count( + PsGrid g, PsCoarse q, int n1, const u32* __restrict__ cells1, const u8* __restrict__ lact, u32* __restrict__ counts) +{ + int p = blockIdx.x * blockDim.x + threadIdx.x; + if (p >= n1) return; + long long c = cells1[p]; + int ic = (int) (c % q.nxc); + int jc = (int) ((c / q.nxc) % q.nyc); + int kc = (int) (c / ((long long) q.nxc * q.nyc)); + int i0, i1, j0, j1, k0, k1; + ps_box(g, q, ic, jc, kc, &i0, &i1, &j0, &j1, &k0, &k1); + long long nxy = (long long) g.nx * g.ny; + u32 m = 0; + for (int k = k0; k <= k1; ++k) + for (int j = j0; j <= j1; ++j) + for (int i = i0; i <= i1; ++i) + if (lact[k * nxy + (long long) j * g.nx + i]) ++m; + counts[p] = m; +} + +/* ... and the fill (ascending fine index within each list). */ +extern "C" __global__ void e3_ps_children_fill( + PsGrid g, PsCoarse q, int n1, const u32* __restrict__ cells1, const u8* __restrict__ lact, + const u32* __restrict__ ptr, u32* __restrict__ out) +{ + int p = blockIdx.x * blockDim.x + threadIdx.x; + if (p >= n1) return; + long long c = cells1[p]; + int ic = (int) (c % q.nxc); + int jc = (int) ((c / q.nxc) % q.nyc); + int kc = (int) (c / ((long long) q.nxc * q.nyc)); + int i0, i1, j0, j1, k0, k1; + ps_box(g, q, ic, jc, kc, &i0, &i1, &j0, &j1, &k0, &k1); + long long nxy = (long long) g.nx * g.ny; + u32 w = ptr[p]; + for (int k = k0; k <= k1; ++k) + for (int j = j0; j <= j1; ++j) + for (int i = i0; i <= i1; ++i) { + long long idx = k * nxy + (long long) j * g.nx + i; + if (lact[idx]) out[w++] = (u32) idx; + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device.rs index 928d6a4..35acd3c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device.rs @@ -119,51 +119,57 @@ pub(crate) fn cfg(n_items: usize) -> LaunchConfig { } } +/// One exported level uploaded (its work vectors zeroed). +pub(crate) fn upload_level(l: &LevelExport) -> DevLevel { + let rt = runtime(); + let up_u = |v: &[u32]| -> CudaSlice { + rt.stream + .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v }) + .expect("upload") + }; + let up_f = |v: &[f32]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; + let n = l.nx * l.ny * l.nz; + DevLevel { + n, + nx: l.nx as i32, + n_cells: l.cells.len(), + n_red: l.red.len(), + n_black: l.black.len(), + cells: up_u(&l.cells), + red: up_u(&l.red), + black: up_u(&l.black), + top: up_u(&l.top), + bot: up_u(&l.bot), + coarse_of: up_u(&l.coarse_of), + children_ptr: up_u(&l.children_ptr), + children_idx: up_u(&l.children_idx), + ae: up_f(&l.ae), + aw: up_f(&l.aw), + an: up_f(&l.an), + as_: up_f(&l.as_), + at: up_f(&l.at), + ab: up_f(&l.ab), + ap: up_f(&l.ap), + b: rt.stream.alloc_zeros::(n).expect("alloc"), + x: rt.stream.alloc_zeros::(n).expect("alloc"), + r: rt.stream.alloc_zeros::(n).expect("alloc"), + } +} + impl DeviceVcycle { pub fn new(levels: &[LevelExport], sweeps: usize) -> Self { - let rt = runtime(); - let up_u = |v: &[u32]| -> CudaSlice { - rt.stream - .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v }) - .expect("upload") - }; - let up_f = |v: &[f32]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; - let dev: Vec = levels - .iter() - .map(|l| { - let n = l.nx * l.ny * l.nz; - DevLevel { - n, - nx: l.nx as i32, - n_cells: l.cells.len(), - n_red: l.red.len(), - n_black: l.black.len(), - cells: up_u(&l.cells), - red: up_u(&l.red), - black: up_u(&l.black), - top: up_u(&l.top), - bot: up_u(&l.bot), - coarse_of: up_u(&l.coarse_of), - children_ptr: up_u(&l.children_ptr), - children_idx: up_u(&l.children_idx), - ae: up_f(&l.ae), - aw: up_f(&l.aw), - an: up_f(&l.an), - as_: up_f(&l.as_), - at: up_f(&l.at), - ab: up_f(&l.ab), - ap: up_f(&l.ap), - b: rt.stream.alloc_zeros::(n).expect("alloc"), - x: rt.stream.alloc_zeros::(n).expect("alloc"), - r: rt.stream.alloc_zeros::(n).expect("alloc"), - } - }) - .collect(); + let dev: Vec = levels.iter().map(upload_level).collect(); + Self::from_levels(dev, sweeps, levels[0].cells.clone()) + } + + /// R6-3: from device levels (level 0 built on the device: its cell list + /// is read back on demand when `fine_cells` is empty and stale). + pub(crate) fn from_levels(dev: Vec, sweeps: usize, fine_cells: Vec) -> Self { let n0 = dev[0].n; Self { levels: dev, sweeps, - fine_cells: levels[0].cells.clone(), + fine_cells, fine_cells_stale: false, r_f32: vec![0.0; n0], z_f32: vec![0.0; n0], 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 c9705b8..ba767ed 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 @@ -217,7 +217,9 @@ impl DeviceCg { singular_count ); let levels = export_hierarchy(problem, params); + let l_export = lap.elapsed(); let vcycle = DeviceVcycle::new(&levels, params.smoother_sweeps.max(1)); + let l_vcycle = lap.elapsed(); let n = problem.nx * problem.ny * problem.nz; let n_cells = fine.cells.len(); let n_blocks = n_cells.div_ceil(256).max(1); @@ -245,8 +247,23 @@ impl DeviceCg { } link_ptr.push(link_idx.len() as u32); } + let l_links = lap.elapsed(); + let key = Some(OperatorKey::of(problem, params)); + let l_key = lap.elapsed(); + if profile { + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + eprintln!( + " new laps: level {:.0} ms, components {:.0} ms, export_hierarchy {:.0} ms, vcycle upload {:.0} ms, link lists {:.0} ms, key {:.0} ms", + ms(l_level), + ms(l_components - l_level), + ms(l_export - l_components), + ms(l_vcycle - l_export), + ms(l_links - l_vcycle), + ms(l_key - l_links) + ); + } Self { - key: Some(OperatorKey::of(problem, params)), + key, n, nx: problem.nx as i32, n_cells, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg/setup.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg/setup.rs index 9c0b91f..94d2798 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg/setup.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg/setup.rs @@ -9,9 +9,12 @@ //! V-cycle's persistent buffers; no host level, no host operator. //! `RTX_E3_BAND_CHECK=1`: every structure against the host build, bit for bit. +use super::super::device::{DevLevel, DeviceVcycle, upload_level}; use super::super::device::{cfg, load_module, runtime}; -use super::super::{Components, Level, Problem}; +use super::super::export::{LevelExport, export_hierarchy, export_levels}; +use super::super::{Components, Hierarchy, Level, Problem}; use super::DeviceCg; +use crate::solvers::incompressible::poisson::MultigridParameters; use cudarc::driver::{ CudaFunction, CudaModule, CudaSlice, DevicePtr, DeviceRepr, LaunchConfig, PushKernelArg, ValidAsZeroBits, @@ -41,6 +44,9 @@ struct PsKernels { top: CudaFunction, compact: CudaFunction, values: CudaFunction, + coarsen: CudaFunction, + children_count: CudaFunction, + children_fill: CudaFunction, } static PS: OnceLock = OnceLock::new(); @@ -67,6 +73,9 @@ fn kernels() -> &'static PsKernels { top: f("e3_ps_scan_top"), compact: f("e3_ps_scan_compact"), values: f("e3_ps_scan_values"), + coarsen: f("e3_ps_coarsen"), + children_count: f("e3_ps_children_count"), + children_fill: f("e3_ps_children_fill"), _module: module, } }) @@ -105,6 +114,17 @@ struct PsOp { unsafe impl DeviceRepr for PsOp {} unsafe impl ValidAsZeroBits for PsOp {} +#[repr(C)] +#[derive(Clone, Copy)] +struct PsCoarse { + nxc: i32, + nyc: i32, + nzc: i32, + pad: i32, +} +unsafe impl DeviceRepr for PsCoarse {} +unsafe impl ValidAsZeroBits for PsCoarse {} + #[repr(C)] #[derive(Clone, Copy)] struct PsIn { @@ -176,6 +196,24 @@ impl PsState { } } +fn ps_grid(grid: &OperatorGrid) -> PsGrid { + let (dx, dy, dz, dt) = (grid.dx, grid.dy, grid.dz, grid.dt); + let o = |b: bool| i32::from(b); + PsGrid { + nx: grid.nx as i32, + ny: grid.ny as i32, + nz: grid.nz as i32, + periodic: o(grid.periodic), + out: grid.outlet.map(o), + ae_int: dt * (dy * dz) / dx, + an_int: dt * (dx * dz) / dy, + at_int: dt * (dx * dy) / dz, + ae_out: dt * (dy * dz) / (0.5 * dx), + an_out: dt * (dx * dz) / (0.5 * dy), + at_out: dt * (dx * dy) / (0.5 * dz), + } +} + fn scan_cfg(n: usize) -> LaunchConfig { LaunchConfig { grid_dim: (n.div_ceil(SCAN_BLOCK).max(1) as u32, 1, 1), @@ -274,21 +312,7 @@ impl DeviceCg { let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz); let n = nx * ny * nz; assert_eq!(n, self.n, "refresh_device: grid size"); - let (dx, dy, dz, dt) = (grid.dx, grid.dy, grid.dz, grid.dt); - let o = |b: bool| i32::from(b); - let pg = PsGrid { - nx: nx as i32, - ny: ny as i32, - nz: nz as i32, - periodic: o(grid.periodic), - out: grid.outlet.map(o), - ae_int: dt * (dy * dz) / dx, - an_int: dt * (dx * dz) / dy, - at_int: dt * (dx * dy) / dz, - ae_out: dt * (dy * dz) / (0.5 * dx), - an_out: dt * (dx * dz) / (0.5 * dy), - at_out: dt * (dx * dy) / (0.5 * dz), - }; + let pg = ps_grid(grid); let mut ps = match self.pset.take() { Some(p) => p, None => Box::new(PsState::new(n)), @@ -768,3 +792,324 @@ impl DeviceCg { } } } + +impl DeviceCg { + /// R6-3: the hierarchy build with level 0 on the device — the fine + /// operator, level, lists, links, components and export by + /// [`Self::refresh_device`], level 1's operator coarsened on the device + /// (`Level::coarsen` per coarse cell) and brought down, the levels below + /// built on the host by `Hierarchy::build`'s loop, level 0's children + /// lists on the device. `None` when the hierarchy would stop at level 0 + /// (the caller builds on the host). `reference`: the host operator — + /// every structure against `DeviceCg::new`'s (`RTX_E3_BAND_CHECK=1`). + pub fn new_device( + grid: &OperatorGrid, + inp: &OperatorInputs<'_>, + params: &MultigridParameters, + reference: Option, + t: f64, + ) -> Option { + let rt = runtime(); + let s = &rt.stream; + let k = kernels(); + let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); + let lap = Instant::now(); + let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz); + let n = nx * ny * nz; + let zf = |m: usize| s.alloc_zeros::(m.max(1)).expect("alloc"); + let zf32 = |m: usize| s.alloc_zeros::(m.max(1)).expect("alloc"); + let zu = |m: usize| s.alloc_zeros::(m.max(1)).expect("alloc"); + let level0 = DevLevel { + n, + nx: nx as i32, + n_cells: 0, + n_red: 0, + n_black: 0, + cells: zu(n), + red: zu(n), + black: zu(n), + top: zu(n), + bot: zu(n), + coarse_of: zu(n), + children_ptr: zu(1), + children_idx: zu(1), + ae: zf32(n), + aw: zf32(n), + an: zf32(n), + as_: zf32(n), + at: zf32(n), + ab: zf32(n), + ap: zf32(n), + b: zf32(n), + x: zf32(n), + r: zf32(n), + }; + let sweeps = params.smoother_sweeps.max(1); + let mut cg = Self { + key: None, + n, + nx: nx as i32, + n_cells: 0, + n_blocks: 1, + cells: zu(n), + top: zu(n), + bot: zu(n), + ae: zf(n), + aw: zf(n), + an: zf(n), + as_: zf(n), + at: zf(n), + ab: zf(n), + ap: zf(n), + link_ptr: zu(n + 1), + link_idx: zu(1), + link_coef: zf(1), + b: zf(n), + r: zf(n), + z: zf(n), + d: zf(n), + q: zf(n), + partial: zf(1), + scalar: zf(1), + vcycle: DeviceVcycle::from_levels(vec![level0], sweeps, Vec::new()), + singular: false, + fine: None, + fine_export: None, + active_host: Vec::new(), + pset: None, + max_iterations: params.max_iterations, + scalar_host: vec![0.0], + }; + // The fine structures are checked there (`reference`), the hierarchy below. + cg.refresh_device(grid, inp, reference.clone(), t); + let l_fine = lap.elapsed(); + if cg.n_cells <= params.coarsest_cells.max(1) { + return None; + } + // Level 1's operator on the device, then down. + let pg = ps_grid(grid); + let (nxc, nyc, nzc) = ((nx / 2).max(1), (ny / 2).max(1), (nz / 2).max(1)); + let q = PsCoarse { + nxc: nxc as i32, + nyc: nyc as i32, + nzc: nzc as i32, + pad: 0, + }; + let nc = nxc * nyc * nzc; + let ps = cg.pset.as_ref().expect("device setup"); + let pf = |x: &CudaSlice| x.device_ptr(s).0; + let op = PsOp { + ptrs: [ + pf(&cg.ae), + pf(&cg.aw), + pf(&cg.an), + pf(&cg.as_), + pf(&cg.at), + pf(&cg.ab), + pf(&ps.extra), + ], + }; + let c_arr: Vec> = (0..7).map(|_| zf(nc)).collect(); + let co = PsOp { + ptrs: std::array::from_fn(|i| pf(&c_arr[i])), + }; + let mut cact = s.alloc_zeros::(nc).expect("alloc"); + unsafe { + s.launch_builder(&k.coarsen) + .arg(&pg) + .arg(&q) + .arg(&op) + .arg(&ps.lact) + .arg(&cg.top) + .arg(&cg.bot) + .arg(&co) + .arg(&mut cact) + .launch(cfg(nc)) + .expect("e3_ps_coarsen"); + } + let mut coarse = Problem::new(nxc, nyc, nzc); + coarse.periodic_z = grid.periodic; + coarse.active = s + .memcpy_dtov(&cact) + .expect("cact") + .into_iter() + .map(|v| v != 0) + .collect(); + let down = |x: &CudaSlice| s.memcpy_dtov(&x.slice(0..nc)).expect("coarse"); + coarse.ae = down(&c_arr[0]); + coarse.aw = down(&c_arr[1]); + coarse.an = down(&c_arr[2]); + coarse.as_ = down(&c_arr[3]); + coarse.at = down(&c_arr[4]); + coarse.ab = down(&c_arr[5]); + coarse.extra_diag = down(&c_arr[6]); + drop(c_arr); + let l_coarsen = lap.elapsed(); + // Levels 1.. on the host (Hierarchy::build's loop). + let first = Level::::new(coarse); + if first.cells.len() >= cg.n_cells { + return None; + } + let levels = Hierarchy::::chain_below(first, params); + let exports = export_levels(&levels); + let l_host = lap.elapsed(); + // Level 0's children lists for level 1's cells, on the device. + let n1 = exports[0].cells.len(); + let d_cells1 = s.memcpy_stod(&exports[0].cells).expect("cells1"); + let n1i = n1 as i32; + let mut counts = zu(n1); + let mut sums = zu(n1.div_ceil(SCAN_BLOCK) + 1); + let mut ptr = zu(n1 + 1); + let n1_64 = n1 as i64; + let nb32 = n1.div_ceil(SCAN_BLOCK).max(1) as i32; + unsafe { + s.launch_builder(&k.children_count) + .arg(&pg) + .arg(&q) + .arg(&n1i) + .arg(&d_cells1) + .arg(&ps.lact) + .arg(&mut counts) + .launch(cfg(n1)) + .expect("e3_ps_children_count"); + s.launch_builder(&k.count_vals) + .arg(&n1_64) + .arg(&counts) + .arg(&mut sums) + .launch(scan_cfg(n1)) + .expect("e3_ps_scan_count_vals"); + s.launch_builder(&k.top) + .arg(&nb32) + .arg(&mut sums) + .launch(one_block()) + .expect("e3_ps_scan_top"); + s.launch_builder(&k.values) + .arg(&n1_64) + .arg(&counts) + .arg(&sums) + .arg(&mut ptr) + .launch(scan_cfg(n1)) + .expect("e3_ps_scan_values"); + } + let nb = n1.div_ceil(SCAN_BLOCK).max(1); + let total = s.memcpy_dtov(&sums.slice(nb..nb + 1)).expect("total")[0]; + s.memcpy_htod(&[total][..], &mut ptr.slice_mut(n1..n1 + 1)) + .expect("children_ptr end"); + let mut children = zu(total as usize); + unsafe { + s.launch_builder(&k.children_fill) + .arg(&pg) + .arg(&q) + .arg(&n1i) + .arg(&d_cells1) + .arg(&ps.lact) + .arg(&ptr) + .arg(&mut children) + .launch(cfg(n1)) + .expect("e3_ps_children_fill"); + } + { + let l0 = &mut cg.vcycle.levels[0]; + l0.children_ptr = ptr; + l0.children_idx = children; + } + cg.vcycle.levels.extend(exports.iter().map(upload_level)); + s.synchronize().expect("sync"); + let l_up = lap.elapsed(); + if profile { + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + eprintln!( + " new_device laps: total {:.0} ms (fine {:.0}, coarsen + download {:.0}, host levels {:.0}, children + upload {:.0}); depth {}", + ms(l_up), + ms(l_fine), + ms(l_coarsen - l_fine), + ms(l_host - l_coarsen), + ms(l_up - l_host), + cg.vcycle.levels.len() + ); + } + if let Some(problem) = reference { + cg.check_hierarchy(problem, params, &exports, t); + } + Some(cg) + } + + /// `RTX_E3_BAND_CHECK=1` for [`Self::new_device`]: the fine structures as + /// [`Self::refresh_device`]'s check, `DeviceCg::new`'s components (the + /// raw operator's), and every level of `export_hierarchy` — level 0 read + /// back from the device, the host-built levels below compared as built. + fn check_hierarchy( + &self, + problem: Problem, + params: &MultigridParameters, + coarse: &[LevelExport], + t: f64, + ) { + let s = &runtime().stream; + let reference = export_hierarchy(&problem, params); + let fine_ref = Level::::new(problem.clone()); + let comps = Components::find_planes(&problem, &fine_ref.cells); + let raw_singular = comps.singular.iter().filter(|&&f| f).count(); + drop(fine_ref); + let mut differs: Vec = Vec::new(); + if (raw_singular > 0) != self.singular { + differs.push("singular (raw components)".into()); + } + if reference.len() != coarse.len() + 1 { + differs.push(format!("depth {} vs {}", coarse.len() + 1, reference.len())); + } + let l0 = &self.vcycle.levels[0]; + let r0 = &reference[0]; + let down_u = |x: &CudaSlice, m: usize| -> Vec { + if m == 0 { + Vec::new() + } else { + s.memcpy_dtov(&x.slice(0..m)).expect("down") + } + }; + if down_u(&l0.children_ptr, r0.children_ptr.len()) != r0.children_ptr { + differs.push("level 0 children_ptr".into()); + } + if down_u(&l0.children_idx, r0.children_idx.len()) != r0.children_idx { + differs.push("level 0 children_idx".into()); + } + let bits = |a: &[f32], b: &[f32]| { + a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()) + }; + for (l, (mine, want)) in coarse.iter().zip(reference.iter().skip(1)).enumerate() { + let same = mine.nx == want.nx + && mine.ny == want.ny + && mine.nz == want.nz + && mine.cells == want.cells + && mine.red == want.red + && mine.black == want.black + && mine.top == want.top + && mine.bot == want.bot + && mine.coarse_of == want.coarse_of + && mine.children_ptr == want.children_ptr + && mine.children_idx == want.children_idx + && bits(&mine.ae, &want.ae) + && bits(&mine.aw, &want.aw) + && bits(&mine.an, &want.an) + && bits(&mine.as_, &want.as_) + && bits(&mine.at, &want.at) + && bits(&mine.ab, &want.ab) + && bits(&mine.ap, &want.ap); + if !same { + differs.push(format!("level {}", l + 1)); + } + } + if differs.is_empty() { + eprintln!( + " hierarchy check t {t:.6}: device level 0 + host levels 1..{} IDENTICAL to export_hierarchy (children {} entries)", + coarse.len(), + r0.children_idx.len() + ); + } else { + eprintln!(" hierarchy check t {t:.6}: DIFFERS in {differs:?}"); + if std::env::var("RTX_E3_GEOM_CHECK_SOFT").is_err() { + panic!("R6-3 hierarchy check: {differs:?}"); + } + } + } +} 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 cd22114..e9d57ed 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 @@ -33,7 +33,14 @@ pub struct LevelExport { /// The f32 hierarchy of `problem`, level 0 fine. pub fn export_hierarchy(problem: &Problem, params: &MultigridParameters) -> Vec { let hier = Hierarchy::::build(problem, params); - let depth = hier.levels.len(); + export_levels(&hier.levels) +} + +/// The exports of a chain of levels (each level's children lists for the +/// next one; R6-3 exports the host-built levels below a device level 0). +#[cfg_attr(not(feature = "cuda"), allow(dead_code))] +pub(crate) fn export_levels(levels: &[Level]) -> Vec { + let depth = levels.len(); let to_u32 = |v: &[usize]| { v.iter() .map(|&i| if i == usize::MAX { u32::MAX } else { i as u32 }) @@ -41,9 +48,9 @@ pub fn export_hierarchy(problem: &Problem, params: &MultigridParameters) -> Vec< }; (0..depth) .map(|l| { - let lv = &hier.levels[l]; + let lv = &levels[l]; let (children_ptr, children_idx) = if l + 1 < depth { - let coarse = &hier.levels[l + 1]; + let coarse = &levels[l + 1]; let nc = coarse.problem.nx * coarse.problem.ny * coarse.problem.nz; let mut pos = vec![usize::MAX; nc]; for (k, &c) in coarse.cells.iter().enumerate() { 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 3dfaa6d..3ee9bdb 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 @@ -441,6 +441,29 @@ impl Hierarchy { self.levels.len() } + /// R6-3: the levels below a level 0 that lives elsewhere — `first` is + /// level 1 (already accepted by [`Self::build`]'s rule), coarsened by the + /// same loop with level 0 counted against `MAX_LEVELS`. + #[cfg_attr(not(feature = "cuda"), allow(dead_code))] + pub(crate) fn chain_below(first: Level, params: &MultigridParameters) -> Vec> { + let mut levels = vec![first]; + while levels.len() + 1 < MAX_LEVELS { + let fine = levels.last().expect("at least one level"); + if fine.cells.len() <= params.coarsest_cells.max(1) { + break; + } + let (coarse, coarse_of) = fine.coarsen(); + let coarse = Level::::new(coarse); + if coarse.cells.len() >= fine.cells.len() { + break; + } + let last = levels.len() - 1; + levels[last].coarse_of = coarse_of; + levels.push(coarse); + } + levels + } + /// `z = M⁻¹ r`: one V-cycle from zero (the 2D sequence). pub fn apply_preconditioner(&mut self, r: &[f64], z: &mut [f64]) { let depth = self.levels.len(); 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 0458b17..63685b9 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 @@ -1139,13 +1139,30 @@ impl DeviceStep { } let cptrs = self.cut.as_ref().expect("cut").ptrs(); if self.cg.is_none() || self.cg_dt != dt { - let problem = self.solver.poisson_operator(g, dt); let params = MultigridParameters { precision: self.solver.params.poisson_precision, smoother: self.solver.params.poisson_smoother, ..MultigridParameters::default() }; - self.cg = Some(DeviceCg::new(&problem, ¶ms)); + // R6-3 (`RTX_E3_POISSON_DEVICE=1`): a moving body's hierarchy with + // level 0 built on the device; the host build otherwise (or when + // the hierarchy would stop at level 0). + let device = match (self.solver.is_moving(), self.cut.as_ref()) { + (true, Some(dc)) => super::poisson_setup::device_grid(&self.solver, g, dt) + .and_then(|grid| { + let check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); + let reference = check.then(|| self.solver.poisson_operator(g, dt)); + DeviceCg::new_device(&grid, &dc.poisson_inputs(), ¶ms, reference, self.solver.time() + dt) + }), + _ => None, + }; + self.cg = Some(match device { + Some(cg) => cg, + None => { + let problem = self.solver.poisson_operator(g, dt); + DeviceCg::new(&problem, ¶ms) + } + }); self.cg_dt = dt; } let rebuild = t_rebuild.elapsed();