diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu index be125a7..ac61d02 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/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 +__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; - cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &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; } 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 e28165e..74713c0 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 @@ -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,