From 58a5501cefeee59a9a05bc91be4ceee2aae96f13 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 19:48:15 -0500 Subject: [PATCH] embedded3: the cut mask's surface velocity uses the interpolant's distance and normal (one body call per face instead of eight); the device shares the surface-velocity tables between a step's phases; sphere gate re-read 0.965/0.961 Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/cutwall.rs | 67 +++++++++++++++++-- .../embedded3/step/device/cut.rs | 33 +++++++-- 2 files changed, 89 insertions(+), 11 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs index a066067..a2d96a1 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -361,16 +361,73 @@ impl Mask { } /// The surface velocity component `c` at the foot of the normal from - /// the face centre. + /// the face centre `x`. With a cut geometry the signed distance and + /// the normal come from the geometry's own corner values (the + /// trilinear interpolant at the face centre, its gradient by central + /// differences of the neighbouring face centres) — one call of the + /// body's velocity per face instead of eight of its distance. pub fn surface_velocity_at(&self, body: &Body, x: [f64; 3], c: usize, t: f64) -> f64 { let g = self.grid; - let eps = 1e-6 * g.dx.min(g.dy).min(g.dz); - let s = body.phi(x[0], x[1], x[2], t); - let (n1, n2, n3) = body.normal(x[0], x[1], x[2], t, eps); - let v = body.surface_velocity(x[0] - s * n1, x[1] - s * n2, x[2] - s * n3, t); + let (s, n) = match self.cut.as_ref() { + Some(cut) => { + let (s, n) = self.interpolant_distance_and_normal(cut, x); + (s, n) + } + None => { + let eps = 1e-6 * g.dx.min(g.dy).min(g.dz); + let s = body.phi(x[0], x[1], x[2], t); + let (n1, n2, n3) = body.normal(x[0], x[1], x[2], t, eps); + (s, [n1, n2, n3]) + } + }; + let v = body.surface_velocity(x[0] - s * n[0], x[1] - s * n[1], x[2] - s * n[2], t); [v.0, v.1, v.2][c] } + /// φ and its unit gradient at `x` from the trilinear interpolant of the + /// corner values (the cut geometry's own surface). + fn interpolant_distance_and_normal(&self, cut: &CutGeometry, x: [f64; 3]) -> (f64, [f64; 3]) { + let g = self.grid; + let (nx, ny, nz) = (g.nx as i64, g.ny as i64, g.nz as i64); + let h = [g.dx, g.dy, g.dz]; + let node = |i: i64, j: i64, k: i64| -> f64 { + let i = i.clamp(0, nx); + let j = j.clamp(0, ny); + let k = k.clamp(0, nz); + cut.phi[((k * (ny + 1) + j) * (nx + 1) + i) as usize] + }; + let gx = x[0] / h[0]; + let gy = x[1] / h[1]; + let gz = x[2] / h[2]; + let (i0, j0, k0) = (gx.floor() as i64, gy.floor() as i64, gz.floor() as i64); + let (fx, fy, fz) = (gx - i0 as f64, gy - j0 as f64, gz - k0 as f64); + // Trilinear value and its partial derivatives. + let c = |di: i64, dj: i64, dk: i64| node(i0 + di, j0 + dj, k0 + dk); + let lerp = |a: f64, b: f64, f: f64| a + f * (b - a); + let c00 = lerp(c(0, 0, 0), c(1, 0, 0), fx); + let c10 = lerp(c(0, 1, 0), c(1, 1, 0), fx); + let c01 = lerp(c(0, 0, 1), c(1, 0, 1), fx); + let c11 = lerp(c(0, 1, 1), c(1, 1, 1), fx); + let c0 = lerp(c00, c10, fy); + let c1 = lerp(c01, c11, fy); + let s = lerp(c0, c1, fz); + let dx0 = lerp(c(1, 0, 0) - c(0, 0, 0), c(1, 1, 0) - c(0, 1, 0), fy); + let dx1 = lerp(c(1, 0, 1) - c(0, 0, 1), c(1, 1, 1) - c(0, 1, 1), fy); + let dphi_dx = lerp(dx0, dx1, fz) / h[0]; + let dy0 = lerp(c(0, 1, 0) - c(0, 0, 0), c(1, 1, 0) - c(1, 0, 0), fx); + let dy1 = lerp(c(0, 1, 1) - c(0, 0, 1), c(1, 1, 1) - c(1, 0, 1), fx); + let dphi_dy = lerp(dy0, dy1, fz) / h[1]; + let dz0 = lerp(c(0, 0, 1) - c(0, 0, 0), c(1, 0, 1) - c(1, 0, 0), fx); + let dz1 = lerp(c(0, 1, 1) - c(0, 1, 0), c(1, 1, 1) - c(1, 1, 0), fx); + let dphi_dz = lerp(dz0, dz1, fy) / h[2]; + let norm = (dphi_dx * dphi_dx + dphi_dy * dphi_dy + dphi_dz * dphi_dz).sqrt(); + if norm > 0.0 { + (s, [dphi_dx / norm, dphi_dy / norm, dphi_dz / norm]) + } else { + (s, [1.0, 0.0, 0.0]) + } + } + /// The volume fluxes of the surface velocity through every cell's wall /// into the body, `U_b·W_c` (zero for a body at rest; the porous /// manufactured surface's flux otherwise), made compatible: the net 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 index c7f930b..941d442 100644 --- 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 @@ -62,6 +62,8 @@ unsafe impl ValidAsZeroBits for E3CutPtrs {} /// The static cut-cell mask on the device. pub(super) struct DeviceCut { + /// The surface-velocity tables (kept to share between the phases). + ub_host: [Vec; 3], a: [CudaSlice; 3], d: [CudaSlice; 3], ub: [CudaSlice; 3], @@ -88,6 +90,18 @@ impl DeviceCut { /// The tables of the solver's cut mask (`None` without one) for /// `phase`, the surface velocities at the mask's time `t`. pub(super) fn build(solver: &Solver, g: Grid, phase: Phase, t: f64) -> Option { + Self::build_with(solver, g, phase, t, None) + } + + /// As [`Self::build`], reusing the surface-velocity tables of `shared` + /// (built at the same mask and time) instead of evaluating them. + pub(super) fn build_with( + solver: &Solver, + g: Grid, + phase: Phase, + t: f64, + shared: Option<&[Vec; 3]>, + ) -> Option { let mask = solver.mask()?; let cut = mask.cut()?; let body = solver.body()?; @@ -130,10 +144,12 @@ impl DeviceCut { (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] = if dists[c][idx].abs() <= band { - mask.surface_velocity_at(body, x, c, t) - } else { - 0.0 + ubc[idx] = match shared { + Some(sh) => sh[c][idx], + None if dists[c][idx].abs() <= band => { + mask.surface_velocity_at(body, x, c, t) + } + None => 0.0, }; opc[idx] = i32::from(if projection { match c { @@ -198,10 +214,13 @@ impl DeviceCut { fold_idx.extend_from_slice(list); fold_ptr.push(fold_idx.len() as u32); } + let ub_host = ub; + let ub_dev = [up_f(&ub_host[0]), up_f(&ub_host[1]), up_f(&ub_host[2])]; Some(Self { + ub_host, a: [up_f(ap_u), up_f(ap_v), up_f(ap_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])], + ub: ub_dev, wall_flux: up_f(&wall_flux), open: [up_i(&open[0]), up_i(&open[1]), up_i(&open[2])], active: up_i(&active), @@ -441,7 +460,9 @@ impl DeviceStep { // The prescribed faces take the surface velocity (the host's // end-of-step impose), and the next predictor's tables. if self.solver.is_moving() { - self.cut = DeviceCut::build(&self.solver, g, Phase::Predictor, t_new); + let shared = self.cut.take().map(|c| c.ub_host); + self.cut = + DeviceCut::build_with(&self.solver, g, Phase::Predictor, t_new, shared.as_ref()); let cptrs = self.cut.as_ref().expect("cut").ptrs(); for c in 0..3i32 { unsafe {