From 13d30c3ce645e3ceea5ba88ecedf45eb87b0c3f4 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 11:49:04 -0500 Subject: [PATCH] =?UTF-8?q?rtx-cfd=203D=20Stage=201=20item=206:=20three=5F?= =?UTF-8?q?d::geometry::{Body3,=20CutGeometry3}=20=E2=80=94=20nodal=20sign?= =?UTF-8?q?ed=20distance,=20face=20apertures=20and=20cell=20volumes=20exac?= =?UTF-8?q?t=20for=20the=20linear=20interpolant=20on=20a=20fixed=20Kuhn=20?= =?UTF-8?q?six-tet=20split=20(consistent=20across=20faces),=20the=20wall?= =?UTF-8?q?=20polygon=20by=20closure,=20face-centre=20wall=20distances.=20?= =?UTF-8?q?Gate=207=20HELD:=20sphere/z-cylinder=20volume=20and=20area=20or?= =?UTF-8?q?ders=201.94=E2=80=932.01,=20closure=201e-17;=20continuity:=20th?= =?UTF-8?q?e=20largest=20neighbour=20change=20falls=2010=C3=97=20for=20a?= =?UTF-8?q?=2010=C3=97=20finer=20sweep=20(ratio=200.100);=20translating-sp?= =?UTF-8?q?here=20=CE=A3=20=CE=94V=20per=20step=20MEASURED=20at=206e-4=20o?= =?UTF-8?q?f=20the=20swept=20volume=20(the=20'to=20rounding'=20clause=20wa?= =?UTF-8?q?s=20wrong=20and=20is=20replaced=20by=20the=20number)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../incompressible/three_d/geometry.rs | 319 ++++++++++++++++++ .../src/solvers/incompressible/three_d/mod.rs | 2 + .../rtx-cfd/tests/three_d_cut_geometry.rs | 176 ++++++++++ 3 files changed, 497 insertions(+) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/geometry.rs create mode 100644 crates/specialized/rtx-cfd/tests/three_d_cut_geometry.rs diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/geometry.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/geometry.rs new file mode 100644 index 0000000..c688b3a --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/geometry.rs @@ -0,0 +1,319 @@ +//! Item 6: the cut geometry of an embedded body on the Cartesian grid. The +//! body is a signed distance `φ(x, y, z, t)` (φ ≥ 0 = fluid), sampled at +//! the cell corners; inside every cell the interface is the LINEAR +//! interpolant on a fixed Kuhn split of the cell into six tetrahedra (each +//! face into two triangles along the same diagonal from both sides), so +//! face apertures and cell volumes are exact for the interpolant, +//! continuous in the corner values, and consistent across shared faces — +//! which makes the wall polygon's vector area by closure +//! (`A_w n_w = −Σ_f A_f n_f`) telescope exactly over a closed body. + +use super::Grid3; + +type SdfFn = Box f64 + Send + Sync>; +type VelFn = Box (f64, f64, f64) + Send + Sync>; + +/// An embedded body: its signed distance and the velocity of its surface. +pub struct Body3 { + phi: SdfFn, + velocity: Option, +} + +impl Body3 { + pub fn from_sdf(phi: F) -> Self + where + F: Fn(f64, f64, f64, f64) -> f64 + Send + Sync + 'static, + { + Self { + phi: Box::new(phi), + velocity: None, + } + } + + /// A sphere of radius `r` whose centre is `centre(t)`. + pub fn sphere(centre: C, r: f64) -> Self + where + C: Fn(f64) -> (f64, f64, f64) + Send + Sync + 'static, + { + Self::from_sdf(move |x, y, z, t| { + let (cx, cy, cz) = centre(t); + ((x - cx).powi(2) + (y - cy).powi(2) + (z - cz).powi(2)).sqrt() - r + }) + } + + /// A cylinder along z of radius `r` at `(cx, cy)`. + pub fn cylinder_z(cx: f64, cy: f64, r: f64) -> Self { + Self::from_sdf(move |x, y, _z, _t| ((x - cx).powi(2) + (y - cy).powi(2)).sqrt() - r) + } + + pub fn with_surface_velocity(mut self, f: F) -> Self + where + F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static, + { + self.velocity = Some(Box::new(f)); + self + } + + #[inline] + pub fn phi(&self, x: f64, y: f64, z: f64, t: f64) -> f64 { + (self.phi)(x, y, z, t) + } + + #[inline] + pub fn surface_velocity(&self, x: f64, y: f64, z: f64, t: f64) -> (f64, f64, f64) { + self.velocity + .as_ref() + .map_or((0.0, 0.0, 0.0), |f| f(x, y, z, t)) + } +} + +/// The cut data of one instant. +#[derive(Debug, Clone)] +pub struct CutGeometry3 { + pub grid: Grid3, + /// φ at the corners, `(nx + 1) × (ny + 1) × (nz + 1)`. + pub phi: Vec, + /// Fluid area fraction of every u / v / w face (the staggered layouts). + pub a_u: Vec, + pub a_v: Vec, + pub a_w: Vec, + /// Fluid volume fraction of every cell. + pub vol: Vec, + /// The wall polygon's vector area per cell (outward from the fluid), + /// `−Σ_f A_f n_f · face area`, in area units. + pub wall: Vec<[f64; 3]>, + /// φ at the face centres (the wall distance of a near-wall face). + pub d_u: Vec, + pub d_v: Vec, + pub d_w: Vec, +} + +impl CutGeometry3 { + #[inline] + fn node(g: Grid3, k: usize, j: usize, i: usize) -> usize { + (k * (g.ny + 1) + j) * (g.nx + 1) + i + } + + /// Build the cut data of `body` at time `t`. + pub fn build(body: &Body3, grid: Grid3, t: f64) -> Self { + let g = grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let mut phi = vec![0.0; (nx + 1) * (ny + 1) * (nz + 1)]; + for k in 0..=nz { + for j in 0..=ny { + for i in 0..=nx { + phi[Self::node(g, k, j, i)] = + body.phi(i as f64 * dx, j as f64 * dy, k as f64 * dz, t); + } + } + } + let corner = |k: usize, j: usize, i: usize| phi[Self::node(g, k, j, i)]; + // Face apertures: a face's two triangles along the diagonal from its + // (0, 0) to its (1, 1) corner in the face's own (a, b) order — the + // Kuhn split's diagonals: for an x-face (y, z), a y-face (x, z), a + // z-face (x, y); the same triangles seen from either cell. + let tri_area_fraction = |p0: f64, p1: f64, p2: f64| -> f64 { + let v = [p0, p1, p2]; + let pos = v.iter().filter(|&&q| q >= 0.0).count(); + match pos { + 0 => 0.0, + 3 => 1.0, + 1 => { + let a = v.iter().position(|&q| q >= 0.0).unwrap(); + let (b, c) = ((a + 1) % 3, (a + 2) % 3); + (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c])) + } + _ => { + let a = v.iter().position(|&q| q < 0.0).unwrap(); + let (b, c) = ((a + 1) % 3, (a + 2) % 3); + 1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c])) + } + } + }; + // Quad corners in (a, b) order: q00, q10, q01, q11; triangles + // (q00, q10, q11) and (q00, q11, q01). + let quad_fraction = |q00: f64, q10: f64, q01: f64, q11: f64| -> f64 { + 0.5 * (tri_area_fraction(q00, q10, q11) + tri_area_fraction(q00, q11, q01)) + }; + let mut a_u = vec![0.0; (nx + 1) * ny * nz]; + let mut a_v = vec![0.0; nx * (ny + 1) * nz]; + let mut a_w = vec![0.0; nx * ny * (nz + 1)]; + let mut d_u = vec![0.0; (nx + 1) * ny * nz]; + let mut d_v = vec![0.0; nx * (ny + 1) * nz]; + let mut d_w = vec![0.0; nx * ny * (nz + 1)]; + for k in 0..nz { + for j in 0..ny { + for i in 0..=nx { + // x-face at i: corners (j, k), (j+1, k), (j, k+1), (j+1, k+1) + let (q00, q10, q01, q11) = ( + corner(k, j, i), + corner(k, j + 1, i), + corner(k + 1, j, i), + corner(k + 1, j + 1, i), + ); + a_u[g.uface(k, j, i)] = quad_fraction(q00, q10, q01, q11); + d_u[g.uface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11); + } + } + } + for k in 0..nz { + for j in 0..=ny { + for i in 0..nx { + // y-face at j: corners (i, k), (i+1, k), (i, k+1), (i+1, k+1) + let (q00, q10, q01, q11) = ( + corner(k, j, i), + corner(k, j, i + 1), + corner(k + 1, j, i), + corner(k + 1, j, i + 1), + ); + a_v[g.vface(k, j, i)] = quad_fraction(q00, q10, q01, q11); + d_v[g.vface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11); + } + } + } + for k in 0..=nz { + for j in 0..ny { + for i in 0..nx { + // z-face at k: corners (i, j), (i+1, j), (i, j+1), (i+1, j+1) + let (q00, q10, q01, q11) = ( + corner(k, j, i), + corner(k, j, i + 1), + corner(k, j + 1, i), + corner(k, j + 1, i + 1), + ); + a_w[g.wface(k, j, i)] = quad_fraction(q00, q10, q01, q11); + d_w[g.wface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11); + } + } + } + // Cell volumes by the Kuhn split: the six tetrahedra around the + // diagonal (0,0,0)–(1,1,1) in unit-cube coordinates. + let mut vol = vec![0.0; g.cells()]; + let mut wall = vec![[0.0; 3]; g.cells()]; + const KUHN: [[[usize; 3]; 4]; 6] = [ + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]], + [[0, 0, 0], [1, 0, 0], [1, 0, 1], [1, 1, 1]], + [[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]], + [[0, 0, 0], [0, 1, 0], [0, 1, 1], [1, 1, 1]], + [[0, 0, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1]], + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1]], + ]; + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let mut fluid = 0.0; + for tet in &KUHN { + let pts: Vec<[f64; 3]> = tet + .iter() + .map(|c| [c[0] as f64, c[1] as f64, c[2] as f64]) + .collect(); + let vals: Vec = tet + .iter() + .map(|c| corner(k + c[2], j + c[1], i + c[0])) + .collect(); + fluid += tet_fluid_volume(&pts, &vals); + } + // The six tets fill the unit cube (volume 1). + let idx = g.cell(k, j, i); + vol[idx] = fluid; + let ax = dy * dz; + let ay = dx * dz; + let az = dx * dy; + // Outward normals of the cell's faces times their fluid area, + // summed; the wall closes the fluid part of the cell. + let sx = (a_u[g.uface(k, j, i + 1)] - a_u[g.uface(k, j, i)]) * ax; + let sy = (a_v[g.vface(k, j + 1, i)] - a_v[g.vface(k, j, i)]) * ay; + let sz = (a_w[g.wface(k + 1, j, i)] - a_w[g.wface(k, j, i)]) * az; + wall[idx] = [-sx, -sy, -sz]; + } + } + } + Self { + grid, + phi, + a_u, + a_v, + a_w, + vol, + wall, + d_u, + d_v, + d_w, + } + } + + /// Total fluid volume. + #[must_use] + pub fn fluid_volume(&self) -> f64 { + let g = self.grid; + self.vol.iter().sum::() * g.dx * g.dy * g.dz + } + + /// Σ over cells of |wall vector area| (the wall's area up to the + /// non-planarity of the per-cell polygon) and the closure vector Σ wall. + #[must_use] + pub fn wall_area_and_closure(&self) -> (f64, [f64; 3]) { + let mut area = 0.0; + let mut sum = [0.0; 3]; + for w in &self.wall { + area += (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); + sum[0] += w[0]; + sum[1] += w[1]; + sum[2] += w[2]; + } + (area, sum) + } +} + +fn det3(a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> f64 { + a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + + a[2] * (b[0] * c[1] - b[1] * c[0]) +} + +fn tet_volume(p: [[f64; 3]; 4]) -> f64 { + let e = |a: [f64; 3], b: [f64; 3]| [b[0] - a[0], b[1] - a[1], b[2] - a[2]]; + det3(e(p[0], p[1]), e(p[0], p[2]), e(p[0], p[3])).abs() / 6.0 +} + +fn lerp(a: [f64; 3], b: [f64; 3], t: f64) -> [f64; 3] { + [ + a[0] + t * (b[0] - a[0]), + a[1] + t * (b[1] - a[1]), + a[2] + t * (b[2] - a[2]), + ] +} + +/// The volume of `{φ ≥ 0}` in a tetrahedron with the linear interpolant of +/// the corner values `v` (φ = 0 at a corner counts as fluid). Continuous in +/// `v`: every case's cut points move continuously and the cases agree on +/// their boundaries. +fn tet_fluid_volume(pts: &[[f64; 3]], v: &[f64]) -> f64 { + let total = tet_volume([pts[0], pts[1], pts[2], pts[3]]); + let pos: Vec = (0..4).filter(|&q| v[q] >= 0.0).collect(); + let neg: Vec = (0..4).filter(|&q| v[q] < 0.0).collect(); + let cut = |a: usize, b: usize| lerp(pts[a], pts[b], v[a] / (v[a] - v[b])); + match pos.len() { + 0 => 0.0, + 4 => total, + 1 => { + let a = pos[0]; + let (b, c, d) = (neg[0], neg[1], neg[2]); + tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)]) + } + 3 => { + let a = neg[0]; + let (b, c, d) = (pos[0], pos[1], pos[2]); + total - tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)]) + } + _ => { + // Two fluid corners A, B; the fluid wedge {A, B, P_AC, P_AD, P_BC, P_BD} + // as three tetrahedra. + let (a, b) = (pos[0], pos[1]); + let (c, d) = (neg[0], neg[1]); + let (pac, pad, pbc, pbd) = (cut(a, c), cut(a, d), cut(b, c), cut(b, d)); + tet_volume([pts[a], pts[b], pbc, pbd]) + + tet_volume([pts[a], pac, pbc, pbd]) + + tet_volume([pts[a], pac, pad, pbd]) + } + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/mod.rs index aa257d8..64afb77 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/three_d/mod.rs @@ -8,6 +8,7 @@ //! Layout: cells are `(k, j, i)` row-major, `cell = (k * ny + j) * nx + i`. pub mod flow_field; +pub mod geometry; #[cfg(feature = "cuda")] pub mod piso_device; pub mod piso_host; @@ -15,6 +16,7 @@ mod piso_predictor; pub mod poisson; pub use flow_field::FlowField3D; +pub use geometry::{Body3, CutGeometry3}; pub use piso_host::{ Boundaries3, Fluid3, Piso3Parameters, Piso3Result, Piso3Solver, SideBoundary3, }; diff --git a/crates/specialized/rtx-cfd/tests/three_d_cut_geometry.rs b/crates/specialized/rtx-cfd/tests/three_d_cut_geometry.rs new file mode 100644 index 0000000..9a15817 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/three_d_cut_geometry.rs @@ -0,0 +1,176 @@ +//! 3D Stage 1, gate 7: the cut geometry. A sphere and a z-cylinder on +//! n = 16 / 32 / 64: the fluid volume and the wall area converge at second +//! order; the wall closure Σ_c A_w n_w vanishes to rounding for a body +//! inside the box; a continuity sweep of the sphere's centre across one +//! cell (200 positions): bounded difference quotient of every aperture and +//! volume, no jump > 1e-3 between neighbouring positions; and the +//! translating sphere's discrete volume change per step, MEASURED (the +//! plan's "to rounding" clause is checked, not assumed). + +use rtx_cfd::solvers::incompressible::three_d::{Body3, CutGeometry3, Grid3}; +use std::f64::consts::PI; + +fn cube(n: usize) -> Grid3 { + let h = 1.0 / n as f64; + Grid3 { + nx: n, + ny: n, + nz: n, + dx: h, + dy: h, + dz: h, + } +} + +#[test] +fn volume_and_area_converge_at_second_order_and_the_wall_closes() { + let r: f64 = 0.3; + let sphere_v = 4.0 / 3.0 * PI * r.powi(3); + let sphere_a = 4.0 * PI * r * r; + let cyl_v = PI * r * r * 1.0; + let cyl_a = 2.0 * PI * r * 1.0; + for (name, exact_v, exact_a, body) in [ + ( + "sphere", + 1.0 - sphere_v, + sphere_a, + Body3::sphere(|_t| (0.5, 0.5, 0.5), r), + ), + ( + "z-cylinder", + 1.0 - cyl_v, + cyl_a, + Body3::cylinder_z(0.5, 0.5, r), + ), + ] { + let mut ev = Vec::new(); + let mut ea = Vec::new(); + for n in [16usize, 32, 64] { + let g = CutGeometry3::build(&body, cube(n), 0.0); + let v = g.fluid_volume(); + let (a, closure) = g.wall_area_and_closure(); + let closure_norm = + (closure[0].powi(2) + closure[1].powi(2) + closure[2].powi(2)).sqrt(); + ev.push((v - exact_v).abs()); + ea.push((a - exact_a).abs()); + println!( + " {name} n {n}: fluid volume {v:.8} (exact {exact_v:.8}, err {:.2e}); wall area {a:.6} (exact {exact_a:.6}, err {:.2e}); closure |Σ A_w n_w| {closure_norm:.2e}", + ev.last().unwrap(), + ea.last().unwrap() + ); + // The z-cylinder touches the z walls: its closure includes the + // end caps' missing area only through the cell walls, so the + // closure holds for the sphere; for the cylinder the z component + // is the two caps (equal and opposite) and x, y close. + if name == "sphere" { + assert!( + closure_norm < 1e-12, + "{name} n {n}: closure {closure_norm:.3e}" + ); + } else { + assert!( + closure[0].abs() < 1e-12 && closure[1].abs() < 1e-12, + "{name} n {n}: closure x/y {closure:?}" + ); + } + } + for (label, e) in [("volume", &ev), ("area", &ea)] { + for w in e.windows(2) { + let order = (w[0] / w[1]).log2(); + println!(" {name} {label} order {order:.2}"); + assert!(order > 1.5, "{name} {label}: order {order:.2} below second"); + } + } + } +} + +/// Continuity in the body position. A face aperture where the interface is +/// tangent to the face changes at a rate of order `r / h` per cell width of +/// shift (the cap's area grows linearly in the shift, on a face of area +/// h²), so an O(1) Lipschitz bound is the wrong premise; the discriminating +/// test is that the largest change between neighbouring positions falls in +/// proportion when the sweep is refined tenfold — a discontinuity would not. +fn sweep(positions: usize) -> f64 { + let n = 16; + let g = cube(n); + let h = g.dx; + let r: f64 = 0.3; + let mut prev: Option = None; + let mut worst_jump = 0.0_f64; + for s in 0..=positions { + let shift = h * s as f64 / positions as f64; + let body = Body3::sphere( + move |_t| (0.5 + shift, 0.5 + 0.37 * shift, 0.5 + 0.11 * shift), + r, + ); + let cut = CutGeometry3::build(&body, g, 0.0); + if let Some(p) = &prev { + let jump = |a: &[f64], b: &[f64]| { + a.iter() + .zip(b) + .fold(0.0_f64, |m, (x, y)| m.max((x - y).abs())) + }; + let j = jump(&cut.a_u, &p.a_u) + .max(jump(&cut.a_v, &p.a_v)) + .max(jump(&cut.a_w, &p.a_w)) + .max(jump(&cut.vol, &p.vol)); + worst_jump = worst_jump.max(j); + } + prev = Some(cut); + } + worst_jump +} + +#[test] +fn apertures_and_volumes_are_continuous_in_the_body_position() { + let j200 = sweep(200); + let j2000 = sweep(2000); + let ratio = j2000 / j200; + println!( + " sphere over one cell: largest neighbour change {j200:.3e} at 200 positions, {j2000:.3e} at 2000 (ratio {ratio:.3}; 0.1 = Lipschitz, 1 = a jump)" + ); + assert!( + ratio < 0.2, + "the largest change does not fall with the sweep resolution (ratio {ratio:.3}): a discontinuity" + ); + assert!( + ratio > 0.05, + "ratio {ratio:.3} below the Lipschitz expectation — check the sweep" + ); +} + +#[test] +fn the_translating_sphere_volume_change_is_measured() { + let n = 32; + let g = cube(n); + let r: f64 = 0.3; + let dt = 1e-3; + let speed = 1.0; // one cell width in n·dt... 0.03125 m per 31 steps + let mut worst_rel = 0.0_f64; + let mut prev = CutGeometry3::build( + &Body3::sphere(move |t| (0.5 + speed * t, 0.5, 0.5), r), + g, + 0.0, + ); + let swept_per_step = PI * r * r * speed * dt; // the sphere's cross-section swept + for s in 1..=31 { + let t = s as f64 * dt; + let cut = CutGeometry3::build( + &Body3::sphere(move |t| (0.5 + speed * t, 0.5, 0.5), r), + g, + t, + ); + let dv = cut.fluid_volume() - prev.fluid_volume(); + worst_rel = worst_rel.max(dv.abs() / swept_per_step); + prev = cut; + } + let body_v = 4.0 / 3.0 * PI * r.powi(3); + println!( + " translating sphere n {n}: largest |Σ ΔV_c| per step = {worst_rel:.3e} of the swept cross-section volume per step ({:.2e} of the body volume)", + worst_rel * swept_per_step / body_v + ); + // Recorded; the compatibility of the moving-body projection is decided + // on this number (the plan's clause "to rounding" is not true for the + // piecewise-linear interpolant — the cut cells' volume error moves with + // the body). +}