From 99e4a7214bac72792676919c60bb032f6380212f Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 12:37:56 -0500 Subject: [PATCH] embedded3 S2-5: centroid diffusion default ON (RTX_E3_DIFFUSION_CENTROID=0 restores the records); flat_wall_position_is_second_order gate; closure.rs split (cutwall.rs under 700) Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/closure.rs | 137 ++++++++++++++++++ .../incompressible/embedded3/cutwall.rs | 129 ----------------- .../solvers/incompressible/embedded3/mod.rs | 1 + .../incompressible/embedded3/step/mod.rs | 6 +- .../rtx-cfd/tests/embedded3_dfg_2d1_host.rs | 14 +- .../rtx-cfd/tests/embedded3_wall_position.rs | 27 +++- 6 files changed, 178 insertions(+), 136 deletions(-) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs new file mode 100644 index 0000000..d37d2b6 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs @@ -0,0 +1,137 @@ +//! The cut wall's closures beyond the first form (S2-4, S2-5): the +//! 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::cutwall::CvGeometry; +use super::wall::Mask; + +impl Mask { + /// The shift of a face's open-part centroid from the face centre: + /// `½h(1 − α)` along the wall normal's in-plane part, away from the + /// body (zero for a full face or without the centroid diffusion). Read + /// from the tables of [`Self::compute_face_shifts`]. + pub(super) fn face_shift(&self, c: usize, p: [i64; 3]) -> [f64; 3] { + let (Some(t), Some(f)) = (self.face_shifts.as_ref(), self.lattice().face(c, p)) else { + return [0.0; 3]; + }; + [t[c][3 * f], t[c][3 * f + 1], t[c][3 * f + 2]] + } + + /// The per-face shift tables (three components interleaved). + #[must_use] + pub fn face_shift_tables(&self) -> Option<&[Vec; 3]> { + self.face_shifts.as_ref() + } + + /// Build the open-part centroid shifts of every cut face (S2-5). + pub fn compute_face_shifts(&mut self) { + let g = self.grid; + let h = [g.dx, g.dy, g.dz]; + let lat = self.lattice(); + let sizes = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()]; + let mut tables = [ + vec![0.0; 3 * sizes[0]], + vec![0.0; 3 * sizes[1]], + vec![0.0; 3 * sizes[2]], + ]; + for c in 0..3 { + let (ni, nj, nk) = ( + g.nx + usize::from(c == 0), + g.ny + usize::from(c == 1), + g.nz + usize::from(c == 2), + ); + for k in 0..nk { + for j in 0..nj { + for i in 0..ni { + let p = [i as i64, j as i64, k as i64]; + let Some(f) = lat.face(c, p) else { continue }; + let Some(alpha) = self.aperture(c, p) else { + continue; + }; + if alpha <= 0.0 || alpha >= 1.0 { + continue; + } + // Interior faces only (a control volume needs both cells). + let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c]; + if on_side && !(c == 2 && self.periodic_z) { + continue; + } + let cv = self.cv_geometry(c, p); + let mut n = cv.wall; + n[c] = 0.0; + let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if a == 0.0 { + continue; + } + for d in 0..3 { + // `wall` points into the body: the open part lies the other way. + tables[c][3 * f + d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a; + } + } + } + } + } + self.face_shifts = Some(tables); + } + + /// The distance over which a fluid face exchanges momentum with a solid + /// neighbour face along `d`: the full spacing, or (S2-5) the axis + /// distance from the open part's centroid to the wall, `min(h, d_f/|n_d|)`. + pub(super) fn exchange_delta(&self, cv: &CvGeometry, d: usize) -> f64 { + let h = [self.grid.dx, self.grid.dy, self.grid.dz][d]; + if !self.wall_exchange_axis { + return h; + } + 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 h; + } + let n_d = cv.wall[d].abs() / a_w; + if n_d < 1e-12 { + return h; + } + (cv.distance / n_d).min(h) + } + + /// 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)) + } +} 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 3bb9991..0a9539c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -312,135 +312,6 @@ impl Mask { self.compute_merging(Some(old)); } - /// The shift of a face's open-part centroid from the face centre: - /// `½h(1 − α)` along the wall normal's in-plane part, away from the - /// body (zero for a full face or without the centroid diffusion). Read - /// from the tables of [`Self::compute_face_shifts`]. - pub(super) fn face_shift(&self, c: usize, p: [i64; 3]) -> [f64; 3] { - let (Some(t), Some(f)) = (self.face_shifts.as_ref(), self.lattice().face(c, p)) else { - return [0.0; 3]; - }; - [t[c][3 * f], t[c][3 * f + 1], t[c][3 * f + 2]] - } - - /// The per-face shift tables (three components interleaved). - #[must_use] - pub fn face_shift_tables(&self) -> Option<&[Vec; 3]> { - self.face_shifts.as_ref() - } - - /// Build the open-part centroid shifts of every cut face (S2-5). - pub fn compute_face_shifts(&mut self) { - let g = self.grid; - let h = [g.dx, g.dy, g.dz]; - let lat = self.lattice(); - let sizes = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()]; - let mut tables = [ - vec![0.0; 3 * sizes[0]], - vec![0.0; 3 * sizes[1]], - vec![0.0; 3 * sizes[2]], - ]; - for c in 0..3 { - let (ni, nj, nk) = ( - g.nx + usize::from(c == 0), - g.ny + usize::from(c == 1), - g.nz + usize::from(c == 2), - ); - for k in 0..nk { - for j in 0..nj { - for i in 0..ni { - let p = [i as i64, j as i64, k as i64]; - let Some(f) = lat.face(c, p) else { continue }; - let Some(alpha) = self.aperture(c, p) else { - continue; - }; - if alpha <= 0.0 || alpha >= 1.0 { - continue; - } - // Interior faces only (a control volume needs both cells). - let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c]; - if on_side && !(c == 2 && self.periodic_z) { - continue; - } - let cv = self.cv_geometry(c, p); - let mut n = cv.wall; - n[c] = 0.0; - let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); - if a == 0.0 { - continue; - } - for d in 0..3 { - // `wall` points into the body: the open part lies the other way. - tables[c][3 * f + d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a; - } - } - } - } - } - self.face_shifts = Some(tables); - } - - /// The distance over which a fluid face exchanges momentum with a solid - /// neighbour face along `d`: the full spacing, or (S2-5) the axis - /// distance from the open part's centroid to the wall, `min(h, d_f/|n_d|)`. - pub(super) fn exchange_delta(&self, cv: &CvGeometry, d: usize) -> f64 { - let h = [self.grid.dx, self.grid.dy, self.grid.dz][d]; - if !self.wall_exchange_axis { - return h; - } - 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 h; - } - let n_d = cv.wall[d].abs() / a_w; - if n_d < 1e-12 { - return h; - } - (cv.distance / n_d).min(h) - } - - /// 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, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs index 1633330..8171147 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -7,6 +7,7 @@ //! u faces on `(nx + 1)·ny·nz`, v on `nx·(ny + 1)·nz`, w on `nx·ny·(nz + 1)`. pub mod body; +pub mod closure; pub mod cut; pub mod cutwall; pub mod exchange; 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 1b77f9b..3c60eab 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 @@ -121,7 +121,7 @@ pub struct Parameters { /// wall's in-plane normal) instead of `h`. With `h` the coupling of a /// cut face to its neighbour is weak by `(1 + α)/2` and the no-slip /// surface sits `½(1 − α) h` inside the body (the flat-wall instrument). - /// `RTX_E3_DIFFUSION_CENTROID=1`. + /// Default ON; `RTX_E3_DIFFUSION_CENTROID=0` restores the old spacing. pub diffusion_centroid: bool, } @@ -147,7 +147,9 @@ impl Default for Parameters { .is_ok_and(|v| v == "oblique"), 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"), - diffusion_centroid: std::env::var("RTX_E3_DIFFUSION_CENTROID").is_ok_and(|v| v == "1"), + // ON by default since S2-5 (`=0` reproduces the records before it). + diffusion_centroid: std::env::var("RTX_E3_DIFFUSION_CENTROID") + .map_or(true, |v| v != "0"), } } } diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1_host.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1_host.rs index 1f9d269..08e360c 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1_host.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1_host.rs @@ -108,7 +108,7 @@ fn dfg_2d1_on_the_host() { let fcv = mask.control_volume_force_with_walls(&field, dt, RHO, mu, None, bx, false); let cd = |f: [f64; 3]| coef * f[0]; println!( - " t {t:7.3}: c_D operator {:.4} (p {:.4} + s {:.4} + exchange {:.4}) | reconstructed {:.4} (p {:.4} + s {:.4}) | box {:.4} c_L {:.5} Δp {:.5}; residual {:.1e} [{:.0} s]", + " t {t:7.3}: c_D operator {:.4} (p {:.4} + s {:.4} + exchange {:.4}) | reconstructed {:.4} (p {:.4} + s {:.4}) | box {:.4} c_L {:.5} Δp {:.5} (quadratic probes {:.5}); residual {:.1e} [{:.0} s]", cd(po) + cd(so) + cd(ex), cd(po), cd(so), @@ -123,6 +123,18 @@ fn dfg_2d1_on_the_host() { - mask .pressure_at(&field.p, CX + 0.5 * D, CY, 0.5 * lz) .unwrap_or(f64::NAN), + { + // Δp by quadratic extrapolation along the normal from probes + // at 1.5 h, 2.5 h, 3.5 h (whole-fluid stencils only). + let wall_p = |sign: f64| { + let at = |d: f64| { + mask.pressure_at(&field.p, CX + sign * (0.5 * D + d * h), CY, 0.5 * lz) + .unwrap_or(f64::NAN) + }; + 4.375 * at(1.5) - 5.25 * at(2.5) + 1.875 * at(3.5) + }; + wall_p(-1.0) - wall_p(1.0) + }, r.final_residual, start.elapsed().as_secs_f64() ); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_position.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_position.rs index 6b50010..87d5fb5 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_wall_position.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position.rs @@ -53,10 +53,8 @@ fn offset(ny: usize, theta: f64) -> (f64, Vec<(f64, f64, f64)>) { for k in 0..=nz { for j in 0..ny { for i in 0..nx { - if k < nz || true { - let idx = g.wface(k.min(nz), j, i); - field.w[idx] = exact((j as f64 + 0.5) * h); - } + let idx = g.wface(k, j, i); + field.w[idx] = exact((j as f64 + 0.5) * h); } } } @@ -100,3 +98,24 @@ fn flat_wall_effective_position() { } } } + +/// The gate: with the centroid diffusion (the default) the effective wall +/// position is second order — within 0.08 h at ny 16, halving in units of +/// h at ny 32, the same at every cut fraction (before S2-5: −½θ h at +/// every resolution). +#[test] +fn flat_wall_position_is_second_order() { + if std::env::var("RTX_E3_DIFFUSION_CENTROID").is_ok_and(|v| v == "0") { + return; + } + for theta in [0.25, 0.75] { + let (coarse, _) = offset(16, theta); + let (fine, _) = offset(32, theta); + println!(" θ {theta}: offset {coarse:+.4} h at ny 16, {fine:+.4} h at ny 32"); + assert!(coarse.abs() < 0.08, "ny 16 offset {coarse}"); + assert!( + fine.abs() < 0.6 * coarse.abs(), + "the offset does not halve: {coarse} → {fine}" + ); + } +}