From d7250bc4cf54102471d647324a0335a6ee040eef Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sat, 19 Sep 2026 19:19:47 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20A3-i:=20the=20advancing-wall=20fric?= =?UTF-8?q?tion=20closure=20(RTX=5FE3=5FWALL=5FADVANCING=3D1,=20host=20+?= =?UTF-8?q?=20device,=20per-face=20normal-velocity=20table;=20f(0)=20=3D?= =?UTF-8?q?=201=20exactly)=20=E2=80=94=20REFUTED=20as=20built=20on=20the?= =?UTF-8?q?=20slab=20ladder=20(flag=20shear=2034.2=20/=2033.6=20/=2031.1?= =?UTF-8?q?=20vs=2028.9=20/=2032.3=20/=2030.1:=20spread=203.4=20->=203.1,?= =?UTF-8?q?=20now=20from=20above;=20worst=20mass=20residual=201e-8=20->=20?= =?UTF-8?q?1e-4);=20knob=20stays=20off,=20default=20digit-identical?= 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_cut.cu | 12 ++++++++ .../incompressible/embedded3/closure.rs | 29 ++++++++++++++++++- .../incompressible/embedded3/cutwall.rs | 12 ++++++-- .../embedded3/step/cut_predictor.rs | 7 ++++- .../incompressible/embedded3/step/device.rs | 4 ++- .../embedded3/step/device/cut.rs | 29 ++++++++++++++++++- .../incompressible/embedded3/step/mod.rs | 8 +++++ .../solvers/incompressible/embedded3/wall.rs | 2 ++ 8 files changed, 97 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 68f7146..bb71d10 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -29,8 +29,15 @@ struct E3Cut { const unsigned int *fold_ptr, *fold_idx; /* CSR: the slaves of every cell */ double *cell_flux; /* scratch per cell */ const double *s_u, *s_v, *s_w; /* open-part centroid shifts per face, 3 interleaved (S2-5) */ + const double *vn_u, *vn_v, *vn_w; /* the wall's normal velocity into the fluid per face (A3-i) */ }; +/* f(xi) = xi / (1 - exp(-xi)), f(0) = 1 exactly (closure.rs advancing_factor). */ +__device__ __forceinline__ double adv_factor(double xi) +{ + return fabs(xi) < 1e-8 ? 1.0 + 0.5 * xi : xi / (-expm1(-xi)); +} + /* 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) { @@ -311,6 +318,11 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu open face away from the body along the wall normal's dominant axis (order 2; the neighbour's old value explicit). */ double c1 = 1.0 / distance, shear_explicit = 0.0; + /* bit 9 of wall_order: the advancing-wall closure on the one-point coefficient (A3-i) */ + if (g.wall_order & 512) { + const double* vnt = c == 0 ? m.vn_u : (c == 1 ? m.vn_v : m.vn_w); + c1 = adv_factor(vnt[fidx] * distance / g.nu) / distance; + } if ((g.wall_order & 15) >= 2 && a_w > 0.0) { double nw[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w }; int d = 0; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs index dd82190..ca39135 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs @@ -2,6 +2,7 @@ //! open-part centroid shifts of the cut faces and the spacing they give //! the cross-direction diffusion (the default since S2-5), the exchange //! distance toward solid neighbours, and the quadratic wall gradient. +use super::body::Body; use super::cutwall::CvGeometry; use super::wall::Mask; @@ -107,18 +108,33 @@ impl Mask { (cv.distance / n_d).min(h) } + /// The wall's velocity along the interpolant's normal at the foot of + /// `x` (positive INTO the fluid: the wall advancing on it), A3-i's v_n. + pub fn surface_normal_velocity_at(&self, body: &Body, x: [f64; 3], t: f64) -> f64 { + let Some(cut) = self.cut.as_ref() else { + return 0.0; + }; + let (s, n) = self.interpolant_distance_and_normal(cut, x); + let v = body.surface_velocity(x[0] - s * n[0], x[1] - s * n[1], x[2] - s * n[2], t); + v.0 * n[0] + v.1 * n[1] + v.2 * n[2] + } + /// The wall-gradient coefficients of the unknown face of component /// `c` at `p` with control volume `cv`: `u'(0) = c_1 (u_f − U_b) + c_2 /// (u_n − U_b)` with `u_n` the face returned (one lattice step away /// from the body along the wall normal's dominant axis). Order 1, or /// no open neighbour: `(1/d_f, 0, None)`. + /// `xi` = v_n d_f / ν of the advancing-wall closure (0 when it is off): + /// the one-point coefficient is multiplied by `advancing_factor(xi)`; the + /// quadratic closure (order 2) is left as it is. pub(super) fn wall_gradient( &self, c: usize, p: [i64; 3], cv: &CvGeometry, + xi: f64, ) -> (f64, f64, Option) { - let linear = (1.0 / cv.distance, 0.0, None); + let linear = (advancing_factor(xi) / cv.distance, 0.0, None); if self.wall_order < 2 { return linear; } @@ -148,3 +164,14 @@ impl Mask { (d2 / (d1 * (d2 - d1)), -d1 / (d2 * (d2 - d1)), Some(f)) } } + +/// `f(ξ) = ξ / (1 − e^{−ξ})`: the ratio of the asymptotic-suction layer's wall +/// gradient to the linear one at the same value and distance; f(0) = 1 exactly. +#[must_use] +pub fn advancing_factor(xi: f64) -> f64 { + if xi.abs() < 1e-8 { + 1.0 + 0.5 * xi + } else { + xi / (-(-xi).exp_m1()) + } +} 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 ca420ff..8d4f60a 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -216,6 +216,7 @@ impl Mask { wall_distance_oblique: false, diffusion_transverse: false, distance_floor_fine: false, + wall_advancing: false, wall_exchange_axis: false, grad_weights: None, diffusion_centroid: false, @@ -658,8 +659,15 @@ impl Mask { if a_w == 0.0 { continue; } - let ub = self.surface_velocity_at(body, lat.face_position(c, p), c, t); - let (c1, c2, nb) = self.wall_gradient(c, p, &cv); + let x = lat.face_position(c, p); + let ub = self.surface_velocity_at(body, x, c, t); + let xi = if self.wall_advancing { + self.surface_normal_velocity_at(body, x, t) * cv.distance * self.density + / mu + } else { + 0.0 + }; + let (c1, c2, nb) = self.wall_gradient(c, p, &cv, xi); let un = nb.map_or(ub, |f| values[c][f]); force[c] += mu * a_w * (c1 * (values[c][idx] - ub) + c2 * (un - ub)); } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs index 8e4e4ca..1209abd 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs @@ -339,7 +339,12 @@ impl Solver { }); let a_w = (cv.wall[0] * cv.wall[0] + cv.wall[1] * cv.wall[1] + cv.wall[2] * cv.wall[2]).sqrt(); - let (c1, c2, nb) = mask.wall_gradient(c, p, &cv); + let xi = if mask.wall_advancing { + mask.surface_normal_velocity_at(body, x, t_old) * cv.distance * rho / mu + } else { + 0.0 + }; + let (c1, c2, nb) = mask.wall_gradient(c, p, &cv, xi); let shear = mu * a_w * c1; // The explicit part of the quadratic wall gradient (order 2). let shear_explicit = nb.map_or(0.0, |f| mu * a_w * c2 * (old[c][f] - ub)); 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 72ee676..cd1ba8e 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 @@ -302,7 +302,9 @@ impl DeviceStep { self.solver.params.diffusion_transverse && self.solver.params.diffusion_centroid, ) - + 256 * i32::from(self.solver.params.distance_floor_fine), + + 256 * i32::from(self.solver.params.distance_floor_fine) + // bit 9: the advancing-wall friction closure (A3-i). + + 512 * i32::from(self.solver.params.wall_advancing), dx: g.dx, dy: g.dy, dz: g.dz, 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 4014371..b8c04b4 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 @@ -56,7 +56,7 @@ fn cut_kernels() -> &'static CutKernels { #[repr(C)] #[derive(Clone, Copy)] struct E3CutPtrs { - ptrs: [u64; 21], + ptrs: [u64; 24], } unsafe impl DeviceRepr for E3CutPtrs {} unsafe impl ValidAsZeroBits for E3CutPtrs {} @@ -77,6 +77,8 @@ pub(super) struct DeviceCut { cell_flux: CudaSlice, /// The open-part centroid shifts per face (S2-5; one dummy entry when off). shift: [CudaSlice; 3], + /// The wall's normal velocity into the fluid per face (A3-i; one dummy entry when off). + vn: [CudaSlice; 3], pub(super) merged: usize, } @@ -132,6 +134,8 @@ impl DeviceCut { 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()]; + let advancing = mask.wall_advancing; + let mut vn: [Vec; 3] = [vec![0.0], vec![0.0], vec![0.0]]; for c in 0..3 { let (ni, nj, nk) = match c { 0 => (nx + 1, ny, nz), @@ -178,6 +182,25 @@ impl DeviceCut { }); ub[c] = ubc; open[c] = opc; + // A3-i: the wall's normal velocity into the fluid per face within the band. + if advancing { + vn[c] = (0..counts[c]) + .into_par_iter() + .map(|idx| { + if dists[c][idx].abs() <= band { + let (k, j, i) = (idx / (nj * ni), (idx / ni) % nj, idx % ni); + 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], + ]; + mask.surface_normal_velocity_at(body, x, t) + } else { + 0.0 + } + }) + .collect(); + } } let l_faces = lap.elapsed(); // The solver's current table (the GCL table on a moving body) when @@ -250,6 +273,7 @@ impl DeviceCut { Some(t) => [up_f(&t[0]), up_f(&t[1]), up_f(&t[2])], None => [up_f(&[0.0]), up_f(&[0.0]), up_f(&[0.0])], }, + vn: [up_f(&vn[0]), up_f(&vn[1]), up_f(&vn[2])], merged, }; if profile { @@ -339,6 +363,9 @@ impl DeviceCut { pf(&self.shift[0]), pf(&self.shift[1]), pf(&self.shift[2]), + pf(&self.vn[0]), + pf(&self.vn[1]), + pf(&self.vn[2]), ], } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index 73d626f..53caba5 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -114,6 +114,12 @@ pub struct Parameters { /// 0.05: the coarse floor doubles the distance of faces with α < 0.1 /// (the flat wall's θ 0.95 excess). ON since S2-6 (`RTX_E3_DISTANCE_FLOOR=coarse`). pub distance_floor_fine: bool, + /// The advancing-wall friction closure (A3-i): the one-point wall gradient + /// times `f(ξ) = ξ / (1 − e^{−ξ})`, `ξ = v_n d_f / ν`, with `v_n` the wall's + /// velocity into the fluid — the asymptotic-suction layer a wall moving + /// normal to itself sets up (thinner than a cell on the fast flag). f(0) = 1 + /// exactly: static cases unchanged. `RTX_E3_WALL_ADVANCING=1`. + pub wall_advancing: bool, /// The diffusive exchange of a fluid face with a SOLID neighbour face /// over the axis distance to the wall, `δ = min(h, d_f/|n_d|)`, and /// implicit — instead of the full `h`, which places the no-slip value @@ -164,6 +170,7 @@ impl Default for Parameters { .map_or(true, |v| v != "0"), distance_floor_fine: std::env::var("RTX_E3_DISTANCE_FLOOR") .map_or(true, |v| v != "coarse"), + wall_advancing: std::env::var("RTX_E3_WALL_ADVANCING").is_ok_and(|v| v == "1"), wall_exchange_axis: std::env::var("RTX_E3_WALL_EXCHANGE").is_ok_and(|v| v == "axis"), pressure_centroid: std::env::var("RTX_E3_PRESSURE_CENTROID").is_ok_and(|v| v == "1"), // ON by default since S2-5 (`=0` reproduces the records before it). @@ -303,6 +310,7 @@ impl Solver { m.diffusion_transverse = self.params.diffusion_transverse && self.params.diffusion_centroid; m.distance_floor_fine = self.params.distance_floor_fine; + m.wall_advancing = self.params.wall_advancing; m.diffusion_centroid = self.params.diffusion_centroid; if self.params.diffusion_centroid { m.compute_face_shifts(); diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs index a14bdd6..6ef8c11 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -104,6 +104,7 @@ pub struct Mask { /// The transverse centroid correction and the fine distance floor (S2-6). pub(super) diffusion_transverse: bool, pub(super) distance_floor_fine: bool, + pub(super) wall_advancing: bool, /// The axis-distance implicit wall exchange (S2-5). pub(super) wall_exchange_axis: bool, /// The centroid prototype's pressure-gradient weights per u / v / w face. @@ -529,6 +530,7 @@ impl Mask { wall_distance_oblique: false, diffusion_transverse: false, distance_floor_fine: false, + wall_advancing: false, wall_exchange_axis: false, grad_weights: None, diffusion_centroid: false,