From 62b46194ddd551c7aa2fc524896458ebf0f605a9 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 04:34:54 -0500 Subject: [PATCH] embedded3: reconstructed route moved to reconstruct.rs (cutwall.rs back under 700 lines) Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/cutwall.rs | 104 +---------------- .../solvers/incompressible/embedded3/mod.rs | 1 + .../incompressible/embedded3/reconstruct.rs | 108 ++++++++++++++++++ 3 files changed, 114 insertions(+), 99 deletions(-) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs 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 12fcf8b..62497b6 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -413,7 +413,11 @@ impl Mask { /// φ and its unit gradient at `x` from the trilinear interpolant of the /// corner values (the cut geometry's own surface). - fn interpolant_distance_and_normal(&self, cut: &CutGeometry, x: [f64; 3]) -> (f64, [f64; 3]) { + pub(super) fn interpolant_distance_and_normal( + &self, + cut: &CutGeometry, + x: [f64; 3], + ) -> (f64, [f64; 3]) { let g = self.grid; let (nx, ny, nz) = (g.nx as i64, g.ny as i64, g.nz as i64); let h = [g.dx, g.dy, g.dz]; @@ -570,104 +574,6 @@ impl Mask { Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]) } - /// The reconstructed wall route (S2-1 remedy): on every wall polygon - /// (the cell's closure `W_c`, its centroid taken as the cell centre's - /// foot on the interpolant surface) the traction from two probes at - /// `h` and `2h` along the interpolant normal — the wall pressure by - /// linear extrapolation, the wall shear from the quadratic fit of the - /// tangential velocity through the probes and the wall velocity (the - /// binary wall's validated sampler, on the cut geometry's own - /// polygons; the probes are a cell away, past the merged slivers). - /// Force on the body: `Σ_c (p_w W_c + μ ∂ₙu_t |W_c| t)` over the - /// planes `k0..k1`; `None` without a cut geometry. - pub fn cut_wall_force_reconstructed( - &self, - body: &Body, - f: &Field, - mu: f64, - t: f64, - planes: Option<(usize, usize)>, - ) -> Option<[f64; 3]> { - let (p, s) = self.cut_wall_force_reconstructed_parts(body, f, mu, t, planes)?; - Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]]) - } - - /// The reconstructed route split into its pressure and shear parts. - pub fn cut_wall_force_reconstructed_parts( - &self, - body: &Body, - f: &Field, - mu: f64, - t: f64, - planes: Option<(usize, usize)>, - ) -> Option<([f64; 3], [f64; 3])> { - let cut = self.cut.as_ref()?; - let g = self.grid; - let (k0, k1) = planes.unwrap_or((0, g.nz)); - let h = g.dx.min(g.dy).min(g.dz); - let (d1, d2) = (h, 2.0 * h); - let mut force = [0.0; 3]; - let mut shear = [0.0; 3]; - for (idx, w) in cut.wall.iter().enumerate() { - let area = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); - if area == 0.0 || !self.cell_fluid[idx] { - continue; - } - let (k, j, i) = g.kji(idx); - if k < k0 || k >= k1 { - continue; - } - let xc = [ - (i as f64 + 0.5) * g.dx, - (j as f64 + 0.5) * g.dy, - (k as f64 + 0.5) * g.dz, - ]; - let (s, n) = self.interpolant_distance_and_normal(cut, xc); - // The wall point and the outward (into the fluid) normal. - let x = [xc[0] - s * n[0], xc[1] - s * n[1], xc[2] - s * n[2]]; - let at = |d: f64| [x[0] + d * n[0], x[1] + d * n[1], x[2] + d * n[2]]; - let (x1, x2) = (at(d1), at(d2)); - let (Some(p1), Some(p2)) = ( - self.pressure_at(&f.p, x1[0], x1[1], x1[2]), - self.pressure_at(&f.p, x2[0], x2[1], x2[2]), - ) else { - // No fit: the operator's own pressure on this polygon. - for c in 0..3 { - force[c] += f.p[idx] * w[c]; - } - continue; - }; - let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1); - for c in 0..3 { - force[c] += p_wall * w[c]; - } - let (Some(u1), Some(u2)) = ( - self.velocity_at(body, f, x1[0], x1[1], x1[2], t), - self.velocity_at(body, f, x2[0], x2[1], x2[2], t), - ) else { - continue; - }; - let us = body.surface_velocity(x[0], x[1], x[2], t); - let us = [us.0, us.1, us.2]; - // Tangential components (the normal removed) and the wall - // gradient of the quadratic through 0, d1, d2. - let tang = |v: [f64; 3]| { - let vn = v[0] * n[0] + v[1] * n[1] + v[2] * n[2]; - [v[0] - vn * n[0], v[1] - vn * n[1], v[2] - vn * n[2]] - }; - let (t1, t2, ts) = (tang(u1), tang(u2), tang(us)); - let wall_gradient = - |f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1)); - for c in 0..3 { - let dn = wall_gradient(t1[c] - ts[c], t2[c] - ts[c]); - // Traction on the body = −(fluid stress on the fluid side): - // the shear the fluid exerts on the wall along +t. - shear[c] += mu * dn * area; - } - } - Some((force, shear)) - } - /// The cut-cell load route restricted to the cells (and faces) of the /// planes `k0..k1`, divided by the slab's thickness: the load per unit /// span on a body's mid-section. diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs index d060461..1633330 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -16,6 +16,7 @@ pub mod grid; pub mod impose; pub mod loads; pub mod poisson; +pub mod reconstruct; pub mod step; pub mod wall; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs new file mode 100644 index 0000000..2a70c17 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs @@ -0,0 +1,108 @@ +//! The reconstructed wall route (S2-1 remedy): two probes along the +//! interpolant normal on every wall polygon. Its pressure part converges +//! with the operator route's; its shear part over-reads on this scheme's +//! near-wall profile (ny 31 / 62: 1.99 / 2.32 against the discrete +//! 1.77 / 1.81), so it is kept for pressure readings, not as a drag route. +use super::body::Body; +use super::field::Field; +use super::wall::Mask; + +impl Mask { + /// The reconstructed wall route (S2-1 remedy): on every wall polygon + /// (the cell's closure `W_c`, its centroid taken as the cell centre's + /// foot on the interpolant surface) the traction from two probes at + /// `h` and `2h` along the interpolant normal — the wall pressure by + /// linear extrapolation, the wall shear from the quadratic fit of the + /// tangential velocity through the probes and the wall velocity (the + /// binary wall's validated sampler, on the cut geometry's own + /// polygons; the probes are a cell away, past the merged slivers). + /// Force on the body: `Σ_c (p_w W_c + μ ∂ₙu_t |W_c| t)` over the + /// planes `k0..k1`; `None` without a cut geometry. + pub fn cut_wall_force_reconstructed( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + planes: Option<(usize, usize)>, + ) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_reconstructed_parts(body, f, mu, t, planes)?; + Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]]) + } + + /// The reconstructed route split into its pressure and shear parts. + pub fn cut_wall_force_reconstructed_parts( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + planes: Option<(usize, usize)>, + ) -> Option<([f64; 3], [f64; 3])> { + let cut = self.cut.as_ref()?; + let g = self.grid; + let (k0, k1) = planes.unwrap_or((0, g.nz)); + let h = g.dx.min(g.dy).min(g.dz); + let (d1, d2) = (h, 2.0 * h); + let mut force = [0.0; 3]; + let mut shear = [0.0; 3]; + for (idx, w) in cut.wall.iter().enumerate() { + let area = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); + if area == 0.0 || !self.cell_fluid[idx] { + continue; + } + let (k, j, i) = g.kji(idx); + if k < k0 || k >= k1 { + continue; + } + let xc = [ + (i as f64 + 0.5) * g.dx, + (j as f64 + 0.5) * g.dy, + (k as f64 + 0.5) * g.dz, + ]; + let (s, n) = self.interpolant_distance_and_normal(cut, xc); + // The wall point and the outward (into the fluid) normal. + let x = [xc[0] - s * n[0], xc[1] - s * n[1], xc[2] - s * n[2]]; + let at = |d: f64| [x[0] + d * n[0], x[1] + d * n[1], x[2] + d * n[2]]; + let (x1, x2) = (at(d1), at(d2)); + let (Some(p1), Some(p2)) = ( + self.pressure_at(&f.p, x1[0], x1[1], x1[2]), + self.pressure_at(&f.p, x2[0], x2[1], x2[2]), + ) else { + // No fit: the operator's own pressure on this polygon. + for c in 0..3 { + force[c] += f.p[idx] * w[c]; + } + continue; + }; + let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1); + for c in 0..3 { + force[c] += p_wall * w[c]; + } + let (Some(u1), Some(u2)) = ( + self.velocity_at(body, f, x1[0], x1[1], x1[2], t), + self.velocity_at(body, f, x2[0], x2[1], x2[2], t), + ) else { + continue; + }; + let us = body.surface_velocity(x[0], x[1], x[2], t); + let us = [us.0, us.1, us.2]; + // Tangential components (the normal removed) and the wall + // gradient of the quadratic through 0, d1, d2. + let tang = |v: [f64; 3]| { + let vn = v[0] * n[0] + v[1] * n[1] + v[2] * n[2]; + [v[0] - vn * n[0], v[1] - vn * n[1], v[2] - vn * n[2]] + }; + let (t1, t2, ts) = (tang(u1), tang(u2), tang(us)); + let wall_gradient = + |f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1)); + for c in 0..3 { + let dn = wall_gradient(t1[c] - ts[c], t2[c] - ts[c]); + // Traction on the body = −(fluid stress on the fluid side): + // the shear the fluid exerts on the wall along +t. + shear[c] += mu * dn * area; + } + } + Some((force, shear)) + } +}