From 9ae0e42dcd3ebf01162a4af8bc31b109481bb09f Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 12:58:30 -0500 Subject: [PATCH] embedded3 S2-3d: face-momentum unsteady term for the box (Mask::box_momentum_rates, Solver::previous_apertures); falsifier prints the residual against it Co-Authored-By: Claude Fable 5.1 --- .../solvers/incompressible/embedded3/loads.rs | 83 +++++++++++++++++++ .../incompressible/embedded3/step/mod.rs | 10 +++ .../incompressible/embedded3/step/moving.rs | 3 + .../rtx-cfd/tests/embedded3_falsifier.rs | 16 +++- 4 files changed, 111 insertions(+), 1 deletion(-) 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 11bdac4..0f7e5b1 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs @@ -485,3 +485,86 @@ impl Mask { force } } + +impl Mask { + /// S2-3d: the rate of change of the momentum the SCHEME carries in a + /// box — on the face control volumes, `Σ_f ρ h A (α_f u_f − α_fⁿ u_fⁿ)/dt` + /// (faces on the box's own surfaces at half weight) — against the + /// cell-volume form of the box route's unsteady term, returned second. + /// Near an accelerating wall the two differ at O(h) in the cut cells. + #[must_use] + pub fn box_momentum_rates( + &self, + f: &Field, + dt: f64, + rho: f64, + (i0, i1, j0, j1, k0, k1): (usize, usize, usize, usize, usize, usize), + apertures_old: Option<&[Vec; 3]>, + vol_old: Option<&[f64]>, + ) -> ([f64; 3], [f64; 3]) { + let g = self.grid(); + let dv = g.dx * g.dy * g.dz; + let mut faces = [0.0; 3]; + let mut cells = [0.0; 3]; + let old = |c: usize, idx: usize, new: f64| apertures_old.map_or(new, |a| a[c][idx]); + for k in k0..k1 { + for j in j0..j1 { + for i in i0..=i1 { + let idx = g.uface(k, j, i); + let w = if i == i0 || i == i1 { 0.5 } else { 1.0 }; + let a = self.a_u(idx); + faces[0] += w * rho * dv * (a * f.u[idx] - old(0, idx, a) * f.u_old[idx]) / dt; + } + } + for j in j0..=j1 { + for i in i0..i1 { + let idx = g.vface(k, j, i); + let w = if j == j0 || j == j1 { 0.5 } else { 1.0 }; + let a = self.a_v(idx); + faces[1] += w * rho * dv * (a * f.v[idx] - old(1, idx, a) * f.v_old[idx]) / dt; + } + } + } + for k in k0..=k1 { + for j in j0..j1 { + for i in i0..i1 { + let kk = if k == g.nz && self.periodic_z() { 0 } else { k }; + let idx = g.wface(kk, j, i); + let w = if (k == k0 || k == k1) && !self.periodic_z() { + 0.5 + } else if k == k1 { + 0.0 + } else { + 1.0 + }; + let a = self.a_w(idx); + faces[2] += w * rho * dv * (a * f.w[idx] - old(2, idx, a) * f.w_old[idx]) / dt; + } + } + } + for k in k0..k1 { + for j in j0..j1 { + for i in i0..i1 { + let idx = g.cell(k, j, i); + let v_new = if self.is_fluid_cell(idx) { + self.vol(idx) + } else { + 0.0 + }; + let v_old = vol_old.map_or(v_new, |vo| vo[idx]); + let c = |a: &[f64], f0: usize, f1: usize| 0.5 * (a[f0] + a[f1]); + let (fu0, fu1) = (g.uface(k, j, i), g.uface(k, j, i + 1)); + let (fv0, fv1) = (g.vface(k, j, i), g.vface(k, j + 1, i)); + let (fw0, fw1) = (g.wface(k, j, i), g.wface(k + 1, j, i)); + cells[0] += + rho * dv * (v_new * c(&f.u, fu0, fu1) - v_old * c(&f.u_old, fu0, fu1)) / dt; + cells[1] += + rho * dv * (v_new * c(&f.v, fv0, fv1) - v_old * c(&f.v_old, fv0, fv1)) / dt; + cells[2] += + rho * dv * (v_new * c(&f.w, fw0, fw1) - v_old * c(&f.w_old, fw0, fw1)) / dt; + } + } + } + (faces, cells) + } +} 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 3e3b151..954935c 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 @@ -183,6 +183,8 @@ pub struct Solver { /// The previous mask's cell volumes (moving bodies; the box route's /// unsteady term). pub(super) vol_old: Vec, + /// The previous mask's face apertures (moving bodies; S2-3d). + pub(super) apertures_old: Option<[Vec; 3]>, boundary_velocity: Option, body: Option, /// The body moves: the mask is rebuilt at every step's new time. @@ -216,6 +218,7 @@ impl Solver { floor_source: std::cell::Cell::new([0.0; 3]), pressure_lag: std::cell::Cell::new([0.0; 3]), vol_old: Vec::new(), + apertures_old: None, boundary_velocity: None, body: None, moving: false, @@ -310,6 +313,13 @@ impl Solver { self.pressure_lag.get() } + /// The previous step's cell volumes (fraction of the cell; 0 for a + /// solid cell), once a moving body's mask has been rebuilt. + #[must_use] + pub fn previous_apertures(&self) -> Option<&[Vec; 3]> { + self.apertures_old.as_ref() + } + /// The previous step's cell volumes (fraction of the cell; 0 for a /// solid cell), once a moving body's mask has been rebuilt. #[must_use] diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs index 9acfefd..dec66d8 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs @@ -70,6 +70,9 @@ impl Solver { self.last_ghost_correction = correction; } if let Some(old) = &self.mask { + self.apertures_old = old + .cut() + .map(|c| [c.a_u.clone(), c.a_v.clone(), c.a_w.clone()]); self.vol_old = (0..field.grid.cells()) .map(|i| { if old.is_fluid_cell(i) { diff --git a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs index 5404855..fe1f3c5 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs @@ -210,6 +210,8 @@ struct Record { fy_lag: f64, /// The force the momentum equation applied (shear on u*, exchange on uⁿ). fy_applied: f64, + /// The box route with the scheme's own (face control volume) momentum rate. + fy_cv_faces: f64, fresh: usize, skipped: usize, p_far: f64, @@ -320,6 +322,15 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { false, solver.previous_volumes(), )[1] / lz; + let (rate_faces, rate_cells) = mask.box_momentum_rates( + &field, + dt, + RHO, + (margin, N - margin, margin, N - margin, 0, nz), + solver.previous_apertures(), + solver.previous_volumes(), + ); + let fy_cv_faces = fy_cv + (rate_cells[1] - rate_faces[1]) / lz; let p_far = field.p[g.cell(kp, jp, ip)]; let mut ke = 0.0; for k in 0..nz { @@ -357,6 +368,7 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { fy_floor, fy_lag, fy_applied, + fy_cv_faces, fresh: result.fresh_cells, skipped, p_far, @@ -370,13 +382,15 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { }; let mean = |f: &dyn Fn(&Record) -> f64| records.iter().map(f).sum::() / n; println!( - " route residual (wall − box) RMS {:.3e} → with the floor source {:.3e}, with the closure lag {:.3e}, both {:.3e}, APPLIED form {:.3e} (mean {:+.3e}); means wall {:+.3e} box {:+.3e} floor {:+.3e} lag {:+.3e}; box RMS {:.3e}", + " route residual (wall − box) RMS {:.3e} → with the floor source {:.3e}, with the closure lag {:.3e}, both {:.3e}, APPLIED form {:.3e} (mean {:+.3e}), against the FACE-momentum box {:.3e} (its mean {:+.3e}); means wall {:+.3e} box {:+.3e} floor {:+.3e} lag {:+.3e}; box RMS {:.3e}", rms(&|r| r.fy - r.fy_cv), rms(&|r| r.fy + r.fy_floor - r.fy_cv), rms(&|r| r.fy + r.fy_lag - r.fy_cv), rms(&|r| r.fy + r.fy_floor + r.fy_lag - r.fy_cv), rms(&|r| r.fy_applied - r.fy_cv), mean(&|r| r.fy_applied), + rms(&|r| r.fy - r.fy_cv_faces), + mean(&|r| r.fy_cv_faces), mean(&|r| r.fy), mean(&|r| r.fy_cv), mean(&|r| r.fy_floor),