embedded3 item 9b: the device step carries a static cut-cell mask (e3_cut.cu: cut predictor, apertured merged continuity with the fold, owner-read corrections; device CG off-stencil links) — host = device to 2e-10 (CFD1 cylinder nz 4) and 4e-14 (sphere) under tight tolerances

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

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

@@ -7,9 +7,12 @@
  */
 #define NONE 0xFFFFFFFFu
 
+/* The seven-point neighbour sum plus the cell's off-stencil links (CSR:
+ * link_ptr[g] .. link_ptr[g + 1]; the virtually merged small cells). */
 __device__ __forceinline__ double nb_sum3d(
     int g, int nx, const double* ae, const double* aw, const double* an, const double* as_,
     const double* at, const double* ab, const unsigned int* top, const unsigned int* bot,
+    const unsigned int* link_ptr, const unsigned int* link_idx, const double* link_coef,
     const double* x)
 {
     double s = 0.0;
@@ -19,6 +22,7 @@ __device__ __forceinline__ double nb_sum3d(
     double ss = as_[g]; if (ss != 0.0) s += ss * x[g - nx];
     double t = at[g]; if (t != 0.0) s += t * x[top[g]];
     double b = ab[g]; if (b != 0.0) s += b * x[bot[g]];
+    for (unsigned int l = link_ptr[g]; l < link_ptr[g + 1]; ++l) s += link_coef[l] * x[link_idx[l]];
     return s;
 }
 
@@ -42,12 +46,14 @@ extern "C" __global__ void e3_cg_spmv(
     const double* __restrict__ an, const double* __restrict__ as_,
     const double* __restrict__ at, const double* __restrict__ ab,
     const unsigned int* __restrict__ top, const unsigned int* __restrict__ bot,
+    const unsigned int* __restrict__ link_ptr, const unsigned int* __restrict__ link_idx,
+    const double* __restrict__ link_coef,
     const double* __restrict__ ap, const double* __restrict__ d, double* __restrict__ q, int nx)
 {
     int t = blockIdx.x * blockDim.x + threadIdx.x;
     if (t >= n_cells) return;
     int g = cells[t];
-    q[g] = ap[g] * d[g] - nb_sum3d(g, nx, ae, aw, an, as_, at, ab, top, bot, d);
+    q[g] = ap[g] * d[g] - nb_sum3d(g, nx, ae, aw, an, as_, at, ab, top, bot, link_ptr, link_idx, link_coef, d);
 }
 
 /* r = b − A p on the cells; partial[block] = Σ |r| over the block's cells. */
@@ -57,6 +63,8 @@ extern "C" __global__ void e3_cg_residual(
     const double* __restrict__ an, const double* __restrict__ as_,
     const double* __restrict__ at, const double* __restrict__ ab,
     const unsigned int* __restrict__ top, const unsigned int* __restrict__ bot,
+    const unsigned int* __restrict__ link_ptr, const unsigned int* __restrict__ link_idx,
+    const double* __restrict__ link_coef,
     const double* __restrict__ ap, const double* __restrict__ b,
     const double* __restrict__ p, double* __restrict__ r, double* __restrict__ partial, int nx)
 {
@@ -64,7 +72,7 @@ extern "C" __global__ void e3_cg_residual(
     double v = 0.0;
     if (t < n_cells) {
         int g = cells[t];
-        v = b[g] - (ap[g] * p[g] - nb_sum3d(g, nx, ae, aw, an, as_, at, ab, top, bot, p));
+        v = b[g] - (ap[g] * p[g] - nb_sum3d(g, nx, ae, aw, an, as_, at, ab, top, bot, link_ptr, link_idx, link_coef, p));
         r[g] = v;
         v = fabs(v);
     }

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

@@ -0,0 +1,309 @@
+/**
+ * embedded3 item 9b: the apertured cut-cell wall's maps on the device — the
+ * host `step/cut_predictor.rs` and `cutwall.rs` expression for expression
+ * (FMA contraction off): the predictor on the open faces of one component
+ * (the momentum volume `V_u = α h A`, mass fluxes averaged from the two
+ * cells, upwind + TVD, apertured diffusion, the projection's pressure
+ * force, the net mass flux times the face's value, the implicit wall
+ * shear, the inertia floor), the apertured continuity with the wall flux
+ * and the fold of the virtually merged small cells into their masters,
+ * the corrections on the open faces reading the owner's p', the pressure
+ * update. Static bodies (the tables are built once).
+ *
+ * Appended after `e3_step.cu` at load (shares E3Params / E3Ptrs and the
+ * helpers there).
+ */
+#define CUT_INERTIA_FLOOR 0.1
+#define CUT_DISTANCE_FLOOR 0.05
+
+struct E3Cut {
+    const double *a_u, *a_v, *a_w;      /* apertures per face */
+    const double *d_u, *d_v, *d_w;      /* φ at the face centres */
+    const double *ub_u, *ub_v, *ub_w;   /* surface velocity component at the foot per face */
+    const double *wall_flux;            /* per cell, compatible */
+    const int *open_u, *open_v, *open_w;/* unknown faces */
+    const int *active;                  /* cells with an equation */
+    const unsigned int *owner;          /* the master of a merged cell (itself otherwise) */
+    const unsigned int *fold_ptr, *fold_idx; /* CSR: the slaves of every cell */
+    double *cell_flux;                  /* scratch per cell */
+};
+
+/* Face index of component c at lattice (i, j, k); −1 outside (z wraps when periodic). */
+__device__ __forceinline__ int cut_face(const E3Params& g, int c, int i, int j, int k)
+{
+    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 (i < 0 || i >= ni || j < 0 || j >= nj) return -1;
+    if (g.periodic_z) { k = ((k % g.nz) + g.nz) % g.nz; }
+    else if (k < 0 || k >= nk) return -1;
+    if (c == 0) return (k * g.ny + j) * (g.nx + 1) + i;
+    if (c == 1) return (k * (g.ny + 1) + j) * g.nx + i;
+    return (k * g.ny + j) * g.nx + i;
+}
+
+__device__ __forceinline__ int cut_cell(const E3Params& g, int i, int j, int k)
+{
+    if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return -1;
+    if (g.periodic_z) { k = ((k % g.nz) + g.nz) % g.nz; }
+    else if (k < 0 || k >= g.nz) return -1;
+    return (k * g.ny + j) * g.nx + i;
+}
+
+__device__ __forceinline__ const double* cut_ap(const E3Cut& m, int c) { return c == 0 ? m.a_u : (c == 1 ? m.a_v : m.a_w); }
+__device__ __forceinline__ const double* cut_old(const E3Ptrs& f, int c) { return c == 0 ? f.uo : (c == 1 ? f.vo : f.wo); }
+
+/* Aperture of component `cc` at the lattice point, or `dflt` outside. */
+__device__ __forceinline__ double cut_ap_at(const E3Params& g, const E3Cut& m, int cc, int i, int j, int k, double dflt)
+{
+    int fidx = cut_face(g, cc, i, j, k);
+    return fidx < 0 ? dflt : cut_ap(m, cc)[fidx];
+}
+
+/* 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)
+{
+    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 };
+    double rho = g.rho, mu = g.nu * g.rho;
+    int sides[3][2] = { { g.bx0, g.bx1 }, { g.by0, g.by1 }, { g.bz0, g.bz1 } };
+    int p[3] = { i, j, k };
+    const double* old_c = cut_old(f, c);
+    double u0 = old_c[fidx];
+    /* the control volume's geometry (cutwall.rs cv_geometry) */
+    double alpha = cut_ap(m, c)[fidx];
+    double apm[3], app[3];
+    int cm[3] = { i, j, k }; cm[c] -= 1;   /* cell minus */
+    int cp[3] = { i, j, k };               /* cell plus */
+    for (int d = 0; d < 3; ++d) {
+        if (d == c) {
+            int qm[3] = { i, j, k }; qm[c] -= 1;
+            int qp[3] = { i, j, k }; qp[c] += 1;
+            double am = cut_ap_at(g, m, c, qm[0], qm[1], qm[2], alpha);
+            double apl = cut_ap_at(g, m, c, qp[0], qp[1], qp[2], alpha);
+            apm[d] = 0.5 * (am + alpha);
+            app[d] = 0.5 * (alpha + apl);
+        } else {
+            int q1[3] = { cm[0], cm[1], cm[2] }; q1[d] += 1;
+            int q2[3] = { cp[0], cp[1], cp[2] }; q2[d] += 1;
+            apm[d] = 0.5 * (cut_ap_at(g, m, d, cm[0], cm[1], cm[2], 1.0) + cut_ap_at(g, m, d, cp[0], cp[1], cp[2], 1.0));
+            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));
+        }
+    }
+    double wall[3];
+    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);
+    double distance = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha), CUT_DISTANCE_FLOOR * h_min);
+    const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
+    double ub = ubt[fidx];
+    /* face position */
+    double x[3];
+    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;
+    for (int d = 0; d < 3; ++d) {
+        double a_d = area[d];
+        int q[3];
+        /* neighbouring faces of this component along d */
+        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] += 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] -= 1;  int f_dn1 = 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 up1 = f_up1 >= 0 ? old_c[f_up1] : 0.0;
+        double up2 = f_up2 >= 0 ? old_c[f_up2] : 0.0;
+        double dn1 = f_dn1 >= 0 ? old_c[f_dn1] : 0.0;
+        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 f0 = alpha * u0;
+            m_plus = 0.5 * (f0 + f_up) * a_d;
+            m_minus = 0.5 * (f_dn + f0) * a_d;
+        } else {
+            const double* old_d = cut_old(f, 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 ? 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;
+            m_plus = 0.5 * (fla + flb) * a_d;
+            m_minus = 0.5 * (flc + fld) * a_d;
+        }
+        mass_out += m_plus - m_minus;
+        /* beyond a domain side along d */
+        double beyond_p = u0, beyond_m = u0;
+        if (sides[d][1] == SIDE_VELOCITY || sides[d][0] == SIDE_VELOCITY) {
+            /* the boundary tables: [side][component] at the face's transverse position */
+            const double* tp = 0; const double* tm = 0; int idx_p = 0, idx_m = 0;
+            if (d == 0) { tp = c == 0 ? f.bx1u : (c == 1 ? f.bx1v : f.bx1w); tm = c == 0 ? f.bx0u : (c == 1 ? f.bx0v : f.bx0w);
+                          int nj = c == 1 ? g.ny + 1 : g.ny; idx_p = idx_m = k * nj + j; }
+            else if (d == 1) { tp = c == 0 ? f.by1u : (c == 1 ? f.by1v : f.by1w); tm = c == 0 ? f.by0u : (c == 1 ? f.by0v : f.by0w);
+                          int ni = c == 0 ? g.nx + 1 : g.nx; idx_p = idx_m = k * ni + i; }
+            else { tp = c == 0 ? f.bz1u : (c == 1 ? f.bz1v : f.bz1w); tm = c == 0 ? f.bz0u : (c == 1 ? f.bz0v : f.bz0w);
+                          int ni = c == 0 ? g.nx + 1 : g.nx; idx_p = idx_m = j * ni + i; }
+            if (sides[d][1] == SIDE_VELOCITY) beyond_p = tp[idx_p];
+            if (sides[d][0] == SIDE_VELOCITY) beyond_m = tm[idx_m];
+        }
+        (void)n; (void)x;
+        /* convection through the plus face */
+        double u_plus, delta_plus;
+        if (f_up1 >= 0) {
+            if (g.scheme == SCHEME_UPWIND) delta_plus = 0.0;
+            else if (m_plus >= 0.0) delta_plus = face_corr3(g.scheme, f_dn1 >= 0, dn1, u0, up1);
+            else delta_plus = face_corr3(g.scheme, f_up2 >= 0, up2, up1, u0);
+            u_plus = upwind3(m_plus, u0, up1);
+        } else { u_plus = upwind3(m_plus, u0, beyond_p); delta_plus = 0.0; }
+        double u_minus, delta_minus;
+        if (f_dn1 >= 0) {
+            if (g.scheme == SCHEME_UPWIND) delta_minus = 0.0;
+            else if (m_minus >= 0.0) delta_minus = face_corr3(g.scheme, f_dn2 >= 0, dn2, dn1, u0);
+            else delta_minus = face_corr3(g.scheme, f_up1 >= 0, up1, u0, dn1);
+            u_minus = upwind3(m_minus, dn1, u0);
+        } else { u_minus = upwind3(m_minus, beyond_m, u0); delta_minus = 0.0; }
+        conv += m_plus * (u_plus + delta_plus) - m_minus * (u_minus + delta_minus);
+        /* diffusion */
+        double g_minus = apm[d], g_plus = app[d];
+        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 (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d];
+        else if (sides[d][0] == SIDE_VELOCITY) diff -= mu * g_minus * a_d * (u0 - beyond_m) / (0.5 * h[d]);
+    }
+    conv -= mass_out * u0;
+    int cpi = cut_cell(g, cp[0], cp[1], cp[2]);
+    int cmi = cut_cell(g, cm[0], cm[1], cm[2]);
+    double p_plus = cpi >= 0 ? f.p[cpi] : 0.0;
+    double p_minus = cmi >= 0 ? f.p[cmi] : 0.0;
+    double pressure = -(p_plus - p_minus) * alpha * area[c];
+    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]);
+    double shear = mu * a_w / distance;
+    double v_eff = fmax(alpha, CUT_INERTIA_FLOOR) * h[c] * area[c];
+    double inertia = rho * v_eff / g.dt;
+    return (inertia * u0 - conv + diff + pressure + source + shear * ub) / (inertia + shear);
+}
+
+/* The predictor on the open interior faces of component c. */
+extern "C" __global__ void e3_cut_predict(E3Params g, E3Ptrs f, E3Cut m, int c)
+{
+    int t = 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 >= ni * nj * nk) return;
+    int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj);
+    if (c == 0 && (i == 0 || i == g.nx)) return;
+    if (c == 1 && (j == 0 || j == g.ny)) return;
+    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* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
+    out[t] = val;
+}
+
+/* cell_flux = ρ (Σ apertured flux out + wall flux) on the active cells (0 elsewhere). */
+extern "C" __global__ void e3_cut_divergence(E3Params g, E3Ptrs f, E3Cut m, int starred)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    int total = g.nx * g.ny * g.nz;
+    if (t >= total) return;
+    if (!m.active[t]) { m.cell_flux[t] = 0.0; return; }
+    int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
+    const double* u = starred ? f.us : f.u;
+    const double* v = starred ? f.vs : f.v;
+    const double* w = starred ? f.ws : f.w;
+    int ue = uf3(g, k, j, i + 1), uw = uf3(g, k, j, i);
+    int vn = vf3(g, k, j + 1, i), vs = vf3(g, k, j, i);
+    int wt = wf3(g, k + 1, j, i), wb = wf3(g, k, j, i);
+    double divergence_flux = g.rho
+        * ((m.a_u[ue] * u[ue] - m.a_u[uw] * u[uw]) * (g.dy * g.dz)
+           + (m.a_v[vn] * v[vn] - m.a_v[vs] * v[vs]) * (g.dx * g.dz)
+           + (m.a_w[wt] * w[wt] - m.a_w[wb] * w[wb]) * (g.dx * g.dy));
+    divergence_flux += g.rho * m.wall_flux[t];
+    m.cell_flux[t] = divergence_flux;
+}
+
+/* The merged fluxes: sp = −(own + slaves) on the owning cells when `write_sp`;
+ * partial[block] = Σ |merged flux| (the source scale / the mass imbalance). */
+extern "C" __global__ void e3_cut_fold(E3Params g, E3Ptrs f, E3Cut m, int write_sp, 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) {
+        if (m.active[t] && m.owner[t] == (unsigned int)t) {
+            double s = m.cell_flux[t];
+            for (unsigned int l = m.fold_ptr[t]; l < m.fold_ptr[t + 1]; ++l) s += m.cell_flux[m.fold_idx[l]];
+            if (write_sp) f.sp[t] = -s;
+            v = fabs(s);
+        } else if (write_sp) {
+            f.sp[t] = 0.0;
+        }
+    }
+    __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 on the open interior faces of component c from the owners' p'. */
+extern "C" __global__ void e3_cut_correct(E3Params g, E3Ptrs f, E3Cut m, int c)
+{
+    int t = 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 >= ni * nj * nk) return;
+    int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj);
+    double cc = g.dt / g.rho;
+    const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
+    double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
+    const double* star = c == 0 ? f.us : (c == 1 ? f.vs : f.ws);
+    int interior = c == 0 ? (i >= 1 && i < g.nx) : (c == 1 ? (j >= 1 && j < g.ny)
+                 : (g.periodic_z ? (k < g.nz) : (k >= 1 && k < g.nz)));
+    if (interior) {
+        if (!open[t]) return;
+        int cp[3] = { i, j, k }, cm[3] = { i, j, k }; cm[c] -= 1;
+        int a = cut_cell(g, cp[0], cp[1], cp[2]);
+        int b = cut_cell(g, cm[0], cm[1], cm[2]);
+        double h = c == 0 ? g.dx : (c == 1 ? g.dy : g.dz);
+        double dp = (f.pp[m.owner[a]] - f.pp[m.owner[b]]) / h;
+        out[t] = star[t] - cc * dp;
+        return;
+    }
+    /* outlet faces against 0 outside (no body reaches an outlet) */
+    int s0 = c == 0 ? g.bx0 : (c == 1 ? g.by0 : g.bz0);
+    int s1 = c == 0 ? g.bx1 : (c == 1 ? g.by1 : g.bz1);
+    int own = c == 0 ? i : (c == 1 ? j : k);
+    int nn = c == 0 ? g.nx : (c == 1 ? g.ny : g.nz);
+    double h = c == 0 ? g.dx : (c == 1 ? g.dy : g.dz);
+    if (own == 0 && s0 == SIDE_OUTLET) {
+        int cp[3] = { i, j, k };
+        int a = cut_cell(g, cp[0], cp[1], cp[2]);
+        out[t] = star[t] - cc * (f.pp[a] - 0.0) / (0.5 * h);
+    } else if (own == nn && s1 == SIDE_OUTLET) {
+        int cm[3] = { i, j, k }; cm[c] -= 1;
+        int b = cut_cell(g, cm[0], cm[1], cm[2]);
+        out[t] = star[t] - cc * (0.0 - f.pp[b]) / (0.5 * h);
+    }
+}
+
+/* p += p'(owner) on the active cells. */
+extern "C" __global__ void e3_cut_add_p(E3Params g, E3Ptrs f, E3Cut m)
+{
+    int t = blockIdx.x * blockDim.x + threadIdx.x;
+    if (t >= g.nx * g.ny * g.nz) return;
+    if (!m.active[t]) return;
+    f.p[t] += f.pp[m.owner[t]];
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs

@@ -70,6 +70,10 @@ pub struct DeviceCg {
     at: CudaSlice<f64>,
     ab: CudaSlice<f64>,
     ap: CudaSlice<f64>,
+    /// Off-stencil links per cell (CSR; the merged small cells).
+    link_ptr: CudaSlice<u32>,
+    link_idx: CudaSlice<u32>,
+    link_coef: CudaSlice<f64>,
     b: CudaSlice<f64>,
     r: CudaSlice<f64>,
     z: CudaSlice<f64>,
@@ -118,6 +122,18 @@ impl DeviceCg {
         };
         let up_f = |v: &[f64]| -> CudaSlice<f64> { rt.stream.memcpy_stod(v).expect("upload") };
         let zeros = || rt.stream.alloc_zeros::<f64>(n).expect("alloc");
+        let lists = problem.link_lists();
+        let mut link_ptr = Vec::with_capacity(n + 1);
+        let mut link_idx = Vec::new();
+        let mut link_coef = Vec::new();
+        link_ptr.push(0u32);
+        for list in &lists {
+            for &(other, c) in list {
+                link_idx.push(other as u32);
+                link_coef.push(c);
+            }
+            link_ptr.push(link_idx.len() as u32);
+        }
         Self {
             key: OperatorKey::of(problem, params),
             n,
@@ -134,6 +150,13 @@ impl DeviceCg {
             at: up_f(&fine.at),
             ab: up_f(&fine.ab),
             ap: up_f(&fine.ap),
+            link_ptr: up_u(&link_ptr),
+            link_idx: up_u(&link_idx),
+            link_coef: up_f(if link_coef.is_empty() {
+                &[0.0]
+            } else {
+                &link_coef
+            }),
             b: zeros(),
             r: zeros(),
             z: zeros(),
@@ -324,6 +347,9 @@ impl DeviceCg {
                 .arg(&self.ab)
                 .arg(&self.top)
                 .arg(&self.bot)
+                .arg(&self.link_ptr)
+                .arg(&self.link_idx)
+                .arg(&self.link_coef)
                 .arg(&self.ap)
                 .arg(&self.b)
                 .arg(p)
@@ -481,6 +507,9 @@ impl DeviceCg {
                     .arg(&self.ab)
                     .arg(&self.top)
                     .arg(&self.bot)
+                    .arg(&self.link_ptr)
+                    .arg(&self.link_idx)
+                    .arg(&self.link_coef)
                     .arg(&self.ap)
                     .arg(&self.d)
                     .arg(&mut self.q)

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs

@@ -208,7 +208,14 @@ impl Problem {
             }
         }
         for &(a, b, c) in &self.links {
-            if a >= n || b >= n || a == b || !self.active[a] || !self.active[b] || c <= 0.0 || c.is_nan() {
+            if a >= n
+                || b >= n
+                || a == b
+                || !self.active[a]
+                || !self.active[b]
+                || c <= 0.0
+                || c.is_nan()
+            {
                 return Err(format!("invalid link ({a}, {b}, {c})"));
             }
         }

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

@@ -5,6 +5,8 @@
 //! fields into a `Field` at instants. Compiled with FMA contraction
 //! off so the predictors are the host's arithmetic to the bit.
 
+mod cut;
+
 use super::{Side, Solver, StepResult};
 use crate::solvers::incompressible::embedded3::Grid;
 use crate::solvers::incompressible::embedded3::field::Field;
@@ -150,6 +152,8 @@ pub struct DeviceStep {
     cg_dt: f64,
     timers: Option<StepTimers>,
     initialized: bool,
+    /// A static cut-cell mask's tables (item 9b), when the solver has one.
+    cut: Option<cut::DeviceCut>,
 }
 
 impl DeviceStep {
@@ -177,6 +181,7 @@ impl DeviceStep {
         let timers = std::env::var("RTX_PROFILE")
             .is_ok()
             .then(StepTimers::default);
+        let cut = cut::DeviceCut::build(&solver, grid);
         Self {
             solver,
             grid,
@@ -203,9 +208,15 @@ impl DeviceStep {
             cg_dt: 0.0,
             timers,
             initialized: false,
+            cut,
         }
     }
 
+    /// The number of virtually merged cells on the device mask (0 without one).
+    pub fn merged_cells(&self) -> usize {
+        self.cut.as_ref().map_or(0, |c| c.merged)
+    }
+
     pub fn timers(&self) -> Option<StepTimers> {
         self.timers
     }
@@ -385,6 +396,9 @@ impl DeviceStep {
         if !self.initialized {
             self.initialize();
         }
+        if self.cut.is_some() {
+            return self.advance_cut(dt);
+        }
         let rt = runtime();
         let k = kernels();
         let g = self.grid;

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

@@ -0,0 +1,380 @@
+//! Item 9b: the device step carrying a static cut-cell mask — the mask's
+//! tables uploaded once (apertures, face distances, surface velocities at
+//! the feet, the compatible wall flux, the open faces, the active cells,
+//! the merged cells' owners and the fold lists) and the step's cut
+//! kernels (`e3_cut.cu`, appended to `e3_step.cu` at load).
+
+use super::{DeviceStep, E3Params, E3Ptrs, StepResult};
+use crate::solvers::incompressible::embedded3::Grid;
+use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime};
+use crate::solvers::incompressible::embedded3::poisson::device_cg::DeviceCg;
+use crate::solvers::incompressible::embedded3::step::Solver;
+use crate::solvers::incompressible::poisson::MultigridParameters;
+use cudarc::driver::{
+    CudaFunction, CudaModule, CudaSlice, DevicePtr, DeviceRepr, PushKernelArg, ValidAsZeroBits,
+};
+use std::sync::{Arc, OnceLock};
+use std::time::Instant;
+
+const CUT_KERNELS: &str = concat!(
+    include_str!("../../../../../kernels/cuda/e3_step.cu"),
+    include_str!("../../../../../kernels/cuda/e3_cut.cu")
+);
+
+struct CutKernels {
+    _module: Arc<CudaModule>,
+    predict: CudaFunction,
+    divergence: CudaFunction,
+    fold: CudaFunction,
+    correct: CudaFunction,
+    add_p: CudaFunction,
+}
+
+static CUT_KERNELS_ONCE: OnceLock<CutKernels> = OnceLock::new();
+
+fn cut_kernels() -> &'static CutKernels {
+    CUT_KERNELS_ONCE.get_or_init(|| {
+        let module = load_module(CUT_KERNELS, "e3_cut.cu", true);
+        let f = |name: &str| module.load_function(name).expect(name);
+        CutKernels {
+            predict: f("e3_cut_predict"),
+            divergence: f("e3_cut_divergence"),
+            fold: f("e3_cut_fold"),
+            correct: f("e3_cut_correct"),
+            add_p: f("e3_cut_add_p"),
+            _module: module,
+        }
+    })
+}
+
+/// `struct E3Cut` in e3_cut.cu: 18 device pointers.
+#[repr(C)]
+#[derive(Clone, Copy)]
+struct E3CutPtrs {
+    ptrs: [u64; 18],
+}
+unsafe impl DeviceRepr for E3CutPtrs {}
+unsafe impl ValidAsZeroBits for E3CutPtrs {}
+
+/// The static cut-cell mask on the device.
+pub(super) struct DeviceCut {
+    a: [CudaSlice<f64>; 3],
+    d: [CudaSlice<f64>; 3],
+    ub: [CudaSlice<f64>; 3],
+    wall_flux: CudaSlice<f64>,
+    open: [CudaSlice<i32>; 3],
+    active: CudaSlice<i32>,
+    owner: CudaSlice<u32>,
+    fold_ptr: CudaSlice<u32>,
+    fold_idx: CudaSlice<u32>,
+    cell_flux: CudaSlice<f64>,
+    pub(super) merged: usize,
+}
+
+impl DeviceCut {
+    /// The tables of the solver's cut mask (`None` without one).
+    pub(super) fn build(solver: &Solver, g: Grid) -> Option<Self> {
+        let mask = solver.mask()?;
+        let cut = mask.cut()?;
+        let body = solver.body()?;
+        let rt = runtime();
+        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
+        let h = [g.dx, g.dy, g.dz];
+        let counts = [(nx + 1) * ny * nz, nx * (ny + 1) * nz, nx * ny * (nz + 1)];
+        let up_f = |v: &[f64]| -> CudaSlice<f64> {
+            rt.stream
+                .memcpy_stod(if v.is_empty() { &[0.0f64][..] } else { v })
+                .expect("upload")
+        };
+        let up_i = |v: &[i32]| -> CudaSlice<i32> { rt.stream.memcpy_stod(v).expect("upload") };
+        let up_u = |v: &[u32]| -> CudaSlice<u32> {
+            rt.stream
+                .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v })
+                .expect("upload")
+        };
+        // Surface velocity at the foot per face, and the open flags.
+        let mut ub: [Vec<f64>; 3] = [Vec::new(), Vec::new(), Vec::new()];
+        let mut open: [Vec<i32>; 3] = [Vec::new(), Vec::new(), Vec::new()];
+        for c in 0..3 {
+            let (ni, nj, nk) = match c {
+                0 => (nx + 1, ny, nz),
+                1 => (nx, ny + 1, nz),
+                _ => (nx, ny, nz + 1),
+            };
+            let mut ubc = vec![0.0; counts[c]];
+            let mut opc = vec![0i32; counts[c]];
+            for k in 0..nk {
+                for j in 0..nj {
+                    for i in 0..ni {
+                        let idx = (k * nj + j) * ni + i;
+                        let x = [
+                            (i as f64 + if c == 0 { 0.0 } else { 0.5 }) * h[0],
+                            (j as f64 + if c == 1 { 0.0 } else { 0.5 }) * h[1],
+                            (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2],
+                        ];
+                        ubc[idx] = mask.surface_velocity_at(body, x, c, 0.0);
+                        opc[idx] = i32::from(match c {
+                            0 => mask.u_open(idx),
+                            1 => mask.v_open(idx),
+                            _ => mask.w_open(idx),
+                        });
+                    }
+                }
+            }
+            ub[c] = ubc;
+            open[c] = opc;
+        }
+        let (wall_flux, _) = mask.wall_flux_table(body, 0.0);
+        let nc = g.cells();
+        let active: Vec<i32> = (0..nc).map(|i| i32::from(mask.cell_active(i))).collect();
+        let owner: Vec<u32> = (0..nc)
+            .map(|i| mask.master(i).unwrap_or(i) as u32)
+            .collect();
+        let mut fold_ptr = Vec::with_capacity(nc + 1);
+        let mut fold_idx = Vec::new();
+        let mut slaves_of: Vec<Vec<u32>> = vec![Vec::new(); nc];
+        let mut merged = 0;
+        for i in 0..nc {
+            if let Some(m) = mask.master(i) {
+                slaves_of[m].push(i as u32);
+                merged += 1;
+            }
+        }
+        fold_ptr.push(0u32);
+        for list in &slaves_of {
+            fold_idx.extend_from_slice(list);
+            fold_ptr.push(fold_idx.len() as u32);
+        }
+        Some(Self {
+            a: [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)],
+            d: [up_f(&cut.d_u), up_f(&cut.d_v), up_f(&cut.d_w)],
+            ub: [up_f(&ub[0]), up_f(&ub[1]), up_f(&ub[2])],
+            wall_flux: up_f(&wall_flux),
+            open: [up_i(&open[0]), up_i(&open[1]), up_i(&open[2])],
+            active: up_i(&active),
+            owner: up_u(&owner),
+            fold_ptr: up_u(&fold_ptr),
+            fold_idx: up_u(&fold_idx),
+            cell_flux: rt.stream.alloc_zeros::<f64>(nc).expect("alloc"),
+            merged,
+        })
+    }
+
+    fn ptrs(&self) -> E3CutPtrs {
+        let rt = runtime();
+        let s = &rt.stream;
+        let pf = |x: &CudaSlice<f64>| x.device_ptr(s).0;
+        let pi = |x: &CudaSlice<i32>| x.device_ptr(s).0;
+        let pu = |x: &CudaSlice<u32>| x.device_ptr(s).0;
+        E3CutPtrs {
+            ptrs: [
+                pf(&self.a[0]),
+                pf(&self.a[1]),
+                pf(&self.a[2]),
+                pf(&self.d[0]),
+                pf(&self.d[1]),
+                pf(&self.d[2]),
+                pf(&self.ub[0]),
+                pf(&self.ub[1]),
+                pf(&self.ub[2]),
+                pf(&self.wall_flux),
+                pi(&self.open[0]),
+                pi(&self.open[1]),
+                pi(&self.open[2]),
+                pi(&self.active),
+                pu(&self.owner),
+                pu(&self.fold_ptr),
+                pu(&self.fold_idx),
+                pf(&self.cell_flux),
+            ],
+        }
+    }
+}
+
+impl DeviceStep {
+    /// One step on the cut-cell mask (the host `advance` with the cut
+    /// predictor, the apertured merged continuity and the owner-read
+    /// corrections).
+    pub(super) fn advance_cut(&mut self, dt: f64) -> StepResult {
+        let rt = runtime();
+        let k = cut_kernels();
+        let g = self.grid;
+        let t_old = self.solver.time();
+        let t_new = t_old + dt;
+        let t0 = Instant::now();
+        rt.stream
+            .memcpy_dtod(&self.u, &mut self.u_old)
+            .expect("u_old");
+        rt.stream
+            .memcpy_dtod(&self.v, &mut self.v_old)
+            .expect("v_old");
+        rt.stream
+            .memcpy_dtod(&self.w, &mut self.w_old)
+            .expect("w_old");
+        self.upload_tables(t_old);
+        let prm: E3Params = self.params(dt);
+        let ptrs: E3Ptrs = 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),
+        ];
+        for c in 0..3i32 {
+            unsafe {
+                rt.stream
+                    .launch_builder(&k.predict)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&c)
+                    .launch(cfg(counts[c as usize]))
+                    .expect("e3_cut_predict");
+            }
+        }
+        self.launch_sides(prm, &ptrs, 0);
+        self.upload_tables(t_new);
+        self.launch_sides(prm, &ptrs, 1);
+        rt.stream
+            .memcpy_dtod(&self.u, &mut self.u_star)
+            .expect("u*");
+        rt.stream
+            .memcpy_dtod(&self.v, &mut self.v_star)
+            .expect("v*");
+        rt.stream
+            .memcpy_dtod(&self.w, &mut self.w_star)
+            .expect("w*");
+        rt.stream.synchronize().expect("sync");
+        let t_pred = t0.elapsed();
+        if self.cg.is_none() || self.cg_dt != dt {
+            let problem = self.solver.poisson_operator(g, dt);
+            let params = MultigridParameters {
+                precision: self.solver.params.poisson_precision,
+                smoother: self.solver.params.poisson_smoother,
+                ..MultigridParameters::default()
+            };
+            self.cg = Some(DeviceCg::new(&problem, &params));
+            self.cg_dt = dt;
+        }
+        let anchor = self.solver.anchor_cell(g);
+        let mut total = 0;
+        let mut final_residual = f64::INFINITY;
+        let mut cg_iterations = 0;
+        let mut t_poisson = std::time::Duration::ZERO;
+        let mut t_apply = std::time::Duration::ZERO;
+        let one = 1i32;
+        let zero = 0i32;
+        for corrector in 0..self.solver.params.corrector_steps.max(1) {
+            let tp = Instant::now();
+            unsafe {
+                rt.stream
+                    .launch_builder(&k.divergence)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&one)
+                    .launch(cfg(g.cells()))
+                    .expect("e3_cut_divergence");
+                rt.stream
+                    .launch_builder(&k.fold)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&one)
+                    .arg(&mut self.partial)
+                    .launch(cfg(g.cells()))
+                    .expect("e3_cut_fold");
+            }
+            let source_scale = self.reduce();
+            let inner_stop = self.solver.inner_stop(g, source_scale);
+            if corrector > 0 {
+                rt.stream.memset_zeros(&mut self.p_prime).expect("p' = 0");
+            }
+            let sol = {
+                let cg = self.cg.as_mut().expect("cg");
+                cg.solve_device(&self.sp, &mut self.p_prime, inner_stop, anchor, 0)
+            };
+            cg_iterations += sol.iterations;
+            t_poisson += tp.elapsed();
+            let ta = Instant::now();
+            for c in 0..3i32 {
+                unsafe {
+                    rt.stream
+                        .launch_builder(&k.correct)
+                        .arg(&prm)
+                        .arg(&ptrs)
+                        .arg(&cptrs)
+                        .arg(&c)
+                        .launch(cfg(counts[c as usize]))
+                        .expect("e3_cut_correct");
+                }
+            }
+            if prm.periodic_z != 0 {
+                self.launch_sides(prm, &ptrs, 0);
+            }
+            unsafe {
+                rt.stream
+                    .launch_builder(&k.add_p)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .launch(cfg(g.cells()))
+                    .expect("e3_cut_add_p");
+                rt.stream
+                    .launch_builder(&k.divergence)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&zero)
+                    .launch(cfg(g.cells()))
+                    .expect("e3_cut_divergence");
+                rt.stream
+                    .launch_builder(&k.fold)
+                    .arg(&prm)
+                    .arg(&ptrs)
+                    .arg(&cptrs)
+                    .arg(&zero)
+                    .arg(&mut self.partial)
+                    .launch(cfg(g.cells()))
+                    .expect("e3_cut_fold");
+            }
+            let imbalance = self.reduce();
+            let reference_flux = self.solver.reference_flux(g);
+            let mass_residual = if reference_flux > 0.0 {
+                imbalance / reference_flux
+            } else {
+                imbalance
+            };
+            final_residual = mass_residual;
+            total += 1;
+            t_apply += ta.elapsed();
+            if mass_residual < self.solver.params.tolerance {
+                break;
+            }
+            rt.stream
+                .memcpy_dtod(&self.u, &mut self.u_star)
+                .expect("u*");
+            rt.stream
+                .memcpy_dtod(&self.v, &mut self.v_star)
+                .expect("v*");
+            rt.stream
+                .memcpy_dtod(&self.w, &mut self.w_star)
+                .expect("w*");
+        }
+        self.solver.set_time(t_new);
+        if let Some(tm) = self.timers.as_mut() {
+            tm.predictor_ns += t_pred.as_nanos() as u64;
+            tm.poisson_ns += t_poisson.as_nanos() as u64;
+            tm.apply_ns += t_apply.as_nanos() as u64;
+            tm.steps += 1;
+            tm.cg_iterations += cg_iterations as u64;
+        }
+        StepResult {
+            fresh_cells: 0,
+            converged: final_residual < self.solver.params.tolerance,
+            corrector_steps_performed: total,
+            final_residual,
+            poisson_iterations: cg_iterations,
+        }
+    }
+}