From 05bf4c74a06d24edf1e69a051959ccaf5df61445 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 25 Sep 2026 20:45:04 -0500 Subject: [PATCH] R8-c: drop zero summands from cut_wall_loads; extrapolation diagnostics beyond FAR_OUTSIDE with the worst point Co-Authored-By: Claude Opus 5.5 (1M context) --- .../incompressible/embedded3/interface.rs | 40 ++++++++++++------- .../solvers/incompressible/embedded3/mod.rs | 2 +- .../rtx-cfd/tests/embedded3_flag_wake.rs | 9 +++-- 3 files changed, 33 insertions(+), 18 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs index 4b595d6..a398422 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs @@ -79,18 +79,17 @@ impl Mask { if self.grad_weights.is_some() { return None; } + // The loops visit every fluid cell and face; only the non-zero + // summands are loads (the zero ones change no sum). let mut raw: Vec<(LoadKind, [f64; 3], [f64; 3])> = Vec::new(); - let (p, s) = self - .cut_wall_force_parts_sink(body, f, mu, t, None, &mut |k, x, v| raw.push((k, x, v)))?; - let (d, c) = self.cut_wall_exchange_parts_sink( - body, - f, - mu, - self.density, - t, - None, - &mut |k, x, v| raw.push((k, x, v)), - )?; + let mut keep = |k: LoadKind, x: [f64; 3], v: [f64; 3]| { + if v != [0.0; 3] { + raw.push((k, x, v)); + } + }; + let (p, s) = self.cut_wall_force_parts_sink(body, f, mu, t, None, &mut keep)?; + let (d, c) = + self.cut_wall_exchange_parts_sink(body, f, mu, self.density, t, None, &mut keep)?; let xch = [d[0] + c[0], d[1] + c[1], d[2] + c[2]]; let total = [ p[0] + s[0] + xch[0], @@ -209,6 +208,10 @@ pub struct Location { pub outside: f64, } +/// The natural-coordinate distance outside an element from which a load +/// counts as extrapolated in [`PlateTransfer`] (0.02 of the half-element). +pub const FAR_OUTSIDE: f64 = 0.02; + /// The result of one load transfer (see [`HexPlate::transfer`]). #[derive(Debug, Clone)] pub struct PlateTransfer { @@ -223,9 +226,13 @@ pub struct PlateTransfer { /// The largest Newton residual (m) and the largest outside-ness. pub max_residual: f64, pub max_outside: f64, - /// Loads located outside every element (extrapolated), and their Σ|F|. + /// Loads located outside every element by more than [`FAR_OUTSIDE`] + /// in natural coordinates (a foot off the plate by more than a few % + /// of an element — the surfaces' round-off-level mismatch is excluded), + /// their Σ|F|, and the worst one's point. pub extrapolated: usize, pub extrapolated_load: f64, + pub worst_point: [f64; 3], /// The share of Σ|f_a| on nodes off the wetted surface (the interior /// layers and the clamped root face). pub interior_share: f64, @@ -481,6 +488,7 @@ impl HexPlate { let (mut force_in, mut moment_in) = ([0.0f64; 3], [0.0f64; 3]); let (mut max_residual, mut max_outside) = (0.0f64, f64::NEG_INFINITY); let (mut extrapolated, mut extrapolated_load) = (0usize, 0.0f64); + let mut worst_point = [0.0f64; 3]; for (&(p, f), loc) in loads.iter().zip(&located) { for (&n, &w) in loc.nodes.iter().zip(&loc.weights) { for c in 0..3 { @@ -493,8 +501,11 @@ impl HexPlate { moment_in[c] += m[c]; } max_residual = max_residual.max(loc.residual); - max_outside = max_outside.max(loc.outside); - if loc.outside > 1e-9 { + if loc.outside > max_outside { + max_outside = loc.outside; + worst_point = p; + } + if loc.outside > FAR_OUTSIDE { extrapolated += 1; extrapolated_load += norm(f); } @@ -524,6 +535,7 @@ impl HexPlate { max_outside, extrapolated, extrapolated_load, + worst_point, interior_share: if total_abs > 0.0 { interior_abs / total_abs } else { 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 b516907..10f7273 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -29,7 +29,7 @@ pub use cut::CutGeometry; pub use export_vtk::write_vtk; pub use field::Field; pub use grid::Grid; -pub use interface::{HexPlate, LoadKind, PlateTransfer, WallLoad}; +pub use interface::{FAR_OUTSIDE, HexPlate, LoadKind, PlateTransfer, WallLoad}; pub use loads::SurfaceForce; pub use plate::PlateSurface; pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult}; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index a971530..7c08033 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -26,8 +26,8 @@ use embedded3_flag_kinematics::{Recorded, recorded}; use rtx_cfd::solvers::incompressible::ConvectionScheme; use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep; use rtx_cfd::solvers::incompressible::embedded3::{ - Body, Boundaries, DeviceSdf, Field, Fluid, Grid, HexPlate, Parameters, PlateSurface, Side, - Solver, WallScheme, write_vtk, + Body, Boundaries, DeviceSdf, FAR_OUTSIDE, Field, Fluid, Grid, HexPlate, Parameters, + PlateSurface, Side, Solver, WallScheme, write_vtk, }; use std::io::Write as _; @@ -735,7 +735,7 @@ fn flag_wake_on_the_device() { let extrap_share = tr.extrapolated_load / abs_f.max(1e-300); let ms = lap.elapsed().as_secs_f64() * 1e3; println!( - " transfer t {t:.4}: {} loads ({} flag); route − Σ {:.1e} N; flag F {:+.4} {:+.4} {:+.4} N, M {:+.5} {:+.5} {:+.5} N m; |ΔF|/Σ|F| {rel_f:.1e}, |ΔM|/(Σ|F| L) {rel_m:.1e}; Newton ≤ {:.1e} m, outside ≤ {:.2e} ({} extrapolated, {:.1e} of Σ|F|), interior share {:.3}; {ms:.0} ms", + " transfer t {t:.4}: {} loads ({} flag); route − Σ {:.1e} N; flag F {:+.4} {:+.4} {:+.4} N, M {:+.5} {:+.5} {:+.5} N m; |ΔF|/Σ|F| {rel_f:.1e}, |ΔM|/(Σ|F| L) {rel_m:.1e}; Newton ≤ {:.1e} m, outside ≤ {:.2e} at ({:.4}, {:.4}, {:.4}) ({} beyond {FAR_OUTSIDE}, {:.1e} of Σ|F|), interior share {:.3}; {ms:.0} ms", loads.len(), flag.len(), nrm([route[0] - sum[0], route[1] - sum[1], route[2] - sum[2]]), @@ -747,6 +747,9 @@ fn flag_wake_on_the_device() { tr.moment_in[2], tr.max_residual, tr.max_outside, + tr.worst_point[0], + tr.worst_point[1], + tr.worst_point[2], tr.extrapolated, extrap_share, tr.interior_share