rtx-cfd embedded3 item 7: e3_step.cu + step::device::{DeviceStep, StepTimers} on the shared runtime (FMA off); gate 7 HELD: device = host ≤ 7e-12 under tight tolerances on MMS/Beltrami/Poiseuille, equal CG counts on every step, periodic planes within 2e-16; the default-tolerance differences are the projection's inner stop

Co-Authored-By: Claude Fable 5.1 <[email protected]>
2026-09-17 14:56:46 -05:00
parent 8821e18520
commit 54911b4db3
4 changed files with 1555 additions and 0 deletions
@@ -0,0 +1,449 @@
+/**
+ * embedded3 item 7: the PISO step's maps on the device — the three
+ * predictors (the host `step/predictor.rs` expression for expression; the
+ * module is compiled with FMA contraction OFF so the f64 arithmetic is the
+ * host's), the normal-velocity stamping from per-side tables, the
+ * continuity source, the corrections, the pressure update and the mass
+ * imbalance partials. One thread per face / cell.
+ */
+#define SIDE_VELOCITY 0
+#define SIDE_SLIP 1
+#define SIDE_OUTLET 2
+#define SIDE_PERIODIC 3
+#define SCHEME_UPWIND 0
+#define SCHEME_VAN_ALBADA 1
+#define SCHEME_VAN_LEER 2
+
+struct E3Params {
+    int nx, ny, nz;
+    int periodic_z;
+    int bx0, bx1, by0, by1, bz0, bz1;
+    int scheme;
+    int pad;
+    double dx, dy, dz, dt, rho, nu;
+};
+
+/* Per-side tables: [side][component]; see `BoundaryTables` in step/device.rs. */
+struct E3Ptrs {
+    double *u, *v, *w, *p, *uo, *vo, *wo, *us, *vs, *ws, *pp, *sp;
+    const double *su, *sv, *sw;
+    const double *bx0u, *bx0v, *bx0w, *bx1u, *bx1v, *bx1w;
+    const double *by0u, *by0v, *by0w, *by1u, *by1v, *by1w;
+    const double *bz0u, *bz0v, *bz0w, *bz1u, *bz1v, *bz1w;
+};
+
+__device__ __forceinline__ int cell3(const E3Params& g, int k, int j, int i) { return (k * g.ny + j) * g.nx + i; }
+__device__ __forceinline__ int uf3(const E3Params& g, int k, int j, int i) { return (k * g.ny + j) * (g.nx + 1) + i; }
+__device__ __forceinline__ int vf3(const E3Params& g, int k, int j, int i) { return (k * (g.ny + 1) + j) * g.nx + i; }
+__device__ __forceinline__ int wf3(const E3Params& g, int k, int j, int i) { return (k * g.ny + j) * g.nx + i; }
+
+__device__ __forceinline__ double upwind3(double face_velocity, double upstream, double downstream) {
+    return face_velocity >= 0.0 ? upstream : downstream;
+}
+
+__device__ __forceinline__ double limiter3(int scheme, double r) {
+    if (scheme == SCHEME_VAN_ALBADA) return r > 0.0 ? (r * r + r) / (r * r + 1.0) : 0.0;
+    if (scheme == SCHEME_VAN_LEER) return (r + fabs(r)) / (1.0 + fabs(r));
+    return 0.0;
+}
+
+/* The limited correction; `has_far` = the far-upwind node exists. */
+__device__ __forceinline__ double face_corr3(int scheme, int has_far, double far, double up, double down) {
+    if (!has_far) return 0.0;
+    double denominator = down - up;
+    if (fabs(denominator) < 1e-300) return 0.0;
+    double r = (up - far) / denominator;
+    return 0.5 * limiter3(scheme, r) * denominator;
+}
+
+/* k above / below with the periodic wrap; −1 = wall. */
+__device__ __forceinline__ int k_up3(const E3Params& g, int k) { return k + 1 < g.nz ? k + 1 : (g.periodic_z ? 0 : -1); }
+__device__ __forceinline__ int k_dn3(const E3Params& g, int k) { return k > 0 ? k - 1 : (g.periodic_z ? g.nz - 1 : -1); }
+
+extern "C" __global__ void e3_step_predict_u(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int nxp = g.nx + 1;
+    int total = nxp * g.ny * g.nz;
+    if (t >= total) return;
+    int i = t % nxp; int j = (t / nxp) % g.ny; int k = t / (nxp * g.ny);
+    if (i == 0 || i == g.nx) return;
+    const double *uo = f.uo, *vo = f.vo, *wo = f.wo;
+    double dx = g.dx, dy = g.dy, dz = g.dz, rho = g.rho, nu = g.nu;
+    int scheme = g.scheme;
+    double u_p = uo[uf3(g, k, j, i)];
+    double ue_face = 0.5 * (uo[uf3(g, k, j, i)] + uo[uf3(g, k, j, i + 1)]);
+    double uw_face = 0.5 * (uo[uf3(g, k, j, i - 1)] + uo[uf3(g, k, j, i)]);
+    int south_is_wall = (j == 0);
+    int north_is_wall = (j + 1 == g.ny);
+    double vn_face = 0.5 * (vo[vf3(g, k, j + 1, i - 1)] + vo[vf3(g, k, j + 1, i)]);
+    double vs_face = 0.5 * (vo[vf3(g, k, j, i - 1)] + vo[vf3(g, k, j, i)]);
+    double beyond_north = (g.by1 == SIDE_VELOCITY) ? f.by1u[k * nxp + i] : u_p;
+    double beyond_south = (g.by0 == SIDE_VELOCITY) ? f.by0u[k * nxp + i] : u_p;
+    double conv_x = (ue_face * upwind3(ue_face, uo[uf3(g, k, j, i)], uo[uf3(g, k, j, i + 1)])
+        - uw_face * upwind3(uw_face, uo[uf3(g, k, j, i - 1)], uo[uf3(g, k, j, i)])) / dx;
+    double conv_y = (vn_face * (north_is_wall ? upwind3(vn_face, u_p, beyond_north)
+                                              : upwind3(vn_face, uo[uf3(g, k, j, i)], uo[uf3(g, k, j + 1, i)]))
+        - vs_face * (south_is_wall ? upwind3(vs_face, beyond_south, u_p)
+                                   : upwind3(vs_face, uo[uf3(g, k, j - 1, i)], uo[uf3(g, k, j, i)]))) / dy;
+    if (scheme != SCHEME_UPWIND) {
+        double delta_e, delta_w, delta_n, delta_s;
+        if (ue_face >= 0.0) delta_e = face_corr3(scheme, 1, uo[uf3(g, k, j, i - 1)], uo[uf3(g, k, j, i)], uo[uf3(g, k, j, i + 1)]);
+        else { int has = (i + 2 <= g.nx); delta_e = face_corr3(scheme, has, has ? uo[uf3(g, k, j, i + 2)] : 0.0, uo[uf3(g, k, j, i + 1)], uo[uf3(g, k, j, i)]); }
+        if (uw_face >= 0.0) { int has = (i >= 2); delta_w = face_corr3(scheme, has, has ? uo[uf3(g, k, j, i - 2)] : 0.0, uo[uf3(g, k, j, i - 1)], uo[uf3(g, k, j, i)]); }
+        else delta_w = face_corr3(scheme, 1, uo[uf3(g, k, j, i + 1)], uo[uf3(g, k, j, i)], uo[uf3(g, k, j, i - 1)]);
+        if (north_is_wall) delta_n = 0.0;
+        else if (vn_face >= 0.0) { int has = (j >= 1); delta_n = face_corr3(scheme, has, has ? uo[uf3(g, k, j - 1, i)] : 0.0, uo[uf3(g, k, j, i)], uo[uf3(g, k, j + 1, i)]); }
+        else { int has = (j + 2 < g.ny); delta_n = face_corr3(scheme, has, has ? uo[uf3(g, k, j + 2, i)] : 0.0, uo[uf3(g, k, j + 1, i)], uo[uf3(g, k, j, i)]); }
+        if (south_is_wall) delta_s = 0.0;
+        else if (vs_face >= 0.0) { int has = (j >= 2); delta_s = face_corr3(scheme, has, has ? uo[uf3(g, k, j - 2, i)] : 0.0, uo[uf3(g, k, j - 1, i)], uo[uf3(g, k, j, i)]); }
+        else { int has = (j + 1 < g.ny); delta_s = face_corr3(scheme, has, has ? uo[uf3(g, k, j + 1, i)] : 0.0, uo[uf3(g, k, j, i)], uo[uf3(g, k, j - 1, i)]); }
+        conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+        conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+    }
+    double diff_x = nu * (uo[uf3(g, k, j, i + 1)] - 2.0 * u_p + uo[uf3(g, k, j, i - 1)]) / (dx * dx);
+    double flux_north = north_is_wall ? (g.by1 == SIDE_VELOCITY ? nu * (f.by1u[k * nxp + i] - u_p) / (0.5 * dy) : 0.0)
+                                      : nu * (uo[uf3(g, k, j + 1, i)] - u_p) / dy;
+    double flux_south = south_is_wall ? (g.by0 == SIDE_VELOCITY ? nu * (u_p - f.by0u[k * nxp + i]) / (0.5 * dy) : 0.0)
+                                      : nu * (u_p - uo[uf3(g, k, j - 1, i)]) / dy;
+    double diff_y = (flux_north - flux_south) / dy;
+    double pressure_gradient = -(f.p[cell3(g, k, j, i)] - f.p[cell3(g, k, j, i - 1)]) / (rho * dx);
+    double body_force = f.su[uf3(g, k, j, i)] / rho;
+    double rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force;
+    /* z terms */
+    int ku = k_up3(g, k), kd = k_dn3(g, k);
+    int top_is_wall = (ku < 0), bottom_is_wall = (kd < 0);
+    double wt_face = 0.5 * (wo[wf3(g, k + 1, j, i - 1)] + wo[wf3(g, k + 1, j, i)]);
+    double wb_face = 0.5 * (wo[wf3(g, k, j, i - 1)] + wo[wf3(g, k, j, i)]);
+    double beyond_top = (g.bz1 == SIDE_VELOCITY) ? f.bz1u[j * nxp + i] : u_p;
+    double beyond_bottom = (g.bz0 == SIDE_VELOCITY) ? f.bz0u[j * nxp + i] : u_p;
+    double u_up = top_is_wall ? u_p : uo[uf3(g, ku, j, i)];
+    double u_dn = bottom_is_wall ? u_p : uo[uf3(g, kd, j, i)];
+    double conv_z = (wt_face * (top_is_wall ? upwind3(wt_face, u_p, beyond_top) : upwind3(wt_face, u_p, u_up))
+        - wb_face * (bottom_is_wall ? upwind3(wb_face, beyond_bottom, u_p) : upwind3(wb_face, u_dn, u_p))) / dz;
+    if (scheme != SCHEME_UPWIND) {
+        int ku2 = top_is_wall ? -1 : k_up3(g, ku);
+        int kd2 = bottom_is_wall ? -1 : k_dn3(g, kd);
+        double far_up2 = ku2 >= 0 ? uo[uf3(g, ku2, j, i)] : 0.0;
+        double far_dn2 = kd2 >= 0 ? uo[uf3(g, kd2, j, i)] : 0.0;
+        double delta_t, delta_b;
+        if (top_is_wall) delta_t = 0.0;
+        else if (wt_face >= 0.0) delta_t = face_corr3(scheme, !bottom_is_wall, u_dn, u_p, u_up);
+        else delta_t = face_corr3(scheme, ku2 >= 0, far_up2, u_up, u_p);
+        if (bottom_is_wall) delta_b = 0.0;
+        else if (wb_face >= 0.0) delta_b = face_corr3(scheme, kd2 >= 0, far_dn2, u_dn, u_p);
+        else delta_b = face_corr3(scheme, !top_is_wall, u_up, u_p, u_dn);
+        conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+    }
+    double flux_top = top_is_wall ? (g.bz1 == SIDE_VELOCITY ? nu * (beyond_top - u_p) / (0.5 * dz) : 0.0) : nu * (u_up - u_p) / dz;
+    double flux_bottom = bottom_is_wall ? (g.bz0 == SIDE_VELOCITY ? nu * (u_p - beyond_bottom) / (0.5 * dz) : 0.0) : nu * (u_p - u_dn) / dz;
+    double diff_z = (flux_top - flux_bottom) / dz;
+    double rhs = rhs_2d - conv_z + diff_z;
+    f.u[uf3(g, k, j, i)] = uo[uf3(g, k, j, i)] + g.dt * rhs;
+}
+
+extern "C" __global__ void e3_step_predict_v(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int nyp = g.ny + 1;
+    int total = g.nx * nyp * g.nz;
+    if (t >= total) return;
+    int i = t % g.nx; int j = (t / g.nx) % nyp; int k = t / (g.nx * nyp);
+    if (j == 0 || j == g.ny) return;
+    const double *uo = f.uo, *vo = f.vo, *wo = f.wo;
+    double dx = g.dx, dy = g.dy, dz = g.dz, rho = g.rho, nu = g.nu;
+    int scheme = g.scheme;
+    double v_p = vo[vf3(g, k, j, i)];
+    double vn_face = 0.5 * (vo[vf3(g, k, j, i)] + vo[vf3(g, k, j + 1, i)]);
+    double vs_face = 0.5 * (vo[vf3(g, k, j - 1, i)] + vo[vf3(g, k, j, i)]);
+    int west_is_wall = (i == 0), east_is_wall = (i + 1 == g.nx);
+    double ue_face = 0.5 * (uo[uf3(g, k, j - 1, i + 1)] + uo[uf3(g, k, j, i + 1)]);
+    double uw_face = 0.5 * (uo[uf3(g, k, j - 1, i)] + uo[uf3(g, k, j, i)]);
+    double beyond_east = (g.bx1 == SIDE_VELOCITY) ? f.bx1v[k * nyp + j] : v_p;
+    double beyond_west = (g.bx0 == SIDE_VELOCITY) ? f.bx0v[k * nyp + j] : v_p;
+    double conv_y = (vn_face * upwind3(vn_face, vo[vf3(g, k, j, i)], vo[vf3(g, k, j + 1, i)])
+        - vs_face * upwind3(vs_face, vo[vf3(g, k, j - 1, i)], vo[vf3(g, k, j, i)])) / dy;
+    double conv_x = (ue_face * (east_is_wall ? upwind3(ue_face, v_p, beyond_east)
+                                             : upwind3(ue_face, vo[vf3(g, k, j, i)], vo[vf3(g, k, j, i + 1)]))
+        - uw_face * (west_is_wall ? upwind3(uw_face, beyond_west, v_p)
+                                  : upwind3(uw_face, vo[vf3(g, k, j, i - 1)], vo[vf3(g, k, j, i)]))) / dx;
+    if (scheme != SCHEME_UPWIND) {
+        double delta_n, delta_s, delta_e, delta_w;
+        if (vn_face >= 0.0) delta_n = face_corr3(scheme, 1, vo[vf3(g, k, j - 1, i)], vo[vf3(g, k, j, i)], vo[vf3(g, k, j + 1, i)]);
+        else { int has = (j + 2 <= g.ny); delta_n = face_corr3(scheme, has, has ? vo[vf3(g, k, j + 2, i)] : 0.0, vo[vf3(g, k, j + 1, i)], vo[vf3(g, k, j, i)]); }
+        if (vs_face >= 0.0) { int has = (j >= 2); delta_s = face_corr3(scheme, has, has ? vo[vf3(g, k, j - 2, i)] : 0.0, vo[vf3(g, k, j - 1, i)], vo[vf3(g, k, j, i)]); }
+        else delta_s = face_corr3(scheme, 1, vo[vf3(g, k, j + 1, i)], vo[vf3(g, k, j, i)], vo[vf3(g, k, j - 1, i)]);
+        if (east_is_wall) delta_e = 0.0;
+        else if (ue_face >= 0.0) { int has = (i >= 1); delta_e = face_corr3(scheme, has, has ? vo[vf3(g, k, j, i - 1)] : 0.0, vo[vf3(g, k, j, i)], vo[vf3(g, k, j, i + 1)]); }
+        else { int has = (i + 2 < g.nx); delta_e = face_corr3(scheme, has, has ? vo[vf3(g, k, j, i + 2)] : 0.0, vo[vf3(g, k, j, i + 1)], vo[vf3(g, k, j, i)]); }
+        if (west_is_wall) delta_w = 0.0;
+        else if (uw_face >= 0.0) { int has = (i >= 2); delta_w = face_corr3(scheme, has, has ? vo[vf3(g, k, j, i - 2)] : 0.0, vo[vf3(g, k, j, i - 1)], vo[vf3(g, k, j, i)]); }
+        else { int has = (i + 1 < g.nx); delta_w = face_corr3(scheme, has, has ? vo[vf3(g, k, j, i + 1)] : 0.0, vo[vf3(g, k, j, i)], vo[vf3(g, k, j, i - 1)]); }
+        conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+        conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+    }
+    double diff_y = nu * (vo[vf3(g, k, j + 1, i)] - 2.0 * v_p + vo[vf3(g, k, j - 1, i)]) / (dy * dy);
+    double flux_east = east_is_wall ? (g.bx1 == SIDE_VELOCITY ? nu * (f.bx1v[k * nyp + j] - v_p) / (0.5 * dx) : 0.0)
+                                    : nu * (vo[vf3(g, k, j, i + 1)] - v_p) / dx;
+    double flux_west = west_is_wall ? (g.bx0 == SIDE_VELOCITY ? nu * (v_p - f.bx0v[k * nyp + j]) / (0.5 * dx) : 0.0)
+                                    : nu * (v_p - vo[vf3(g, k, j, i - 1)]) / dx;
+    double diff_x = (flux_east - flux_west) / dx;
+    double pressure_gradient = -(f.p[cell3(g, k, j, i)] - f.p[cell3(g, k, j - 1, i)]) / (rho * dy);
+    double body_force = f.sv[vf3(g, k, j, i)] / rho;
+    double rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force;
+    /* z terms */
+    int ku = k_up3(g, k), kd = k_dn3(g, k);
+    int top_is_wall = (ku < 0), bottom_is_wall = (kd < 0);
+    double wt_face = 0.5 * (wo[wf3(g, k + 1, j - 1, i)] + wo[wf3(g, k + 1, j, i)]);
+    double wb_face = 0.5 * (wo[wf3(g, k, j - 1, i)] + wo[wf3(g, k, j, i)]);
+    double beyond_top = (g.bz1 == SIDE_VELOCITY) ? f.bz1v[j * g.nx + i] : v_p;
+    double beyond_bottom = (g.bz0 == SIDE_VELOCITY) ? f.bz0v[j * g.nx + i] : v_p;
+    double v_up = top_is_wall ? v_p : vo[vf3(g, ku, j, i)];
+    double v_dn = bottom_is_wall ? v_p : vo[vf3(g, kd, j, i)];
+    double conv_z = (wt_face * (top_is_wall ? upwind3(wt_face, v_p, beyond_top) : upwind3(wt_face, v_p, v_up))
+        - wb_face * (bottom_is_wall ? upwind3(wb_face, beyond_bottom, v_p) : upwind3(wb_face, v_dn, v_p))) / dz;
+    if (scheme != SCHEME_UPWIND) {
+        int ku2 = top_is_wall ? -1 : k_up3(g, ku);
+        int kd2 = bottom_is_wall ? -1 : k_dn3(g, kd);
+        double far_up2 = ku2 >= 0 ? vo[vf3(g, ku2, j, i)] : 0.0;
+        double far_dn2 = kd2 >= 0 ? vo[vf3(g, kd2, j, i)] : 0.0;
+        double delta_t, delta_b;
+        if (top_is_wall) delta_t = 0.0;
+        else if (wt_face >= 0.0) delta_t = face_corr3(scheme, !bottom_is_wall, v_dn, v_p, v_up);
+        else delta_t = face_corr3(scheme, ku2 >= 0, far_up2, v_up, v_p);
+        if (bottom_is_wall) delta_b = 0.0;
+        else if (wb_face >= 0.0) delta_b = face_corr3(scheme, kd2 >= 0, far_dn2, v_dn, v_p);
+        else delta_b = face_corr3(scheme, !top_is_wall, v_up, v_p, v_dn);
+        conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+    }
+    double flux_top = top_is_wall ? (g.bz1 == SIDE_VELOCITY ? nu * (beyond_top - v_p) / (0.5 * dz) : 0.0) : nu * (v_up - v_p) / dz;
+    double flux_bottom = bottom_is_wall ? (g.bz0 == SIDE_VELOCITY ? nu * (v_p - beyond_bottom) / (0.5 * dz) : 0.0) : nu * (v_p - v_dn) / dz;
+    double diff_z = (flux_top - flux_bottom) / dz;
+    double rhs = rhs_2d - conv_z + diff_z;
+    f.v[vf3(g, k, j, i)] = vo[vf3(g, k, j, i)] + g.dt * rhs;
+}
+
+extern "C" __global__ void e3_step_predict_w(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int total = g.nx * g.ny * (g.nz + 1);
+    if (t >= total) return;
+    int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
+    int periodic = g.periodic_z;
+    if (periodic) { if (k == g.nz) return; } else { if (k == 0 || k == g.nz) return; }
+    const double *uo = f.uo, *vo = f.vo, *wo = f.wo;
+    double dx = g.dx, dy = g.dy, dz = g.dz, rho = g.rho, nu = g.nu;
+    int scheme = g.scheme;
+    int nz = g.nz;
+    int k_below = k > 0 ? k - 1 : nz - 1;
+    int k_above = k % nz;
+    int w_dn_idx = k > 0 ? wf3(g, k - 1, j, i) : wf3(g, nz - 1, j, i);
+    int w_up_idx = (k + 1 == nz && periodic) ? wf3(g, 0, j, i) : wf3(g, k + 1, j, i);
+    double w_p = wo[wf3(g, k, j, i)];
+    double wt_face = 0.5 * (wo[wf3(g, k, j, i)] + wo[w_up_idx]);
+    double wb_face = 0.5 * (wo[w_dn_idx] + wo[wf3(g, k, j, i)]);
+    int west_is_wall = (i == 0), east_is_wall = (i + 1 == g.nx);
+    int south_is_wall = (j == 0), north_is_wall = (j + 1 == g.ny);
+    double ue_face = 0.5 * (uo[uf3(g, k_below, j, i + 1)] + uo[uf3(g, k_above, j, i + 1)]);
+    double uw_face = 0.5 * (uo[uf3(g, k_below, j, i)] + uo[uf3(g, k_above, j, i)]);
+    double vn_face = 0.5 * (vo[vf3(g, k_below, j + 1, i)] + vo[vf3(g, k_above, j + 1, i)]);
+    double vs_face = 0.5 * (vo[vf3(g, k_below, j, i)] + vo[vf3(g, k_above, j, i)]);
+    double beyond_east = (g.bx1 == SIDE_VELOCITY) ? f.bx1w[k * g.ny + j] : w_p;
+    double beyond_west = (g.bx0 == SIDE_VELOCITY) ? f.bx0w[k * g.ny + j] : w_p;
+    double beyond_north = (g.by1 == SIDE_VELOCITY) ? f.by1w[k * g.nx + i] : w_p;
+    double beyond_south = (g.by0 == SIDE_VELOCITY) ? f.by0w[k * g.nx + i] : w_p;
+    double conv_z = (wt_face * upwind3(wt_face, w_p, wo[w_up_idx]) - wb_face * upwind3(wb_face, wo[w_dn_idx], w_p)) / dz;
+    double conv_x = (ue_face * (east_is_wall ? upwind3(ue_face, w_p, beyond_east) : upwind3(ue_face, w_p, wo[wf3(g, k, j, i + 1)]))
+        - uw_face * (west_is_wall ? upwind3(uw_face, beyond_west, w_p) : upwind3(uw_face, wo[wf3(g, k, j, i - 1)], w_p))) / dx;
+    double conv_y = (vn_face * (north_is_wall ? upwind3(vn_face, w_p, beyond_north) : upwind3(vn_face, w_p, wo[wf3(g, k, j + 1, i)]))
+        - vs_face * (south_is_wall ? upwind3(vs_face, beyond_south, w_p) : upwind3(vs_face, wo[wf3(g, k, j - 1, i)], w_p))) / dy;
+    if (scheme != SCHEME_UPWIND) {
+        int has_up2, has_dn2; double far_up2 = 0.0, far_dn2 = 0.0;
+        if (periodic) { has_up2 = 1; far_up2 = wo[wf3(g, (k + 2) % nz, j, i)]; has_dn2 = 1; far_dn2 = wo[wf3(g, (k + nz - 2) % nz, j, i)]; }
+        else { has_up2 = (k + 2 <= nz); if (has_up2) far_up2 = wo[wf3(g, k + 2, j, i)]; has_dn2 = (k >= 2); if (has_dn2) far_dn2 = wo[wf3(g, k - 2, j, i)]; }
+        double delta_t = wt_face >= 0.0 ? face_corr3(scheme, 1, wo[w_dn_idx], w_p, wo[w_up_idx]) : face_corr3(scheme, has_up2, far_up2, wo[w_up_idx], w_p);
+        double delta_b = wb_face >= 0.0 ? face_corr3(scheme, has_dn2, far_dn2, wo[w_dn_idx], w_p) : face_corr3(scheme, 1, wo[w_up_idx], w_p, wo[w_dn_idx]);
+        double delta_e, delta_w, delta_n, delta_s;
+        if (east_is_wall) delta_e = 0.0;
+        else if (ue_face >= 0.0) { int has = (i >= 1); delta_e = face_corr3(scheme, has, has ? wo[wf3(g, k, j, i - 1)] : 0.0, w_p, wo[wf3(g, k, j, i + 1)]); }
+        else { int has = (i + 2 < g.nx); delta_e = face_corr3(scheme, has, has ? wo[wf3(g, k, j, i + 2)] : 0.0, wo[wf3(g, k, j, i + 1)], w_p); }
+        if (west_is_wall) delta_w = 0.0;
+        else if (uw_face >= 0.0) { int has = (i >= 2); delta_w = face_corr3(scheme, has, has ? wo[wf3(g, k, j, i - 2)] : 0.0, wo[wf3(g, k, j, i - 1)], w_p); }
+        else { int has = (i + 1 < g.nx); delta_w = face_corr3(scheme, has, has ? wo[wf3(g, k, j, i + 1)] : 0.0, w_p, wo[wf3(g, k, j, i - 1)]); }
+        if (north_is_wall) delta_n = 0.0;
+        else if (vn_face >= 0.0) { int has = (j >= 1); delta_n = face_corr3(scheme, has, has ? wo[wf3(g, k, j - 1, i)] : 0.0, w_p, wo[wf3(g, k, j + 1, i)]); }
+        else { int has = (j + 2 < g.ny); delta_n = face_corr3(scheme, has, has ? wo[wf3(g, k, j + 2, i)] : 0.0, wo[wf3(g, k, j + 1, i)], w_p); }
+        if (south_is_wall) delta_s = 0.0;
+        else if (vs_face >= 0.0) { int has = (j >= 2); delta_s = face_corr3(scheme, has, has ? wo[wf3(g, k, j - 2, i)] : 0.0, wo[wf3(g, k, j - 1, i)], w_p); }
+        else { int has = (j + 1 < g.ny); delta_s = face_corr3(scheme, has, has ? wo[wf3(g, k, j + 1, i)] : 0.0, w_p, wo[wf3(g, k, j - 1, i)]); }
+        conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+        conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+        conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+    }
+    double diff_z = nu * (wo[w_up_idx] - 2.0 * w_p + wo[w_dn_idx]) / (dz * dz);
+    double flux_east = east_is_wall ? (g.bx1 == SIDE_VELOCITY ? nu * (beyond_east - w_p) / (0.5 * dx) : 0.0) : nu * (wo[wf3(g, k, j, i + 1)] - w_p) / dx;
+    double flux_west = west_is_wall ? (g.bx0 == SIDE_VELOCITY ? nu * (w_p - beyond_west) / (0.5 * dx) : 0.0) : nu * (w_p - wo[wf3(g, k, j, i - 1)]) / dx;
+    double diff_x = (flux_east - flux_west) / dx;
+    double flux_north = north_is_wall ? (g.by1 == SIDE_VELOCITY ? nu * (beyond_north - w_p) / (0.5 * dy) : 0.0) : nu * (wo[wf3(g, k, j + 1, i)] - w_p) / dy;
+    double flux_south = south_is_wall ? (g.by0 == SIDE_VELOCITY ? nu * (w_p - beyond_south) / (0.5 * dy) : 0.0) : nu * (w_p - wo[wf3(g, k, j - 1, i)]) / dy;
+    double diff_y = (flux_north - flux_south) / dy;
+    double pressure_gradient = -(f.p[cell3(g, k_above, j, i)] - f.p[cell3(g, k_below, j, i)]) / (rho * dz);
+    double body_force = f.sw[wf3(g, k, j, i)] / rho;
+    double rhs = -conv_x - conv_y - conv_z + diff_x + diff_y + diff_z + pressure_gradient + body_force;
+    f.w[wf3(g, k, j, i)] = wo[wf3(g, k, j, i)] + g.dt * rhs;
+}
+
+/* After the predictor: the periodic copy w[nz] = w[0], the outlet
+ * zero-gradient faces, and the normal stamping from the tables. One thread
+ * per (k, j) for the x sides, (k, i) for the y sides, (j, i) for the z
+ * sides — three kernels. */
+extern "C" __global__ void e3_step_sides_x(E3Params g, E3Ptrs f, int stamp)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    if (t >= g.ny * g.nz) return;
+    int j = t % g.ny, k = t / g.ny;
+    if (g.bx0 == SIDE_OUTLET) f.u[uf3(g, k, j, 0)] = f.u[uf3(g, k, j, 1)];
+    else if (stamp) f.u[uf3(g, k, j, 0)] = f.bx0u[k * g.ny + j];
+    if (g.bx1 == SIDE_OUTLET) f.u[uf3(g, k, j, g.nx)] = f.u[uf3(g, k, j, g.nx - 1)];
+    else if (stamp) f.u[uf3(g, k, j, g.nx)] = f.bx1u[k * g.ny + j];
+}
+
+extern "C" __global__ void e3_step_sides_y(E3Params g, E3Ptrs f, int stamp)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    if (t >= g.nx * g.nz) return;
+    int i = t % g.nx, k = t / g.nx;
+    if (g.by0 == SIDE_OUTLET) f.v[vf3(g, k, 0, i)] = f.v[vf3(g, k, 1, i)];
+    else if (stamp) f.v[vf3(g, k, 0, i)] = f.by0v[k * g.nx + i];
+    if (g.by1 == SIDE_OUTLET) f.v[vf3(g, k, g.ny, i)] = f.v[vf3(g, k, g.ny - 1, i)];
+    else if (stamp) f.v[vf3(g, k, g.ny, i)] = f.by1v[k * g.nx + i];
+}
+
+extern "C" __global__ void e3_step_sides_z(E3Params g, E3Ptrs f, int stamp)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    if (t >= g.nx * g.ny) return;
+    int i = t % g.nx, j = t / g.nx;
+    if (g.periodic_z) { f.w[wf3(g, g.nz, j, i)] = f.w[wf3(g, 0, j, i)]; return; }
+    if (g.bz0 == SIDE_OUTLET) f.w[wf3(g, 0, j, i)] = f.w[wf3(g, 1, j, i)];
+    else if (stamp) f.w[wf3(g, 0, j, i)] = f.bz0w[j * g.nx + i];
+    if (g.bz1 == SIDE_OUTLET) f.w[wf3(g, g.nz, j, i)] = f.w[wf3(g, g.nz - 1, j, i)];
+    else if (stamp) f.w[wf3(g, g.nz, j, i)] = f.bz1w[j * g.nx + i];
+}
+
+/* sp = −ρ Σ (flux out) over the cells; partial[block] = Σ |flux| (source scale). */
+extern "C" __global__ void e3_step_divergence(E3Params g, E3Ptrs f, double* __restrict__ partial)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int total = g.nx * g.ny * g.nz;
+    double v = 0.0;
+    if (t < total) {
+        int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
+        double divergence_flux = g.rho
+            * ((f.us[uf3(g, k, j, i + 1)] - f.us[uf3(g, k, j, i)]) * (g.dy * g.dz)
+               + (f.vs[vf3(g, k, j + 1, i)] - f.vs[vf3(g, k, j, i)]) * (g.dx * g.dz)
+               + (f.ws[wf3(g, k + 1, j, i)] - f.ws[wf3(g, k, j, i)]) * (g.dx * g.dy));
+        f.sp[t] = -divergence_flux;
+        v = fabs(divergence_flux);
+    }
+    __shared__ double sh[256];
+    sh[threadIdx.x] = v;
+    __syncthreads();
+    for (int s = 128; s > 0; s >>= 1) { if (threadIdx.x < s) sh[threadIdx.x] += sh[threadIdx.x + s]; __syncthreads(); }
+    if (threadIdx.x == 0) partial[blockIdx.x] = sh[0];
+}
+
+/* The corrections from p' (interior faces, outlet faces against 0 outside). */
+extern "C" __global__ void e3_step_correct_u(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int nxp = g.nx + 1;
+    if (t >= nxp * g.ny * g.nz) return;
+    int i = t % nxp; int j = (t / nxp) % g.ny; int k = t / (nxp * g.ny);
+    double c = g.dt / g.rho;
+    if (i >= 1 && i < g.nx) {
+        double dp_dx = (f.pp[cell3(g, k, j, i)] - f.pp[cell3(g, k, j, i - 1)]) / g.dx;
+        f.u[t] = f.us[t] - c * dp_dx;
+    } else if (i == 0 && g.bx0 == SIDE_OUTLET) {
+        double dp_dx = (f.pp[cell3(g, k, j, 0)] - 0.0) / (0.5 * g.dx);
+        f.u[t] = f.us[t] - c * dp_dx;
+    } else if (i == g.nx && g.bx1 == SIDE_OUTLET) {
+        double dp_dx = (0.0 - f.pp[cell3(g, k, j, g.nx - 1)]) / (0.5 * g.dx);
+        f.u[t] = f.us[t] - c * dp_dx;
+    }
+}
+
+extern "C" __global__ void e3_step_correct_v(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int nyp = g.ny + 1;
+    if (t >= g.nx * nyp * g.nz) return;
+    int i = t % g.nx; int j = (t / g.nx) % nyp; int k = t / (g.nx * nyp);
+    double c = g.dt / g.rho;
+    if (j >= 1 && j < g.ny) {
+        double dp_dy = (f.pp[cell3(g, k, j, i)] - f.pp[cell3(g, k, j - 1, i)]) / g.dy;
+        f.v[t] = f.vs[t] - c * dp_dy;
+    } else if (j == 0 && g.by0 == SIDE_OUTLET) {
+        double dp_dy = (f.pp[cell3(g, k, 0, i)] - 0.0) / (0.5 * g.dy);
+        f.v[t] = f.vs[t] - c * dp_dy;
+    } else if (j == g.ny && g.by1 == SIDE_OUTLET) {
+        double dp_dy = (0.0 - f.pp[cell3(g, k, g.ny - 1, i)]) / (0.5 * g.dy);
+        f.v[t] = f.vs[t] - c * dp_dy;
+    }
+}
+
+extern "C" __global__ void e3_step_correct_w(E3Params g, E3Ptrs f)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    if (t >= g.nx * g.ny * (g.nz + 1)) return;
+    int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
+    double c = g.dt / g.rho;
+    if (g.periodic_z) {
+        if (k == g.nz) return;
+        int below = k > 0 ? k - 1 : g.nz - 1;
+        double dp_dz = (f.pp[cell3(g, k, j, i)] - f.pp[cell3(g, below, j, i)]) / g.dz;
+        f.w[t] = f.ws[t] - c * dp_dz;
+        return;
+    }
+    if (k >= 1 && k < g.nz) {
+        double dp_dz = (f.pp[cell3(g, k, j, i)] - f.pp[cell3(g, k - 1, j, i)]) / g.dz;
+        f.w[t] = f.ws[t] - c * dp_dz;
+    } else if (k == 0 && g.bz0 == SIDE_OUTLET) {
+        double dp_dz = (f.pp[cell3(g, 0, j, i)] - 0.0) / (0.5 * g.dz);
+        f.w[t] = f.ws[t] - c * dp_dz;
+    } else if (k == g.nz && g.bz1 == SIDE_OUTLET) {
+        double dp_dz = (0.0 - f.pp[cell3(g, g.nz - 1, j, i)]) / (0.5 * g.dz);
+        f.w[t] = f.ws[t] - c * dp_dz;
+    }
+}
+
+/* p += p'; partial[block] = Σ |mass imbalance| of the corrected field. */
+extern "C" __global__ void e3_step_add_p_and_imbalance(E3Params g, E3Ptrs f, double* __restrict__ partial)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int total = g.nx * g.ny * g.nz;
+    double v = 0.0;
+    if (t < total) {
+        int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
+        f.p[t] += f.pp[t];
+        double divergence_flux = g.rho
+            * ((f.u[uf3(g, k, j, i + 1)] - f.u[uf3(g, k, j, i)]) * (g.dy * g.dz)
+               + (f.v[vf3(g, k, j + 1, i)] - f.v[vf3(g, k, j, i)]) * (g.dx * g.dz)
+               + (f.w[wf3(g, k + 1, j, i)] - f.w[wf3(g, k, j, i)]) * (g.dx * g.dy));
+        v = fabs(divergence_flux);
+    }
+    __shared__ double sh[256];
+    sh[threadIdx.x] = v;
+    __syncthreads();
+    for (int s = 128; s > 0; s >>= 1) { if (threadIdx.x < s) sh[threadIdx.x] += sh[threadIdx.x + s]; __syncthreads(); }
+    if (threadIdx.x == 0) partial[blockIdx.x] = sh[0];
+}
+
+/* Σ partial in index order (one thread) — the fixed-order final sum. */
+extern "C" __global__ void e3_step_reduce(int n, const double* __restrict__ partial, double* __restrict__ out)
+{
+    if (blockIdx.x * blockDim.x + threadIdx.x != 0) return;
+    double s = 0.0;
+    for (int i = 0; i < n; ++i) s += partial[i];
+    out[0] = s;
+}