diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_loads.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_loads.cu
new file mode 100644
index 0000000..9eea1d5
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_loads.cu
@@ -0,0 +1,237 @@
+/**
+ * R8-f: the operator-route wall load on the device — `Mask::cut_wall_force`
+ * (the cells' p W and the fluid faces' implicit wall shear, cutwall.rs) and
+ * `Mask::cut_wall_exchange_parts` (the diffusive and convective exchange of
+ * the fluid faces with their prescribed neighbours, exchange.rs), summand
+ * for summand in the host's expression order (FMA contraction off, so each
+ * summand is the host's to the bit). The kernels do not sum: every summand
+ * is appended as a record (key, index, value) to one buffer; the host sorts
+ * the records into its loop order and sums them sequentially, so the totals
+ * are the host's to the bit as well.
+ *
+ * Appended after `e3_step.cu` and `e3_cut.cu` at load (shares E3Params,
+ * E3Ptrs, E3Cut and `cut_face` / `cut_cv` / `face_corr3` / `upwind3`).
+ * The E3Cut pointers are the PREDICTOR set: the instantaneous apertures, the
+ * fluid-kind open flags and the fluid cells of the mask at the solver's time.
+ */
+
+/* key = kind << 28 | c << 24 | slot; kind 0 the cell's p W (v[3]), 1 the
+   face's wall shear, 2 an exchange summand (slot = 4 d + 2 side + sub; side
+   0 plus, 1 minus; sub 0 convective, 1 diffusive). */
+struct LoadRec {
+    unsigned int key;
+    unsigned int idx;
+    double v[3];
+};
+
+struct LoadArgs {
+    double mu, rho;
+    int k0, k1;          /* the z planes */
+    int shifts;          /* the centroid-shift tables are the mask's (face_shifts Some) */
+    int centroid;        /* diffusion_centroid */
+    int axis;            /* wall_exchange_axis */
+    int conv_off;        /* exchange_convection_off */
+    int order;           /* the shear closure's order (wall_order) */
+    unsigned int cap;    /* the record buffer's capacity */
+};
+
+__device__ __forceinline__ void load_push(LoadRec* out, unsigned int* count, unsigned int cap,
+                                          unsigned int key, unsigned int idx, double v0, double v1, double v2)
+{
+    unsigned int at = atomicAdd(count, 1u);
+    if (at < cap) {
+        out[at].key = key;
+        out[at].idx = idx;
+        out[at].v[0] = v0;
+        out[at].v[1] = v1;
+        out[at].v[2] = v2;
+    }
+}
+
+/* The cells' p W over the planes k0..k1 (thread per cell of those planes):
+   W = −Σ A n from the apertures (cut.rs, `[-sx, -sy, -sz]`); a fluid cell
+   with W = 0 adds exact zeros to every sum and is not recorded. */
+extern "C" __global__ void e3_loads_cells(E3Params g, E3Ptrs f, E3Cut m, LoadArgs a,
+                                          LoadRec* __restrict__ out, unsigned int* __restrict__ count)
+{
+    long long t = (long long)blockIdx.x * blockDim.x + threadIdx.x;
+    long long plane = (long long)g.nx * g.ny;
+    if (t >= plane * (a.k1 - a.k0)) return;
+    long long idx = (long long)a.k0 * plane + t;
+    if (!m.active[idx]) return;
+    int i = (int)(idx % g.nx), j = (int)((idx / g.nx) % g.ny), k = (int)(idx / plane);
+    double ax = g.dy * g.dz, ay = g.dx * g.dz, az = g.dx * g.dy;
+    double sx = (m.a_u[uf3(g, k, j, i + 1)] - m.a_u[uf3(g, k, j, i)]) * ax;
+    double sy = (m.a_v[vf3(g, k, j + 1, i)] - m.a_v[vf3(g, k, j, i)]) * ay;
+    double sz = (m.a_w[wf3(g, k + 1, j, i)] - m.a_w[wf3(g, k, j, i)]) * az;
+    double w0 = -sx, w1 = -sy, w2 = -sz;
+    if (w0 == 0.0 && w1 == 0.0 && w2 == 0.0) return;
+    double p = f.p[idx];
+    load_push(out, count, a.cap, 0u, (unsigned int)idx, p * w0, p * w1, p * w2);
+}
+
+__device__ __forceinline__ const double* load_val(const E3Ptrs& f, int c) { return c == 0 ? f.u : (c == 1 ? f.v : f.w); }
+__device__ __forceinline__ const int* load_open(const E3Cut& m, int c) { return c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w); }
+
+/* The fluid faces of component c (thread per face of the component's
+   lattice): the implicit wall shear (cutwall.rs) and the exchange with the
+   prescribed neighbours (exchange.rs). */
+extern "C" __global__ void e3_loads_faces(E3Params g, E3Ptrs f, E3Cut m, LoadArgs a, int c,
+                                          LoadRec* __restrict__ out, unsigned int* __restrict__ count)
+{
+    long long t = (long long)blockIdx.x * blockDim.x + threadIdx.x;
+    int ni = c == 0 ? g.nx + 1 : g.nx;
+    int nj = c == 1 ? g.ny + 1 : g.ny;
+    int nk = c == 2 ? g.nz + 1 : g.nz;
+    if (t >= (long long)ni * nj * nk) return;
+    int i = (int)(t % ni), j = (int)((t / ni) % nj), k = (int)(t / ((long long)ni * nj));
+    /* the host's loop ranges */
+    if (c == 0 && (i < 1 || i >= g.nx)) return;
+    if (c == 1 && (j < 1 || j >= g.ny)) return;
+    int ks = a.k0, ke = a.k1;
+    if (c == 2) {
+        int ws = g.periodic_z ? 0 : 1;
+        if (ws > ks) ks = ws;
+        if (g.nz < ke) ke = g.nz;
+    }
+    if (k < ks || k >= ke) return;
+    int fidx = (int)t;
+    const int* open = load_open(m, c);
+    if (!open[fidx]) return;
+    const double* vals = load_val(f, c);
+    double h[3] = { g.dx, g.dy, g.dz };
+    double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy };
+    double alpha, apm[3], app[3], wall[3], distance;
+    cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &distance);
+    double u0 = vals[fidx];
+    unsigned int kc = (unsigned int)c << 24;
+
+    /* The implicit wall shear: mu A_w (c1 (u − U_b) + c2 (u_n − U_b)). */
+    double a_w = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
+    if (a_w != 0.0) {
+        const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
+        double ub = ubt[fidx];
+        /* advancing_factor(0) / d = (1 + 0.5 * 0) / d */
+        double c1 = (1.0 + 0.5 * 0.0) / distance, c2 = 0.0;
+        int nb = -1;
+        if (a.order >= 2) {
+            double n[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w };
+            int d = 0;
+            for (int kk = 1; kk < 3; ++kk) if (fabs(n[kk]) > fabs(n[d])) d = kk;
+            int q[3] = { i, j, k };
+            q[d] -= n[d] > 0.0 ? 1 : -1;
+            int fn = cut_face(g, c, q[0], q[1], q[2]);
+            if (fn >= 0 && cut_ap(m, c)[fn] > 0.0) {
+                double d1 = distance;
+                double d2 = d1 + h[d] * fabs(n[d]);
+                if (!(d2 <= d1)) {
+                    c1 = d2 / (d1 * (d2 - d1));
+                    c2 = -d1 / (d2 * (d2 - d1));
+                    nb = fn;
+                }
+            }
+        }
+        double un = nb >= 0 ? vals[nb] : ub;
+        double v = a.mu * a_w * (c1 * (u0 - ub) + c2 * (un - ub));
+        load_push(out, count, a.cap, (1u << 28) | kc, (unsigned int)fidx,
+                  c == 0 ? v : 0.0, c == 1 ? v : 0.0, c == 2 ? v : 0.0);
+    }
+
+    /* The exchange with prescribed neighbours. */
+    const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w);
+    double shift0[3] = { 0.0, 0.0, 0.0 };
+    if (a.shifts) { shift0[0] = sh[3 * fidx]; shift0[1] = sh[3 * fidx + 1]; shift0[2] = sh[3 * fidx + 2]; }
+    int cm[3] = { i, j, k }; cm[c] -= 1;   /* cell minus */
+    int cp[3] = { i, j, k };               /* cell plus */
+    const double* apc = cut_ap(m, c);
+    for (int d = 0; d < 3; ++d) {
+        double a_d = area[d];
+        int q[3];
+        q[0] = i; q[1] = j; q[2] = k; q[d] += 1; int f_up1 = cut_face(g, c, q[0], q[1], q[2]);
+        q[0] = i; q[1] = j; q[2] = k; q[d] -= 1; int f_dn1 = cut_face(g, c, q[0], q[1], q[2]);
+        int presc_up = f_up1 >= 0 && !open[f_up1];
+        int presc_dn = f_dn1 >= 0 && !open[f_dn1];
+        if (!presc_up && !presc_dn) continue;
+        q[0] = i; q[1] = j; q[2] = k; q[d] += 2; int f_up2 = cut_face(g, c, q[0], q[1], q[2]);
+        q[0] = i; q[1] = j; q[2] = k; q[d] -= 2; int f_dn2 = cut_face(g, c, q[0], q[1], q[2]);
+        double m_plus, m_minus;
+        if (d == c) {
+            double f_up = (f_up1 >= 0 ? apc[f_up1] : alpha) * (f_up1 >= 0 ? vals[f_up1] : u0);
+            double f_dn = (f_dn1 >= 0 ? apc[f_dn1] : alpha) * (f_dn1 >= 0 ? vals[f_dn1] : u0);
+            double f0 = alpha * u0;
+            m_plus = 0.5 * (f0 + f_up) * a_d;
+            m_minus = 0.5 * (f_dn + f0) * a_d;
+        } else {
+            const double* vd = load_val(f, d);
+            const double* apd = cut_ap(m, d);
+            int q1[3] = { cm[0], cm[1], cm[2] }; q1[d] += 1;
+            int q2[3] = { cp[0], cp[1], cp[2] }; q2[d] += 1;
+            int fa = cut_face(g, d, q1[0], q1[1], q1[2]);
+            int fb = cut_face(g, d, q2[0], q2[1], q2[2]);
+            int fc = cut_face(g, d, cm[0], cm[1], cm[2]);
+            int fd = cut_face(g, d, cp[0], cp[1], cp[2]);
+            double fla = (fa >= 0 ? apd[fa] : 1.0) * (fa >= 0 ? vd[fa] : 0.0);
+            double flb = (fb >= 0 ? apd[fb] : 1.0) * (fb >= 0 ? vd[fb] : 0.0);
+            double flc = (fc >= 0 ? apd[fc] : 1.0) * (fc >= 0 ? vd[fc] : 0.0);
+            double fld = (fd >= 0 ? apd[fd] : 1.0) * (fd >= 0 ? vd[fd] : 0.0);
+            m_plus = 0.5 * (fla + flb) * a_d;
+            m_minus = 0.5 * (flc + fld) * a_d;
+        }
+        /* exchange.rs `solid_spacing`: the centroid spacing in a cross
+           direction without the axis exchange, else `exchange_delta`. */
+        double delta_x = h[d];
+        if (a.axis) {
+            double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
+            if (a_w0 != 0.0) {
+                double n_d = fabs(wall[d]) / a_w0;
+                if (!(n_d < 1e-12)) delta_x = fmin(distance / n_d, h[d]);
+            }
+        }
+        int centroid_spacing = a.centroid && d != c && !a.axis;
+        unsigned int kd = (2u << 28) | kc | (unsigned int)(4 * d);
+        if (presc_up) {
+            double un = vals[f_up1];
+            double up1 = un;
+            double dn1 = f_dn1 >= 0 ? vals[f_dn1] : 0.0;
+            double up2 = f_up2 >= 0 ? vals[f_up2] : 0.0;
+            (void)up1;
+            double delta;
+            if (g.scheme == SCHEME_UPWIND) delta = 0.0;
+            else if (m_plus >= 0.0) delta = face_corr3(g.scheme, f_dn1 >= 0, dn1, u0, un);
+            else delta = face_corr3(g.scheme, f_up2 >= 0, up2, un, u0);
+            double u_face = upwind3(m_plus, u0, un) + delta;
+            if (!a.conv_off) {
+                double v = -a.rho * m_plus * (u_face - u0);
+                double r = -v;
+                load_push(out, count, a.cap, kd | 0u, (unsigned int)fidx,
+                          c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
+            }
+            double s = centroid_spacing ? fmin(fmax(h[d] - 1.0 * shift0[d], 0.25 * h[d]), 2.0 * h[d]) : delta_x;
+            double v = a.mu * app[d] * a_d * (un - u0) / s;
+            double r = -v;
+            load_push(out, count, a.cap, kd | 1u, (unsigned int)fidx,
+                      c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
+        }
+        if (presc_dn) {
+            double ud = vals[f_dn1];
+            double dn2 = f_dn2 >= 0 ? vals[f_dn2] : 0.0;
+            double up1 = f_up1 >= 0 ? vals[f_up1] : 0.0;
+            double delta;
+            if (g.scheme == SCHEME_UPWIND) delta = 0.0;
+            else if (m_minus >= 0.0) delta = face_corr3(g.scheme, f_dn2 >= 0, dn2, ud, u0);
+            else delta = face_corr3(g.scheme, f_up1 >= 0, up1, u0, ud);
+            double u_face = upwind3(m_minus, ud, u0) + delta;
+            if (!a.conv_off) {
+                double v = a.rho * m_minus * (u_face - u0);
+                double r = -v;
+                load_push(out, count, a.cap, kd | 2u, (unsigned int)fidx,
+                          c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
+            }
+            double s = centroid_spacing ? fmin(fmax(h[d] - (-1.0) * shift0[d], 0.25 * h[d]), 2.0 * h[d]) : delta_x;
+            double v = a.mu * apm[d] * a_d * (ud - u0) / s;
+            double r = -v;
+            load_push(out, count, a.cap, kd | 3u, (unsigned int)fidx,
+                      c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
+        }
+    }
+}
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
index 999412c..7a0828d 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
@@ -7,10 +7,12 @@
 
 mod cut;
 mod geom;
+mod loads;
 mod mask;
 mod poisson_setup;
 mod snapshot;
 
+pub use loads::{loads_device_enabled, DeviceWallLoads, WallLoadRecord, WallLoadsCheck};
 pub use snapshot::DeviceSnapshot;
 
 use super::{Side, Solver, StepResult};
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 c72cc6a..92dcf54 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
@@ -55,7 +55,7 @@ fn cut_kernels() -> &'static CutKernels {
 /// `struct E3Cut` in e3_cut.cu: 27 device pointers.
 #[repr(C)]
 #[derive(Clone, Copy)]
-struct E3CutPtrs {
+pub(super) struct E3CutPtrs {
     ptrs: [u64; 27],
 }
 unsafe impl DeviceRepr for E3CutPtrs {}
@@ -751,7 +751,7 @@ impl DeviceCut {
         }
     }
 
-    fn ptrs(&self) -> E3CutPtrs {
+    pub(super) fn ptrs(&self) -> E3CutPtrs {
         let rt = runtime();
         let s = &rt.stream;
         let pf = |x: &CudaSlice<f64>| x.device_ptr(s).0;
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/loads.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/loads.rs
new file mode 100644
index 0000000..c3372c9
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/loads.rs
@@ -0,0 +1,627 @@
+//! R8-f: the operator-route wall load on the device (`RTX_E3_LOADS_DEVICE=1`,
+//! default off). The host route — `Mask::cut_wall_force` (the cells' `p W`
+//! and the fluid faces' implicit wall shear, `cutwall.rs`) plus
+//! `Mask::cut_wall_exchange_force` (`exchange.rs`) — loops over every cell
+//! and face of the grid on the host, after the whole field comes down; in the
+//! 3D coupled step (1.45 M cells) that was 75 % of the step.
+//!
+//! Here the kernels (`e3_loads.cu`) evaluate every summand on the device from
+//! the predictor set of the device mask (the instantaneous apertures, the
+//! fluid-kind open flags, the distances, the surface velocities and the
+//! centroid shifts of the mask at the solver's time — the tables the host
+//! route reads) and the device fields, in the host's expression order with
+//! FMA contraction off, so each summand is the host's to the bit. They do
+//! not sum: each summand is appended as a record to one buffer, which comes
+//! down once per call; the host sorts the records into its own loop order
+//! (cells ascending; shear by component then face; exchange by component,
+//! face, direction, side, convective before diffusive) and sums them
+//! sequentially — the totals are therefore the host route's TO THE BIT, not
+//! merely within a tolerance (a device reduction would reorder the sums).
+//!
+//! What is not recorded: summands that are exactly zero because their cell
+//! has no wall (`W = 0`: the host adds `p · 0 = ±0` for every fluid cell).
+//! Adding `±0` to a sequential sum that starts at `+0` never changes its bits
+//! (a round-to-nearest sum of non-zero terms is never `−0`, and `x ± 0 = x`
+//! otherwise), so the totals and every consumer's sums are unchanged; the
+//! record lists replayed into the sinks drop exact zeros too (R8-c's own
+//! `cut_wall_loads` already drops all-zero loads).
+//!
+//! Supported: the cut-cell wall with the device's own closures (order 1 or
+//! 2, oblique distance, fine floor, axis exchange, centroid diffusion,
+//! exchange convection on/off). Not supported (`None`, the caller keeps the
+//! host route; logged once): the diagnostic load window, the gradient
+//! weights, the advancing-wall closure, the centroid foot, the S2-7 host
+//! prototypes, or a device mask not in its predictor phase.
+
+use super::DeviceStep;
+use super::cut::Phase;
+use crate::solvers::incompressible::embedded3::Grid;
+use crate::solvers::incompressible::embedded3::body::Body;
+use crate::solvers::incompressible::embedded3::exchange::{in_load_window, to_load_sink};
+use crate::solvers::incompressible::embedded3::interface::{LoadKind, WallLoad};
+use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime};
+use cudarc::driver::{
+    CudaFunction, CudaModule, CudaSlice, DeviceRepr, PushKernelArg, ValidAsZeroBits,
+};
+use std::sync::{Arc, OnceLock};
+
+const LOAD_KERNELS: &str = concat!(
+    include_str!("../../../../../kernels/cuda/e3_step.cu"),
+    include_str!("../../../../../kernels/cuda/e3_cut.cu"),
+    include_str!("../../../../../kernels/cuda/e3_loads.cu")
+);
+
+struct LoadKernels {
+    _module: Arc<CudaModule>,
+    cells: CudaFunction,
+    faces: CudaFunction,
+}
+
+static LOAD_ONCE: OnceLock<LoadKernels> = OnceLock::new();
+
+/// The record buffer's capacity (grown when a pass overflows it).
+static LOADS_CAP: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0);
+
+fn kernels() -> &'static LoadKernels {
+    LOAD_ONCE.get_or_init(|| {
+        // FMA contraction off: the host's roundings.
+        let module = load_module(LOAD_KERNELS, "e3_loads.cu", true);
+        let f = |name: &str| module.load_function(name).expect(name);
+        LoadKernels {
+            cells: f("e3_loads_cells"),
+            faces: f("e3_loads_faces"),
+            _module: module,
+        }
+    })
+}
+
+/// `RTX_E3_LOADS_DEVICE=1` (read once): callers take the device route.
+pub fn loads_device_enabled() -> bool {
+    static ON: OnceLock<bool> = OnceLock::new();
+    *ON.get_or_init(|| std::env::var("RTX_E3_LOADS_DEVICE").is_ok_and(|v| v == "1"))
+}
+
+fn log_fallback(reason: &str) {
+    static LOGGED: std::sync::Once = std::sync::Once::new();
+    LOGGED.call_once(|| {
+        eprintln!("    R8-f device loads: host route ({reason})");
+    });
+}
+
+/// `struct LoadRec` in e3_loads.cu.
+#[repr(C)]
+#[derive(Clone, Copy, Debug, Default)]
+struct LoadRec {
+    key: u32,
+    idx: u32,
+    v: [f64; 3],
+}
+unsafe impl DeviceRepr for LoadRec {}
+unsafe impl ValidAsZeroBits for LoadRec {}
+
+/// `struct LoadArgs` in e3_loads.cu.
+#[repr(C)]
+#[derive(Clone, Copy)]
+struct LoadArgs {
+    mu: f64,
+    rho: f64,
+    k0: i32,
+    k1: i32,
+    shifts: i32,
+    centroid: i32,
+    axis: i32,
+    conv_off: i32,
+    order: i32,
+    cap: u32,
+}
+unsafe impl DeviceRepr for LoadArgs {}
+unsafe impl ValidAsZeroBits for LoadArgs {}
+
+/// One summand of the operator route (a force ON the body), in the host's
+/// loop order.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct WallLoadRecord {
+    /// R8-a's part numbering: 0 the cell's `p W`, 1 the face's wall shear,
+    /// 2 the diffusive and 3 the convective exchange.
+    pub part: usize,
+    /// The face's component (0 for a cell).
+    pub c: usize,
+    /// The cell or face index (in its own lattice).
+    pub index: usize,
+    /// The cell centre or the face position.
+    pub pos: [f64; 3],
+    /// The force; a face summand has only component `c`.
+    pub value: [f64; 3],
+}
+
+/// The device route's parts (the host's sums, to the bit) and its records.
+#[derive(Clone, Debug, Default)]
+pub struct DeviceWallLoads {
+    pub pressure: [f64; 3],
+    pub shear: [f64; 3],
+    pub diffusive: [f64; 3],
+    pub convective: [f64; 3],
+    pub records: Vec<WallLoadRecord>,
+    /// `Some((k0, k1))` for a plane-restricted route.
+    pub planes: Option<(usize, usize)>,
+}
+
+impl DeviceWallLoads {
+    /// The route's total as `Mask::cut_wall_force` forms it.
+    #[must_use]
+    pub fn total(&self) -> [f64; 3] {
+        let (p, s) = (self.pressure, self.shear);
+        let x = [
+            self.diffusive[0] + self.convective[0],
+            self.diffusive[1] + self.convective[1],
+            self.diffusive[2] + self.convective[2],
+        ];
+        [p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]
+    }
+
+    /// The summands as the R8-a load sink receives them (position,
+    /// component, part, value), exact zeros dropped.
+    #[must_use]
+    pub fn sink_items(&self) -> Vec<([f64; 3], usize, usize, f64)> {
+        let mut out = Vec::with_capacity(self.records.len() + self.records.len() / 4);
+        for r in &self.records {
+            if r.part == 0 {
+                for c in 0..3 {
+                    if r.value[c] != 0.0 {
+                        out.push((r.pos, c, 0, r.value[c]));
+                    }
+                }
+            } else if r.value[r.c] != 0.0 {
+                out.push((r.pos, r.c, r.part, r.value[r.c]));
+            }
+        }
+        out
+    }
+
+    /// Hand the summands to the process-wide R8-a load sink (`exchange::
+    /// set_load_sink`), in the host route's order, exact zeros dropped.
+    pub fn replay_load_sink(&self) {
+        for (pos, c, part, v) in self.sink_items() {
+            to_load_sink(pos, c, part, v);
+        }
+    }
+
+    /// R8-c's summands (`Mask::cut_wall_loads`: kind, position, foot, force;
+    /// all-zero loads dropped), the foot by the body's own φ as there.
+    #[must_use]
+    pub fn wall_loads(&self, body: &Body, grid: Grid, t: f64) -> Vec<WallLoad> {
+        let raw: Vec<(LoadKind, [f64; 3], [f64; 3])> = self
+            .records
+            .iter()
+            .filter(|r| r.value != [0.0; 3])
+            .map(|r| {
+                let kind = match r.part {
+                    0 => LoadKind::Pressure,
+                    1 => LoadKind::Shear,
+                    2 => LoadKind::ExchangeDiffusive,
+                    _ => LoadKind::ExchangeConvective,
+                };
+                (kind, r.pos, r.value)
+            })
+            .collect();
+        // interface.rs `cut_wall_loads`: two projections onto φ = 0.
+        let eps = 1e-6 * grid.dx.min(grid.dy).min(grid.dz);
+        let foot = |x: [f64; 3]| {
+            let mut q = x;
+            for _ in 0..2 {
+                let s = body.phi(q[0], q[1], q[2], t);
+                let n = body.normal(q[0], q[1], q[2], t, eps);
+                q = [q[0] - s * n.0, q[1] - s * n.1, q[2] - s * n.2];
+            }
+            q
+        };
+        use rayon::prelude::*;
+        raw.par_iter()
+            .map(|&(kind, x, f)| WallLoad {
+                kind,
+                x,
+                foot: foot(x),
+                f,
+            })
+            .collect()
+    }
+}
+
+/// The G1 instrument's verdict for one call ([`DeviceStep::check_wall_loads`]).
+#[derive(Clone, Debug)]
+pub struct WallLoadsCheck {
+    /// Host and device totals (per span when plane-restricted).
+    pub host: [f64; 3],
+    pub device: [f64; 3],
+    /// The totals bit-identical.
+    pub totals_identical: bool,
+    /// The host sink's summands (incl. exact zeros), the non-zero ones, and
+    /// the device's replayed ones.
+    pub host_items: usize,
+    pub host_nonzero: usize,
+    pub device_items: usize,
+    /// Non-zero host summands and device summands identical item by item
+    /// (position, component, part, value bits).
+    pub items_identical: bool,
+    /// The sums of the two item lists (component-wise, sequential).
+    pub host_item_sum: [f64; 3],
+    pub device_item_sum: [f64; 3],
+    /// The records summed in the device's (atomic, unsorted) order: the
+    /// relative difference a device-order reduction would have made.
+    pub device_order_rel: f64,
+    /// R8-c's summands (whole body, no planes): identical to the host's
+    /// `cut_wall_loads` (None when not checked).
+    pub r8c_identical: Option<bool>,
+    /// Wall time of the host route (on the downloaded field) and of the
+    /// device route (kernels, download of the records, sort, sums), ms.
+    pub host_ms: f64,
+    pub device_ms: f64,
+}
+
+impl DeviceStep {
+    /// Why the device route cannot serve this solver (None = it can).
+    fn loads_unsupported(&self) -> Option<&'static str> {
+        let mask = self.solver.mask()?;
+        if mask.cut.is_none() {
+            return Some("no cut geometry");
+        }
+        let Some(dc) = self.cut.as_ref() else {
+            return Some("no device mask");
+        };
+        if dc.phase != Phase::Predictor {
+            return Some("the device mask is not in its predictor phase");
+        }
+        if !in_load_window(f64::NAN) {
+            return Some("the diagnostic load window is set");
+        }
+        if mask.grad_weights.is_some() {
+            return Some("the gradient weights (S2-5 host prototype)");
+        }
+        if mask.wall_advancing {
+            return Some("the advancing-wall closure");
+        }
+        if mask.wall_foot_centroid {
+            return Some("the centroid foot");
+        }
+        if mask.cv_sides_exact || mask.wall_order2_centroid {
+            return Some("an S2-7 host prototype");
+        }
+        if mask.periodic_z
+            != (self.solver.params.boundaries.z0
+                == crate::solvers::incompressible::embedded3::step::Side::Periodic)
+        {
+            return Some("the mask's z periodicity differs from the device's");
+        }
+        None
+    }
+
+    /// The operator route's parts and records on the device over the planes
+    /// `planes` (all when `None`), with `mu` the dynamic viscosity (the host
+    /// route's argument). `None` when the device cannot serve it (the caller
+    /// keeps the host route; the reason is logged once).
+    pub fn cut_wall_loads_device(
+        &self,
+        mu: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<DeviceWallLoads> {
+        if let Some(reason) = self.loads_unsupported() {
+            log_fallback(reason);
+            return None;
+        }
+        self.solver.mask()?;
+        let (recs, _) = self.load_records(mu, planes);
+        let mut records = recs;
+        records.sort_unstable_by_key(|r| (r.key >> 24, r.idx, r.key & 0x00ff_ffff));
+        Some(self.assemble(&records, planes))
+    }
+
+    /// The route's total like `Mask::cut_wall_force` (whole body, `None`
+    /// planes) or `Mask::cut_wall_force_per_span` (the planes, divided by
+    /// their thickness), with the summands handed to the R8-a load sink as
+    /// the host route hands them (exact zeros dropped). `None`: the caller
+    /// keeps the host route.
+    pub fn cut_wall_force_device(
+        &self,
+        mu: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<[f64; 3]> {
+        let loads = self.cut_wall_loads_device(mu, planes)?;
+        loads.replay_load_sink();
+        let t = loads.total();
+        Some(match planes {
+            None => t,
+            Some((k0, k1)) => {
+                let lz = (k1 - k0) as f64 * self.grid.dz;
+                [t[0] / lz, t[1] / lz, t[2] / lz]
+            }
+        })
+    }
+
+    /// G1 instrument (`RTX_E3_LOADS_CHECK=1` in the drivers): the device
+    /// route against the host route on the downloaded field — totals, the
+    /// R8-a sink's items and (whole body) R8-c's summands, bit for bit.
+    /// Installs its own capture sink: no load sink may be installed.
+    pub fn check_wall_loads(
+        &self,
+        mu: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<WallLoadsCheck> {
+        use crate::solvers::incompressible::embedded3::exchange::set_load_sink;
+        use crate::solvers::incompressible::embedded3::field::Field;
+        type Items = Vec<([f64; 3], usize, usize, f64)>;
+        if self.loads_unsupported().is_some() {
+            return None;
+        }
+        let mut field = Field::new(self.grid);
+        self.download(&mut field);
+        let t = self.solver.time();
+        let capture = || {
+            let store = Arc::new(std::sync::Mutex::new(Items::new()));
+            let s = store.clone();
+            let prev = set_load_sink(Some(Box::new(move |pos, c, part, v| {
+                s.lock().expect("capture").push((pos, c, part, v));
+            })));
+            assert!(
+                prev.is_none(),
+                "check_wall_loads: a load sink was already installed"
+            );
+            store
+        };
+        let take = |store: Arc<std::sync::Mutex<Items>>| -> Items {
+            set_load_sink(None);
+            std::mem::take(&mut *store.lock().expect("capture"))
+        };
+        // Host.
+        let (mask, body) = (self.solver.mask()?, self.solver.body()?);
+        let lap = std::time::Instant::now();
+        let store = capture();
+        let host = match planes {
+            None => mask.cut_wall_force(body, &field, mu, t),
+            Some(p) => mask.cut_wall_force_per_span(body, &field, mu, t, p),
+        };
+        let host_items = take(store);
+        let host_ms = lap.elapsed().as_secs_f64() * 1e3;
+        let host = host?;
+        // The device route's own cost (what a caller pays).
+        let lap = std::time::Instant::now();
+        let timed = self.cut_wall_loads_device(mu, planes);
+        let device_ms = lap.elapsed().as_secs_f64() * 1e3;
+        drop(timed);
+        // Device (records in both orders).
+        let (raw, _) = self.load_records(mu, planes);
+        let mut device_order = [0.0f64; 3];
+        for r in &raw {
+            let ll = self.assemble(std::slice::from_ref(r), planes);
+            let v = ll.total();
+            for c in 0..3 {
+                device_order[c] += v[c];
+            }
+        }
+        let mut sorted = raw;
+        sorted.sort_unstable_by_key(|r| (r.key >> 24, r.idx, r.key & 0x00ff_ffff));
+        let loads = self.assemble(&sorted, planes);
+        let store = capture();
+        loads.replay_load_sink();
+        let device_items = take(store);
+        let device = match planes {
+            None => loads.total(),
+            Some((k0, k1)) => {
+                let lz = (k1 - k0) as f64 * self.grid.dz;
+                let v = loads.total();
+                [v[0] / lz, v[1] / lz, v[2] / lz]
+            }
+        };
+        let bits = |v: [f64; 3]| v.map(f64::to_bits);
+        let nonzero: Items = host_items.iter().copied().filter(|x| x.3 != 0.0).collect();
+        let same_item = |a: &([f64; 3], usize, usize, f64), b: &([f64; 3], usize, usize, f64)| {
+            bits(a.0) == bits(b.0) && a.1 == b.1 && a.2 == b.2 && a.3.to_bits() == b.3.to_bits()
+        };
+        let items_identical = nonzero.len() == device_items.len()
+            && nonzero
+                .iter()
+                .zip(&device_items)
+                .all(|(a, b)| same_item(a, b));
+        let item_sum = |v: &Items| {
+            let mut s = [0.0f64; 3];
+            for x in v {
+                s[x.1] += x.3;
+            }
+            s
+        };
+        let scale = host
+            .iter()
+            .fold(0.0f64, |m, v| m.max(v.abs()))
+            .max(f64::MIN_POSITIVE);
+        let lz = planes.map_or(1.0, |(k0, k1)| (k1 - k0) as f64 * self.grid.dz);
+        let device_order_rel = (0..3)
+            .map(|c| (device_order[c] / lz - host[c]).abs() / scale)
+            .fold(0.0f64, f64::max);
+        let r8c_identical = if planes.is_none() {
+            let mask = self.solver.mask()?;
+            let body = self.solver.body()?;
+            mask.cut_wall_loads(body, &field, mu, t).map(|(h, route)| {
+                let d = loads.wall_loads(body, self.grid, t);
+                bits(route) == bits(loads.total())
+                    && h.len() == d.len()
+                    && h.iter().zip(&d).all(|(a, b)| {
+                        a.kind == b.kind
+                            && bits(a.x) == bits(b.x)
+                            && bits(a.foot) == bits(b.foot)
+                            && bits(a.f) == bits(b.f)
+                    })
+            })
+        } else {
+            None
+        };
+        Some(WallLoadsCheck {
+            host,
+            device,
+            totals_identical: bits(host) == bits(device),
+            host_items: host_items.len(),
+            host_nonzero: nonzero.len(),
+            device_items: device_items.len(),
+            items_identical,
+            host_item_sum: item_sum(&nonzero),
+            device_item_sum: item_sum(&device_items),
+            device_order_rel,
+            r8c_identical,
+            host_ms,
+            device_ms,
+        })
+    }
+
+    /// The records of one pass (unsorted: the device's append order) and
+    /// the number the kernels produced.
+    fn load_records(&self, mu: f64, planes: Option<(usize, usize)>) -> (Vec<LoadRec>, usize) {
+        let rt = runtime();
+        let k = kernels();
+        let g = self.grid;
+        let mask = self.solver.mask().expect("mask");
+        let (k0, k1) = planes.unwrap_or((0, g.nz));
+        let mut prm = self.params(1.0);
+        // cut_cv reads the order's low bits, the oblique distance (16) and
+        // the fine floor (256): the MASK's settings (the host route's).
+        prm.wall_order = i32::from(mask.wall_order)
+            + 16 * i32::from(mask.wall_distance_oblique)
+            + 256 * i32::from(mask.distance_floor_fine);
+        prm.scheme = super::scheme_code(mask.scheme);
+        let mut args = LoadArgs {
+            mu,
+            rho: mask.density,
+            k0: k0 as i32,
+            k1: k1 as i32,
+            shifts: i32::from(mask.face_shifts.is_some()),
+            centroid: i32::from(mask.diffusion_centroid),
+            axis: i32::from(mask.wall_exchange_axis),
+            conv_off: i32::from(mask.exchange_convection_off),
+            order: i32::from(mask.wall_order),
+            cap: 0,
+        };
+        let ptrs = self.ptrs();
+        let cptrs = self.cut.as_ref().expect("cut").ptrs();
+        let counts = [
+            (g.nx + 1) * g.ny * g.nz,
+            g.nx * (g.ny + 1) * g.nz,
+            g.nx * g.ny * (g.nz + 1),
+        ];
+        let cells = g.nx * g.ny * (k1 - k0);
+        let mut cap = LOADS_CAP
+            .load(std::sync::atomic::Ordering::Relaxed)
+            .max(1 << 16);
+        loop {
+            args.cap = cap as u32;
+            let mut out = rt
+                .stream
+                .alloc_zeros::<LoadRec>(cap)
+                .expect("alloc records");
+            let mut count = rt.stream.alloc_zeros::<u32>(1).expect("alloc count");
+            unsafe {
+                rt.stream
+                    .launch_builder(&k.cells)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&args)
+                    .arg(&mut out)
+                    .arg(&mut count)
+                    .launch(cfg(cells))
+                    .expect("e3_loads_cells");
+                for c in 0..3i32 {
+                    rt.stream
+                        .launch_builder(&k.faces)
+                        .arg(&prm)
+                        .arg(&ptrs)
+                        .arg(&cptrs)
+                        .arg(&args)
+                        .arg(&c)
+                        .arg(&mut out)
+                        .arg(&mut count)
+                        .launch(cfg(counts[c as usize]))
+                        .expect("e3_loads_faces");
+                }
+            }
+            let n: Vec<u32> = rt.stream.memcpy_dtov(&count).expect("count");
+            let n = n[0] as usize;
+            if n <= cap {
+                let recs: Vec<LoadRec> = if n == 0 {
+                    Vec::new()
+                } else {
+                    rt.stream.memcpy_dtov(&out.slice(0..n)).expect("records")
+                };
+                LOADS_CAP.store(cap, std::sync::atomic::Ordering::Relaxed);
+                return (recs, n);
+            }
+            cap = n + n / 4;
+        }
+    }
+
+    /// The parts summed in the order of `records` (sequentially, as the host
+    /// loops do) and the records with their positions.
+    fn assemble(&self, records: &[LoadRec], planes: Option<(usize, usize)>) -> DeviceWallLoads {
+        let g = self.grid;
+        let h = [g.dx, g.dy, g.dz];
+        let (nx, ny) = (g.nx, g.ny);
+        let mut out = DeviceWallLoads {
+            planes,
+            records: Vec::with_capacity(records.len()),
+            ..DeviceWallLoads::default()
+        };
+        for r in records {
+            let kind = r.key >> 28;
+            let c = ((r.key >> 24) & 0xf) as usize;
+            let slot = (r.key & 0x00ff_ffff) as usize;
+            let idx = r.idx as usize;
+            let (part, pos) = if kind == 0 {
+                let (k, j, i) = (idx / (nx * ny), (idx / nx) % ny, idx % nx);
+                (
+                    0,
+                    [
+                        (i as f64 + 0.5) * g.dx,
+                        (j as f64 + 0.5) * g.dy,
+                        (k as f64 + 0.5) * g.dz,
+                    ],
+                )
+            } else {
+                let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1));
+                let p = [idx % ni, (idx / ni) % nj, idx / (ni * nj)];
+                let mut x = [0.0; 3];
+                for d in 0..3 {
+                    // cutwall.rs `Lattice::face_position`.
+                    let off = if d == c { 0.0 } else { 0.5 };
+                    x[d] = (p[d] as f64 + off) * h[d];
+                }
+                let part = if kind == 1 {
+                    1
+                } else if slot & 1 == 0 {
+                    3
+                } else {
+                    2
+                };
+                (part, x)
+            };
+            match part {
+                0 => {
+                    for cc in 0..3 {
+                        out.pressure[cc] += r.v[cc];
+                    }
+                }
+                1 => out.shear[c] += r.v[c],
+                2 => out.diffusive[c] += r.v[c],
+                _ => out.convective[c] += r.v[c],
+            }
+            out.records.push(WallLoadRecord {
+                part,
+                c,
+                index: idx,
+                pos,
+                value: r.v,
+            });
+        }
+        // cutwall.rs adds the gradient-weight force (zero: the weights are
+        // off on this route) to the pressure part.
+        for cc in 0..3 {
+            out.pressure[cc] += 0.0;
+        }
+        out
+    }
+}
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs
index 7c08033..859bea3 100644
--- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs
+++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs
@@ -650,23 +650,80 @@ fn flag_wake_on_the_device() {
     let (mut drag_sum, mut lift_min, mut lift_max, mut samples) =
         (0.0, f64::INFINITY, f64::NEG_INFINITY, 0usize);
     let mut worst_residual = 0.0_f64;
+    // R8-f: `RTX_E3_LOADS_DEVICE=1` takes the operator route on the device
+    // (the per-span and whole-body totals below); `RTX_E3_LOADS_CHECK=1`
+    // checks the device route against the host on EVERY step (both the
+    // slab's planes and the whole body), bit for bit.
+    let loads_device =
+        rtx_cfd::solvers::incompressible::embedded3::step::device::loads_device_enabled();
+    let loads_check = std::env::var("RTX_E3_LOADS_CHECK").is_ok_and(|v| v == "1");
+    let mut loads_bad = 0usize;
     for step in 0..steps {
         let r = device.advance(dt);
         worst_residual = worst_residual.max(r.final_residual);
         assert!(r.final_residual.is_finite(), "death at step {step}");
         let t = device.solver.time();
+        if loads_check {
+            for (name, planes) in [("slab", Some(slab)), ("whole", None)] {
+                let ck = device
+                    .check_wall_loads(RHO * NU, planes)
+                    .expect("R8-f: the device route is unsupported here");
+                let ok =
+                    ck.totals_identical && ck.items_identical && ck.r8c_identical.unwrap_or(true);
+                loads_bad += usize::from(!ok);
+                println!(
+                    "  loads check step {step} t {t:.5} {name}: {} — totals host ({:+.12e}, {:+.12e}, {:+.12e}) device ({:+.12e}, {:+.12e}, {:+.12e}); items host {} ({} non-zero) device {} {}; item sums x {:+.12e}/{:+.12e} y {:+.12e}/{:+.12e}; R8-c {}; device-order rel {:.2e}; host {:.1} ms device {:.1} ms",
+                    if ok { "IDENTICAL" } else { "DIFFERS" },
+                    ck.host[0],
+                    ck.host[1],
+                    ck.host[2],
+                    ck.device[0],
+                    ck.device[1],
+                    ck.device[2],
+                    ck.host_items,
+                    ck.host_nonzero,
+                    ck.device_items,
+                    if ck.items_identical {
+                        "identical"
+                    } else {
+                        "DIFFER"
+                    },
+                    ck.host_item_sum[0],
+                    ck.device_item_sum[0],
+                    ck.host_item_sum[1],
+                    ck.device_item_sum[1],
+                    match ck.r8c_identical {
+                        Some(true) => "identical",
+                        Some(false) => "DIFFERS",
+                        None => "n/a",
+                    },
+                    ck.device_order_rel,
+                    ck.host_ms,
+                    ck.device_ms
+                );
+            }
+        }
         let sample = (step + 1) % 10 == 0 || step + 1 == steps;
         let phase_due = vtk_dir.is_some()
             && t >= last_period_start + next_phase as f64 * period / phases as f64
             && next_phase < total_phases;
         if sample || phase_due {
+            let dev = loads_device.then(|| {
+                (
+                    device.cut_wall_force_device(RHO * NU, Some(slab)),
+                    device.cut_wall_force_device(RHO * NU, None),
+                )
+            });
             device.download(&mut field);
             let solver = &device.solver;
             let mask = solver.mask().expect("mask");
             let body = solver.body().expect("body");
-            let fs = mask
-                .cut_wall_force_per_span(body, &field, RHO * NU, t, slab)
-                .expect("wall");
+            let fs = match dev {
+                Some((Some(v), _)) => v,
+                _ => mask
+                    .cut_wall_force_per_span(body, &field, RHO * NU, t, slab)
+                    .expect("wall"),
+            };
             // The reconstructed wall route on the same slab, per span.
             let fr = mask
                 .cut_wall_force_reconstructed(body, &field, RHO * NU, t, Some(slab))
@@ -675,9 +732,12 @@ fn flag_wake_on_the_device() {
                     [v[0] / lz, v[1] / lz, v[2] / lz]
                 })
                 .expect("reconstructed");
-            let ft = mask
-                .cut_wall_force(body, &field, RHO * NU, t)
-                .expect("wall");
+            let ft = match dev {
+                Some((_, Some(v))) => v,
+                _ => mask
+                    .cut_wall_force(body, &field, RHO * NU, t)
+                    .expect("wall"),
+            };
             if sample {
                 samples_seen += 1;
             }
@@ -910,4 +970,8 @@ fn flag_wake_on_the_device() {
     if let Some(t) = device.timers() {
         println!("  timers: {t:?}");
     }
+    if loads_check {
+        println!("  R8-f LOADS CHECK: {loads_bad} calls differ");
+        assert_eq!(loads_bad, 0, "R8-f: the device route differs from the host");
+    }
 }
diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs
index 4ba4ef8..29af92b 100644
--- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs
+++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs
@@ -20,7 +20,9 @@ use std::sync::{Arc, Mutex, RwLock};
 
 use rtx_cfd::solvers::incompressible::ConvectionScheme;
 use rtx_cfd::solvers::incompressible::embedded3::exchange::set_load_sink;
-use rtx_cfd::solvers::incompressible::embedded3::step::device::{DeviceSnapshot, DeviceStep};
+use rtx_cfd::solvers::incompressible::embedded3::step::device::{
+    DeviceSnapshot, DeviceStep, loads_device_enabled,
+};
 use rtx_cfd::solvers::incompressible::embedded3::{
     Body, Boundaries, DeviceSdf, Field, Fluid, Grid, Parameters, Side, Solver, StepResult,
     WallScheme,
@@ -334,7 +336,24 @@ impl E3Fluid {
         // flow; the whole span by default.
         let nz = self.grid.nz;
         let n = (super::env_f("RTX_E3FSI_LOAD_PLANES", 0.0) as usize).min(nz);
-        let f = if n > 0 && n < nz {
+        // R8-f (`RTX_E3_LOADS_DEVICE=1`): the same route on the device, its
+        // summands replayed into the sink in the host's order (exact zeros
+        // dropped); the host route when the device cannot serve it.
+        let planes = (n > 0 && n < nz).then(|| ((nz - n) / 2, (nz - n) / 2 + n));
+        let on_device = if loads_device_enabled() {
+            self.device.cut_wall_force_device(RHO * NU, planes)
+        } else {
+            None
+        };
+        let f = if let Some(f) = on_device {
+            if let Some((k0, k1)) = planes {
+                self.load_width = (k1 - k0) as f64 * self.h;
+                f
+            } else {
+                self.load_width = self.width;
+                [f[0] / self.width, f[1] / self.width, f[2] / self.width]
+            }
+        } else if n > 0 && n < nz {
             let k0 = (nz - n) / 2;
             self.load_width = n as f64 * self.h;
             mask.cut_wall_force_per_span(body, &self.field, RHO * NU, t, (k0, k0 + n))