From aca77fe1ac229ae3356f18d1372e059109f19bc9 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 19:33:33 -0500 Subject: [PATCH] embedded3: the cut mask's imposition and the device surface-velocity tables restricted to a 4-cell band (deeper prescribed faces are never read; the surface velocity was the moving step's cost); the operator rebuild counted in the rebuild share Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/src/kernels/cuda/e3_cut.cu | 4 +++ .../incompressible/embedded3/impose.rs | 28 +++++++++++++++++-- .../embedded3/step/device/cut.rs | 13 +++++++-- 3 files changed, 39 insertions(+), 6 deletions(-) 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 3796086..c1698e8 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -323,6 +323,10 @@ extern "C" __global__ void e3_cut_impose(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; + /* only within the imposition band (host `impose`: 4 cells) */ + const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w); + double h_min = fmin(fmin(g.dx, g.dy), g.dz); + if (fabs(dist[t]) > 4.0 * h_min) return; const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w); double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w); out[t] = ubt[t]; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs index 6b0a176..1cda0a7 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs @@ -6,7 +6,19 @@ use super::body::Body; use super::wall::{FaceKind, Mask}; +/// The imposition band of a cut mask, in cells. +pub(crate) const IMPOSE_BAND_CELLS: f64 = 4.0; + impl Mask { + /// The band (in length units) within which a cut mask's prescribed + /// faces are imposed; `None` on the binary wall (every solid face). + #[must_use] + pub fn impose_band(&self) -> Option { + self.cut + .as_ref() + .map(|_| IMPOSE_BAND_CELLS * self.grid.dx.min(self.grid.dy).min(self.grid.dz)) + } + /// Impose the wall on `(u, v, w)` from the same field. pub fn impose(&self, body: &Body, u: &mut [f64], v: &mut [f64], w: &mut [f64], t: f64) -> f64 { let (us, vs, ws) = (u.to_vec(), v.to_vec(), w.to_vec()); @@ -30,11 +42,21 @@ impl Mask { ) -> f64 { let g = self.grid; let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + // On a cut mask only the prescribed faces within the band are + // imposed: deeper ones are never read (the predictor reaches two + // faces past an open one) and the surface velocity is the costly + // part of a moving body's step. + let band = self.impose_band(); + let near = |d: &[f64], idx: usize| band.is_none_or(|b| d[idx].abs() <= b); + let (d_u, d_v, d_w): (&[f64], &[f64], &[f64]) = match &self.cut { + Some(c) => (&c.d_u, &c.d_v, &c.d_w), + None => (&[], &[], &[]), + }; for k in 0..nz { for j in 0..ny { for i in 1..nx { let idx = g.uface(k, j, i); - if self.u_kind[idx] == FaceKind::Solid { + if self.u_kind[idx] == FaceKind::Solid && near(d_u, idx) { u[idx] = body .surface_velocity( i as f64 * dx, @@ -49,7 +71,7 @@ impl Mask { for j in 1..ny { for i in 0..nx { let idx = g.vface(k, j, i); - if self.v_kind[idx] == FaceKind::Solid { + if self.v_kind[idx] == FaceKind::Solid && near(d_v, idx) { v[idx] = body .surface_velocity( (i as f64 + 0.5) * dx, @@ -66,7 +88,7 @@ impl Mask { for j in 0..ny { for i in 0..nx { let idx = g.wface(k, j, i); - if self.w_kind[idx] == FaceKind::Solid { + if self.w_kind[idx] == FaceKind::Solid && near(d_w, idx) { w[idx] = body .surface_velocity( (i as f64 + 0.5) * dx, 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 1b92ab3..c7f930b 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 @@ -107,7 +107,10 @@ impl DeviceCut { .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v }) .expect("upload") }; - // Surface velocity at the foot per face, and the open flags. + // Surface velocity at the foot per face within the imposition band + // (zero beyond it: never read), and the open flags. + let band = mask.impose_band().unwrap_or(f64::INFINITY); + let dists: [&[f64]; 3] = [&cut.d_u, &cut.d_v, &cut.d_w]; 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 { @@ -127,7 +130,11 @@ 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] = mask.surface_velocity_at(body, x, c, t); + ubc[idx] = if dists[c][idx].abs() <= band { + mask.surface_velocity_at(body, x, c, t) + } else { + 0.0 + }; opc[idx] = i32::from(if projection { match c { 0 => mask.u_open(idx), @@ -315,7 +322,6 @@ impl DeviceStep { .expect("w*"); } let cptrs = self.cut.as_ref().expect("cut").ptrs(); - let rebuild = t_rebuild.elapsed(); if self.cg.is_none() || self.cg_dt != dt { let problem = self.solver.poisson_operator(g, dt); let params = MultigridParameters { @@ -326,6 +332,7 @@ impl DeviceStep { self.cg = Some(DeviceCg::new(&problem, ¶ms)); self.cg_dt = dt; } + let rebuild = t_rebuild.elapsed(); let anchor = self.solver.anchor_cell(g); let mut total = 0; let mut final_residual = f64::INFINITY;