From e9b2e7887b428a7fbc58fa39b4d567fefe2f27b3 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 16:22:38 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20S2-6:=20the=20oblique-wall=20instru?= =?UTF-8?q?ment=20(tests/embedded3=5Fwall=5Fposition=5Foblique.rs:=20z-flo?= =?UTF-8?q?w=20/=20in-plane=20Poiseuille=20+=20Couette=20linear=20exactnes?= =?UTF-8?q?s;=20gate=20oblique=5Fwall=5Fposition=5Fis=5Fsecond=5Forder)=20?= =?UTF-8?q?and=20the=20closures=20it=20found,=20host=20+=20device,=20defau?= =?UTF-8?q?lt=20OFF=20=E2=80=94=20the=20transverse=20centroid=20correction?= =?UTF-8?q?=20(RTX=5FE3=5FDIFFUSION=5FTRANSVERSE=3D1,=20wall=5Forder=20bit?= =?UTF-8?q?=207),=20the=20fine=20distance=20floor=20(RTX=5FE3=5FDISTANCE?= =?UTF-8?q?=5FFLOOR=3Dfine,=20bit=208);=20kernel=20geometry=20factored=20i?= =?UTF-8?q?nto=20cut=5Fcv;=20DFG=202D-1=20x-shift=20knob=20(RTX=5FE3=5FDFG?= =?UTF-8?q?=5FSHIFT=5FX);=20FLAG=5FX0=200.6=20->=200.25=20in=20the=20three?= =?UTF-8?q?=20flag=20tests=20(the=20flag=20of=20every=20record=20was=20det?= =?UTF-8?q?ached)?= 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 | 91 ++++- .../incompressible/embedded3/closure.rs | 9 + .../incompressible/embedded3/cutwall.rs | 7 +- .../embedded3/step/cut_predictor.rs | 48 +++ .../incompressible/embedded3/step/device.rs | 9 +- .../incompressible/embedded3/step/mod.rs | 18 + .../solvers/incompressible/embedded3/wall.rs | 5 + .../rtx-cfd/tests/embedded3_dfg_2d1.rs | 20 +- .../rtx-cfd/tests/embedded3_flag_geometry.rs | 6 +- .../tests/embedded3_flag_reference_2d.rs | 6 +- .../rtx-cfd/tests/embedded3_flag_wake.rs | 6 +- .../tests/embedded3_wall_position_oblique.rs | 378 ++++++++++++++++++ 12 files changed, 579 insertions(+), 24 deletions(-) create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs 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 1504303..43b16d0 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -15,6 +15,7 @@ */ #define CUT_INERTIA_FLOOR 0.1 #define CUT_DISTANCE_FLOOR 0.05 +#define CUT_DISTANCE_FLOOR_FINE 0.01 struct E3Cut { const double *a_u, *a_v, *a_w; /* apertures per face */ @@ -61,20 +62,15 @@ __device__ __forceinline__ double cut_ap_at(const E3Params& g, const E3Cut& m, i return fidx < 0 ? dflt : cut_ap(m, cc)[fidx]; } -/* The predicted value of the open face of component c at (i, j, k). */ -__device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cut& m, int c, int i, int j, int k, int fidx) +/* The control volume's geometry of the face of component c at (i, j, k) + (cutwall.rs cv_geometry): side apertures, the closing wall vector, the + floored wall distance. */ +__device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, int k, int fidx, + double* alpha_out, double* apm, double* app, double* wall, double* distance_out) { double h[3] = { g.dx, g.dy, g.dz }; - double n[3] = { (double)g.nx, (double)g.ny, (double)g.nz }; double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy }; - double rho = g.rho, mu = g.nu * g.rho; - int sides[3][2] = { { g.bx0, g.bx1 }, { g.by0, g.by1 }, { g.bz0, g.bz1 } }; - int p[3] = { i, j, k }; - const double* old_c = cut_old(f, c); - double u0 = old_c[fidx]; - /* the control volume's geometry (cutwall.rs cv_geometry) */ double alpha = cut_ap(m, c)[fidx]; - double apm[3], app[3]; int cm[3] = { i, j, k }; cm[c] -= 1; /* cell minus */ int cp[3] = { i, j, k }; /* cell plus */ for (int d = 0; d < 3; ++d) { @@ -92,7 +88,6 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu app[d] = 0.5 * (cut_ap_at(g, m, d, q1[0], q1[1], q1[2], 1.0) + cut_ap_at(g, m, d, q2[0], q2[1], q2[2], 1.0)); } } - double wall[3]; for (int d = 0; d < 3; ++d) wall[d] = -(app[d] - apm[d]) * area[d]; double h_min = fmin(fmin(g.dx, g.dy), g.dz); const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w); @@ -102,7 +97,64 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]); if (a_w0 > 0.0) { double n_c = wall[c] / a_w0; n_t = sqrt(fmax(1.0 - n_c * n_c, 0.0)); } } - double distance = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha) * n_t, CUT_DISTANCE_FLOOR * h_min); + /* bit 8 of wall_order: the fine distance floor (S2-6) */ + double floor_h = (g.wall_order & 256) ? CUT_DISTANCE_FLOOR_FINE : CUT_DISTANCE_FLOOR; + *alpha_out = alpha; + *distance_out = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha) * n_t, floor_h * h_min); +} + +/* S2-6 (bit 7 of wall_order): the transverse part of the centroid diffusion + between the face (fidx, its geometry given) and its open neighbour fq at + q along d — grad(u) . ds_perp with the cut faces' own wall-normal + gradients (u - U_b)/d_f (cut_predictor.rs `transverse`). */ +__device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d, int fidx, int fq, const int* q, + double sign, double u0, double ub, double alpha, const double* wall, double distance, + double uq, const double* sh, const double* ubt) +{ + double ds[3] = { 0.0, 0.0, 0.0 }; + int any = 0; + for (int e = 0; e < 3; ++e) { + if (e == d) continue; + ds[e] = sign * (sh[3 * fq + e] - sh[3 * fidx + e]); + if (ds[e] != 0.0) any = 1; + } + if (!any) return 0.0; + double alpha_q, apm_q[3], app_q[3], wall_q[3], distance_q; + cut_cv(g, m, c, q[0], q[1], q[2], fq, &alpha_q, apm_q, app_q, wall_q, &distance_q); + double gr[3] = { 0.0, 0.0, 0.0 }; + double count = 0.0; + double a0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]); + if (a0 > 0.0 && alpha < 1.0) { + double slope = (u0 - ub) / (distance * a0); + for (int e = 0; e < 3; ++e) gr[e] -= slope * wall[e]; + count += 1.0; + } + double aq = sqrt(wall_q[0] * wall_q[0] + wall_q[1] * wall_q[1] + wall_q[2] * wall_q[2]); + if (aq > 0.0 && alpha_q < 1.0) { + double slope = (uq - ubt[fq]) / (distance_q * aq); + for (int e = 0; e < 3; ++e) gr[e] -= slope * wall_q[e]; + count += 1.0; + } + if (count == 0.0) return 0.0; + return (gr[0] * ds[0] + gr[1] * ds[1] + gr[2] * ds[2]) / count; +} + +/* The predicted value of the open face of component c at (i, j, k). */ +__device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cut& m, int c, int i, int j, int k, int fidx) +{ + double h[3] = { g.dx, g.dy, g.dz }; + double n[3] = { (double)g.nx, (double)g.ny, (double)g.nz }; + double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy }; + double rho = g.rho, mu = g.nu * g.rho; + int sides[3][2] = { { g.bx0, g.bx1 }, { g.by0, g.by1 }, { g.bz0, g.bz1 } }; + int p[3] = { i, j, k }; + const double* old_c = cut_old(f, c); + double u0 = old_c[fidx]; + /* the control volume's geometry (cutwall.rs cv_geometry) */ + double alpha, apm[3], app[3], wall[3], distance; + cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &distance); + int cm[3] = { i, j, k }; cm[c] -= 1; /* cell minus */ + int cp[3] = { i, j, k }; /* cell plus */ const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w); double ub = ubt[fidx]; /* face position */ @@ -194,6 +246,21 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu the part beyond 1/h point-implicit (S2-5) */ int centroid = (g.wall_order & 64) && d != c; const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w); + /* bit 7 of wall_order: the transverse centroid correction, cross and own direction (S2-6) */ + if ((g.wall_order & 128) && (g.wall_order & 64)) { + if (f_up1 >= 0 && !solid_up && cut_ap(m, c)[f_up1] > 0.0) { + int qn[3] = { i, j, k }; qn[d] += 1; + double dl = h[d]; + if (d != c) { dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]); dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]); } + diff -= mu * g_plus * a_d * cut_transverse(g, m, c, d, fidx, f_up1, qn, 1.0, u0, ub, alpha, wall, distance, up1, sh, ubt) / dl; + } + if (f_dn1 >= 0 && !solid_dn && cut_ap(m, c)[f_dn1] > 0.0) { + int qn[3] = { i, j, k }; qn[d] -= 1; + double dl = h[d]; + if (d != c) { dl = h[d] - (sh[3 * f_dn1 + d] - sh[3 * fidx + d]); dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]); } + diff += mu * g_minus * a_d * cut_transverse(g, m, c, d, fidx, f_dn1, qn, -1.0, u0, ub, alpha, wall, distance, dn1, sh, ubt) / dl; + } + } if (solid_up) { double kx = mu * g_plus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * up1; } else if (f_up1 >= 0 && centroid) { double dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]); 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 d37d2b6..4e6e23a 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs @@ -17,6 +17,15 @@ impl Mask { [t[c][3 * f], t[c][3 * f + 1], t[c][3 * f + 2]] } + /// The wall-distance floor in units of the smallest spacing (S2-6). + pub(super) fn distance_floor(&self) -> f64 { + if self.distance_floor_fine { + super::cutwall::DISTANCE_FLOOR_FINE + } else { + super::cutwall::DISTANCE_FLOOR + } + } + /// The per-face shift tables (three components interleaved). #[must_use] pub fn face_shift_tables(&self) -> Option<&[Vec; 3]> { 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 0a9539c..52f66dd 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -25,6 +25,8 @@ use super::wall::{FaceKind, Mask}; pub(super) const INERTIA_FLOOR: f64 = 0.1; /// The wall-distance floor of a face, in units of the smallest spacing. pub(super) const DISTANCE_FLOOR: f64 = 0.05; +/// The fine floor (S2-6, `Parameters::distance_floor_fine`). +pub(super) const DISTANCE_FLOOR_FINE: f64 = 0.01; /// Virtual merging: a cell whose fluid fraction (at either end of the /// step) stays below this shares its pressure unknown with a neighbour. pub(super) const MERGE_FRACTION: f64 = 0.1; @@ -211,6 +213,8 @@ impl Mask { density: 1.0, wall_order: 1, wall_distance_oblique: false, + diffusion_transverse: false, + distance_floor_fine: false, wall_exchange_axis: false, grad_weights: None, diffusion_centroid: false, @@ -398,7 +402,8 @@ impl Mask { } else { 1.0 }; - let distance = (phi_face + 0.5 * h[c] * (1.0 - alpha) * n_t).max(DISTANCE_FLOOR * h_min); + let distance = + (phi_face + 0.5 * h[c] * (1.0 - alpha) * n_t).max(self.distance_floor() * h_min); CvGeometry { alpha, ap, 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 07ff6d3..c37b000 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 @@ -197,6 +197,54 @@ impl Solver { delta.clamp(0.25 * h[d], 2.0 * h[d]) }; let solid = |q: [i64; 3]| mask.wall_exchange_axis && ap(c, q) == Some(0.0); + // S2-6: the centroids of two neighbours with different apertures + // are offset ACROSS the difference's axis (in the cross and in + // the own direction); the two-point difference then carries + // `∇u · Δs_⊥`, removed with the cut faces' own wall-normal + // gradients `(u − U_b)/d_f` (old values: explicit). + let transverse = |q: [i64; 3], uq: f64, sign: f64| -> f64 { + let sq = mask.face_shift(c, q); + let mut ds = [0.0; 3]; + for e2 in 0..3 { + if e2 != d { + ds[e2] = sign * (sq[e2] - shift0[e2]); + } + } + if ds == [0.0; 3] { + return 0.0; + } + let cvq = mask.cv_geometry(c, q); + let ubq = mask.surface_velocity_at(body, lat.face_position(c, q), c, t_old); + let mut g = [0.0; 3]; + let mut count = 0.0; + for (cvk, uk, ubk) in [(&cv, u0, ub), (&cvq, uq, ubq)] { + let a = (cvk.wall[0] * cvk.wall[0] + + cvk.wall[1] * cvk.wall[1] + + cvk.wall[2] * cvk.wall[2]) + .sqrt(); + if a > 0.0 && cvk.alpha < 1.0 { + let slope = (uk - ubk) / (cvk.distance * a); + for e2 in 0..3 { + g[e2] -= slope * cvk.wall[e2]; + } + count += 1.0; + } + } + if count == 0.0 { + return 0.0; + } + (g[0] * ds[0] + g[1] * ds[1] + g[2] * ds[2]) / count + }; + if mask.diffusion_transverse { + for (nb, sign, gap) in [(up1, 1.0, g_plus), (dn1, -1.0, g_minus)] { + let q = add(p, ed, sign as i64); + let Some(uq) = nb else { continue }; + if solid(q) || !ap(c, q).is_some_and(|a| a > 0.0) { + continue; + } + diff -= sign * mu * gap * a_d * transverse(q, uq, sign) / spacing(q, sign); + } + } diff += match up1 { Some(un) if solid(add(p, ed, 1)) => { let k = mu * g_plus * a_d / mask.exchange_delta(&cv, d); 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 671a0d7..d82f2a3 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 @@ -292,7 +292,14 @@ impl DeviceStep { wall_order: i32::from(self.solver.params.wall_order) + 16 * i32::from(self.solver.params.wall_distance_oblique) + 32 * i32::from(self.solver.params.wall_exchange_axis) - + 64 * i32::from(self.solver.params.diffusion_centroid), + + 64 * i32::from(self.solver.params.diffusion_centroid) + // bit 7: the transverse centroid correction; bit 8: the fine floor (S2-6). + + 128 + * i32::from( + self.solver.params.diffusion_transverse + && self.solver.params.diffusion_centroid, + ) + + 256 * i32::from(self.solver.params.distance_floor_fine), 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 954935c..4ca8d88 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 @@ -102,6 +102,18 @@ pub struct Parameters { /// by `½h(1 − α)(1 − |n_t|)` at every h (S2-5). `Parameters::default()` /// reads `RTX_E3_WALL_DISTANCE=oblique` (default: the recorded form). pub wall_distance_oblique: bool, + /// The transverse part of the centroid diffusion (S2-6): two neighbour + /// faces with different apertures have open-part centroids offset + /// ACROSS the axis of their difference, so the two-point difference + /// carries `∇u · Δs_⊥` — a first-order wall position on oblique walls + /// (≈ 0.1 h inside the body). Removed with the faces' own wall-normal + /// gradients `(u − U_b)/d_f`, in the cross AND the own direction. + /// Needs `diffusion_centroid`. `RTX_E3_DIFFUSION_TRANSVERSE=1`. + pub diffusion_transverse: bool, + /// The wall-distance floor at 0.01 of the smallest spacing instead of + /// 0.05: the coarse floor doubles the distance of faces with α < 0.1 + /// (the flat wall's θ 0.95 excess, S2-6). `RTX_E3_DISTANCE_FLOOR=fine`. + pub distance_floor_fine: 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 @@ -145,6 +157,9 @@ impl Default for Parameters { .unwrap_or(1), wall_distance_oblique: std::env::var("RTX_E3_WALL_DISTANCE") .is_ok_and(|v| v == "oblique"), + diffusion_transverse: std::env::var("RTX_E3_DIFFUSION_TRANSVERSE") + .is_ok_and(|v| v == "1"), + distance_floor_fine: std::env::var("RTX_E3_DISTANCE_FLOOR").is_ok_and(|v| v == "fine"), 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). @@ -281,6 +296,9 @@ impl Solver { m.wall_order = self.params.wall_order; m.wall_distance_oblique = self.params.wall_distance_oblique; m.wall_exchange_axis = self.params.wall_exchange_axis; + m.diffusion_transverse = + self.params.diffusion_transverse && self.params.diffusion_centroid; + m.distance_floor_fine = self.params.distance_floor_fine; 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 80869b6..a14bdd6 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -101,6 +101,9 @@ pub struct Mask { pub(super) wall_order: u8, /// The oblique wall distance of the cut faces (S2-5). pub(super) wall_distance_oblique: bool, + /// The transverse centroid correction and the fine distance floor (S2-6). + pub(super) diffusion_transverse: bool, + pub(super) distance_floor_fine: 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. @@ -524,6 +527,8 @@ impl Mask { density: 1.0, wall_order: 1, wall_distance_oblique: false, + diffusion_transverse: false, + distance_floor_fine: false, wall_exchange_axis: false, grad_weights: None, diffusion_centroid: false, diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs index 3b82694..b072fd0 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs @@ -47,6 +47,10 @@ fn dfg_2d_1_on_the_device() { let nz = NZ; let lz = nz as f64 * h; let t_end = env_f("RTX_E3_DFG_T", 10.0); + // S2-6: the body (and the Δp probes) shifted in x by a fraction of h — + // the lift's grid-registration scatter (a y-shift would change the + // physical lift: the cylinder is 5 mm off the centreline). + let cx = CX + env_f("RTX_E3_DFG_SHIFT_X", 0.0) * h; // Explicit stability: CFL 0.3 on U_m and half the viscous limit. let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU)); let mut solver = Solver::new( @@ -81,7 +85,7 @@ fn dfg_2d_1_on_the_device() { // traction route (S2-1). solver.set_body(Body::extruded( rtx_cfd::solvers::incompressible::EmbeddedBody::circle( - CX, + cx, CY, 0.5 * D + env_f("RTX_E3_DFG_DR", 0.0) * h, ), @@ -101,8 +105,10 @@ fn dfg_2d_1_on_the_device() { solver.initialize(&mut field); let mask_cells = solver.mask().map_or(0, |m| m.fluid_cells()); println!( - " DFG 2D-1 ny {ny}: {nx}×{ny}×{nz} = {} cells ({mask_cells} fluid), h {h:.4e}, dt {dt:.3e}, t_end {t_end}", - g.cells() + " DFG 2D-1 ny {ny}: {nx}×{ny}×{nz} = {} cells ({mask_cells} fluid), h {h:.4e}, dt {dt:.3e}, t_end {t_end}, centre x/h {:.3} y/h {:.3}", + g.cells(), + cx / h, + CY / h ); let mut device = DeviceStep::new(solver, g); device.upload(&field); @@ -138,8 +144,8 @@ fn dfg_2d_1_on_the_device() { let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2); let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2); let bx = ( - ci(CX - margin), - ci(CX + margin), + ci(cx - margin), + ci(cx + margin), cj(CY - 0.15), cj(CY + 0.15), 0, @@ -159,10 +165,10 @@ fn dfg_2d_1_on_the_device() { let (cd_s, cl_s) = (coef * fs.f[0], coef * fs.f[1]); let zc = 0.5 * lz; let p_front = mask - .pressure_at(&field.p, CX - 0.5 * D, CY, zc) + .pressure_at(&field.p, cx - 0.5 * D, CY, zc) .unwrap_or(f64::NAN); let p_back = mask - .pressure_at(&field.p, CX + 0.5 * D, CY, zc) + .pressure_at(&field.p, cx + 0.5 * D, CY, zc) .unwrap_or(f64::NAN); let dp = p_front - p_back; let (cd, cl, cd_cv, cl_cv) = (coef * fw[0], coef * fw[1], coef * fcv[0], coef * fcv[1]); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs index c777f89..6449b48 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs @@ -4,7 +4,11 @@ use rtx_cfd::solvers::incompressible::embedded3::{Body, Boundaries, Grid, Mask, Side}; const H: f64 = 0.41; -const FLAG_X0: f64 = 0.6; +/// The flag's ROOT: the cylinder's rear (the Turek–Hron flag runs from the +/// cylinder to its tip A at x = 0.6). Until 2026-09-18 this was 0.6 — the +/// flag sat DETACHED, its root where the benchmark's tip is; every flag +/// record before that date is of that geometry. +const FLAG_X0: f64 = 0.25; const FLAG_LEN: f64 = 0.35; const FLAG_HALF: f64 = 0.01; const FLAG_SPAN: f64 = 0.2; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs index 24c3ef2..c7bcee5 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs @@ -19,7 +19,11 @@ const L: f64 = 2.5; const CX: f64 = 0.2; const CY: f64 = 0.2; const R_CYL: f64 = 0.05; -const FLAG_X0: f64 = 0.6; +/// The flag's ROOT: the cylinder's rear (the Turek–Hron flag runs from the +/// cylinder to its tip A at x = 0.6). Until 2026-09-18 this was 0.6 — the +/// flag sat DETACHED, its root where the benchmark's tip is; every flag +/// record before that date is of that geometry. +const FLAG_X0: f64 = 0.25; const FLAG_LEN: f64 = 0.35; const FLAG_HALF: f64 = 0.01; const AMP: f64 = 0.084; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index 3590779..9adba35 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -31,7 +31,11 @@ const L: f64 = 2.5; const CX: f64 = 0.2; const CY: f64 = 0.2; const R_CYL: f64 = 0.05; -const FLAG_X0: f64 = 0.6; +/// The flag's ROOT: the cylinder's rear (the Turek–Hron flag runs from the +/// cylinder to its tip A at x = 0.6). Until 2026-09-18 this was 0.6 — the +/// flag sat DETACHED, its root where the benchmark's tip is; every flag +/// record before that date is of that geometry. +const FLAG_X0: f64 = 0.25; const FLAG_LEN: f64 = 0.35; const FLAG_HALF: f64 = 0.01; const FLAG_SPAN: f64 = 0.2; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs new file mode 100644 index 0000000..5ae1227 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs @@ -0,0 +1,378 @@ +//! S2-6 instrument: the cut wall's effective position on an OBLIQUE wall +//! with the flow IN the plane of the cut. Poiseuille flow (body force `F` +//! along the tangent) in a channel between two embedded parallel planes +//! `y = s x + c0` and `y = s x + c0 + w`; the sides carry the exact +//! solution, z is periodic. The exact field is `U(r) t`, `U = F/(2μ) +//! (G²/4 − r²)`, constant pressure, and the 5-point Laplacian is exact on +//! it, so on full faces of a central window `u/t_x + F r²/(2μ)` is the +//! constant `F G_eff²/(8μ)`: the effective gap from the u faces and from +//! the v faces separately, each an effective wall offset per wall in +//! units of h. Unlike the z-directed flat-wall test this one exercises the +//! own-direction coupling of cut faces, the convective terms' cut-face +//! values and the projection next to the wall (the spurious pressure is +//! printed). +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; + +const MU: f64 = 0.1; +const LX: f64 = 2.0; +const LY: f64 = 2.0; +const W: f64 = 0.5; + +/// The cut wall's parameters: the environment's (`None`), or the S2-6 +/// closures forced on / off for the gate. +fn parameters(s26: Option) -> Parameters { + let mut p = Parameters { + corrector_steps: 2, + tolerance: 1e-10, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + ..Parameters::default() + }; + if let Some(on) = s26 { + p.diffusion_centroid = true; + p.wall_distance_oblique = on; + p.diffusion_transverse = on; + p.distance_floor_fine = on; + } + p +} + +struct Reading { + /// Effective wall offset per wall from the fitted profile of the u / v faces. + off_u: f64, + off_v: f64, + /// The driving force the profile's curvature implies, over F. + force_u: f64, + /// 1 − (the fitted streamwise pressure slope)/F: must equal `force_u`. + force_p: f64, + /// RMS of the pressure about its linear fit, over F·G: full cells, cut cells. + p_full: f64, + p_cut: f64, +} + +/// Least squares of `y = a − b x`: returns (a, b). +fn fit(points: &[(f64, f64)]) -> (f64, f64) { + let n = points.len() as f64; + let (sx, sy) = points + .iter() + .fold((0.0, 0.0), |s, p| (s.0 + p.0, s.1 + p.1)); + let (mx, my) = (sx / n, sy / n); + let (sxx, sxy) = points.iter().fold((0.0, 0.0), |s, p| { + (s.0 + (p.0 - mx) * (p.0 - mx), s.1 + (p.0 - mx) * (p.1 - my)) + }); + let slope = sxy / sxx; + (my - slope * mx, -slope) +} + +fn reading(n: usize, slope: f64, c0: f64, along_z: bool, s26: Option) -> Reading { + // The driving force (`RTX_E3_OBLIQUE_F`): the problem is linear in it + // but for the convective terms, so a small value switches them off. + #[allow(non_snake_case)] + let F: f64 = std::env::var("RTX_E3_OBLIQUE_F") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(1.0); + let h = 1.0 / n as f64; + let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2); + let norm = (1.0 + slope * slope).sqrt(); + let (tx, ty) = (1.0 / norm, slope / norm); + let gap = W / norm; + let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm; + let speed = move |x: f64, y: f64| { + let r = r_of(x, y); + if r.abs() < 0.5 * gap { + F / (2.0 * MU) * (0.25 * gap * gap - r * r) + } else { + 0.0 + } + }; + let mut solver = Solver::new( + Fluid { + density: 1.0, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + }, + parameters(s26), + ); + // `along_z`: the same channel with the flow along the periodic z (the + // cross-direction diffusion of w alone: no pressure, no convection). + solver.set_boundary_velocity(move |x, y, _z, _t| { + let s = speed(x, y); + if along_z { + (0.0, 0.0, s) + } else { + (s * tx, s * ty, 0.0) + } + }); + solver.set_momentum_source(move |_, _, _, _| { + if along_z { + (0.0, 0.0, F) + } else { + (F * tx, F * ty, 0.0) + } + }); + solver.set_body(Body::from_sdf(move |x, y, _z, _t| { + 0.5 * gap - r_of(x, y).abs() + })); + let g = Grid::cubic(nx, ny, nz, h); + let mut field = Field::new(g); + for k in 0..=nz { + for j in 0..ny { + for i in 0..nx { + if along_z { + field.w[g.wface(k, j, i)] = speed((i as f64 + 0.5) * h, (j as f64 + 0.5) * h); + } + } + } + } + for k in 0..nz { + if along_z { + break; + } + for j in 0..ny { + for i in 0..=nx { + field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h); + } + } + for j in 0..=ny { + for i in 0..nx { + field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h); + } + } + } + solver.initialize(&mut field); + let dt = 0.5 * h * h / (6.0 * MU); + let steps = (2.0 / dt).ceil() as usize; + for _ in 0..steps { + solver.advance(&mut field, dt); + } + let mask = solver.mask().expect("mask"); + // The sides pin the flow RATE (the exact profile), so a displaced wall + // appears as a streamwise pressure slope: F_eff = F − dp/ds, and on the + // full faces of a central window `u/t_x = F_eff/(2μ) (G_eff²/4 − r²)` + // exactly (the 5-point Laplacian is exact on it). Fit both constants. + let window = |x: f64| (x - 0.5 * LX).abs() < 0.3; + let (mut pu, mut pv) = (Vec::new(), Vec::new()); + let (mut pf, mut pc) = (Vec::new(), Vec::new()); + for j in 0..ny { + for i in 0..nx { + let (xu, yu) = (i as f64 * h, (j as f64 + 0.5) * h); + let fu = g.uface(0, j, i); + let (xw, yw) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h); + let fw = g.wface(0, j, i); + if along_z { + if window(xw) && r_of(xw, yw).abs() < 0.4 * gap && mask.a_w(fw) >= 1.0 { + pu.push((r_of(xw, yw).powi(2), field.w[fw])); + } + } else if window(xu) && r_of(xu, yu).abs() < 0.4 * gap && mask.a_u(fu) >= 1.0 { + pu.push((r_of(xu, yu).powi(2), field.u[fu] / tx)); + } + let (xv, yv) = ((i as f64 + 0.5) * h, j as f64 * h); + let fv = g.vface(0, j, i); + if !along_z + && slope > 0.0 + && window(xv) + && r_of(xv, yv).abs() < 0.4 * gap + && mask.a_v(fv) >= 1.0 + { + pv.push((r_of(xv, yv).powi(2), field.v[fv] / ty)); + } + let c = g.cell(0, j, i); + let (xc, yc) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h); + if window(xc) && mask.cell_active(c) && r_of(xc, yc).abs() < 0.5 * gap + h { + let s_along = xc * tx + yc * ty; + let full = r_of(xc, yc).abs() < 0.5 * gap - 1.5 * h; + if full { + pf.push((s_along, field.p[c])); + } else { + pc.push((s_along, field.p[c])); + } + } + } + } + let profile = |points: &[(f64, f64)]| { + if points.is_empty() { + return (f64::NAN, f64::NAN); + } + let (a, b) = fit(points); + let gap_eff = 2.0 * (a / b).sqrt(); + (0.5 * (gap - gap_eff) / h, 2.0 * MU * b / F) + }; + let (off_u, force_u) = profile(&pu); + let (off_v, _) = profile(&pv); + let (p0, minus_slope) = fit(&pf); + let rms = |points: &[(f64, f64)]| { + (points + .iter() + .map(|(s, p)| (p - (p0 - minus_slope * s)).powi(2)) + .sum::() + / points.len().max(1) as f64) + .sqrt() + / (F * gap) + }; + Reading { + off_u, + off_v, + force_u, + force_p: 1.0 + minus_slope / F, + p_full: rms(&pf), + p_cut: rms(&pc), + } +} + +/// The linear-exactness mode: in-plane Couette flow `u = K dist t` over ONE +/// embedded oblique wall (no force, constant pressure, the sides carry the +/// exact field). A scheme exact on linear fields returns the wall position +/// to round-off; the fitted zero of the profile on full faces is the offset. +fn couette(n: usize, slope: f64, c0: f64, s26: Option) -> (f64, f64) { + const K: f64 = 1.0; + let h = 1.0 / n as f64; + let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2); + let norm = (1.0 + slope * slope).sqrt(); + let (tx, ty) = (1.0 / norm, slope / norm); + let dist = move |x: f64, y: f64| (y - slope * x - c0) / norm; + let speed = move |x: f64, y: f64| K * dist(x, y).max(0.0); + let mut solver = Solver::new( + Fluid { + density: 1.0, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + }, + parameters(s26), + ); + solver.set_boundary_velocity(move |x, y, _z, _t| { + let s = speed(x, y); + (s * tx, s * ty, 0.0) + }); + solver.set_body(Body::from_sdf(move |x, y, _z, _t| dist(x, y))); + let g = Grid::cubic(nx, ny, nz, h); + let mut field = Field::new(g); + for k in 0..nz { + for j in 0..ny { + for i in 0..=nx { + field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h); + } + } + for j in 0..=ny { + for i in 0..nx { + field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h); + } + } + } + solver.initialize(&mut field); + let dt = 0.5 * h * h / (6.0 * MU); + let steps = (2.0 / dt).ceil() as usize; + for _ in 0..steps { + solver.advance(&mut field, dt); + } + let mask = solver.mask().expect("mask"); + // u/t_x = K (dist − δ): fit on full faces of the central window, two to + // six cells off the wall. + let mut points = Vec::new(); + let mut worst: f64 = 0.0; + for j in 0..ny { + for i in 0..nx { + let (x, y) = (i as f64 * h, (j as f64 + 0.5) * h); + let f = g.uface(0, j, i); + let dd = dist(x, y); + if (x - 0.5 * LX).abs() < 0.3 && dd > 2.0 * h && dd < 6.0 * h && mask.a_u(f) >= 1.0 { + points.push((dd, field.u[f] / tx)); + worst = worst.max((field.u[f] / tx - K * dd).abs() / (K * h)); + } + } + } + let (a, minus_b) = fit(&points); + // y = a − (−b) x with b the slope: the zero sits at dist = −a / b. + let b = -minus_b; + (-a / b / h, worst) +} + +#[test] +#[ignore = "S2-6 instrument: linear exactness of the cut wall on an oblique wall (in-plane Couette; a minute on the host)"] +fn oblique_wall_linear_exactness() { + for slope in [0.0, 0.25, 0.5, 1.0] { + for n in [16usize, 32] { + let (off, worst) = couette(n, slope, 0.53, None); + println!( + " couette slope {slope:.2} n {n}: wall offset {off:+.4} h (negative = inside the body); worst full-face error {worst:.4} of K·h" + ); + } + } +} + +#[test] +#[ignore = "S2-6 instrument: the oblique cut wall's effective position with in-plane flow (minutes on the host)"] +fn oblique_wall_effective_position() { + for (slope, c0) in [ + (0.0, 0.53), + (0.0, 0.77), + (0.25, 0.53), + (0.5, 0.53), + (1.0, 0.53), + ] { + // `RTX_E3_OBLIQUE_N=64` adds a finer rung to the in-plane mode. + let extra: Option = std::env::var("RTX_E3_OBLIQUE_N") + .ok() + .and_then(|v| v.parse().ok()); + let mut runs = vec![(16usize, false), (32, false), (16, true), (32, true)]; + if let Some(n) = extra { + runs = vec![(n, false)]; + } + for (n, along_z) in runs { + let r = reading(n, slope, c0, along_z, None); + println!( + " slope {slope:.2} c0 {c0} n {n} {}: wall offset {:+.4} h (u faces) {:+.4} h (v faces), positive = inside the fluid; F_eff/F {:.5} (profile) {:.5} (pressure slope); pressure about its fit: {:.2e} full cells, {:.2e} near-wall cells (of F·G)", + if along_z { + "z-flow (w faces)" + } else { + "in-plane" + }, + r.off_u, + r.off_v, + r.force_u, + r.force_p, + r.p_full, + r.p_cut + ); + } + } +} + +/// The gate (S2-6): with the oblique distance, the transverse centroid +/// correction and the fine floor the cut wall is linear-exact to 0.02 h on +/// an oblique wall (without them it sits 0.05–0.07 h inside the body at +/// every h), and the z-directed Poiseuille offset halves per rung at +/// slope ½ (without them: −0.087 → −0.083 h). +#[test] +fn oblique_wall_position_is_second_order() { + for slope in [0.5, 1.0] { + let (fixed, _) = couette(32, slope, 0.53, Some(true)); + let (before, _) = couette(32, slope, 0.53, Some(false)); + println!( + " couette slope {slope}: offset {fixed:+.4} h with the S2-6 closures, {before:+.4} h without" + ); + assert!(fixed.abs() < 0.02, "slope {slope}: {fixed}"); + assert!( + before.abs() > 2.0 * fixed.abs(), + "the instrument lost its contrast" + ); + } + let coarse = reading(16, 0.5, 0.53, true, Some(true)).off_u; + let fine = reading(32, 0.5, 0.53, true, Some(true)).off_u; + println!(" z-flow slope 0.5: offset {coarse:+.4} h at n 16, {fine:+.4} h at n 32"); + assert!(coarse.abs() < 0.04, "n 16 offset {coarse}"); + assert!( + fine.abs() < 0.65 * coarse.abs(), + "the offset does not halve: {coarse} → {fine}" + ); +}