From cda973df26e3e69eac949327cd6b0ccaa90dec04 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 07:21:57 -0500 Subject: [PATCH] embedded3 S2-4: quadratic wall gradient behind Parameters::wall_order (host predictor, operator route, e3_cut.cu; RTX_E3_WALL_ORDER); density pin holds at order 2 Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/src/kernels/cuda/e3_cut.cu | 22 ++++++++- .../rtx-cfd/src/kernels/cuda/e3_step.cu | 2 +- .../incompressible/embedded3/cutwall.rs | 46 ++++++++++++++++++- .../embedded3/step/cut_predictor.rs | 9 ++-- .../incompressible/embedded3/step/device.rs | 4 +- .../incompressible/embedded3/step/mod.rs | 10 ++++ .../solvers/incompressible/embedded3/wall.rs | 3 ++ 7 files changed, 87 insertions(+), 9 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 be389ba..b45478d 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -187,10 +187,28 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu 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; + /* The wall gradient: one-point (order 1) or quadratic through the next + 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; + if (g.wall_order >= 2 && a_w > 0.0) { + double nw[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w }; + int d = 0; + for (int kk = 1; kk < 3; ++kk) if (fabs(nw[kk]) > fabs(nw[d])) d = kk; + int q[3] = { i, j, k }; + q[d] -= nw[d] > 0.0 ? 1 : -1; + int fn = cut_face(g, c, q[0], q[1], q[2]); + if (fn >= 0 && cut_ap(m, c)[fn] > 0.0) { + double d1 = distance, d2 = d1 + h[d] * fabs(nw[d]); + c1 = d2 / (d1 * (d2 - d1)); + double c2 = -d1 / (d2 * (d2 - d1)); + shear_explicit = mu * a_w * c2 * (old_c[fn] - ub); + } + } + double shear = mu * a_w * c1; 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); + return (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit) / (inertia + shear); } /* The predictor on the open interior faces of component c. */ diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_step.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_step.cu index 687116d..6e10fa4 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_step.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_step.cu @@ -19,7 +19,7 @@ struct E3Params { int periodic_z; int bx0, bx1, by0, by1, bz0, bz1; int scheme; - int pad; + int wall_order; double dx, dy, dz, dt, rho, nu; }; 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 62497b6..4ec986b 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -209,6 +209,7 @@ impl Mask { merge_master: Vec::new(), scheme: crate::solvers::incompressible::ConvectionScheme::Upwind, density: 1.0, + wall_order: 1, }; mask.compute_merging(None); Ok(mask) @@ -306,6 +307,47 @@ impl Mask { self.compute_merging(Some(old)); } + /// 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)`. + pub(super) fn wall_gradient( + &self, + c: usize, + p: [i64; 3], + cv: &CvGeometry, + ) -> (f64, f64, Option) { + let linear = (1.0 / cv.distance, 0.0, None); + if self.wall_order < 2 { + return linear; + } + let a_w = + (cv.wall[0] * cv.wall[0] + cv.wall[1] * cv.wall[1] + cv.wall[2] * cv.wall[2]).sqrt(); + if a_w == 0.0 { + return linear; + } + let n = [cv.wall[0] / a_w, cv.wall[1] / a_w, cv.wall[2] / a_w]; + let mut d = 0; + for k in 1..3 { + if n[k].abs() > n[d].abs() { + d = k; + } + } + // `n` points into the body: step the other way. + let mut q = p; + q[d] -= if n[d] > 0.0 { 1 } else { -1 }; + let open = self.aperture(c, q).is_some_and(|a| a > 0.0); + if !open { + return linear; + } + let f = self.lattice().face(c, q).expect("open face"); + let h = [self.grid.dx, self.grid.dy, self.grid.dz]; + let d1 = cv.distance; + let d2 = d1 + h[d] * n[d].abs(); + (d2 / (d1 * (d2 - d1)), -d1 / (d2 * (d2 - d1)), Some(f)) + } + pub(super) fn lattice(&self) -> Lattice { Lattice { g: self.grid, @@ -659,7 +701,9 @@ impl Mask { continue; } let ub = self.surface_velocity_at(body, lat.face_position(c, p), c, t); - force[c] += mu * a_w * (values[c][idx] - ub) / cv.distance; + let (c1, c2, nb) = self.wall_gradient(c, p, &cv); + 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 f984203..e3ced72 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 @@ -226,7 +226,10 @@ 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 shear = mu * a_w / cv.distance; + let (c1, c2, nb) = mask.wall_gradient(c, p, &cv); + 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)); let fraction = if self.params.momentum_volume_cell_mean { let vol = |q: [i64; 3]| lat.cell(q).map_or(cv.alpha, |ci| mask.vol(ci)); 0.5 * (vol(cell_minus) + vol(cell_plus)) @@ -235,8 +238,8 @@ impl Solver { }; let v_eff = fraction.max(INERTIA_FLOOR) * h[c] * area[c]; let inertia = rho * v_eff / dt; - let u_star = - (inertia * u0 - conv + diff + pressure + source + shear * ub) / (inertia + shear); + let u_star = (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit) + / (inertia + shear); let v_alpha = fraction * h[c] * area[c]; (u_star, rho * (v_eff - v_alpha) * (u_star - u0) / dt) } 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 4cb3cf9..002570d 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 @@ -78,7 +78,7 @@ struct E3Params { bz0: i32, bz1: i32, scheme: i32, - pad: i32, + wall_order: i32, dx: f64, dy: f64, dz: f64, @@ -284,7 +284,7 @@ impl DeviceStep { bz0: side_code(b.z0), bz1: side_code(b.z1), scheme: scheme_code(self.solver.params.convection_scheme), - pad: 0, + wall_order: i32::from(self.solver.params.wall_order), dx: g.dx, dy: g.dy, dz: g.dz, 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 e8299ce..d190b05 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 @@ -91,6 +91,11 @@ pub struct Parameters { /// inertia of the fluid half of its volume instead of a tenth of a /// cell against O(1) fluxes). pub momentum_volume_cell_mean: bool, + /// The cut wall's shear closure: 1 = one-point `μ A_w (u_f − U_b)/d_f` + /// (the recorded form), 2 = the quadratic wall gradient through the + /// next open face along the wall normal's dominant axis (S2-4). + /// `Parameters::default()` reads `RTX_E3_WALL_ORDER` (default 1). + pub wall_order: u8, } impl Default for Parameters { @@ -107,6 +112,10 @@ impl Default for Parameters { max_surface_speed: None, aperture_substeps: 0, momentum_volume_cell_mean: false, + wall_order: std::env::var("RTX_E3_WALL_ORDER") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(1), } } } @@ -231,6 +240,7 @@ impl Solver { .map(|mut m| { m.scheme = self.params.convection_scheme; m.density = self.fluid.density; + m.wall_order = self.params.wall_order; m }) .expect("embedded mask") 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 2cd3881..443bc81 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -97,6 +97,8 @@ pub struct Mask { pub(super) scheme: crate::solvers::incompressible::ConvectionScheme, /// The fluid's density (the exchange route's convective flux). pub(super) density: f64, + /// The cut wall's shear closure order (S2-4). + pub(super) wall_order: u8, } /// The z lattice position of a query: the lower plane index, the upper @@ -509,6 +511,7 @@ impl Mask { merge_master: Vec::new(), scheme: crate::solvers::incompressible::ConvectionScheme::Upwind, density: 1.0, + wall_order: 1, }) }