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) <[email protected]>
2026-09-25 04:09:51 -05:00
parent b00f297ef7
commit da2945fcf5
7 changed files with 581 additions and 61 deletions
@@ -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;
+            }
+}
@@ -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<u32> {
+        rt.stream
+            .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v })
+            .expect("upload")
+    };
+    let up_f = |v: &[f32]| -> CudaSlice<f32> { 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::<f32>(n).expect("alloc"),
+        x: rt.stream.alloc_zeros::<f32>(n).expect("alloc"),
+        r: rt.stream.alloc_zeros::<f32>(n).expect("alloc"),
+    }
+}
+
 impl DeviceVcycle {
    pub fn new(levels: &[LevelExport], sweeps: usize) -> Self {
-        let rt = runtime();
-        let up_u = |v: &[u32]| -> CudaSlice<u32> {
-            rt.stream
-                .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v })
-                .expect("upload")
-        };
-        let up_f = |v: &[f32]| -> CudaSlice<f32> { rt.stream.memcpy_stod(v).expect("upload") };
-        let dev: Vec<DevLevel> = 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::<f32>(n).expect("alloc"),
-                    x: rt.stream.alloc_zeros::<f32>(n).expect("alloc"),
-                    r: rt.stream.alloc_zeros::<f32>(n).expect("alloc"),
-                }
-            })
-            .collect();
+        let dev: Vec<DevLevel> = 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<DevLevel>, sweeps: usize, fine_cells: Vec<u32>) -> 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],
@@ -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,
@@ -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<PsKernels> = 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<Problem>,
+        t: f64,
+    ) -> Option<Self> {
+        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::<f64>(m.max(1)).expect("alloc");
+        let zf32 = |m: usize| s.alloc_zeros::<f32>(m.max(1)).expect("alloc");
+        let zu = |m: usize| s.alloc_zeros::<u32>(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<f64>| 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<CudaSlice<f64>> = (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::<u8>(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<f64>| 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::<f32>::new(coarse);
+        if first.cells.len() >= cg.n_cells {
+            return None;
+        }
+        let levels = Hierarchy::<f32>::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::<f64>::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<String> = 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<u32>, m: usize| -> Vec<u32> {
+            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:?}");
+            }
+        }
+    }
+}
@@ -33,7 +33,14 @@ pub struct LevelExport {
 /// The f32 hierarchy of `problem`, level 0 fine.
 pub fn export_hierarchy(problem: &Problem, params: &MultigridParameters) -> Vec<LevelExport> {
    let hier = Hierarchy::<f32>::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<f32>]) -> Vec<LevelExport> {
+    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() {
@@ -441,6 +441,29 @@ impl<T: MgScalar> Hierarchy<T> {
        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<T>, params: &MultigridParameters) -> Vec<Level<T>> {
+        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::<T>::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();
@@ -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, &params));
+            // 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(), &params, 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, &params)
+                }
+            });
            self.cg_dt = dt;
        }
        let rebuild = t_rebuild.elapsed();