From c430510802697fdb56cd2708396d4f174758bff6 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sun, 20 Sep 2026 07:34:40 -0500 Subject: [PATCH] embedded3 S2-7b instruments: the quadratic / full-cells-only pressure probe (Mask::pressure_at_quadratic[_from]) and the DFG 2D-1 test's RTX_E3_DFG_DP_PROBE line; the manufactured sphere's pressure-error read by cell class + six signed wall-point reads, RTX_E3_MMS_SCHEME=tvd, the cut_cell_pressure_ladder and sphere_operator_probe tests (defaults untouched; MMS gates green) Co-Authored-By: Claude Fable 5.1 --- .../solvers/incompressible/embedded3/loads.rs | 81 +++++ .../rtx-cfd/tests/embedded3_dfg_2d1.rs | 17 + .../rtx-cfd/tests/embedded3_embedded_mms.rs | 30 +- .../rtx-cfd/tests/embedded3_sphere/mod.rs | 295 +++++++++++++++++- 4 files changed, 421 insertions(+), 2 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs index 0f7e5b1..eaa73eb 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs @@ -124,6 +124,87 @@ impl Mask { linear_fit(&pts, (x, y, z)) } + /// S2-7b: pressure at a point by a QUADRATIC least-squares fit (six + /// coefficients in the plane) through the fluid, unmerged cells whose + /// centres lie within `radius` of the point on the plane nearest `z`; + /// `None` with fewer than eight cells or a singular normal matrix. + pub fn pressure_at_quadratic(&self, p: &[f64], x: f64, y: f64, z: f64, radius: f64) -> Option { + self.pressure_at_quadratic_from(p, x, y, z, radius, false) + } + + /// The quadratic fit through FULL fluid cells only (`full_only`): the + /// cut cells' pressures left out of the read (S2-7b's discriminator). + pub fn pressure_at_quadratic_from( + &self, + p: &[f64], + x: f64, + y: f64, + z: f64, + radius: f64, + full_only: bool, + ) -> Option { + let g = self.grid(); + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let k = ((z / dz - 0.5).round().max(0.0) as usize).min(nz - 1); + let usable = |idx: usize| { + self.is_fluid_cell(idx) + && self.master(idx).is_none() + && (!full_only || self.vol(idx) >= 1.0 - 1e-9) + }; + let (ri, rj) = ((radius / dx).ceil() as i64 + 1, (radius / dy).ceil() as i64 + 1); + let (ic, jc) = ((x / dx - 0.5).round() as i64, (y / dy - 0.5).round() as i64); + // Rows: [1, ξ, η, ξ², ξη, η²] with ξ, η in units of h about the point. + let mut ata = [[0.0f64; 6]; 6]; + let mut atb = [0.0f64; 6]; + let mut count = 0; + for j in (jc - rj).max(0)..=(jc + rj).min(ny as i64 - 1) { + for i in (ic - ri).max(0)..=(ic + ri).min(nx as i64 - 1) { + let idx = g.cell(k, j as usize, i as usize); + if !usable(idx) { + continue; + } + let (xi, eta) = (((i as f64 + 0.5) * dx - x) / dx, ((j as f64 + 0.5) * dy - y) / dy); + if (xi * dx).powi(2) + (eta * dy).powi(2) > radius * radius { + continue; + } + let row = [1.0, xi, eta, xi * xi, xi * eta, eta * eta]; + for a in 0..6 { + for b in 0..6 { + ata[a][b] += row[a] * row[b]; + } + atb[a] += row[a] * p[idx]; + } + count += 1; + } + } + if count < 8 { + return None; + } + // Gaussian elimination with partial pivoting; the value at the point + // is the constant coefficient. + let mut m = [[0.0f64; 7]; 6]; + for a in 0..6 { + m[a][..6].copy_from_slice(&ata[a]); + m[a][6] = atb[a]; + } + for c in 0..6 { + let piv = (c..6).max_by(|&a, &b| m[a][c].abs().partial_cmp(&m[b][c].abs()).unwrap())?; + if m[piv][c].abs() < 1e-12 * count as f64 { + return None; + } + m.swap(c, piv); + for r in 0..6 { + if r != c { + let f = m[r][c] / m[c][c]; + for cc in c..7 { + m[r][cc] -= f * m[c][cc]; + } + } + } + } + Some(m[0][6] / m[0][0]) + } + /// Velocity at a point: trilinear over a component's nodes when all are /// fluid faces, else the fit through the fluid ones and the point's own /// boundary intercept with its surface velocity. diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs index b072fd0..1095d86 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2d1.rs @@ -171,6 +171,23 @@ fn dfg_2d_1_on_the_device() { .pressure_at(&field.p, cx + 0.5 * D, CY, zc) .unwrap_or(f64::NAN); let dp = p_front - p_back; + // S2-7b: the stagnation-point reads by the quadratic fit too. + if std::env::var_os("RTX_E3_DFG_DP_PROBE").is_some() { + let h = field.grid.dx; + let mut line = format!(" Δp probes: linear {dp:.5} (front {p_front:.5} back {p_back:.5})"); + for (r, full) in [(2.5, false), (3.5, false), (4.5, false), (2.5, true), (3.5, true), (4.5, true)] { + let pf = mask.pressure_at_quadratic_from(&field.p, cx - 0.5 * D, CY, zc, r * h, full); + let pb = mask.pressure_at_quadratic_from(&field.p, cx + 0.5 * D, CY, zc, r * h, full); + let tag = if full { "full-only quad" } else { "quad" }; + match (pf, pb) { + (Some(pf), Some(pb)) => { + line += &format!("; {tag} r{r:.1}h {:.5} (front {pf:.5} back {pb:.5})", pf - pb); + } + _ => line += &format!("; {tag} r{r:.1}h n/a"), + } + } + println!("{line}"); + } let (cd, cl, cd_cv, cl_cv) = (coef * fw[0], coef * fw[1], coef * fcv[0], coef * fcv[1]); println!( " t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}, reconstructed {cd_s:.4} skipped {}) c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {cl_s:.5}) Δp {dp:.4} residual {:.1e} CG {} [{:.0} s]", diff --git a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs index 364fc2c..0ea72ba 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs @@ -13,7 +13,7 @@ mod embedded3_sphere; -use embedded3_sphere::{C, Measurement, exact_force_and_flux, measure}; +use embedded3_sphere::{C, Measurement, exact_force_and_flux, measure, operator_probe}; use rtx_cfd::solvers::incompressible::embedded3::WallScheme; fn norm(a: [f64; 3]) -> f64 { @@ -156,3 +156,31 @@ fn embedded_sphere_three_rungs() { fn cut_cell_three_rungs() { compare(&[12, 24, 48], 0.1); } + +/// S2-7b: the cut wall's pressure error by cell class and the wall-point +/// reads on the manufactured sphere (`RTX_E3_MMS_NS=12,24,48`). +#[test] +#[ignore = "S2-7b instrument: pressure error by cell class on the manufactured sphere (minutes to an hour on the host)"] +fn cut_cell_pressure_ladder() { + let ns: Vec = std::env::var("RTX_E3_MMS_NS") + .ok() + .map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect()) + .unwrap_or_else(|| vec![12, 24]); + for n in ns { + let m = measure(n, WallScheme::CutCell, C); + println!(" n {n}: l2 velocity {:.3e}", m.l2_velocity); + } +} + +/// S2-7b: the operator probe on the manufactured sphere (`RTX_E3_MMS_NS`). +#[test] +#[ignore = "S2-7b probe: the discrete operator on the exact manufactured field (seconds per rung)"] +fn sphere_operator_probe() { + let ns: Vec = std::env::var("RTX_E3_MMS_NS") + .ok() + .map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect()) + .unwrap_or_else(|| vec![12, 24, 48]); + for n in ns { + operator_probe(n, C); + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs index 3ef6af6..f10cdc7 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs @@ -131,7 +131,12 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement Parameters { corrector_steps: 2, tolerance: 1e-8, - convection_scheme: ConvectionScheme::Upwind, + // S2-7b: `RTX_E3_MMS_SCHEME=tvd` for a second-order interior. + convection_scheme: if std::env::var("RTX_E3_MMS_SCHEME").is_ok_and(|v| v == "tvd") { + ConvectionScheme::TvdVanAlbada + } else { + ConvectionScheme::Upwind + }, wall_scheme: scheme, momentum_volume_cell_mean: std::env::var("RTX_E3_CELL_MEAN").is_ok_and(|v| v == "1"), ..Parameters::default() @@ -245,6 +250,100 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement " [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}", last.final_residual ); + // S2-7b: the pressure's error against the manufactured p (mean-free + // over the full interior cells) in the full cells and in the cut cells, + // and the six axis wall points' reads by the linear box probe and the + // quadratic fit (r 2.5 h) — in units of the pressure's scale (1). + { + let (mut sum_full, mut n_full) = (0.0, 0usize); + let exact = |idx: usize| { + let (k, j, i) = g.kji(idx); + p3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h) + }; + let is_cut = |idx: usize| mask.vol(idx) < 1.0 - 1e-9; + for idx in 0..g.cells() { + if mask.is_fluid_cell(idx) && !is_cut(idx) && mask.master(idx).is_none() { + sum_full += f.p[idx] - exact(idx); + n_full += 1; + } + } + let level = sum_full / n_full.max(1) as f64; + let (mut sq_full, mut sq_cut, mut n_cut) = (0.0, 0.0, 0usize); + // By fluid fraction (small < 0.5 ≤ large) and the cut cells' mean + // error (a level) against their scatter about it. + let (mut sq_small, mut n_small, mut sq_large, mut n_large, mut sum_cut) = (0.0, 0usize, 0.0, 0usize, 0.0); + let (mut sum_small, mut sum_large) = (0.0, 0.0); + for idx in 0..g.cells() { + if !mask.is_fluid_cell(idx) || mask.master(idx).is_some() { + continue; + } + let e = f.p[idx] - level - exact(idx); + if is_cut(idx) { + sq_cut += e * e; + n_cut += 1; + sum_cut += e; + if mask.vol(idx) < 0.5 { + sq_small += e * e; + sum_small += e; + n_small += 1; + } else { + sq_large += e * e; + sum_large += e; + n_large += 1; + } + } else { + sq_full += e * e; + } + } + let mean_cut = sum_cut / n_cut.max(1) as f64; + println!( + " [{scheme:?} n {n}] cut cells' pressure error: mean {mean_cut:+.3e}, scatter about it {:.3e}; fraction < 0.5: rms {:.3e} mean {:+.3e} ({n_small}), ≥ 0.5: rms {:.3e} mean {:+.3e} ({n_large})", + ((sq_cut / n_cut.max(1) as f64) - mean_cut * mean_cut).max(0.0).sqrt(), + (sq_small / n_small.max(1) as f64).sqrt(), + sum_small / n_small.max(1) as f64, + (sq_large / n_large.max(1) as f64).sqrt(), + sum_large / n_large.max(1) as f64 + ); + // The six axis wall points' signed errors: linear box / quadratic + // r 2.5 h / quadratic through full cells only (r 2.5 h and 3.5 h). + let mut wl = String::new(); + for (dx, dy, dz) in [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0), (0.0, 0.0, -1.0)] { + let (x, y, z) = (c.0 + R * dx, c.1 + R * dy, c.2 + R * dz); + let pe = p3(x, y, z) + level; + let f1 = |v: Option| v.map_or("n/a".to_string(), |v| format!("{:+.4}", v - pe)); + wl += &format!( + " [{}]", + [ + f1(mask.pressure_at(&f.p, x, y, z)), + f1(mask.pressure_at_quadratic(&f.p, x, y, z, 2.5 * h)), + f1(mask.pressure_at_quadratic_from(&f.p, x, y, z, 2.5 * h, true)), + f1(mask.pressure_at_quadratic_from(&f.p, x, y, z, 3.5 * h, true)), + ] + .join(" ") + ); + } + println!(" [{scheme:?} n {n}] wall-point signed errors (linear, quad, full-only quad r2.5, r3.5):{wl}"); + let (mut sq_lin, mut sq_quad, mut n_lin, mut n_quad) = (0.0, 0.0, 0usize, 0usize); + for (dx, dy, dz) in [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0), (0.0, 0.0, -1.0)] { + let (x, y, z) = (c.0 + R * dx, c.1 + R * dy, c.2 + R * dz); + let pe = p3(x, y, z); + if let Some(pl) = mask.pressure_at(&f.p, x, y, z) { + sq_lin += (pl - level - pe).powi(2); + n_lin += 1; + } + if let Some(pq) = mask.pressure_at_quadratic(&f.p, x, y, z, 2.5 * h) { + sq_quad += (pq - level - pe).powi(2); + n_quad += 1; + } + } + println!( + " [{scheme:?} n {n}] pressure error rms: full cells {:.3e} ({n_full}), cut cells {:.3e} ({n_cut}); wall-point reads rms: linear box {:.3e} ({n_lin}/6), quadratic r2.5h {:.3e} ({n_quad}/6)", + (sq_full / n_full.max(1) as f64).sqrt(), + (sq_cut / n_cut.max(1) as f64).sqrt(), + (sq_lin / n_lin.max(1) as f64).sqrt(), + (sq_quad / n_quad.max(1) as f64).sqrt() + ); + } let surface = match scheme { WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h), WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce { @@ -272,3 +371,197 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement force_sampler, } } + +/// S2-7b: the discrete operator applied to the EXACT manufactured field +/// valued at the faces' open-part centroids — one predictor, one +/// corrector; the predictor's acceleration per fluid face (recovered from +/// the one correction) in units of the source's largest acceleration, by +/// aperture band, plus the full faces next to cut cells and the interior; +/// the corrected field's divergence per cut cell over its largest face +/// flux; the correction's pressure at cut cells over the pressure scale. +pub fn operator_probe(n: usize, c: (f64, f64, f64)) { + let h = 1.0 / n as f64; + let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h); + let mut params = Parameters { + corrector_steps: 1, + tolerance: 1e-10, + convection_scheme: if std::env::var("RTX_E3_MMS_SCHEME").is_ok_and(|v| v == "tvd") { + ConvectionScheme::TvdVanAlbada + } else { + ConvectionScheme::Upwind + }, + wall_scheme: WallScheme::CutCell, + ..Parameters::default() + }; + params.momentum_volume_cell_mean = false; + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + }, + params, + ); + solver.set_momentum_source(|x, y, z, _t| source3(x, y, z)); + solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z)); + solver.set_body( + Body::sphere(move |_t| c, R) + .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))), + ); + let g = Grid::cubic(n, n, n, h); + let mut f = Field::new(g); + solver.initialize(&mut f); + let tables = solver + .mask() + .expect("mask") + .face_shift_tables() + .expect("shift tables") + .clone(); + for k in 0..n { + for j in 0..n { + for i in 0..=n { + let idx = g.uface(k, j, i); + let t = &tables[0][3 * idx..3 * idx + 3]; + f.u[idx] = u3(i as f64 * h + t[0], (j as f64 + 0.5) * h + t[1], (k as f64 + 0.5) * h + t[2]); + } + } + for j in 0..=n { + for i in 0..n { + let idx = g.vface(k, j, i); + let t = &tables[1][3 * idx..3 * idx + 3]; + f.v[idx] = v3((i as f64 + 0.5) * h + t[0], j as f64 * h + t[1], (k as f64 + 0.5) * h + t[2]); + } + } + } + for k in 0..=n { + for j in 0..n { + for i in 0..n { + let idx = g.wface(k, j, i); + let t = &tables[2][3 * idx..3 * idx + 3]; + f.w[idx] = w3((i as f64 + 0.5) * h + t[0], (j as f64 + 0.5) * h + t[1], k as f64 * h + t[2]); + } + } + } + for idx in 0..g.cells() { + let (k, j, i) = g.kji(idx); + f.p[idx] = p3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h); + } + { + let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask")); + mask.impose(body, &mut f.u, &mut f.v, &mut f.w, solver.time()); + } + solver.advance(&mut f, dt); + let mask = solver.mask().expect("mask"); + let pp = &f.p_prime; + // The source's largest acceleration (the unit of the read). + let mut a_max: f64 = 0.0; + for k in 0..n { + for j in 0..n { + for i in 0..n { + let s = source3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h); + a_max = a_max.max((s.0 * s.0 + s.1 * s.1 + s.2 * s.2).sqrt() / RHO); + } + } + } + let is_cut = |idx: usize| mask.vol(idx) < 1.0 - 1e-9; + // Bands: α < 1/16, 1/16–1/8, 1/8–1/4, ¼–½, ½–¾, ¾–1, full next to a cut + // cell, interior. Each entry: the acceleration and the FORCE per unit + // `ρ h A` (the acceleration times the floored fraction max(α, 0.1)). + let bin_of = |a: f64, near: bool| -> usize { + if a < 1.0 / 16.0 { + 0 + } else if a < 1.0 / 8.0 { + 1 + } else if a < 0.25 { + 2 + } else if a < 1.0 { + 2 + ((a * 4.0).floor() as usize).min(3) + } else if near { + 6 + } else { + 7 + } + }; + let mut acc: [Vec<(f64, f64)>; 8] = Default::default(); + // u faces + for k in 0..n { + for j in 0..n { + for i in 1..n { + let idx = g.uface(k, j, i); + if mask.u_kind(idx) != FaceKind::Fluid { + continue; + } + let (cm, cp) = (g.cell(k, j, i - 1), g.cell(k, j, i)); + let star = f.u[idx] + (dt / RHO) * mask.grad_weight(0, idx) * (pp[cp] - pp[cm]) / h; + let a = (star - f.u_old[idx]) / dt / a_max; + acc[bin_of(mask.a_u(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_u(idx).max(0.1))); + } + } + } + for k in 0..n { + for j in 1..n { + for i in 0..n { + let idx = g.vface(k, j, i); + if mask.v_kind(idx) != FaceKind::Fluid { + continue; + } + let (cm, cp) = (g.cell(k, j - 1, i), g.cell(k, j, i)); + let star = f.v[idx] + (dt / RHO) * mask.grad_weight(1, idx) * (pp[cp] - pp[cm]) / h; + let a = (star - f.v_old[idx]) / dt / a_max; + acc[bin_of(mask.a_v(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_v(idx).max(0.1))); + } + } + } + for k in 1..n { + for j in 0..n { + for i in 0..n { + let idx = g.wface(k, j, i); + if mask.w_kind(idx) != FaceKind::Fluid { + continue; + } + let (cm, cp) = (g.cell(k - 1, j, i), g.cell(k, j, i)); + let star = f.w[idx] + (dt / RHO) * mask.grad_weight(2, idx) * (pp[cp] - pp[cm]) / h; + let a = (star - f.w_old[idx]) / dt / a_max; + acc[bin_of(mask.a_w(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_w(idx).max(0.1))); + } + } + } + // The correction's pressure at cut cells and interior cells, over the + // pressure scale (1), mean-free over the interior. + let (mut sum_int, mut n_int) = (0.0, 0usize); + for idx in 0..g.cells() { + if mask.is_fluid_cell(idx) && !is_cut(idx) && mask.master(idx).is_none() { + sum_int += pp[idx]; + n_int += 1; + } + } + let lvl = sum_int / n_int.max(1) as f64; + let (mut sq_int, mut sq_cut, mut n_cut) = (0.0, 0.0, 0usize); + for idx in 0..g.cells() { + if !mask.is_fluid_cell(idx) || mask.master(idx).is_some() { + continue; + } + let e = pp[idx] - lvl; + if is_cut(idx) { + sq_cut += e * e; + n_cut += 1; + } else { + sq_int += e * e; + } + } + let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::() / v.len().max(1) as f64).sqrt(); + let names = ["α<1/16", "1/16–1/8", "1/8–1/4", "¼–½", "½–¾", "¾–1", "full next to cut", "interior"]; + let mut line = format!(" probe n {n} (source accel {a_max:.3e}):"); + for (b, name) in names.iter().enumerate() { + let a: Vec = acc[b].iter().map(|x| x.0).collect(); + let fo: Vec = acc[b].iter().map(|x| x.1).collect(); + line += &format!(" {name}: {} acc {:.3e} force {:.3e};", acc[b].len(), rms(&a), rms(&fo)); + } + line += &format!( + " p' rms interior {:.3e} cut {:.3e} ({n_cut})", + (sq_int / n_int.max(1) as f64).sqrt(), + (sq_cut / n_cut.max(1) as f64).sqrt() + ); + println!("{line}"); +}