embedded3 PERF-3 P3-1: the device predictor's interior fast path (bit 11, default ON; RTX_E3_PREDICT_FAST=0 = uniform path), the face update templated on the component + unrolled (392-byte local frame -> 64), exact-zero divisions/sqrt skipped on the fast path; all byte-identical (DFG 2D ny 122, slab flag ny 62, DFG 3D ny 62; device tests); predictor 22.4 -> 12.5 ms/step at DFG 3D ny 62

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu

@@ -73,7 +73,7 @@ __device__ __forceinline__ double cut_ap_at(const E3Params& g, const E3Cut& m, i
 /* The control volume's geometry of the face of component c at (i, j, k)
    (cutwall.rs cv_geometry): side apertures, the closing wall vector, the
    floored wall distance. */
-__device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, int k, int fidx,
+__device__ __forceinline__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, int k, int fidx,
                        double* alpha_out, double* apm, double* app, double* wall, double* distance_out)
 {
     double h[3] = { g.dx, g.dy, g.dz };
@@ -81,6 +81,7 @@ __device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, i
     double alpha = cut_ap(m, c)[fidx];
     int cm[3] = { i, j, k }; cm[c] -= 1;   /* cell minus */
     int cp[3] = { i, j, k };               /* cell plus */
+    #pragma unroll
     for (int d = 0; d < 3; ++d) {
         if (d == c) {
             int qm[3] = { i, j, k }; qm[c] -= 1;
@@ -96,6 +97,7 @@ __device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, i
             app[d] = 0.5 * (cut_ap_at(g, m, d, q1[0], q1[1], q1[2], 1.0) + cut_ap_at(g, m, d, q2[0], q2[1], q2[2], 1.0));
         }
     }
+    #pragma unroll
     for (int d = 0; d < 3; ++d) wall[d] = -(app[d] - apm[d]) * area[d];
     double h_min = fmin(fmin(g.dx, g.dy), g.dz);
     const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w);
@@ -115,12 +117,13 @@ __device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, i
    between the face (fidx, its geometry given) and its open neighbour fq at
    q along d — grad(u) . ds_perp with the cut faces' own wall-normal
    gradients (u - U_b)/d_f (cut_predictor.rs `transverse`). */
-__device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d, int fidx, int fq, const int* q,
+__device__ __forceinline__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d, int fidx, int fq, const int* q,
                                  double sign, double u0, double ub, double alpha, const double* wall, double distance,
                                  double uq, const double* sh, const double* ubt)
 {
     double ds[3] = { 0.0, 0.0, 0.0 };
     int any = 0;
+    #pragma unroll
     for (int e = 0; e < 3; ++e) {
         if (e == d) continue;
         ds[e] = sign * (sh[3 * fq + e] - sh[3 * fidx + e]);
@@ -140,19 +143,28 @@ __device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d
     int cut0 = a0 > 0.0 && alpha < 1.0, cutq = aq > 0.0 && alpha_q < 1.0;
     if (!cut0 && !cutq) return 0.0;
     double n0[3], nq[3];
+    #pragma unroll
     for (int e = 0; e < 3; ++e) { n0[e] = cut0 ? wall[e] / a0 : 0.0; nq[e] = cutq ? wall_q[e] / aq : 0.0; }
-    if (!cut0) for (int e = 0; e < 3; ++e) n0[e] = nq[e];
-    if (!cutq) for (int e = 0; e < 3; ++e) nq[e] = n0[e];
+    if (!cut0) {
+        #pragma unroll
+        for (int e = 0; e < 3; ++e) n0[e] = nq[e];
+    }
+    if (!cutq) {
+        #pragma unroll
+        for (int e = 0; e < 3; ++e) nq[e] = n0[e];
+    }
     int far = distance >= d_min || distance_q >= d_min;
     double gr[3] = { 0.0, 0.0, 0.0 };
     double count = 0.0;
     if (!((far && distance < d_min) || (!far && !cut0))) {
         double slope = (u0 - ub) / fmax(distance, d_min);
+        #pragma unroll
         for (int e = 0; e < 3; ++e) gr[e] -= slope * n0[e];
         count += 1.0;
     }
     if (!((far && distance_q < d_min) || (!far && !cutq))) {
         double slope = (uq - ubt[fq]) / fmax(distance_q, d_min);
+        #pragma unroll
         for (int e = 0; e < 3; ++e) gr[e] -= slope * nq[e];
         count += 1.0;
     }
@@ -161,8 +173,10 @@ __device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d
 }
 
 /* The predicted value of the open face of component c at (i, j, k). */
-__device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cut& m, int c, int i, int j, int k, int fidx)
+template <int C>
+__device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cut& m, int i, int j, int k, int fidx)
 {
+    const int c = C;   /* compile-time: every [c] index is a register, not local memory (P3-1) */
     double h[3] = { g.dx, g.dy, g.dz };
     double n[3] = { (double)g.nx, (double)g.ny, (double)g.nz };
     double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy };
@@ -171,18 +185,34 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
     int p[3] = { i, j, k };
     const double* old_c = cut_old(f, c);
     double u0 = old_c[fidx];
+    /* PERF-3 P3-1 (bit 11 of wall_order): an INTERIOR face — further than the imposition band
+       (4 h_min) from the wall, so its own and its neighbours' apertures are 1, the wall vector 0,
+       the centroid shifts 0 and the surface velocity 0 — skips the geometry loads; every skipped
+       term is an exact 0.0 or a factor 1.0 in the general path, so the result is bit-identical. */
+    const double* dist_c = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w);
+    double h_min_f = fmin(fmin(g.dx, g.dy), g.dz);
+    int fast = (g.wall_order & 2048) && fabs(dist_c[fidx]) > 4.0 * h_min_f;
     /* the control volume's geometry (cutwall.rs cv_geometry) */
     double alpha, apm[3], app[3], wall[3], distance;
+    if (fast) {
+        alpha = 1.0;
+        #pragma unroll
+        for (int d = 0; d < 3; ++d) { apm[d] = 1.0; app[d] = 1.0; wall[d] = 0.0; }
+        distance = dist_c[fidx];
+    } else {
         cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &distance);
+    }
     int cm[3] = { i, j, k }; cm[c] -= 1;   /* cell minus */
     int cp[3] = { i, j, k };               /* cell plus */
     const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
-    double ub = ubt[fidx];
+    double ub = fast ? 0.0 : ubt[fidx];
     /* face position */
     double x[3];
+    #pragma unroll
     for (int d = 0; d < 3; ++d) x[d] = (p[d] + (d == c ? 0.0 : 0.5)) * h[d];
 
     double mass_out = 0.0, conv = 0.0, diff = 0.0, wall_implicit = 0.0, wall_rhs = 0.0;
+    #pragma unroll
     for (int d = 0; d < 3; ++d) {
         double a_d = area[d];
         int q[3];
@@ -197,8 +227,8 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         double dn2 = f_dn2 >= 0 ? old_c[f_dn2] : 0.0;
         double m_plus, m_minus;
         if (d == c) {
-            double f_up = (f_up1 >= 0 ? cut_ap(m, c)[f_up1] : alpha) * (f_up1 >= 0 ? up1 : u0);
-            double f_dn = (f_dn1 >= 0 ? cut_ap(m, c)[f_dn1] : alpha) * (f_dn1 >= 0 ? dn1 : u0);
+            double f_up = (f_up1 >= 0 ? (fast ? 1.0 : cut_ap(m, c)[f_up1]) : alpha) * (f_up1 >= 0 ? up1 : u0);
+            double f_dn = (f_dn1 >= 0 ? (fast ? 1.0 : cut_ap(m, c)[f_dn1]) : alpha) * (f_dn1 >= 0 ? dn1 : u0);
             double f0 = alpha * u0;
             m_plus = 0.5 * (f0 + f_up) * a_d;
             m_minus = 0.5 * (f_dn + f0) * a_d;
@@ -210,10 +240,10 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
             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 ? cut_ap(m, d)[fa] * old_d[fa] : 0.0;
-            double flb = fb >= 0 ? cut_ap(m, d)[fb] * old_d[fb] : 0.0;
-            double flc = fc >= 0 ? cut_ap(m, d)[fc] * old_d[fc] : 0.0;
-            double fld = fd >= 0 ? cut_ap(m, d)[fd] * old_d[fd] : 0.0;
+            double fla = fa >= 0 ? (fast ? 1.0 : cut_ap(m, d)[fa]) * old_d[fa] : 0.0;
+            double flb = fb >= 0 ? (fast ? 1.0 : cut_ap(m, d)[fb]) * old_d[fb] : 0.0;
+            double flc = fc >= 0 ? (fast ? 1.0 : cut_ap(m, d)[fc]) * old_d[fc] : 0.0;
+            double fld = fd >= 0 ? (fast ? 1.0 : cut_ap(m, d)[fd]) * old_d[fd] : 0.0;
             m_plus = 0.5 * (fla + flb) * a_d;
             m_minus = 0.5 * (flc + fld) * a_d;
         }
@@ -237,8 +267,8 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         double u_plus, delta_plus;
         /* bit 10 of wall_order (A2): a prescribed neighbour face exchanges no convective momentum */
         const int* opn = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
-        int presc_up = (g.wall_order & 1024) && f_up1 >= 0 && !opn[f_up1];
-        int presc_dn = (g.wall_order & 1024) && f_dn1 >= 0 && !opn[f_dn1];
+        int presc_up = !fast && (g.wall_order & 1024) && f_up1 >= 0 && !opn[f_up1];
+        int presc_dn = !fast && (g.wall_order & 1024) && f_dn1 >= 0 && !opn[f_dn1];
         if (presc_up) { u_plus = u0; delta_plus = 0.0; }
         else if (f_up1 >= 0) {
             if (g.scheme == SCHEME_UPWIND) delta_plus = 0.0;
@@ -260,21 +290,21 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         /* bit 5 of wall_order: the exchange with a SOLID neighbour over the axis
            distance to the wall, implicit (S2-5) */
         double delta_x = h[d];
-        if (g.wall_order & 32) {
+        if (!fast && (g.wall_order & 32)) {
             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 solid_up = (g.wall_order & 32) && f_up1 >= 0 && cut_ap(m, c)[f_up1] == 0.0;
-        int solid_dn = (g.wall_order & 32) && f_dn1 >= 0 && cut_ap(m, c)[f_dn1] == 0.0;
+        int solid_up = !fast && (g.wall_order & 32) && f_up1 >= 0 && cut_ap(m, c)[f_up1] == 0.0;
+        int solid_dn = !fast && (g.wall_order & 32) && f_dn1 >= 0 && cut_ap(m, c)[f_dn1] == 0.0;
         /* bit 6 of wall_order: cross diffusion over the open-part centroid spacing,
            the part beyond 1/h point-implicit (S2-5) */
         int centroid = (g.wall_order & 64) && d != c;
         const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w);
         /* bit 7 of wall_order: the transverse centroid correction, cross and own direction (S2-6) */
-        if ((g.wall_order & 128) && (g.wall_order & 64)) {
+        if (!fast && (g.wall_order & 128) && (g.wall_order & 64)) {
             if (f_up1 >= 0 && !solid_up && cut_ap(m, c)[f_up1] > 0.0) {
                 int qn[3] = { i, j, k }; qn[d] += 1;
                 double dl = h[d];
@@ -289,6 +319,8 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
             }
         }
         if (solid_up) { double kx = mu * g_plus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * up1; }
+        /* fast: dl == h exactly, kx == +0.0 exactly, so the else branch below is (X + 0.0) * (up1 - u0) == X * (up1 - u0) */
+        else if (fast && f_up1 >= 0 && centroid) diff += (mu * g_plus * a_d / h[d]) * (up1 - u0);
         else if (f_up1 >= 0 && centroid) {
             double dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]);
             dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
@@ -299,6 +331,7 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         else if (f_up1 >= 0) diff += mu * g_plus * a_d * (up1 - u0) / h[d];
         else if (sides[d][1] == SIDE_VELOCITY) diff += mu * g_plus * a_d * (beyond_p - u0) / (0.5 * h[d]);
         if (solid_dn) { double kx = mu * g_minus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * dn1; }
+        else if (fast && f_dn1 >= 0 && centroid) diff -= (mu * g_minus * a_d / h[d]) * (u0 - dn1);
         else if (f_dn1 >= 0 && centroid) {
             double dl = h[d] - (sh[3 * f_dn1 + d] - sh[3 * fidx + d]);
             dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
@@ -319,13 +352,14 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
     double v_u = alpha * h[c] * area[c];
     const double* src = c == 0 ? f.su : (c == 1 ? f.sv : f.sw);
     double source = src[fidx] * v_u;
-    double a_w = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
+    /* fast: no wall vector, so a_w == 0.0 and shear == mu * 0.0 * c1 == 0.0 for any finite c1 */
+    double a_w = fast ? 0.0 : sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
     /* The wall gradient: one-point (order 1) or quadratic through the next
        open face away from the body along the wall normal's dominant axis
        (order 2; the neighbour's old value explicit). */
-    double c1 = 1.0 / distance, shear_explicit = 0.0;
+    double c1 = fast ? 0.0 : 1.0 / distance, shear_explicit = 0.0;
     /* bit 9 of wall_order: the advancing-wall closure on the one-point coefficient (A3-i) */
-    if (g.wall_order & 512) {
+    if (!fast && (g.wall_order & 512)) {
         const double* vnt = c == 0 ? m.vn_u : (c == 1 ? m.vn_v : m.vn_w);
         c1 = adv_factor(vnt[fidx] * distance / g.nu) / distance;
     }
@@ -363,7 +397,8 @@ extern "C" __global__ void e3_cut_predict(E3Params g, E3Ptrs f, E3Cut m, int c)
     if (c == 2) { if (g.periodic_z) { if (k == g.nz) return; } else if (k == 0 || k == g.nz) return; }
     const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
     if (!open[t]) return;
-    double val = cut_face_update(g, f, m, c, i, j, k, t);
+    double val = c == 0 ? cut_face_update<0>(g, f, m, i, j, k, t)
+               : (c == 1 ? cut_face_update<1>(g, f, m, i, j, k, t) : cut_face_update<2>(g, f, m, i, j, k, t));
     double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
     out[t] = val;
 }

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs

@@ -306,7 +306,11 @@ impl DeviceStep {
                 // bit 9: the advancing-wall friction closure (A3-i).
                 + 512 * i32::from(self.solver.params.wall_advancing)
                 // bit 10: no convective exchange with prescribed faces (A2).
-                + 1024 * i32::from(self.solver.params.exchange_convection_off),
+                + 1024 * i32::from(self.solver.params.exchange_convection_off)
+                // bit 11: the predictor's interior fast path (PERF-3 P3-1; bit-identical
+                // by construction and gated so; default ON, `RTX_E3_PREDICT_FAST=0` restores
+                // the uniform path).
+                + 2048 * i32::from(!std::env::var("RTX_E3_PREDICT_FAST").is_ok_and(|v| v == "0")),
             dx: g.dx,
             dy: g.dy,
             dz: g.dz,