diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu index cad4201..cdbd880 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu @@ -25,6 +25,8 @@ struct GeomSdf { int cyl_cut, flag_cut; /* cut to the span */ int npts; /* polyline points (x, y interleaved) */ int nst, ns; /* R8-c plate: stations and points per station (nst 0: the polyline) */ + int flat_tip; /* R8-h: 1 = the flat tip (corners rounded to tip_corner); 0 = the capsule */ + double tip_corner; }; struct GeomGrid { @@ -46,6 +48,24 @@ __device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s return outside + rs_min(rs_max(q1, q2), 0.0) - r; } +/* R8-h: the tip's axial distance along the last segment a→b from b (plate.rs tip_axial). */ +__device__ __forceinline__ double tip_axial(double x, double y, double ax, double ay, double bx, double by) +{ + double ex = bx - ax, ey = by - ay; + double l = sqrt(ex * ex + ey * ey); + return ((x - bx) * ex + (y - by) * ey) / l; +} + +/* R8-h: the flat tip's rounded cap (plate.rs flat_cap). */ +__device__ __forceinline__ double flat_cap(double dl, double a, double r) +{ + double q1 = dl + r; + double q2 = a + r; + double m1 = rs_max(q1, 0.0), m2 = rs_max(q2, 0.0); + double outside = sqrt(m1 * m1 + m2 * m2); + return outside + rs_min(rs_max(q1, q2), 0.0) - r; +} + /* * R8-c: the deformed plate (`plate.rs`, the host twin expression for * expression): `poly` holds the stations' rows (x, y interleaved, row-major @@ -57,7 +77,7 @@ __device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s */ __device__ double plate_dist(double x, double y, double z, const GeomSdf& s, const double* __restrict__ P, const double* __restrict__ V, - double* vx, double* vy) + double* vx, double* vy, double* axial) { int nst = s.nst, ns = s.ns; const double* Z = P + 2 * (long long) nst * ns; @@ -97,14 +117,14 @@ __device__ double plate_dist(double x, double y, double z, const GeomSdf& s, double l2 = ex * ex + ey * ey; double u = ((x - ax) * ex + (y - ay) * ey) / l2; if (u < 0.0) u = 0.0; - if (u > 1.0) u = 1.0; + if (u > 1.0 && !(s.flat_tip && m + 2 == ns)) u = 1.0; double px = ax + u * ex, py = ay + u * ey; double qx = x - px, qy = y - py; double d = sqrt(qx * qx + qy * qy); if (d < best) { best = d; mb = m; - ub = u; + ub = u > 1.0 ? 1.0 : u; qbx = qx; qby = qy; } @@ -118,6 +138,20 @@ __device__ double plate_dist(double x, double y, double z, const GeomSdf& s, double qn = (qbx * cx + qby * cy) / best; best = best / sqrt(1.0 + qn * qn); } + if (s.flat_tip && axial) { + int m = ns - 2; + double ax, ay, bx, by; + if (interp) { + ax = R0[2 * m] + w * (R1[2 * m] - R0[2 * m]); + ay = R0[2 * m + 1] + w * (R1[2 * m + 1] - R0[2 * m + 1]); + bx = R0[2 * m + 2] + w * (R1[2 * m + 2] - R0[2 * m + 2]); + by = R0[2 * m + 3] + w * (R1[2 * m + 3] - R0[2 * m + 3]); + } else { + ax = R0[2 * m]; ay = R0[2 * m + 1]; + bx = R0[2 * m + 2]; by = R0[2 * m + 3]; + } + *axial = tip_axial(x, y, ax, ay, bx, by); + } if (V) { const double* V0 = V + 2 * (long long) k * ns; const double* V1 = interp ? V0 + 2 * ns : V0; @@ -144,22 +178,29 @@ __device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, co double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc; if (s.cyl_cut) dc = span_cut(dc, z, s); /* the capsule: distance to the polyline (or the plate, R8-c) */ - double best = 1.0 / 0.0; - if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr); - else for (int m = 0; m + 1 < s.npts; ++m) { - double ax = poly[2 * m], ay = poly[2 * m + 1]; - double bx = poly[2 * m + 2], by = poly[2 * m + 3]; - double ex = bx - ax, ey = by - ay; - double l2 = ex * ex + ey * ey; - double u = ((x - ax) * ex + (y - ay) * ey) / l2; - if (u < 0.0) u = 0.0; - if (u > 1.0) u = 1.0; - double px = ax + u * ex, py = ay + u * ey; - double qx = x - px, qy = y - py; - double d = sqrt(qx * qx + qy * qy); - if (d < best) best = d; + double best = 1.0 / 0.0, axial = 0.0; + if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr, &axial); + else { + for (int m = 0; m + 1 < s.npts; ++m) { + double ax = poly[2 * m], ay = poly[2 * m + 1]; + double bx = poly[2 * m + 2], by = poly[2 * m + 3]; + double ex = bx - ax, ey = by - ay; + double l2 = ex * ex + ey * ey; + double u = ((x - ax) * ex + (y - ay) * ey) / l2; + if (u < 0.0) u = 0.0; + if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0; + double px = ax + u * ex, py = ay + u * ey; + double qx = x - px, qy = y - py; + double d = sqrt(qx * qx + qy * qy); + if (d < best) best = d; + } + if (s.flat_tip) { + int m = s.npts - 2; + axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]); + } } double df = best - s.half; + if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner); if (s.flag_cut) df = span_cut(df, z, s); double r = s.fillet; if (r > 0.0 && dc < r && df < r) { @@ -391,28 +432,36 @@ extern "C" __global__ void e3_geom_gather( __device__ double body_velocity(double x, double y, double z, int c, const GeomSdf& s, const double* __restrict__ poly, const double* __restrict__ vel) { - double best = 1.0 / 0.0, vx = 0.0, vy = 0.0; - if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy); - else for (int m = 0; m + 1 < s.npts; ++m) { - double ax = poly[2 * m], ay = poly[2 * m + 1]; - double bx = poly[2 * m + 2], by = poly[2 * m + 3]; - double ex = bx - ax, ey = by - ay; - double l2 = ex * ex + ey * ey; - double u = ((x - ax) * ex + (y - ay) * ey) / l2; - if (u < 0.0) u = 0.0; - if (u > 1.0) u = 1.0; - double px = ax + u * ex, py = ay + u * ey; - double qx = x - px, qy = y - py; - double d = sqrt(qx * qx + qy * qy); - if (d < best) { - best = d; - double avx = vel[2 * m], avy = vel[2 * m + 1]; - double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3]; - vx = avx + u * (bvx - avx); - vy = avy + u * (bvy - avy); + double best = 1.0 / 0.0, vx = 0.0, vy = 0.0, axial = 0.0; + if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy, &axial); + else { + for (int m = 0; m + 1 < s.npts; ++m) { + double ax = poly[2 * m], ay = poly[2 * m + 1]; + double bx = poly[2 * m + 2], by = poly[2 * m + 3]; + double ex = bx - ax, ey = by - ay; + double l2 = ex * ex + ey * ey; + double u = ((x - ax) * ex + (y - ay) * ey) / l2; + if (u < 0.0) u = 0.0; + if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0; + double px = ax + u * ex, py = ay + u * ey; + double qx = x - px, qy = y - py; + double d = sqrt(qx * qx + qy * qy); + if (d < best) { + best = d; + double uv = u > 1.0 ? 1.0 : u; + double avx = vel[2 * m], avy = vel[2 * m + 1]; + double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3]; + vx = avx + uv * (bvx - avx); + vy = avy + uv * (bvy - avy); + } + } + if (s.flat_tip) { + int m = s.npts - 2; + axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]); } } double df = best - s.half; + if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner); if (s.flag_cut) df = span_cut(df, z, s); double ex0 = x - s.cx, ey0 = y - s.cy; double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs index ed83221..de3da8f 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs @@ -31,6 +31,12 @@ pub struct DeviceSdf { pub half: f64, /// The root fillet radius (0 = the plain union). pub fillet: f64, + /// R8-h: the tip. `None` = the capsule (the semicircular tip around the + /// last point, as before); `Some(r_c)` = a FLAT tip: the flag ends in + /// the plane through the last point normal to the last segment, its two + /// corners rounded to radius `r_c` (0 ≤ r_c ≤ `half`; `r_c = half` is + /// the capsule's semicircle pulled back by `half`, geometrically). + pub tip_corner: Option, /// The centreline polyline's points (x, y) at `t`. pub poly: Vec<[f64; 2]>, /// R6-2 step 2: the centreline's velocity (vx, vy) per point at `t`, the 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 a398422..d2af901 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs @@ -708,6 +708,7 @@ mod tests { r_edge: 0.0, half, fillet: 0.0, + tip_corner: None, poly: Vec::new(), vel: Vec::new(), plate: Some(surf), diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/plate.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/plate.rs index c908a6a..06363d8 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/plate.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/plate.rs @@ -131,12 +131,16 @@ fn span_cut(d2: f64, z: f64, s: &DeviceSdf) -> f64 { } /// The closest segment of a polyline given by `pt(m)`, `m < n`: the -/// in-plane distance, the segment, its parameter and the offset `(qx, qy)`. +/// in-plane distance, the segment, its parameter (clamped to the segment) +/// and the offset `(qx, qy)`. R8-h: with `ray_last` the last segment is a +/// RAY beyond its end point (its parameter unclamped above for the +/// distance) — the lateral distance the flat tip's cap needs. #[inline] fn closest [f64; 2]>( x: f64, y: f64, n: usize, + ray_last: bool, pt: P, ) -> (f64, usize, f64, [f64; 2]) { let mut best = f64::INFINITY; @@ -151,7 +155,7 @@ fn closest [f64; 2]>( if u < 0.0 { u = 0.0; } - if u > 1.0 { + if u > 1.0 && !(ray_last && m + 2 == n) { u = 1.0; } let px = ax + u * ex; @@ -162,25 +166,63 @@ fn closest [f64; 2]>( if d < best { best = d; mb = m; - ub = u; + ub = if u > 1.0 { 1.0 } else { u }; qb = [qx, qy]; } } (best, mb, ub, qb) } +/// R8-h: the signed distance along the last segment's unit tangent from +/// its end point `b` (the flat tip's plane; > 0 beyond the tip). +#[inline] +fn tip_axial(x: f64, y: f64, a: [f64; 2], b: [f64; 2]) -> f64 { + let ex = b[0] - a[0]; + let ey = b[1] - a[1]; + let l = (ex * ex + ey * ey).sqrt(); + ((x - b[0]) * ex + (y - b[1]) * ey) / l +} + +/// R8-h: the flat tip — the strip of lateral distance `dl` (distance to +/// the centreline with its last segment a ray, minus the half-thickness) +/// cut by the tip's plane (axial distance `a`) with the two corners rounded +/// to radius `r` (the span cut's rounded intersection). Exact near the tip +/// where the centreline is straight over the corner (the last segment). +#[inline] +fn flat_cap(dl: f64, a: f64, r: f64) -> f64 { + let q1 = dl + r; + let q2 = a + r; + let m1 = rs_max(q1, 0.0); + let m2 = rs_max(q2, 0.0); + let outside = (m1 * m1 + m2 * m2).sqrt(); + outside + rs_min(rs_max(q1, q2), 0.0) - r +} + /// The plate's in-plane closest point at `(x, y, z)`, the slope-corrected /// distance to the mid-surface, and the interpolated velocity there /// (zero without velocities). -fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) { +/// +/// R8-h: with `flat` the last segment is a ray (see [`closest`]) and the +/// third value is the tip's axial distance on the interpolated polyline +/// (0 otherwise). +fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64, flat: bool) -> (f64, [f64; 2], f64) { let ns = p.ns; let (k, w) = p.bracket(z); let row = |k: usize, m: usize| p.xy[k * ns + m]; let lerp2 = |a: [f64; 2], b: [f64; 2], w: f64| [a[0] + w * (b[0] - a[0]), a[1] + w * (b[1] - a[1])]; let (best, m, u, q) = match w { - None => closest(x, y, ns, |m| row(k, m)), - Some(w) => closest(x, y, ns, |m| lerp2(row(k, m), row(k + 1, m), w)), + None => closest(x, y, ns, flat, |m| row(k, m)), + Some(w) => closest(x, y, ns, flat, |m| lerp2(row(k, m), row(k + 1, m), w)), + }; + let axial = if flat { + let at = |m: usize| match w { + None => row(k, m), + Some(w) => lerp2(row(k, m), row(k + 1, m), w), + }; + tip_axial(x, y, at(ns - 2), at(ns - 1)) + } else { + 0.0 }; let mut d = best; if let Some(_w) = w { @@ -211,7 +253,7 @@ fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) { let (a, b) = (at(m), at(m + 1)); [a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])] }; - (d, v) + (d, v, axial) } impl DeviceSdf { @@ -231,20 +273,30 @@ impl DeviceSdf { /// and its surface velocity at the closest point. #[must_use] pub fn flag_distance_host(&self, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) { - let (best, v) = match self.plate.as_ref() { - Some(p) => plate_closest(p, x, y, z), + let flat = self.tip_corner.is_some(); + let (best, v, axial) = match self.plate.as_ref() { + Some(p) => plate_closest(p, x, y, z, flat), None => { - let (best, m, u, _) = closest(x, y, self.poly.len(), |m| self.poly[m]); + let n = self.poly.len(); + let (best, m, u, _) = closest(x, y, n, flat, |m| self.poly[m]); let v = if self.vel.len() == self.poly.len() && !self.vel.is_empty() { let (a, b) = (self.vel[m], self.vel[m + 1]); [a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])] } else { [0.0, 0.0] }; - (best, v) + let axial = if flat { + tip_axial(x, y, self.poly[n - 2], self.poly[n - 1]) + } else { + 0.0 + }; + (best, v, axial) } }; let mut df = best - self.half; + if let Some(rc) = self.tip_corner { + df = flat_cap(df, axial, rc); + } if self.flag_cut { df = span_cut(df, z, self); } @@ -344,6 +396,7 @@ mod tests { r_edge: 0.0066, half: 0.01, fillet: 0.0, + tip_corner: None, poly, vel, plate, @@ -406,7 +459,7 @@ mod tests { vel: Vec::new(), }; for &(y, zz) in &[(0.53, 0.05), (0.47, 0.15), (0.6, 0.12)] { - let (d, _) = plate_closest(&p, 0.1, y, zz); + let (d, _, _) = plate_closest(&p, 0.1, y, zz, false); let exact = (y - (0.5 + a * (zz - 0.1))).abs() / (1.0 + a * a).sqrt(); assert!((d - exact).abs() < 1e-15, "{d} vs {exact}"); } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs index cc47da5..caf0783 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs @@ -75,6 +75,9 @@ struct GeomSdf { npts: i32, nst: i32, ns: i32, + /// R8-h: the flat tip (1) and its corner radius. + flat_tip: i32, + tip_corner: f64, } unsafe impl DeviceRepr for GeomSdf {} unsafe impl ValidAsZeroBits for GeomSdf {} @@ -151,6 +154,8 @@ fn geom_sdf(sdf: &DeviceSdf) -> GeomSdf { npts: sdf.poly.len() as i32, nst: sdf.plate.as_ref().map_or(0, |p| p.z.len() as i32), ns: sdf.plate.as_ref().map_or(0, |p| p.ns as i32), + flat_tip: i32::from(sdf.tip_corner.is_some()), + tip_corner: sdf.tip_corner.unwrap_or(0.0), } } diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index 7c08033..ef55561 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -96,8 +96,35 @@ fn root_fillet() -> f64 { /// record before this date had the apex 10 mm beyond A (`=0` restores /// them): at ny 62 on the recorded motion that was drag 253.3 → 240.0 and /// the lift swing 1,005 → 836 (the overset's 867). +/// +/// R8-h: with the flat tip the default inset is 0 (the flat face through A). fn tip_inset() -> f64 { - env_f("RTX_E3_FLAG_TIP_INSET", FLAG_HALF) + env_f( + "RTX_E3_FLAG_TIP_INSET", + if flat_tip().is_some() { 0.0 } else { FLAG_HALF }, + ) +} + +/// R8-h (2026-09-25): the tip's shape. `RTX_E3_FLAG_TIP=flat` gives the +/// flag a FLAT tip through the centreline's last point (A, the inset +/// defaulting to 0), its corners rounded to `RTX_E3_FLAG_TIP_CORNER` +/// metres (default 0.00125, the 2D overset's recommended line; at most +/// `FLAG_HALF`); unset or `capsule` = the capsule (the semicircular tip). +/// The flat tip's host φ and surface velocity are the device form's +/// (`DeviceSdf::phi_host`, the kernel's arithmetic). +fn flat_tip() -> Option { + match std::env::var("RTX_E3_FLAG_TIP").as_deref() { + Err(_) | Ok("capsule") => None, + Ok("flat") => { + let rc = env_f("RTX_E3_FLAG_TIP_CORNER", 0.00125); + assert!( + (0.0..=FLAG_HALF).contains(&rc), + "RTX_E3_FLAG_TIP_CORNER {rc} outside [0, {FLAG_HALF}]" + ); + Some(rc) + } + Ok(v) => panic!("RTX_E3_FLAG_TIP={v}: flat or capsule"), + } } /// Smooth union with a concave fillet of radius `r` (the plain `min` at r = 0). @@ -491,6 +518,7 @@ fn flag_wake_on_the_device() { r_edge, half: FLAG_HALF, fillet: r_fillet, + tip_corner: flat_tip(), poly: match (&plate, recorded()) { (Some(_), _) => Vec::new(), (None, Some(rec)) => recorded_polyline(rec, t) @@ -527,9 +555,9 @@ fn flag_wake_on_the_device() { c.1.clone().expect("sdf") }) }; - let body = if plate_body() { + let body = if plate_body() || flat_tip().is_some() { // R8-c: the host φ and surface velocity ARE the device form's (the - // kernel's arithmetic on the host). + // kernel's arithmetic on the host); R8-h: the flat tip too. Body::from_sdf(move |x, y, z, t| sdf_at(t).phi_host(x, y, z)) .with_surface_velocity(move |x, y, z, t| sdf_at(t).velocity_host(x, y, z)) } else { @@ -587,6 +615,9 @@ fn flag_wake_on_the_device() { g.cells(), (t_end / dt).ceil() as usize ); + if let Some(rc) = flat_tip() { + println!(" R8-h: FLAT tip through the centreline's last point, corner radius {rc:.5} m"); + } unsafe { std::env::set_var("RTX_PROFILE", "1") }; let mut device = DeviceStep::new(solver, g); device.upload(&field); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flat_tip.rs b/crates/specialized/rtx-cfd/tests/embedded3_flat_tip.rs new file mode 100644 index 0000000..6750cce --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_flat_tip.rs @@ -0,0 +1,138 @@ +//! R8-h (2026-09-25): the flat-tipped flag's geometry (`DeviceSdf::tip_corner`, +//! host twin `DeviceSdf::phi_host` = the kernel's `geom_phi_at`). +//! +//! G2: a straight flag (a five-point polyline, root a semicircle, tip flat +//! with corners rounded to `r_c`) at two orientations, cut on the cubic grid +//! at three resolutions: the solid volume and the wall area of the cut +//! geometry against the analytic rounded rectangle; the tip plane through +//! the last point; `r_c = half` = the capsule pulled back by `half`. + +use rtx_cfd::solvers::incompressible::embedded3::{Body, CutGeometry, DeviceSdf, Grid}; + +const HALF: f64 = 0.01; +const LEN: f64 = 0.35; +const P0: [f64; 2] = [0.05, 0.08]; + +fn flag(theta: f64, tip: Option, pull_back: f64) -> DeviceSdf { + let (c, s) = (theta.cos(), theta.sin()); + let len = LEN - pull_back; + let poly = (0..5) + .map(|m| { + let f = len * m as f64 / 4.0; + [P0[0] + f * c, P0[1] + f * s] + }) + .collect(); + DeviceSdf { + // the circle far outside the domain + cyl: [-10.0, -10.0, 0.01], + cyl_cut: false, + flag_cut: false, + zc: 0.0, + span: 1.0, + r_edge: 0.0, + half: HALF, + fillet: 0.0, + tip_corner: tip, + poly, + vel: Vec::new(), + plate: None, + } +} + +/// The analytic cross-section: the strip `LEN × 2 HALF`, the root's half +/// disk, less the two rounded corners' `(1 − π/4) r²` (capsule: a second +/// half disk instead of the flat end). +fn analytic(tip: Option) -> (f64, f64) { + let pi = std::f64::consts::PI; + match tip { + Some(r) => ( + LEN * 2.0 * HALF + 0.5 * pi * HALF * HALF - 2.0 * (1.0 - 0.25 * pi) * r * r, + 2.0 * (LEN - r) + (2.0 * HALF - 2.0 * r) + pi * r + pi * HALF, + ), + None => ( + LEN * 2.0 * HALF + pi * HALF * HALF, + 2.0 * LEN + 2.0 * pi * HALF, + ), + } +} + +fn cut(sdf: DeviceSdf, h: f64) -> (f64, f64) { + let (nx, ny, nz) = ((0.46 / h).ceil() as usize, (0.26 / h).ceil() as usize, 2); + let g = Grid::cubic(nx, ny, nz, h); + let body = Body::from_sdf(move |x, y, z, _t| sdf.phi_host(x, y, z)); + let geo = CutGeometry::build(&body, g, 0.0); + let depth = nz as f64 * h; + let solid = (nx * ny * nz) as f64 * h * h * h - geo.fluid_volume(); + let (wall, _) = geo.wall_area_and_closure(); + (solid / depth, wall / depth) +} + +#[test] +fn flat_tip_plane_passes_through_the_last_point() { + for theta in [0.0, 0.3] { + for rc in [0.0, 0.00125, 0.005, HALF] { + let s = flag(theta, Some(rc), 0.0); + let tip = *s.poly.last().unwrap(); + let (c, sn) = (f64::cos(theta), f64::sin(theta)); + // the face centre on the surface; 1 mm beyond it at 1 mm, inside by 1 mm + assert!(s.phi_host(tip[0], tip[1], 0.0).abs() < 1e-15); + if rc < HALF { + let beyond = s.phi_host(tip[0] + 1e-3 * c, tip[1] + 1e-3 * sn, 0.0); + assert!((beyond - 1e-3).abs() < 1e-12, "{beyond}"); + } + let inside = s.phi_host(tip[0] - 1e-3 * c, tip[1] - 1e-3 * sn, 0.0); + assert!(inside < 0.0); + // the corner's apex region: the rounded corner's arc + let (px, py) = ( + tip[0] - rc * c - (HALF - rc) * sn, + tip[1] - rc * sn + (HALF - rc) * c, + ); + let d = 2e-3; + let (dx, dy) = ((c - sn) / 2f64.sqrt(), (sn + c) / 2f64.sqrt()); + let got = s.phi_host(px + (rc + d) * dx, py + (rc + d) * dy, 0.0); + assert!((got - d).abs() < 1e-12, "corner {rc}: {got} vs {d}"); + } + } +} + +#[test] +fn flat_tip_at_half_is_the_capsule_pulled_back() { + let a = flag(0.2, Some(HALF), 0.0); + let b = flag(0.2, None, HALF); + let mut worst = 0.0f64; + for i in 0..400 { + for j in 0..200 { + let (x, y) = (0.30 + 0.0003 * i as f64, 0.10 + 0.0010 * j as f64); + let (pa, pb) = (a.phi_host(x, y, 0.0), b.phi_host(x, y, 0.0)); + if pa.abs().min(pb.abs()) < 0.02 { + worst = worst.max((pa - pb).abs()); + } + } + } + println!(" r_c = half vs capsule pulled back by half: max |Δφ| {worst:.2e}"); + assert!(worst < 1e-12, "{worst}"); +} + +#[test] +fn flat_tip_cut_volume_and_area_converge_to_the_rounded_rectangle() { + println!(" theta tip h solid err wall err"); + for theta in [0.0, 0.3] { + for tip in [None, Some(0.0), Some(0.00125), Some(0.005)] { + let (va, pa) = analytic(tip); + let mut errs = Vec::new(); + for ny in [62.0, 124.0, 248.0] { + let h = 0.41 / ny; + let (v, p) = cut(flag(theta, tip, 0.0), h); + let (ev, ep) = ((v - va) / va, (p - pa) / pa); + println!( + " {theta:.1} {:>9} {h:.3e} {ev:+.3e} {ep:+.3e}", + tip.map_or("capsule".to_string(), |r| format!("{r:.5}")) + ); + assert!(ev.abs() < 1e-2 && ep.abs() < 2e-2, "{ev} {ep} at h {h}"); + errs.push((ev, ep)); + } + // the volume converges (the rung ny 248 below ny 62's error) + assert!(errs[2].0.abs() < errs[0].0.abs(), "volume {errs:?}"); + } + } +} diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs index 99cf106..38bbc02 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs @@ -56,6 +56,27 @@ const X0: f64 = 0.25; /// last 10 mm is the capsule's cap). const STATIONS: usize = 35; +/// R8-h: the flat tip (`RTX_E3FSI_TIP=flat`, corner radius +/// `RTX_E3FSI_TIP_CORNER`, default 0.00125 m): the centreline gains the +/// tip node A (x 0.60) as its last station and the body's tip is flat +/// through it (`DeviceSdf::tip_corner`); unset = the capsule. +pub fn flat_tip() -> Option { + match std::env::var("RTX_E3FSI_TIP").as_deref() { + Err(_) | Ok("capsule") => None, + Ok("flat") => { + let rc = env_f("RTX_E3FSI_TIP_CORNER", 0.00125); + assert!((0.0..=HALF).contains(&rc), "RTX_E3FSI_TIP_CORNER {rc}"); + Some(rc) + } + Ok(v) => panic!("RTX_E3FSI_TIP={v}: flat or capsule"), + } +} + +/// The centreline's stations (35, or 36 with the flat tip's A). +fn stations() -> usize { + STATIONS + usize::from(flat_tip().is_some()) +} + /// The flag's structure-side bookkeeping: the centreline nodes, the wetted /// edges with their reference coordinates, the tip node A. struct Flag { @@ -82,7 +103,7 @@ impl Flag { .expect("node") .0 }; - let centre = (0..STATIONS) + let centre = (0..stations()) .map(|k| find(X0 + 0.01 * k as f64, 0.2)) .collect(); let interface = Interface::build(mesh); @@ -217,10 +238,10 @@ fn distribute( /// The fluid's centreline from the structure's centreline displacement `c` /// (2 per station) with velocities `(c − c_prev) / dt`. fn line_of(t: f64, c: &[f64], c_prev: &[f64], dt: f64) -> Line { - let pts = (0..STATIONS) + let pts = (0..stations()) .map(|k| [X0 + 0.01 * k as f64 + c[2 * k], 0.2 + c[2 * k + 1]]) .collect(); - let vel = (0..STATIONS) + let vel = (0..stations()) .map(|k| { [ (c[2 * k] - c_prev[2 * k]) / dt, @@ -252,7 +273,7 @@ fn fsi2_on_embedded3() { let mesh = flag_mesh(35, 2); let flag_geo = Flag::build(&mesh); - let zero_c = vec![0.0; 2 * STATIONS]; + let zero_c = vec![0.0; 2 * stations()]; let rest = line_of(0.0, &zero_c, &zero_c, 1.0); // `RTX_E3FSI_LOAD=`: continue from a saved coupled state (extruded // onto the full duct when the saved nz differs); `RTX_E3FSI_SAVE=` @@ -381,7 +402,7 @@ fn fsi2_on_embedded3() { .collect(); let a_dofs = flag.borrow().node_dofs(flag_geo.a_node); let extract = |s: &DynamicState| -> Vec { - let mut c = vec![0.0; 2 * STATIONS]; + let mut c = vec![0.0; 2 * stations()]; for (k, d) in centre_dofs.iter().enumerate() { c[2 * k] = s.displacement[d[0]]; c[2 * k + 1] = s.displacement[d[1]]; @@ -486,8 +507,8 @@ fn fsi2_on_embedded3() { .sqrt(); println!( " step {step} pass: |c_new − c_cand| {res:.3e}, tip cand ({:+.4e}, {:+.4e}), load flag ({:+.3}, {:+.3}) cyl ({:+.3}, {:+.3}) total ({:+.3}, {:+.3}), residual {:.1e}, fresh {}", - cand[2 * STATIONS - 2], - cand[2 * STATIONS - 1], + cand[2 * stations() - 2], + cand[2 * stations() - 1], on_flag[0], on_flag[1], on_cyl[0], diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs index 4ba4ef8..a8349b4 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs @@ -227,17 +227,9 @@ impl E3Fluid { })); VERSION.fetch_add(1, Ordering::AcqRel); let (l1, l2, l3) = (lines.clone(), lines.clone(), lines.clone()); - let body = Body::from_sdf(move |x, y, _z, t| cylinder(x, y).min(capsule(&l1, x, y, t).0)) - .with_surface_velocity(move |x, y, _z, t| { - let (df, (vx, vy)) = capsule(&l2, x, y, t); - if df <= cylinder(x, y) { - (vx, vy, 0.0) - } else { - (0.0, 0.0, 0.0) - } - }); let width = nz as f64 * h; - let body = body.with_device_sdf(move |t| { + let tip_corner = super::flat_tip(); + let device_sdf = move |t: f64| { let line = l3.read().expect("lines").at(t); DeviceSdf { cyl: [CX, CY, R_CYL], @@ -248,12 +240,48 @@ impl E3Fluid { r_edge: h, half: HALF, fillet: 0.0, + tip_corner, poly: line.pts, vel: line.vel, // R8-c's plate body (merged alongside): the span-uniform harness keeps the polyline. plate: None, } - }); + }; + let body = if tip_corner.is_some() { + // R8-h: the flat tip's host φ and velocity are the device form's + // (the kernel's arithmetic), cached per thread, time and lines. + let ds = device_sdf.clone(); + let sdf_at = move |t: f64| -> Arc { + thread_local! { + static SDF: RefCell<(u64, u64, Option>)> = + const { RefCell::new((u64::MAX, u64::MAX, None)) }; + } + let ver = VERSION.load(Ordering::Acquire); + SDF.with(|cell| { + let mut c = cell.borrow_mut(); + if c.0 != t.to_bits() || c.1 != ver || c.2.is_none() { + c.2 = Some(Arc::new(ds(t))); + c.0 = t.to_bits(); + c.1 = ver; + } + c.2.clone().expect("sdf") + }) + }; + let sdf_v = sdf_at.clone(); + Body::from_sdf(move |x, y, z, t| sdf_at(t).phi_host(x, y, z)) + .with_surface_velocity(move |x, y, z, t| sdf_v(t).velocity_host(x, y, z)) + } else { + Body::from_sdf(move |x, y, _z, t| cylinder(x, y).min(capsule(&l1, x, y, t).0)) + .with_surface_velocity(move |x, y, _z, t| { + let (df, (vx, vy)) = capsule(&l2, x, y, t); + if df <= cylinder(x, y) { + (vx, vy, 0.0) + } else { + (0.0, 0.0, 0.0) + } + }) + }; + let body = body.with_device_sdf(device_sdf); solver.set_moving_body(body); let g = Grid::cubic(nx, ny, nz, h); let mut field = Field::new(g);