From 7fffbb1301d04f8cf40a73240e57a832b233178c Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sun, 20 Sep 2026 08:48:41 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20S2-7b:=20the=20curved=20instrument'?= =?UTF-8?q?s=20per-wall=20operator=20probe=20(curved=5Foperator=5Fprobe)?= =?UTF-8?q?=20and=20RTX=5FE3=5FCURVED=5FSCHEME;=20host=20prototype=20RTX?= =?UTF-8?q?=5FE3=5FWALL=5FFOOT=3Dcentroid=20(the=20shear=20closure's=20wal?= =?UTF-8?q?l=20velocity=20at=20the=20foot=20from=20the=20open=20part's=20C?= =?UTF-8?q?ENTROID,=20not=20the=20face=20centre=20=E2=80=94=20on=20a=20rot?= =?UTF-8?q?ating=20wall=20the=20two=20feet=20are=20=C2=BDh(1=E2=88=92?= =?UTF-8?q?=CE=B1)=20apart:=20rigid-mode=20probe=20residual=204.5=C3=97=20?= =?UTF-8?q?smaller,=20the=20rotating=20wall's=20offset=20=E2=88=920.022=20?= =?UTF-8?q?=E2=86=92=20=C2=B10.002=20h,=20its=20cut=20pressure=20zeroth=20?= =?UTF-8?q?=E2=86=92=20first=20order);=20default=20off,=20device=20refuses?= =?UTF-8?q?=20it?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/cutwall.rs | 1 + .../embedded3/step/cut_predictor.rs | 15 +- .../incompressible/embedded3/step/device.rs | 3 +- .../incompressible/embedded3/step/mod.rs | 9 + .../solvers/incompressible/embedded3/wall.rs | 2 + .../tests/embedded3_wall_position_curved.rs | 173 +++++++++++++++++- 6 files changed, 197 insertions(+), 6 deletions(-) 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 e7db24b..0b49b42 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -224,6 +224,7 @@ impl Mask { wall_exchange_foot: false, conv_sides_exact: false, wall_flux_true_normal: false, + wall_foot_centroid: false, grad_weights: None, diffusion_centroid: false, face_shifts: None, 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 4c1c2dd..c6d3cb8 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 @@ -120,12 +120,21 @@ impl Solver { let ec = e(c); let cell_minus = add(p, ec, -1); let cell_plus = p; - let ub = mask.surface_velocity_at(body, x, c, t_old); + let shift0 = mask.face_shift(c, p); + // S2-7b: the wall velocity at the foot from the open part's centroid. + let foot_of = |q: [i64; 3], xq: [f64; 3]| -> [f64; 3] { + if mask.wall_foot_centroid { + let sq = mask.face_shift(c, q); + [xq[0] + sq[0], xq[1] + sq[1], xq[2] + sq[2]] + } else { + xq + } + }; + let ub = mask.surface_velocity_at(body, foot_of(p, x), c, t_old); let mut mass_out = 0.0; let mut conv = 0.0; let mut diff = 0.0; - let shift0 = mask.face_shift(c, p); // The implicit exchange with solid neighbour faces (S2-5). let (mut wall_implicit, mut wall_rhs) = (0.0, 0.0); for d in 0..3 { @@ -251,7 +260,7 @@ impl Solver { 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 ubq = mask.surface_velocity_at(body, foot_of(q, lat.face_position(c, q)), c, t_old); // 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 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 2802033..492fa0b 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 @@ -180,7 +180,8 @@ impl DeviceStep { assert!( !solver.params.wall_exchange_foot && !solver.params.conv_sides_exact - && !solver.params.wall_flux_true_normal, + && !solver.params.wall_flux_true_normal + && !solver.params.wall_foot_centroid, "the axis-foot exchange / exact convective sides (S2-7b) are host prototypes: the device kernels do not carry them" ); assert!( 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 d057df6..d73f582 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 @@ -181,6 +181,13 @@ pub struct Parameters { /// of a rotating or flexing wall. `RTX_E3_WALL_FLUX=true`; the device /// refuses it. pub wall_flux_true_normal: bool, + /// HOST PROTOTYPE (S2-7b, moving walls): the shear closure's wall velocity + /// at the foot of the normal from the face's open-part CENTROID (where + /// its value lives) instead of from the face centre — the two feet are + /// ½h(1 − α) apart along the wall, an O(Ωh) velocity error over d_f on a + /// rotating or flexing wall (zero on a wall at rest or in translation). + /// `RTX_E3_WALL_FOOT=centroid`; the device refuses it. + pub wall_foot_centroid: bool, /// S2-5: the cross-direction diffusion between two faces over the /// distance between their OPEN-PART CENTROIDS (a cut face's velocity /// is its open part's mean, ½h(1 − α) off the face centre along the @@ -226,6 +233,7 @@ impl Default for Parameters { wall_exchange_foot: std::env::var("RTX_E3_WALL_EXCHANGE").is_ok_and(|v| v == "axisfoot"), conv_sides_exact: std::env::var("RTX_E3_CONV_SIDES").is_ok_and(|v| v == "exact"), wall_flux_true_normal: std::env::var("RTX_E3_WALL_FLUX").is_ok_and(|v| v == "true"), + wall_foot_centroid: std::env::var("RTX_E3_WALL_FOOT").is_ok_and(|v| v == "centroid"), pressure_centroid: std::env::var("RTX_E3_PRESSURE_CENTROID").is_ok_and(|v| v == "1"), momentum_volume_tiled: std::env::var("RTX_E3_MOMENTUM_VOLUME") .is_ok_and(|v| v == "tiled"), @@ -375,6 +383,7 @@ impl Solver { m.wall_exchange_foot = self.params.wall_exchange_foot; m.conv_sides_exact = self.params.conv_sides_exact; m.wall_flux_true_normal = self.params.wall_flux_true_normal; + m.wall_foot_centroid = self.params.wall_foot_centroid; 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 0bcdd7c..f0e0aae 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -119,6 +119,7 @@ pub struct Mask { pub(super) wall_exchange_foot: bool, pub(super) conv_sides_exact: bool, pub(super) wall_flux_true_normal: bool, + pub(super) wall_foot_centroid: bool, /// The centroid prototype's pressure-gradient weights per u / v / w face. /// The centroid-distance cross diffusion (S2-5). pub(super) diffusion_centroid: bool, @@ -550,6 +551,7 @@ impl Mask { wall_exchange_foot: false, conv_sides_exact: false, wall_flux_true_normal: false, + wall_foot_centroid: false, grad_weights: None, diffusion_centroid: false, face_shifts: None, diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs index a3ebbaa..83726f8 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs @@ -10,7 +10,7 @@ //! (`u = Ω r e_θ`, `p = ρ Ω² r²/2`). use rtx_cfd::solvers::incompressible::ConvectionScheme; use rtx_cfd::solvers::incompressible::embedded3::{ - Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, + Body, Boundaries, FaceKind, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, }; const MU: f64 = 0.1; @@ -77,7 +77,11 @@ fn parameters() -> Parameters { Parameters { corrector_steps: 2, tolerance: 1e-10, - convection_scheme: ConvectionScheme::TvdVanAlbada, + convection_scheme: if std::env::var("RTX_E3_CURVED_SCHEME").is_ok_and(|v| v == "upwind") { + ConvectionScheme::Upwind + } else { + ConvectionScheme::TvdVanAlbada + }, wall_scheme: WallScheme::CutCell, boundaries: Boundaries { z0: Side::Periodic, @@ -326,3 +330,168 @@ fn curved_wall_effective_position_and_pressure() { reading(n, rigid); } } + +/// S2-7b B1: the operator probe on this instrument — one predictor and one +/// corrector from the exact field at the faces' open-part centroids (exact +/// p at the cells); the predictor's residual per face, recovered from the +/// one correction, in units of the exact field's largest acceleration +/// (Ω² R2), by aperture band and by WALL (inner / outer); the correction's +/// pressure at each wall's cut cells (of ρ(ΩR1)²). +fn probe(n: usize, rigid: bool) { + let h = 1.0 / n as f64; + let (nx, ny, nz) = ((LX * n as f64) as usize, (LX * n as f64) as usize, 2); + let ex = Exact::new(rigid); + let mut params = parameters(); + params.corrector_steps = 1; + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: OMEGA * R1, + reference_length: R2 - R1, + }, + params, + ); + let exb = Exact::new(rigid); + solver.set_boundary_velocity(move |x, y, _z, _t| { + let (u, v) = exb.velocity(x, y, rigid); + (u, v, 0.0) + }); + solver.set_body( + Body::from_sdf(move |x, y, _z, _t| { + let r = ((x - CENTRE.0).powi(2) + (y - CENTRE.1).powi(2)).sqrt(); + (r - R1).min(R2 - r) + }) + .with_surface_velocity(move |x, y, _z, _t| { + let (dx, dy) = (x - CENTRE.0, y - CENTRE.1); + let r = (dx * dx + dy * dy).sqrt().max(1e-12); + let inner_side = r < 0.5 * (R1 + R2); + let moving = rigid || (inner_side != outer_drives()); + let ut = if moving { OMEGA * r } else { 0.0 }; + (-ut * dy / r, ut * dx / r, 0.0) + }), + ); + let g = Grid::cubic(nx, ny, nz, h); + let mut field = Field::new(g); + solver.initialize(&mut field); + let tables = solver.mask().expect("mask").face_shift_tables().expect("shift tables").clone(); + for k in 0..nz { + for j in 0..ny { + for i in 0..=nx { + let f = g.uface(k, j, i); + let t = &tables[0][3 * f..3 * f + 3]; + field.u[f] = ex.velocity(i as f64 * h + t[0], (j as f64 + 0.5) * h + t[1], rigid).0; + } + } + for j in 0..=ny { + for i in 0..nx { + let f = g.vface(k, j, i); + let t = &tables[1][3 * f..3 * f + 3]; + field.v[f] = ex.velocity((i as f64 + 0.5) * h + t[0], j as f64 * h + t[1], rigid).1; + } + } + } + let r_of = |x: f64, y: f64| ((x - CENTRE.0).powi(2) + (y - CENTRE.1).powi(2)).sqrt(); + for idx in 0..g.cells() { + let (_, j, i) = g.kji(idx); + let r = r_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h).clamp(R1, R2); + field.p[idx] = ex.p(r); + } + { + let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask")); + mask.impose(body, &mut field.u, &mut field.v, &mut field.w, solver.time()); + } + let dt = 0.5 * h * h / (6.0 * MU); + solver.advance(&mut field, dt); + let mask = solver.mask().expect("mask"); + let pp = &field.p_prime; + let a_max = OMEGA * OMEGA * R2; + let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9; + // [wall][band]: bands α<¼, ¼–½, ½–¾, ¾–1, full next to cut. + let mut acc: [[Vec; 5]; 2] = Default::default(); + let bin_of = |a: f64, near: bool| -> Option { + if a < 1.0 { + Some(((a * 4.0).floor() as usize).min(3)) + } else if near { + Some(4) + } else { + None + } + }; + let wall_of = |x: f64, y: f64| usize::from(r_of(x, y) >= 0.5 * (R1 + R2)); + for j in 0..ny { + for i in 1..nx { + let f = g.uface(0, j, i); + if mask.u_kind(f) != FaceKind::Fluid { + continue; + } + let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i)); + let Some(b) = bin_of(mask.a_u(f), is_cut(cm) || is_cut(cp)) else { continue }; + let star = field.u[f] + (dt / RHO) * mask.grad_weight(0, f) * (pp[cp] - pp[cm]) / h; + acc[wall_of(i as f64 * h, (j as f64 + 0.5) * h)][b].push((star - field.u_old[f]) / dt / a_max); + } + } + for j in 1..ny { + for i in 0..nx { + let f = g.vface(0, j, i); + if mask.v_kind(f) != FaceKind::Fluid { + continue; + } + let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i)); + let Some(b) = bin_of(mask.a_v(f), is_cut(cm) || is_cut(cp)) else { continue }; + let star = field.v[f] + (dt / RHO) * mask.grad_weight(1, f) * (pp[cp] - pp[cm]) / h; + acc[wall_of((i as f64 + 0.5) * h, j as f64 * h)][b].push((star - field.v_old[f]) / dt / a_max); + } + } + let scale = RHO * (OMEGA * R1).powi(2); + let (mut sum, mut cnt) = (0.0, 0usize); + for j in 0..ny { + for i in 0..nx { + let c = g.cell(0, j, i); + if mask.cell_active(c) && !is_cut(c) && mask.master(c).is_none() { + sum += pp[c]; + cnt += 1; + } + } + } + let lvl = sum / cnt.max(1) as f64; + let mut psq = [0.0f64; 3]; + let mut pn = [0usize; 3]; + for j in 0..ny { + for i in 0..nx { + let c = g.cell(0, j, i); + if !mask.cell_active(c) || mask.master(c).is_some() { + continue; + } + let e = (pp[c] - lvl) / scale; + let w = if is_cut(c) { wall_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h) } else { 2 }; + psq[w] += e * e; + pn[w] += 1; + } + } + let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::() / v.len().max(1) as f64).sqrt(); + let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full next to cut"]; + let mode = if rigid { "rigid" } else if outer_drives() { "outer-driven" } else { "couette" }; + for (w, wname) in ["inner wall", "outer wall"].iter().enumerate() { + let mut line = format!(" probe {mode} n {n} {wname}:"); + for (b, name) in names.iter().enumerate() { + line += &format!(" {name} {} rms {:.3e};", acc[w][b].len(), rms(&acc[w][b])); + } + line += &format!(" p' at its cut cells {:.3e} ({})", (psq[w] / pn[w].max(1) as f64).sqrt(), pn[w]); + println!("{line}"); + } + println!(" probe {mode} n {n} interior p' {:.3e} ({})", (psq[2] / pn[2].max(1) as f64).sqrt(), pn[2]); +} + +#[test] +#[ignore = "S2-7b B1 probe: the discrete operator on the exact Taylor–Couette field (seconds per rung)"] +fn curved_operator_probe() { + let ns: Vec = std::env::var("RTX_E3_CURVED_NS") + .ok() + .map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect()) + .unwrap_or_else(|| vec![16, 32, 64]); + let rigid = std::env::var("RTX_E3_CURVED_MODE").is_ok_and(|v| v == "rigid"); + for n in ns { + probe(n, rigid); + } +}