From 3b3d6c84c0bdc82898bed09687bd061f7141a9d6 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 17:22:52 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20item=209b:=20the=20device=20step=20?= =?UTF-8?q?carries=20a=20static=20cut-cell=20mask=20(e3=5Fcut.cu:=20cut=20?= =?UTF-8?q?predictor,=20apertured=20merged=20continuity=20with=20the=20fol?= =?UTF-8?q?d,=20owner-read=20corrections;=20device=20CG=20off-stencil=20li?= =?UTF-8?q?nks)=20=E2=80=94=20host=20=3D=20device=20to=202e-10=20(CFD1=20c?= =?UTF-8?q?ylinder=20nz=204)=20and=204e-14=20(sphere)=20under=20tight=20to?= =?UTF-8?q?lerances?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/src/kernels/cuda/e3_cg.cu | 12 +- .../rtx-cfd/src/kernels/cuda/e3_cut.cu | 309 ++++++++++++++ .../embedded3/poisson/device_cg.rs | 29 ++ .../embedded3/poisson/problem.rs | 9 +- .../incompressible/embedded3/step/device.rs | 14 + .../embedded3/step/device/cut.rs | 380 ++++++++++++++++++ .../rtx-cfd/tests/embedded3_cg_device.rs | 52 +++ .../rtx-cfd/tests/embedded3_device_cut.rs | 180 +++++++++ 8 files changed, 982 insertions(+), 3 deletions(-) create mode 100644 crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_device_cut.rs diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cg.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cg.cu index 56b5e4f..c72e41b 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cg.cu +++ b/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); } diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu new file mode 100644 index 0000000..d367912 --- /dev/null +++ b/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]]; +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs index c97f0ef..f6aee3d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs @@ -70,6 +70,10 @@ pub struct DeviceCg { at: CudaSlice, ab: CudaSlice, ap: CudaSlice, + /// Off-stencil links per cell (CSR; the merged small cells). + link_ptr: CudaSlice, + link_idx: CudaSlice, + link_coef: CudaSlice, b: CudaSlice, r: CudaSlice, z: CudaSlice, @@ -118,6 +122,18 @@ impl DeviceCg { }; let up_f = |v: &[f64]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; let zeros = || rt.stream.alloc_zeros::(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) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs index e342096..a6a2441 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/problem.rs +++ b/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})")); } } 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 43e49a5..64ede77 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 @@ -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, initialized: bool, + /// A static cut-cell mask's tables (item 9b), when the solver has one. + cut: Option, } 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 { 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; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs new file mode 100644 index 0000000..da82d3c --- /dev/null +++ b/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, + predict: CudaFunction, + divergence: CudaFunction, + fold: CudaFunction, + correct: CudaFunction, + add_p: CudaFunction, +} + +static CUT_KERNELS_ONCE: OnceLock = 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; 3], + d: [CudaSlice; 3], + ub: [CudaSlice; 3], + wall_flux: CudaSlice, + open: [CudaSlice; 3], + active: CudaSlice, + owner: CudaSlice, + fold_ptr: CudaSlice, + fold_idx: CudaSlice, + cell_flux: CudaSlice, + 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 { + 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 { + rt.stream + .memcpy_stod(if v.is_empty() { &[0.0f64][..] } else { v }) + .expect("upload") + }; + let up_i = |v: &[i32]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; + let up_u = |v: &[u32]| -> CudaSlice { + 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; 3] = [Vec::new(), Vec::new(), Vec::new()]; + let mut open: [Vec; 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 = (0..nc).map(|i| i32::from(mask.cell_active(i))).collect(); + let owner: Vec = (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![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::(nc).expect("alloc"), + merged, + }) + } + + fn ptrs(&self) -> E3CutPtrs { + let rt = runtime(); + let s = &rt.stream; + let pf = |x: &CudaSlice| x.device_ptr(s).0; + let pi = |x: &CudaSlice| x.device_ptr(s).0; + let pu = |x: &CudaSlice| 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, ¶ms)); + 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, + } + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_cg_device.rs b/crates/specialized/rtx-cfd/tests/embedded3_cg_device.rs index 44ce289..35c2321 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_cg_device.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_cg_device.rs @@ -176,3 +176,55 @@ fn the_singular_box_honours_the_anchor_and_the_mean() { println!(" singular box without anchor: mean {mean:.3e} of {scale:.3e}"); assert!(mean.abs() <= 1e-13 * scale, "mean {mean:.3e}"); } + +/// Item 9b: off-stencil links (the virtually merged small cells) on the +/// device operator — the channel with 60 symmetric links between active +/// cells two apart; device = host PCG to the solve's accuracy. +#[test] +fn device_cg_carries_the_links() { + let mut problem = channel(96, 40, 8, true, true, 7); + let n = problem.nx * problem.ny * problem.nz; + let mut seed = 12345u64; + let mut links = Vec::new(); + while links.len() < 60 { + seed = seed + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + let a = (seed >> 20) as usize % n; + let b = a + 2; + if b < n && problem.active[a] && problem.active[b] && a % problem.nx + 2 < problem.nx { + links.push((a, b, 0.7 * problem.ae[a].max(problem.aw[a]))); + } + } + problem.links = links; + problem.validate().expect("linked problem"); + let params = MultigridParameters { + smoother: MgSmoother::RedBlack, + ..MultigridParameters::default() + }; + let tol = 1e-9 * problem.rhs.iter().map(|r| r.abs()).sum::(); + let mut p_host = vec![0.0; n]; + let host = solve_pcg(&problem, &mut p_host, ¶ms, tol, None); + let mut p_dev = vec![0.0; n]; + let mut cache = DevicePcgCache::default(); + let dev = solve_pcg_device_cached(&problem, &mut p_dev, ¶ms, tol, None, &mut cache); + let scale = p_host.iter().fold(0.0_f64, |m, &x| m.max(x.abs())); + let diff = p_host + .iter() + .zip(&p_dev) + .fold(0.0_f64, |m, (&a, &b)| m.max((a - b).abs())); + println!( + " linked channel (60 links): host {} it residual {:.3e}; device {} it residual {:.3e}; max |Δp| {:.3e} on {:.3e}", + host.iterations, host.residual, dev.iterations, dev.residual, diff, scale + ); + assert!(host.converged && dev.converged); + assert!( + diff < 1e-8 * scale, + "device with links differs from the host: {diff:.3e} on {scale:.3e}" + ); + assert!( + problem.residual_l1(&p_dev) < 2.0 * tol, + "device residual on the linked operator {:.3e}", + problem.residual_l1(&p_dev) + ); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_device_cut.rs b/crates/specialized/rtx-cfd/tests/embedded3_device_cut.rs new file mode 100644 index 0000000..be37140 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_device_cut.rs @@ -0,0 +1,180 @@ +//! embedded3 item 9b: the device step carrying a static cut-cell mask — +//! host vs device on the manufactured sphere (CutCell, merging in) and on +//! the CFD1 cylinder (nz 4 periodic), fields agreeing to the registered +//! `1e-9·scale` over 100 steps under TIGHT tolerances (inner stop 1e-6, +//! mass 1e-12 — item 7's identity rule; under the default tolerances the +//! two CGs differ at the inner-stop level and the agreement is reported), +//! equal corrector counts; the device step time recorded. +//! +//! `RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_device_cut -- --nocapture` +#![cfg(feature = "cuda")] + +mod embedded3_sphere; + +use embedded3_sphere::{C, MU, R, RHO, boundary3, source3, u3, v3, w3}; +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; + +fn max_diff(a: &[f64], b: &[f64]) -> f64 { + a.iter() + .zip(b) + .fold(0.0_f64, |m, (&x, &y)| m.max((x - y).abs())) +} + +fn scale(a: &[f64]) -> f64 { + a.iter().fold(0.0_f64, |m, &x| m.max(x.abs())) +} + +/// Builds two identical solvers (host and device) and marches both. +fn march( + make: &dyn Fn(bool) -> (Solver, Grid), + dt: f64, + steps: usize, + label: &str, + bound: Option, + tight: bool, +) { + let make = |t: bool| make(t); + let (mut host, g) = make(tight); + let mut fh = Field::new(g); + host.initialize(&mut fh); + let (mut dev_solver, _) = make(tight); + let mut fd = Field::new(g); + dev_solver.initialize(&mut fd); + let mut device = DeviceStep::new(dev_solver, g); + device.upload(&fd); + let merged = device.merged_cells(); + let start = std::time::Instant::now(); + let mut differ = 0; + for _ in 0..steps { + let rh = host.advance(&mut fh, dt); + let rd = device.advance(dt); + if rh.corrector_steps_performed != rd.corrector_steps_performed { + differ += 1; + } + } + let device_seconds = start.elapsed().as_secs_f64(); + device.download(&mut fd); + let du = max_diff(&fh.u, &fd.u) + .max(max_diff(&fh.v, &fd.v)) + .max(max_diff(&fh.w, &fd.w)); + let su = scale(&fh.u).max(scale(&fh.v)).max(scale(&fh.w)); + let dp = max_diff(&fh.p, &fd.p); + let sp = scale(&fh.p).max(RHO); + println!( + " {label} (tight {tight}): {steps} steps, {merged} merged cells; host vs device max |Δu| {du:.3e} on {su:.3e}, max |Δp| {dp:.3e} on {sp:.3e}; corrector counts differ on {differ} steps; {:.1} ms per step (host + device)", + 1e3 * device_seconds / steps as f64 + ); + if let Some(bound) = bound { + assert!(du < bound * su, "velocity differs: {du:.3e} on {su:.3e}"); + assert!(dp < bound * sp, "pressure differs: {dp:.3e} on {sp:.3e}"); + // The corrector counts are reported, not gated: at a mass + // tolerance of 1e-12 the two residuals differ in their last digits + // (the fields agree to rounding regardless). + } +} + +fn tolerances(tight: bool) -> (f64, f64) { + if tight { (1e-12, 1e-6) } else { (1e-8, 1e-2) } +} + +#[test] +fn sphere_cut_cell_host_equals_device() { + let n = 12; + let h = 1.0 / n as f64; + let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h); + let make = |tight: bool| { + let (tolerance, inner_stop_factor) = tolerances(tight); + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + }, + Parameters { + corrector_steps: 2, + tolerance, + inner_stop_factor, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: WallScheme::CutCell, + ..Parameters::default() + }, + ); + solver.set_momentum_source(|x, y, z, _t| source3(x, y, z)); + solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z)); + solver.set_body( + Body::sphere(move |_t| C, R) + .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))), + ); + (solver, Grid::cubic(n, n, n, h)) + }; + march(&make, dt, 100, "sphere MMS n 12 CutCell", None, false); + march(&make, dt, 100, "sphere MMS n 12 CutCell", Some(1e-9), true); +} + +#[test] +fn cylinder_cut_cell_host_equals_device() { + use rtx_cfd::solvers::incompressible::EmbeddedBody; + let ny = 41; + let nz = 4; + let h = 0.41 / ny as f64; + let nx = (2.2 / h).round() as usize; + let dt = 2e-3; + let make = |tight: bool| { + let (tolerance, inner_stop_factor) = tolerances(tight); + let mut solver = Solver::new( + Fluid { + density: 1.0, + viscosity: 1e-3, + reference_velocity: 0.2, + reference_length: 0.1, + }, + Parameters { + corrector_steps: 2, + tolerance, + inner_stop_factor, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + x1: Side::PressureOutlet, + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(move |x, y, _z, _t| { + if x <= 0.0 { + (1.5 * 0.2 * 4.0 * y * (0.41 - y) / (0.41 * 0.41), 0.0, 0.0) + } else { + (0.0, 0.0, 0.0) + } + }); + let circle = EmbeddedBody::from_sdf(|x, y, _t| { + ((x - 0.2_f64).powi(2) + (y - 0.2_f64).powi(2)).sqrt() - 0.05 + }); + solver.set_body(Body::extruded(circle, nz as f64 * h)); + (solver, Grid::cubic(nx, ny, nz, h)) + }; + march( + &make, + dt, + 100, + "CFD1 cylinder ny 41 nz 4 periodic CutCell", + None, + false, + ); + march( + &make, + dt, + 100, + "CFD1 cylinder ny 41 nz 4 periodic CutCell", + Some(1e-9), + true, + ); +}