From 9c3fac0755a80725868c77fbb08b8c17f4e49693 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 16:54:52 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20S2-6:=20the=20transverse=20gradient?= =?UTF-8?q?=20from=20the=20face=20at=20least=200.2=20h=20off=20the=20wall?= =?UTF-8?q?=20(full=20neighbours=20over=20their=20own=20distance=20along?= =?UTF-8?q?=20the=20cut=20partner's=20normal)=20=E2=80=94=20the=20first=20?= =?UTF-8?q?form=20blew=20up=20at=20DFG=20ny=2061=20with=20the=20fine=20flo?= =?UTF-8?q?or=20(coefficient=20~=201/d=5Ff);=20host=20=3D=20device=207.7e-?= =?UTF-8?q?14=20with=20the=20closures=20on;=20Couette=20linear=20exactness?= =?UTF-8?q?=20kept=20(<=3D=200.01=20h);=20RTX=5FE3=5FOBLIQUE=5FN=20rung=20?= =?UTF-8?q?knob?= 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 | 30 +++++++++++---- .../incompressible/embedded3/closure.rs | 4 ++ .../embedded3/step/cut_predictor.rs | 38 ++++++++++++++----- 3 files changed, 55 insertions(+), 17 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 43b16d0..68f7146 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -16,6 +16,7 @@ #define CUT_INERTIA_FLOOR 0.1 #define CUT_DISTANCE_FLOOR 0.05 #define CUT_DISTANCE_FLOOR_FINE 0.01 +#define CUT_TRANSVERSE_FLOOR 0.2 struct E3Cut { const double *a_u, *a_v, *a_w; /* apertures per face */ @@ -121,18 +122,31 @@ __device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d 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); + /* Explicit, coefficient ~ 1/d_f: the gradient from the faces at least + CUT_TRANSVERSE_FLOOR h off the wall — a FULL neighbour too, over its own + distance along the cut partner's normal; only when neither is that far, + from the cut faces over the floored distance (cut_predictor.rs). */ + double h_min = fmin(fmin(g.dx, g.dy), g.dz); + double d_min = CUT_TRANSVERSE_FLOOR * h_min; + double a0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]); + double aq = sqrt(wall_q[0] * wall_q[0] + wall_q[1] * wall_q[1] + wall_q[2] * wall_q[2]); + int cut0 = a0 > 0.0 && alpha < 1.0, cutq = aq > 0.0 && alpha_q < 1.0; + if (!cut0 && !cutq) return 0.0; + double n0[3], nq[3]; + for (int e = 0; e < 3; ++e) { n0[e] = cut0 ? wall[e] / a0 : 0.0; nq[e] = cutq ? wall_q[e] / aq : 0.0; } + if (!cut0) for (int e = 0; e < 3; ++e) n0[e] = nq[e]; + if (!cutq) for (int e = 0; e < 3; ++e) nq[e] = n0[e]; + int far = distance >= d_min || distance_q >= d_min; 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]; + if (!((far && distance < d_min) || (!far && !cut0))) { + double slope = (u0 - ub) / fmax(distance, d_min); + for (int e = 0; e < 3; ++e) gr[e] -= slope * n0[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]; + if (!((far && distance_q < d_min) || (!far && !cutq))) { + double slope = (uq - ubt[fq]) / fmax(distance_q, d_min); + for (int e = 0; e < 3; ++e) gr[e] -= slope * nq[e]; count += 1.0; } if (count == 0.0) return 0.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 4e6e23a..dd82190 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs @@ -5,6 +5,10 @@ use super::cutwall::CvGeometry; use super::wall::Mask; +/// The distance floor inside the explicit transverse gradient (S2-6), in +/// units of the smallest spacing. +pub(in crate::solvers::incompressible::embedded3) const TRANSVERSE_DISTANCE_FLOOR: f64 = 0.2; + 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 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 c37b000..8e4e4ca 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 @@ -11,6 +11,8 @@ //! carries the inertia floor. use super::{Side, Solver}; +use crate::solvers::incompressible::embedded3::closure::TRANSVERSE_DISTANCE_FLOOR; +use crate::solvers::incompressible::embedded3::cutwall::CvGeometry; use crate::solvers::incompressible::embedded3::cutwall::INERTIA_FLOOR; use crate::solvers::incompressible::embedded3::field::Field; use crate::solvers::incompressible::simple::ConvectionScheme; @@ -202,6 +204,7 @@ impl Solver { // 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 h_min = h[0].min(h[1]).min(h[2]); let transverse = |q: [i64; 3], uq: f64, sign: f64| -> f64 { let sq = mask.face_shift(c, q); let mut ds = [0.0; 3]; @@ -215,20 +218,37 @@ impl Solver { } 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)] { + // Explicit, with a coefficient ∝ 1/d_f: take the gradient from + // the faces at least TRANSVERSE_DISTANCE_FLOOR h off the wall — + // a FULL neighbour too, over its own distance along the cut + // partner's normal (a linear field gives the same gradient on + // every one of them); only when neither is that far, from the + // cut faces over the floored distance. + let d_min = TRANSVERSE_DISTANCE_FLOOR * h_min; + let unit = |cvk: &CvGeometry| { 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; + (a > 0.0 && cvk.alpha < 1.0) + .then(|| [cvk.wall[0] / a, cvk.wall[1] / a, cvk.wall[2] / a]) + }; + let (n0, nq) = (unit(&cv), unit(&cvq)); + let faces = [(&cv, u0, ub, n0.or(nq)), (&cvq, uq, ubq, nq.or(n0))]; + let far = faces.iter().any(|f| f.3.is_some() && f.0.distance >= d_min); + let mut g = [0.0; 3]; + let mut count = 0.0; + for (cvk, uk, ubk, normal) in faces { + let Some(nk) = normal else { continue }; + let is_cut = unit(cvk).is_some(); + if far && cvk.distance < d_min || !far && !is_cut { + continue; } + let slope = (uk - ubk) / cvk.distance.max(d_min); + for e2 in 0..3 { + g[e2] -= slope * nk[e2]; + } + count += 1.0; } if count == 0.0 { return 0.0;