diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs
index 7286e06..125dc4a 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs
@@ -69,7 +69,7 @@ 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,
+    PoissonSolution, export_hierarchy, vcycle_f32_reference, vcycle_f32_work,
     PoissonSolverKind, configure_threads, solve_multigrid_pcg, solve_multigrid_pcg_cached,
 };
 pub use polygon_sdf::PolygonSdf;
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs
index bb5a31d..2b111fd 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs
@@ -1094,6 +1094,24 @@ pub fn vcycle_f32_reference(
     hier.apply_preconditioner(r, z);
 }
 
+/// The f32 V-cycle's per-level work vectors after one application on `r`
+/// (for locating where a device V-cycle departs from the CPU one): per
+/// level `(b, x, r)` — the restricted right-hand side, the level's final
+/// correction (after its up-leg smoothing), and its down-leg residual.
+pub fn vcycle_f32_work(
+    problem: &PoissonProblem,
+    params: &MultigridParameters,
+    r: &[f64],
+) -> Vec<(Vec<f32>, Vec<f32>, Vec<f32>)> {
+    let mut hier = Hierarchy::<f32>::build(problem, params);
+    let mut z = vec![0.0; r.len()];
+    hier.apply_preconditioner(r, &mut z);
+    hier.work
+        .iter()
+        .map(|w| (w.b.clone(), w.x.clone(), w.r.clone()))
+        .collect()
+}
+
 /// A reusable prepared solver per V-cycle precision (PERF-2 P1.1): the
 /// operator's hierarchy is rebuilt only when the operator changes.
 #[derive(Default)]
diff --git a/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs b/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs
index 819c962..f17dc1c 100644
--- a/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs
+++ b/crates/specialized/rtx-cfd/tests/gpu_vcycle_bench.rs
@@ -15,7 +15,7 @@ use cudarc::driver::{CudaContext, CudaSlice, LaunchConfig, PushKernelArg};
 use cudarc::nvrtc::{CompileOptions, compile_ptx_with_opts};
 use rtx_cfd::solvers::incompressible::{
     LevelExport, MgSmoother, MultigridParameters, PoissonProblem, export_hierarchy,
-    vcycle_f32_reference,
+    vcycle_f32_reference, vcycle_f32_work,
 };
 use std::sync::Arc;
 use std::time::Instant;
@@ -342,6 +342,33 @@ fn batched_device_vcycle_go_no_go() {
 
         // Correctness at K = 1 (and every batch member at K > 1).
         vcycle(&mut dev, &r_f32, &mut z_out);
+        if k == 1 {
+            // Level by level against the CPU's work vectors (march 0).
+            let work = vcycle_f32_work(&prob, &params, &r_host);
+            for (l, (b_ref, x_ref, r_ref)) in work.iter().enumerate() {
+                let n = dev[l].n;
+                let mut b = vec![0.0f32; n * k];
+                let mut x = vec![0.0f32; n * k];
+                let mut r = vec![0.0f32; n * k];
+                stream.memcpy_dtoh(&dev[l].b, &mut b).unwrap();
+                stream.memcpy_dtoh(&dev[l].x, &mut x).unwrap();
+                stream.memcpy_dtoh(&dev[l].r, &mut r).unwrap();
+                stream.synchronize().unwrap();
+                let cmp = |a: &[f32], c: &[f32]| {
+                    let (mut worst, mut scale) = (0.0f32, 0.0f32);
+                    for &idx in &levels[l].cells {
+                        let i = idx as usize;
+                        worst = worst.max((a[i] - c[i]).abs());
+                        scale = scale.max(c[i].abs());
+                    }
+                    (worst, scale)
+                };
+                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());
+            }
+        }
         let mut worst = 0.0_f64;
         for m in 0..k {
             for idx in 0..n0 {