From 2c94ae0bb14c1b146d48516467fb7e9bc296a14a Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 18:47:46 -0500 Subject: [PATCH] embedded3 S2-1: the traction sampler on the cut field (pressure fit excludes merged cells; sampler route in the sphere MMS and DFG drivers) Co-Authored-By: Claude Fable 5.1 --- .../solvers/incompressible/embedded3/loads.rs | 7 +++-- .../rtx-cfd/tests/embedded3_dfg_2z.rs | 27 +++++++++++++++---- .../rtx-cfd/tests/embedded3_embedded_mms.rs | 16 ++++++++++- .../rtx-cfd/tests/embedded3_sphere/mod.rs | 10 +++++++ 4 files changed, 52 insertions(+), 8 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 5d51829..b75a3b2 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs @@ -83,7 +83,10 @@ impl Mask { } } } - if nodes.iter().all(|&(idx, _)| self.is_fluid_cell(idx)) { + // A virtually merged small cell carries its master's pressure + // unknown, not its own: it is dropped from the fit. + let usable = |idx: usize| self.is_fluid_cell(idx) && self.master(idx).is_none(); + if nodes.iter().all(|&(idx, _)| usable(idx)) { return Some(nodes.iter().map(|&(idx, w)| w * p[idx]).sum()); } // Weighted by the z-direction weight (a z-invariant field then fits @@ -96,7 +99,7 @@ impl Mask { for _ in 0..dirs(ny, fy).len() * dirs(nx, fx).len() { let &(idx, _) = it.next().expect("node"); let _ = dk; - if self.is_fluid_cell(idx) { + if usable(idx) { let (k, j, i) = g.kji(idx); pts.push(( (i as f64 + 0.5) * dx, diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs index a33f82a..18aeaad 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs @@ -75,7 +75,12 @@ fn dfg_3d_2z_on_the_device() { (0.0, 0.0, 0.0) } }); - solver.set_body(Body::cylinder_z(CX, CY, 0.5 * D)); + // The cylinder extruded across the width, with samples for the + // traction route (S2-1). + solver.set_body(Body::extruded( + rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D), + H, + )); let g = Grid::cubic(nx, ny, nz, h); let mut field = Field::new(g); // Start from the inflow profile everywhere (a faster approach to steady). @@ -99,13 +104,18 @@ fn dfg_3d_2z_on_the_device() { let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H); let csv = std::env::var("RTX_E3_DFG_CSV").ok().map(|p| { let mut f = std::fs::File::create(p).expect("csv"); - writeln!(f, "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg").unwrap(); + writeln!( + f, + "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg,cd_sampler,cl_sampler" + ) + .unwrap(); f }); let mut csv = csv; let sample_every = (steps / 100).max(1); let start = std::time::Instant::now(); let mut last: Option<(f64, f64, f64, f64, f64)> = None; + let mut last_sampler = (f64::NAN, f64::NAN); let mut settled = false; for step in 0..steps { let r = device.advance(dt); @@ -131,6 +141,8 @@ fn dfg_3d_2z_on_the_device() { ); let fcv = mask.control_volume_force_with_walls(&field, dt, RHO, RHO * NU, None, bx, true); + let fs = mask.surface_force(body, &field, RHO * NU, t, 0.5 * h); + let (cd_s, cl_s) = (coef * fs.f[0], coef * fs.f[1]); let zc = 0.5 * H; let p_front = mask .pressure_at(&field.p, CX - 0.5 * D, CY, zc) @@ -141,7 +153,8 @@ fn dfg_3d_2z_on_the_device() { 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]); println!( - " t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}) c_L {cl:.5} (CV {cl_cv:.5}) Δp {dp:.4} residual {:.1e} CG {} [{:.0} s]", + " t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}, sampler {cd_s:.4} skipped {}) c_L {cl:.5} (CV {cl_cv:.5}, sampler {cl_s:.5}) Δp {dp:.4} residual {:.1e} CG {} [{:.0} s]", + fs.skipped, r.final_residual, r.poisson_iterations, start.elapsed().as_secs_f64() @@ -149,11 +162,12 @@ fn dfg_3d_2z_on_the_device() { if let Some(f) = csv.as_mut() { writeln!( f, - "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{}", + "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{},{cd_s:.6},{cl_s:.6}", r.final_residual, r.poisson_iterations ) .unwrap(); } + last_sampler = (cd_s, cl_s); if let Some((pcd, pcl, _, _, pdp)) = last { let rel = ((cd - pcd) / cd) .abs() @@ -180,8 +194,11 @@ fn dfg_3d_2z_on_the_device() { } let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples"); println!( - " FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart) c_L {cl:.5} (CV {cl_cv:.5}) Δp {dp:.4} — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s", + " FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; sampler {:.4}, {:.2e} from CV) c_L {cl:.5} (CV {cl_cv:.5}, sampler {:.5}) Δp {dp:.4} — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s", ((cd - cd_cv) / cd).abs(), + last_sampler.0, + ((last_sampler.0 - cd_cv) / cd_cv).abs(), + last_sampler.1, start.elapsed().as_secs_f64() ); if let Some(t) = device.timers() { diff --git a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs index 78d2822..f58feb7 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs @@ -54,8 +54,13 @@ fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder { mm.force_cv[1] - fcv[1], mm.force_cv[2] - fcv[2], ]) / f_scale; + let a = norm([ + mm.force_sampler[0] - fe[0], + mm.force_sampler[1] - fe[1], + mm.force_sampler[2] - fe[2], + ]) / f_scale; println!( - " n = {n:3} L2 u {:.4e} (order {rate}) max div {:.2e} ghost corr {:.2e} F_surface {:.4?} rel {s:.3e} (skipped {}) F_cv {:.4?} rel {c:.3e}", + " n = {n:3} L2 u {:.4e} (order {rate}) max div {:.2e} ghost corr {:.2e} F_surface {:.4?} rel {s:.3e} (skipped {}) F_cv {:.4?} rel {c:.3e} F_sampler rel {a:.3e}", mm.l2_velocity, mm.max_div, mm.ghost_correction, @@ -63,6 +68,15 @@ fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder { mm.skipped, mm.force_cv ); + if scheme == WallScheme::CutCell { + // S2-1 (reported, not gated): the traction sampler on the cut + // field against the box route — 1.15× at n 24 on the first read + // (`docs/embedded3_campaign.md`, S2-1). + println!( + " sampler / box error ratio at n = {n}: {:.3}", + a / c.max(1e-300) + ); + } se.push(s); ce.push(c); } diff --git a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs index 433d593..841f380 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs @@ -107,9 +107,14 @@ pub struct Measurement { pub l2_velocity: f64, pub max_div: f64, pub ghost_correction: f64, + /// The scheme's wall route (the traction sampler on the binary wall, + /// the operator route on the cut wall). pub force_surface: [f64; 3], pub skipped: usize, pub force_cv: [f64; 3], + /// The traction sampler on the cut wall (S2-1); the wall route again + /// on the binary wall. + pub force_sampler: [f64; 3], } /// March the manufactured solution with the sphere at `c` to steady state on grid `n`. @@ -250,6 +255,10 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement let (i0, i1) = (n / 8, n - n / 8); let src = |x: f64, y: f64, z: f64| source3(x, y, z); let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1)); + let force_sampler = match scheme { + WallScheme::GhostBinary => surface.f, + WallScheme::CutCell => mask.surface_force(body, &f, MU, t, 0.5 * h).f, + }; Measurement { l2_velocity: (sq / vol).sqrt(), max_div, @@ -257,5 +266,6 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement force_surface: surface.f, skipped: surface.skipped, force_cv, + force_sampler, } }