overset mesh: the flat-tip outline (P5-3 option B) — cylinder_flag_outline_deformed_tip / cylinder_flag_patch_deformed{,_from}_tip take the tip's corner radius (t = the recorded semicircle, bit for bit; smaller = the benchmark's flat face between corner arcs sampled at the wall spacing, point count fixed across the deformation, the tip's corners kept off the hull source); pins: corner t reproduces the recorded outline, the 2.5 mm outline has the exact area excess with no collapsed spacing, the O-grid builds straight and at ± 80 mm cold and warm at ny 41/62; harness knob RTX_FSI2O_TIP_CORNER (default 0.01) on the cold and warm builds, printed in both headers; pinned-toolchain fmt on touched crates
Performance Benchmarks / Run Benchmarks (push) Failing after 10s
CI / Format Check (push) Failing after 10s
CI / Build (ubuntu-latest) (push) Failing after 10s
CI / Clippy Check (push) Failing after 10s
Documentation / Build User Guide (push) Successful in 9s
CI / Build CPU-Only (Explicit) (push) Failing after 1m13s
Documentation / Build API Documentation (push) Failing after 2m25s
CI / Build (macos-latest) (push) Canceled after 0s
CI / Test (macos-latest) (push) Canceled after 0s
CI / Test (ubuntu-latest) (push) Canceled after 0s
CI / Python Bindings (maturin) (macos-latest) (push) Canceled after 0s
CI / Python Bindings (maturin) (ubuntu-latest) (push) Canceled after 0s
CI / WASM Build + Size Check (push) Canceled after 0s
CI / Distributed Training Tests (push) Canceled after 0s
CI / CI Success (push) Canceled after 0s

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01LzcjQX7tvgn87CQCyg9Cfr
This commit is contained in:
Omar Sobh
2026-09-13 22:22:08 -05:00
co-authored by Claude Fable 5.1
parent 7224521309
commit 415c32ab04
12 changed files with 519 additions and 71 deletions
+153 -10
View File
@@ -667,8 +667,83 @@ pub fn cylinder_flag_patch_deformed_from(
stretch: f64, stretch: f64,
winslow_sweeps: usize, winslow_sweeps: usize,
) -> CfdResult<(PatchMesh, (usize, f64))> { ) -> CfdResult<(PatchMesh, (usize, f64))> {
let (inner, tip) = cylinder_flag_patch_deformed_from_tip(
cylinder_flag_outline_deformed(centre, r, t, edges, x_tip_ref, fillet, 3, 16, h, 1.15); prev,
centre,
r,
t,
edges,
x_tip_ref,
h,
fillet,
offset,
nn,
stretch,
winslow_sweeps,
t,
)
}
/// [`cylinder_flag_patch_deformed`] with the flat tip of
/// [`cylinder_flag_outline_deformed_tip`] (`tip_corner = t` is the recorded
/// patch).
#[allow(clippy::too_many_arguments)]
pub fn cylinder_flag_patch_deformed_tip(
centre: [f64; 2],
r: f64,
t: f64,
edges: &FlagEdges,
x_tip_ref: f64,
h: f64,
fillet: f64,
offset: f64,
nn: usize,
stretch: f64,
winslow_sweeps: usize,
tip_corner: f64,
) -> CfdResult<(PatchMesh, (usize, f64))> {
cylinder_flag_patch_deformed_from_tip(
None,
centre,
r,
t,
edges,
x_tip_ref,
h,
fillet,
offset,
nn,
stretch,
winslow_sweeps,
tip_corner,
)
}
/// [`cylinder_flag_patch_deformed_from`] with the flat tip.
#[allow(clippy::too_many_arguments)]
pub fn cylinder_flag_patch_deformed_from_tip(
prev: Option<&PatchMesh>,
centre: [f64; 2],
r: f64,
t: f64,
edges: &FlagEdges,
x_tip_ref: f64,
h: f64,
fillet: f64,
offset: f64,
nn: usize,
stretch: f64,
winslow_sweeps: usize,
tip_corner: f64,
) -> CfdResult<(PatchMesh, (usize, f64))> {
let (inner, tip) = cylinder_flag_outline_deformed_tip(
centre, r, t, edges, x_tip_ref, fillet, 3, 16, h, 1.15, tip_corner,
);
// With a flat tip the outline's corners would put 90° corners on the
// hull, whose normal offset folds the rays at the tip; the tip disc
// (already a hull source) covers the tip, so the inner points ahead of
// the disc's centre plane are left out of the hull source.
let flat = tip_corner < tip.radius - 1e-12;
// The hull source carries the deformed outline itself: for a straight // The hull source carries the deformed outline itself: for a straight
// flag the cylinder and the tip disc alone give the flag's surfaces as // flag the cylinder and the tip disc alone give the flag's surfaces as
// the hull's tangents, for a BENT flag they give the chord under the // the hull's tangents, for a BENT flag they give the chord under the
@@ -676,7 +751,9 @@ pub fn cylinder_flag_patch_deformed_from(
// where 6 h was asked, the acceptors' donors reaching into the fringe // where 6 h was asked, the acceptors' donors reaching into the fringe
// (P5-2's first death). // (P5-2's first death).
let mut hull_src = hull_source(centre, r, tip.centre, tip.radius, tip.axis); let mut hull_src = hull_source(centre, r, tip.centre, tip.radius, tip.axis);
hull_src.extend(inner.iter().copied()); hull_src.extend(inner.iter().copied().filter(|p| {
!flat || (p[0] - tip.centre[0]) * tip.axis[0] + (p[1] - tip.centre[1]) * tip.axis[1] < 0.0
}));
let start = prev.and_then(|m| { let start = prev.and_then(|m| {
(m.ns() == inner.len() && m.nn() == nn).then(|| { (m.ns() == inner.len() && m.nn() == nn).then(|| {
(0..=nn) (0..=nn)
@@ -767,11 +844,7 @@ fn arclength_at_x(pts: &[[f64; 2]], cum: &[f64], x0: f64, from_end: bool) -> f64
return cum[i] + f * (cum[i + 1] - cum[i]); return cum[i] + f * (cum[i + 1] - cum[i]);
} }
} }
if from_end { if from_end { cum[n - 1] } else { 0.0 }
cum[n - 1]
} else {
0.0
}
} }
/// Points at the arclength fractions `fr` (`0 ..= 1`) along an open /// Points at the arclength fractions `fr` (`0 ..= 1`) along an open
@@ -810,6 +883,33 @@ pub fn cylinder_flag_outline_deformed(
k_tip: usize, k_tip: usize,
d_straight: f64, d_straight: f64,
grade: f64, grade: f64,
) -> (Vec<[f64; 2]>, TipArc) {
cylinder_flag_outline_deformed_tip(
centre, r, t, edges, x_tip_ref, fillet, k_fillet, k_tip, d_straight, grade, t,
)
}
/// [`cylinder_flag_outline_deformed`] with the tip as the benchmark's FLAT
/// face between two corner arcs of radius `tip_corner` (P5-3 option B):
/// `tip_corner = t` is the recorded semicircle, bit for bit; smaller
/// corners keep the along-wall spacing `d0` (the corner arcs and the flat
/// face are sampled at it, so the fluid step does not collapse) and a
/// point count fixed across the deformation (the flat face's count comes
/// from the UNDEFORMED thickness). The returned [`TipArc`] is still the
/// tip disc of radius half the tip edge — the hull source's shape.
#[allow(clippy::too_many_arguments)]
pub fn cylinder_flag_outline_deformed_tip(
centre: [f64; 2],
r: f64,
t: f64,
edges: &FlagEdges,
x_tip_ref: f64,
fillet: f64,
k_fillet: usize,
k_tip: usize,
d_straight: f64,
grade: f64,
tip_corner: f64,
) -> (Vec<[f64; 2]>, TipArc) { ) -> (Vec<[f64; 2]>, TipArc) {
let (cx, cy) = (centre[0], centre[1]); let (cx, cy) = (centre[0], centre[1]);
// Root junction, as the rigid outline. // Root junction, as the rigid outline.
@@ -852,9 +952,31 @@ pub fn cylinder_flag_outline_deformed(
let mut pts: Vec<[f64; 2]> = Vec::new(); let mut pts: Vec<[f64; 2]> = Vec::new();
// 1. Top edge: from the arc's top end toward the root, along the top // 1. Top edge: from the arc's top end toward the root, along the top
// polyline (given tip → root) from arclength `radius` to x = x_f. // polyline (given tip → root) from arclength `radius` to x = x_f.
// The corner radius: `radius` (the semicircle) or the flat tip's.
let rc = tip_corner.min(radius).max(0.0);
let flat = rc < radius - 1e-12;
// Corner arc centres and the tangent points on the faces.
let c_top = [
tc[0] - rc * axis[0] - rc * normal[0],
tc[1] - rc * axis[1] - rc * normal[1],
];
let c_bot = [
b[0] - rc * axis[0] + rc * normal[0],
b[1] - rc * axis[1] + rc * normal[1],
];
let (start_top, start_bot, s_top0, s_bot1) = if flat {
(
[c_top[0] + rc * normal[0], c_top[1] + rc * normal[1]],
[c_bot[0] - rc * normal[0], c_bot[1] - rc * normal[1]],
rc,
rc,
)
} else {
(start_top, start_bot, radius, radius)
};
let cum_t = open_cum(top); let cum_t = open_cum(top);
let s_root_t = arclength_at_x(top, &cum_t, x_f, true); let s_root_t = arclength_at_x(top, &cum_t, x_f, true);
let mut top_run = open_slice(top, &cum_t, radius, s_root_t); let mut top_run = open_slice(top, &cum_t, s_top0, s_root_t);
top_run[0] = start_top; top_run[0] = start_top;
pts.extend(along_fractions(&top_run, &fr)); pts.extend(along_fractions(&top_run, &fr));
// 2. Top fillet, 3. cylinder arc, 4. bottom fillet — the rigid construction. // 2. Top fillet, 3. cylinder arc, 4. bottom fillet — the rigid construction.
@@ -885,13 +1007,34 @@ pub fn cylinder_flag_outline_deformed(
// to `radius` short of the corner, ending at the arc's bottom end. // to `radius` short of the corner, ending at the arc's bottom end.
let cum_b = open_cum(bottom); let cum_b = open_cum(bottom);
let s_root_b = arclength_at_x(bottom, &cum_b, x_f, false); let s_root_b = arclength_at_x(bottom, &cum_b, x_f, false);
let mut bot_run = open_slice(bottom, &cum_b, s_root_b, cum_b[nb - 1] - radius); let mut bot_run = open_slice(bottom, &cum_b, s_root_b, cum_b[nb - 1] - s_bot1);
let last = bot_run.len() - 1; let last = bot_run.len() - 1;
bot_run[last] = start_bot; bot_run[last] = start_bot;
pts.extend(along_fractions(&bot_run, &fr)); pts.extend(along_fractions(&bot_run, &fr));
if flat {
// 6'. Bottom corner arc (from the bottom face to the tip face), the
// flat tip face, the top corner arc (to the top face, whose start
// point opens the outline — excluded here). Counts from the
// UNDEFORMED thickness so the topology is fixed.
let k_c = ((PI * rc / 2.0) / d0).round().max(2.0) as usize;
let k_f = ((2.0 * (t - rc)) / d0).round().max(1.0) as usize;
let a_b0 = (-normal[1]).atan2(-normal[0]);
pts.extend(arc_points(c_bot, rc, a_b0, a_b0 + PI / 2.0, k_c));
let (p0, p1) = (
[c_bot[0] + rc * axis[0], c_bot[1] + rc * axis[1]],
[c_top[0] + rc * axis[0], c_top[1] + rc * axis[1]],
);
for m in 0..k_f {
let f = m as f64 / k_f as f64;
pts.push([p0[0] + f * (p1[0] - p0[0]), p0[1] + f * (p1[1] - p0[1])]);
}
let a_t0 = axis[1].atan2(axis[0]);
pts.extend(arc_points(c_top, rc, a_t0, a_t0 + PI / 2.0, k_c));
} else {
// 6. Tip semicircle from the bottom end through the apex to the top end. // 6. Tip semicircle from the bottom end through the apex to the top end.
let a0 = (start_bot[1] - tip_c[1]).atan2(start_bot[0] - tip_c[0]); let a0 = (start_bot[1] - tip_c[1]).atan2(start_bot[0] - tip_c[0]);
pts.extend(arc_points(tip_c, radius, a0, a0 + PI, k_tip)); pts.extend(arc_points(tip_c, radius, a0, a0 + PI, k_tip));
}
( (
pts, pts,
TipArc { TipArc {
@@ -41,7 +41,7 @@ mod projection;
use super::ale::{AleBoundaries, SideBoundary}; use super::ale::{AleBoundaries, SideBoundary};
use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind}; use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind};
use super::poisson::{ use super::poisson::{
solve_multigrid_pcg, MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg,
}; };
use super::simple::ConvectionScheme; use super::simple::ConvectionScheme;
use super::{FlowField, SolverResult}; use super::{FlowField, SolverResult};
@@ -53,8 +53,8 @@ pub use curvilinear::{
}; };
pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState}; pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
pub use embedded_body::{ pub use embedded_body::{
polygon_interface_velocity, polygon_signed_distance, EmbeddedBody, EmbeddedMask, FaceKind, EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_interface_velocity,
SurfaceForce, SurfaceSample, polygon_signed_distance,
}; };
pub use flow_field::FlowField; pub use flow_field::FlowField;
pub use overset::{ pub use overset::{
@@ -3,8 +3,8 @@
//! S3 harness) — do the Stokes-limit and upwind orders survive the junction //! S3 harness) — do the Stokes-limit and upwind orders survive the junction
//! fillets' skew? Rungs at the benchmark's h = 0.41 / 41, 62, 82. //! fillets' skew? Rungs at the benchmark's h = 0.41 / 41, 62, 82.
use rtx_cfd::mesh::patch_gen::{cylinder_flag_patch, cylinder_flag_patch_deformed, FlagEdges};
use rtx_cfd::mesh::PatchMesh; use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::{FlagEdges, cylinder_flag_patch, cylinder_flag_patch_deformed};
use rtx_cfd::solvers::incompressible::{ use rtx_cfd::solvers::incompressible::{
CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField,
}; };
@@ -7,8 +7,8 @@
//! Reference (FEATFLOW level 6): drag 14.2929, lift 1.11905. The embedded //! Reference (FEATFLOW level 6): drag 14.2929, lift 1.11905. The embedded
//! staircase measured drag 15.71 (surface) / 15.62 (CV) at ny = 41 (+10%). //! staircase measured drag 15.71 (surface) / 15.62 (CV) at ny = 41 (+10%).
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
use rtx_cfd::mesh::PatchSide; use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
use rtx_cfd::solvers::incompressible::{ use rtx_cfd::solvers::incompressible::{
AleBoundaries, CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters, AleBoundaries, CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
EmbeddedPisoSolver, FlowField, MomentumResidual, NormalDiffusion, OversetField, EmbeddedPisoSolver, FlowField, MomentumResidual, NormalDiffusion, OversetField,
@@ -356,13 +356,39 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
let pct = |b: &rtx_cfd::solvers::incompressible::ResidualBucket| 100.0 * b.fx / wall[0]; let pct = |b: &rtx_cfd::solvers::incompressible::ResidualBucket| 100.0 * b.fx / wall[0];
println!( println!(
" momentum residual ny = {ny} [N/m, x / y; % of wall drag; faces evaluated/total]: solved far Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.3e}, {:.3e}) {}/{} | solved near ring Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.4}, {:.4}) max|r| ({:.3e}, {:.3e}) {}/{} | fringefringe ({:+.4}, {:+.4}) {:+.2}% {}/{} | fringehole ({:+.4}, {:+.4}) {:+.2}% {}/{} | holehole skipped {} (ghosts: {} cells, {} faces) | Σ|r| fringefringe ({:.4}, {:.4}) fringehole ({:.4}, {:.4}); ring total ({:+.4}, {:+.4}) {:+.2}% vs routes' ring defect {:+.4} ({:+.2}%)", " momentum residual ny = {ny} [N/m, x / y; % of wall drag; faces evaluated/total]: solved far Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.3e}, {:.3e}) {}/{} | solved near ring Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.4}, {:.4}) max|r| ({:.3e}, {:.3e}) {}/{} | fringefringe ({:+.4}, {:+.4}) {:+.2}% {}/{} | fringehole ({:+.4}, {:+.4}) {:+.2}% {}/{} | holehole skipped {} (ghosts: {} cells, {} faces) | Σ|r| fringefringe ({:.4}, {:.4}) fringehole ({:.4}, {:.4}); ring total ({:+.4}, {:+.4}) {:+.2}% vs routes' ring defect {:+.4} ({:+.2}%)",
mr.solved_far.fx, mr.solved_far.fy, mr.solved_far.abs_x, mr.solved_far.abs_y, mr.solved_far.evaluated, mr.solved_far.total, mr.solved_far.fx,
mr.solved_near.fx, mr.solved_near.fy, mr.solved_near.abs_x, mr.solved_near.abs_y, mr.solved_near.max_abs_x, mr.solved_near.max_abs_y, mr.solved_near.evaluated, mr.solved_near.total, mr.solved_far.fy,
mr.fringe_fringe.fx, mr.fringe_fringe.fy, pct(&mr.fringe_fringe), mr.fringe_fringe.evaluated, mr.fringe_fringe.total, mr.solved_far.abs_x,
mr.fringe_hole.fx, mr.fringe_hole.fy, pct(&mr.fringe_hole), mr.fringe_hole.evaluated, mr.fringe_hole.total, mr.solved_far.abs_y,
mr.hole_hole_skipped, mr.hole_ghosts, mr.ghost_faces, mr.solved_far.evaluated,
mr.fringe_fringe.abs_x, mr.fringe_fringe.abs_y, mr.fringe_hole.abs_x, mr.fringe_hole.abs_y, mr.solved_far.total,
mr.fringe_fringe.fx + mr.fringe_hole.fx, mr.fringe_fringe.fy + mr.fringe_hole.fy, mr.solved_near.fx,
mr.solved_near.fy,
mr.solved_near.abs_x,
mr.solved_near.abs_y,
mr.solved_near.max_abs_x,
mr.solved_near.max_abs_y,
mr.solved_near.evaluated,
mr.solved_near.total,
mr.fringe_fringe.fx,
mr.fringe_fringe.fy,
pct(&mr.fringe_fringe),
mr.fringe_fringe.evaluated,
mr.fringe_fringe.total,
mr.fringe_hole.fx,
mr.fringe_hole.fy,
pct(&mr.fringe_hole),
mr.fringe_hole.evaluated,
mr.fringe_hole.total,
mr.hole_hole_skipped,
mr.hole_ghosts,
mr.ghost_faces,
mr.fringe_fringe.abs_x,
mr.fringe_fringe.abs_y,
mr.fringe_hole.abs_x,
mr.fringe_hole.abs_y,
mr.fringe_fringe.fx + mr.fringe_hole.fx,
mr.fringe_fringe.fy + mr.fringe_hole.fy,
pct(&mr.fringe_fringe) + pct(&mr.fringe_hole), pct(&mr.fringe_fringe) + pct(&mr.fringe_hole),
hole.0 - ring.0, hole.0 - ring.0,
100.0 * (hole.0 - ring.0) / wall[0], 100.0 * (hole.0 - ring.0) / wall[0],
@@ -573,15 +599,36 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
let hole_flux = sf.0 + ring_sum; let hole_flux = sf.0 + ring_sum;
println!( println!(
" patch momentum balance ny = {ny} [N/m x / y; {} solved cells, {} interface faces, {} wall faces]: balance residual ({:+.3e}, {:+.3e}) | flux-form force through the interface ({:.4}, {:.4}) | wall force, scheme fluxes ({:.4}, {:.4}) | wall force, surface formula ({:.4}, {:.4}) | pressure defect δP = p_ls p_face ({:+.4}, {:+.4}) [{:+.2}% of wall drag] | pieces: conv_acc ({:+.4}, {:+.4}) visc_acc ({:+.4}, {:+.4}) p_face_acc ({:+.4}, {:+.4}) visc_wall ({:+.4}, {:+.4}) p_face_wall ({:+.4}, {:+.4}) p_ls ({:+.4}, {:+.4}) | acceptor band: background hole flux {:.4} patch interface {:.4} = {:+.4} ({:+.2}%)", " patch momentum balance ny = {ny} [N/m x / y; {} solved cells, {} interface faces, {} wall faces]: balance residual ({:+.3e}, {:+.3e}) | flux-form force through the interface ({:.4}, {:.4}) | wall force, scheme fluxes ({:.4}, {:.4}) | wall force, surface formula ({:.4}, {:.4}) | pressure defect δP = p_ls p_face ({:+.4}, {:+.4}) [{:+.2}% of wall drag] | pieces: conv_acc ({:+.4}, {:+.4}) visc_acc ({:+.4}, {:+.4}) p_face_acc ({:+.4}, {:+.4}) visc_wall ({:+.4}, {:+.4}) p_face_wall ({:+.4}, {:+.4}) p_ls ({:+.4}, {:+.4}) | acceptor band: background hole flux {:.4} patch interface {:.4} = {:+.4} ({:+.2}%)",
pb.cells, pb.acc_faces, pb.wall_faces, pb.cells,
bal[0], bal[1], pb.acc_faces,
ff[0], ff[1], pb.wall_faces,
fw[0], fw[1], bal[0],
wall[0], wall[1], bal[1],
dp[0], dp[1], 100.0 * dp[0] / wall[0], ff[0],
pb.conv_acc[0], pb.conv_acc[1], pb.visc_acc[0], pb.visc_acc[1], pb.p_face_acc[0], pb.p_face_acc[1], ff[1],
pb.visc_wall[0], pb.visc_wall[1], pb.p_face_wall[0], pb.p_face_wall[1], pb.p_ls[0], pb.p_ls[1], fw[0],
hole_flux, ff[0], hole_flux - ff[0], 100.0 * (hole_flux - ff[0]) / wall[0], fw[1],
wall[0],
wall[1],
dp[0],
dp[1],
100.0 * dp[0] / wall[0],
pb.conv_acc[0],
pb.conv_acc[1],
pb.visc_acc[0],
pb.visc_acc[1],
pb.p_face_acc[0],
pb.p_face_acc[1],
pb.visc_wall[0],
pb.visc_wall[1],
pb.p_face_wall[0],
pb.p_face_wall[1],
pb.p_ls[0],
pb.p_ls[1],
hole_flux,
ff[0],
hole_flux - ff[0],
100.0 * (hole_flux - ff[0]) / wall[0],
); );
// `RTX_OVERSET_CFD1_SAVE=dir`: the settled fields, for offline // `RTX_OVERSET_CFD1_SAVE=dir`: the settled fields, for offline
// diagnostics without the march (background in `FlowField::save`'s // diagnostics without the march (background in `FlowField::save`'s
@@ -758,7 +805,8 @@ async fn momentum_residual_vanishes_on_the_solved_faces() -> CfdResult<()> {
// prescribed velocity) reads zero, so Σ|r| = N_interface · max|r| on // prescribed velocity) reads zero, so Σ|r| = N_interface · max|r| on
// each lattice. // each lattice.
assert!( assert!(
(near.abs_x - mr.interface_u as f64 * near.max_abs_x).abs() <= 1e-6 * near.abs_x.max(1e-300) (near.abs_x - mr.interface_u as f64 * near.max_abs_x).abs()
<= 1e-6 * near.abs_x.max(1e-300)
&& (near.abs_y - mr.interface_v as f64 * near.max_abs_y).abs() && (near.abs_y - mr.interface_v as f64 * near.max_abs_y).abs()
<= 1e-6 * near.abs_y.max(1e-300), <= 1e-6 * near.abs_y.max(1e-300),
"solved near: not a uniform level offset on the interface — Σ|r| ({:.4e}, {:.4e}) vs N·max|r| ({:.4e}, {:.4e}) with N = ({}, {})", "solved near: not a uniform level offset on the interface — Σ|r| ({:.4e}, {:.4e}) vs N·max|r| ({:.4e}, {:.4e}) with N = ({}, {})",
@@ -24,8 +24,8 @@
//! metres instead) and `RTX_OVERSET_PATCH_ROWS` (12; scale it with the //! metres instead) and `RTX_OVERSET_PATCH_ROWS` (12; scale it with the
//! offset to keep the wall spacing). //! offset to keep the wall spacing).
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
use rtx_cfd::mesh::PatchSide; use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
use rtx_cfd::solvers::incompressible::{ use rtx_cfd::solvers::incompressible::{
AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver, AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver,
EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField, EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField,
@@ -134,7 +134,8 @@ async fn pitching_flag_in_still_fluid_extracts_energy_at_every_h() -> CfdResult<
let dt = period / steps_per_period as f64; let dt = period / steps_per_period as f64;
println!( println!(
" ENERGY PIN ny = {ny} (h {h:.4}, nx {nx}, patch ns {} nn {}, inner edge {hs:.4}): tip ± {tip_amp} m at {f} Hz (θ0 {theta0:.4} rad), dt {dt:.3e} ({steps_per_period} per period), {periods} periods", " ENERGY PIN ny = {ny} (h {h:.4}, nx {nx}, patch ns {} nn {}, inner edge {hs:.4}): tip ± {tip_amp} m at {f} Hz (θ0 {theta0:.4} rad), dt {dt:.3e} ({steps_per_period} per period), {periods} periods",
base.ns(), base.nn() base.ns(),
base.nn()
); );
let start = std::time::Instant::now(); let start = std::time::Instant::now();
let mut per_period: Vec<(f64, f64, f64)> = Vec::new(); // (work, +part, part) let mut per_period: Vec<(f64, f64, f64)> = Vec::new(); // (work, +part, part)
@@ -215,7 +216,11 @@ async fn pitching_flag_in_still_fluid_extracts_energy_at_every_h() -> CfdResult<
let last = per_period.last().copied().unwrap_or((f64::NAN, 0.0, 0.0)); let last = per_period.last().copied().unwrap_or((f64::NAN, 0.0, 0.0));
println!( println!(
" ENERGY PIN ny = {ny}: last-period net work {:+.4e} J/m (gross +{:.3e} / {:+.3e}); added inertia {inertia:.4e} kg·m²/m, rotational damping {damping:.4e} N·m·s/m; {} steps in {:.0} s", " ENERGY PIN ny = {ny}: last-period net work {:+.4e} J/m (gross +{:.3e} / {:+.3e}); added inertia {inertia:.4e} kg·m²/m, rotational damping {damping:.4e} N·m·s/m; {} steps in {:.0} s",
last.0, last.1, last.2, periods * steps_per_period, start.elapsed().as_secs_f64() last.0,
last.1,
last.2,
periods * steps_per_period,
start.elapsed().as_secs_f64()
); );
assert!(last.0.is_finite(), "non-finite work"); assert!(last.0.is_finite(), "non-finite work");
Ok(()) Ok(())
@@ -5,14 +5,14 @@
//! the input polygon, and classifies the ny = 41 background (the rigid //! the input polygon, and classifies the ny = 41 background (the rigid
//! patch's own gates: wall row < 0.5 h, along-body > 0.19 h, 80°). //! patch's own gates: wall row < 0.5 h, along-body > 0.19 h, 80°).
use rtx_cfd::mesh::patch_gen::{
cylinder_flag_outline, cylinder_flag_outline_deformed, cylinder_flag_patch,
cylinder_flag_patch_deformed, FlagEdges,
};
use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::solvers::incompressible::overset::overlap::DEFAULT_OVERLAP_ROWS;
use rtx_cfd::solvers::incompressible::OverlapMap;
use rtx_cfd::CfdResult; use rtx_cfd::CfdResult;
use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::{
FlagEdges, cylinder_flag_outline, cylinder_flag_outline_deformed, cylinder_flag_patch,
cylinder_flag_patch_deformed,
};
use rtx_cfd::solvers::incompressible::OverlapMap;
use rtx_cfd::solvers::incompressible::overset::overlap::DEFAULT_OVERLAP_ROWS;
const L: f64 = 2.5; const L: f64 = 2.5;
const H: f64 = 0.41; const H: f64 = 0.41;
@@ -0,0 +1,236 @@
//! P5-3 option B (`docs/overset_metal_campaign.md` §5.12, 2026-09-14): the
//! flag's tip as the benchmark's FLAT face with rounded corners of radius
//! `tip_corner` (the recorded outline is a full semicircle, `tip_corner =
//! t`). Pins: (1) `tip_corner = t` reproduces the recorded outline to
//! rounding, deformed or not; (2) at 2.5 mm the outline is closed, its
//! area is the semicircle's plus the exact flat-tip excess, its spacing
//! never collapses (the dt budget), its point count is fixed across the
//! deformation (the composite's fixed topology); (3) the O-grid builds on
//! it at ny 41 / 62, straight and at ± 80 mm, cold and warm.
use rtx_cfd::CfdResult;
use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::{
FlagEdges, cylinder_flag_outline_deformed, cylinder_flag_outline_deformed_tip,
cylinder_flag_patch_deformed_from_tip, cylinder_flag_patch_deformed_tip,
};
use std::f64::consts::PI;
const CENTRE: [f64; 2] = [0.2, 0.2];
const R: f64 = 0.05;
const T: f64 = 0.01;
const X0: f64 = 0.25;
const X1: f64 = 0.6;
const FILLET: f64 = 0.5 * 0.41 / 41.0;
const CORNER: f64 = 0.0025;
fn edges(a: f64, nx: usize, ny: usize) -> FlagEdges {
let len = X1 - X0;
let centre = |x: f64| {
let xi = (x - X0) / len;
let y = a * xi * xi * (3.0 - xi) / 2.0;
let dy = a * (6.0 * xi - 3.0 * xi * xi) / 2.0 / len;
(y, dy)
};
let edge = |x: f64, side: f64| -> [f64; 2] {
let (y, dy) = centre(x);
let n = (1.0 + dy * dy).sqrt();
[x - side * T * dy / n, CENTRE[1] + y + side * T / n]
};
let m = 2 * nx;
let bottom: Vec<[f64; 2]> = (0..=m)
.map(|i| edge(X0 + len * i as f64 / m as f64, -1.0))
.collect();
let top: Vec<[f64; 2]> = (0..=m)
.rev()
.map(|i| edge(X0 + len * i as f64 / m as f64, 1.0))
.collect();
let (b, t) = (bottom[m], top[0]);
let k = 2 * ny;
let tip: Vec<[f64; 2]> = (0..=k)
.map(|j| {
let f = j as f64 / k as f64;
[b[0] + f * (t[0] - b[0]), b[1] + f * (t[1] - b[1])]
})
.collect();
FlagEdges { bottom, tip, top }
}
fn outline(a: f64, h: f64, corner: f64) -> Vec<[f64; 2]> {
cylinder_flag_outline_deformed_tip(
CENTRE,
R,
T,
&edges(a, 35, 2),
X1,
FILLET,
3,
16,
h,
1.15,
corner,
)
.0
}
fn signed_area(p: &[[f64; 2]]) -> f64 {
let n = p.len();
(0..n)
.map(|i| {
let (a, b) = (p[i], p[(i + 1) % n]);
a[0] * b[1] - b[0] * a[1]
})
.sum::<f64>()
* 0.5
}
#[test]
fn corner_t_reproduces_the_recorded_semicircular_outline() {
for &a in &[0.0, 0.08] {
for &h in &[0.01, 0.41 / 62.0] {
let old = cylinder_flag_outline_deformed(
CENTRE,
R,
T,
&edges(a, 35, 2),
X1,
FILLET,
3,
16,
h,
1.15,
)
.0;
let new = outline(a, h, T);
assert_eq!(old.len(), new.len(), "point count at a = {a}, h = {h}");
let worst = old
.iter()
.zip(&new)
.map(|(p, q)| ((p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2)).sqrt())
.fold(0.0, f64::max);
assert!(
worst < 1e-12,
"a = {a}, h = {h}: worst point move {worst:.2e}"
);
}
}
}
#[test]
fn flat_tip_outline_is_closed_with_the_exact_area_excess_and_no_collapsed_spacing() {
let h = 0.41 / 62.0;
let semi = outline(0.0, h, T);
let flat = outline(0.0, h, CORNER);
// CCW, positive.
let (a_semi, a_flat) = (signed_area(&semi), signed_area(&flat));
assert!(a_semi > 0.0 && a_flat > 0.0, "areas {a_semi} {a_flat}");
// The tip beyond x = x_tip t: semicircle π t²/2 → rectangle 2 t · t minus
// two corner cut-outs (1 π/4) r_c².
let excess = 2.0 * T * T - 2.0 * (1.0 - PI / 4.0) * CORNER * CORNER - PI * T * T / 2.0;
let got = a_flat - a_semi;
assert!(
(got - excess).abs() < 0.03 * excess,
"area excess {got:.4e} vs exact {excess:.4e}"
);
// The flat face sits at x = x_tip; the outline reaches it (within one
// chord of the corner arcs) and never beyond.
let xmax = flat.iter().map(|p| p[0]).fold(f64::MIN, f64::max);
assert!((xmax - X1).abs() < 1e-9, "tip face at x = {xmax}");
let d0 = (PI * T / 16.0).min(FILLET * PI / 2.0 / 3.0);
let n = flat.len();
let (mut smin, mut smax) = (f64::INFINITY, 0.0_f64);
for i in 0..n {
let (p, q) = (flat[i], flat[(i + 1) % n]);
let s = ((p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2)).sqrt();
smin = smin.min(s);
smax = smax.max(s);
}
assert!(
smin > 0.3 * d0,
"spacing collapsed to {smin:.2e} (d0 {d0:.2e})"
);
assert!(smax < 1.5 * h * 1.15_f64.powi(2) + 1e-9 || smax < 3.0 * d0 || smax <= 1.5 * h);
// Fixed topology across the deformation.
assert_eq!(flat.len(), outline(0.08, h, CORNER).len());
assert_eq!(flat.len(), outline(-0.08, h, CORNER).len());
}
fn build(a: f64, h: f64, prev: Option<&PatchMesh>) -> CfdResult<(PatchMesh, (usize, f64))> {
match prev {
None => cylinder_flag_patch_deformed_tip(
CENTRE,
R,
T,
&edges(a, 35, 2),
X1,
h,
FILLET,
6.0 * h,
12,
4.0,
500,
CORNER,
),
Some(m) => cylinder_flag_patch_deformed_from_tip(
Some(m),
CENTRE,
R,
T,
&edges(a, 35, 2),
X1,
h,
FILLET,
6.0 * h,
12,
4.0,
20,
CORNER,
),
}
}
fn min_area(m: &PatchMesh) -> f64 {
(0..m.cell_count())
.map(|c| m.area(c))
.fold(f64::INFINITY, f64::min)
}
/// Cold builds straight and at ± 80 mm (FSI2's amplitude); the warm
/// regeneration over a small increment (the march's per-step use: a jump
/// from straight to 80 mm in one warm start folds for the recorded
/// semicircle too — not the tip's fault, not the march's path).
#[test]
fn flat_tip_patch_builds_straight_and_bent_cold_and_warm() -> CfdResult<()> {
for &h in &[0.01, 0.41 / 62.0] {
let (m0, _) = build(0.0, h, None)?;
assert!(
min_area(&m0) > 0.0,
"h = {h}: min cell area {:.2e}",
min_area(&m0)
);
for &a in &[0.08, -0.08] {
let (mb, _) = build(a, h, None)?;
assert_eq!(mb.ns(), m0.ns(), "topology at a = {a}");
assert!(
min_area(&mb) > 0.0,
"h = {h} a = {a}: min area {:.2e}",
min_area(&mb)
);
}
let (m1, (sweeps, last)) = build(0.005, h, Some(&m0))?;
assert_eq!(m1.ns(), m0.ns());
assert!(
min_area(&m1) > 0.0,
"h = {h} warm: min area {:.2e}",
min_area(&m1)
);
println!(
" flat tip h = {h:.4}: ns {} nn {}, straight min area {:.2e}, ± 80 mm cold OK, warm +5 mm ({sweeps} sweeps, last move {last:.1e}) min area {:.2e}",
m0.ns(),
m0.nn(),
min_area(&m0),
min_area(&m1)
);
}
Ok(())
}
@@ -13,7 +13,7 @@ use std::sync::{Arc, RwLock};
use nalgebra::Vector3; use nalgebra::Vector3;
use rtx_cfd::mesh::PatchSide; use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed; use rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_tip;
use rtx_cfd::solvers::incompressible::{ use rtx_cfd::solvers::incompressible::{
AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver, AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver,
EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField, EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField,
@@ -45,6 +45,17 @@ fn patch_offset_h() -> f64 {
.unwrap_or(6.0) .unwrap_or(6.0)
} }
/// The tip's corner radius (`RTX_FSI2O_TIP_CORNER`, metres; the recorded
/// outline is the full semicircle, corner = t = 0.01): P5-3 option B — the
/// benchmark's flat tip with rounded corners (2.5 mm) against the
/// semicircle, gated on the lift phase and the per-period amplitude.
pub fn tip_corner() -> f64 {
std::env::var("RTX_FSI2O_TIP_CORNER")
.ok()
.and_then(|v| v.parse().ok())
.unwrap_or(FLAG_T)
}
/// The outline's fillet radius (`RTX_FSI2O_FILLET`, metres; 5 mm = t/2 in /// The outline's fillet radius (`RTX_FSI2O_FILLET`, metres; 5 mm = t/2 in
/// every recorded run, so the tip is a full half-round and the root /// every recorded run, so the tip is a full half-round and the root
/// carries 5 mm fillets — the reference's flag is a sharp rectangle). /// carries 5 mm fillets — the reference's flag is a sharp rectangle).
@@ -292,7 +303,7 @@ impl OversetFluid {
} }
}); });
let (cold_mesh, _) = cylinder_flag_patch_deformed( let (cold_mesh, _) = cylinder_flag_patch_deformed_tip(
CYL_CENTRE, CYL_CENTRE,
CYL_R, CYL_R,
FLAG_T, FLAG_T,
@@ -304,6 +315,7 @@ impl OversetFluid {
patch_rows(), patch_rows(),
patch_stretch(), patch_stretch(),
sweeps, sweeps,
tip_corner(),
)?; )?;
// With the warm-started regeneration (`RTX_FSI2O_WARM_SWEEPS`), the // With the warm-started regeneration (`RTX_FSI2O_WARM_SWEEPS`), the
// march starts on the converged fixed point of the warm build's // march starts on the converged fixed point of the warm build's
@@ -314,7 +326,7 @@ impl OversetFluid {
.unwrap_or(0); .unwrap_or(0);
let start_mesh = if warm_sweeps > 0 { let start_mesh = if warm_sweeps > 0 {
let t0 = std::time::Instant::now(); let t0 = std::time::Instant::now();
let (m, (it, last)) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from( let (m, (it, last)) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from_tip(
Some(&cold_mesh), Some(&cold_mesh),
CYL_CENTRE, CYL_CENTRE,
CYL_R, CYL_R,
@@ -327,6 +339,7 @@ impl OversetFluid {
patch_rows(), patch_rows(),
patch_stretch(), patch_stretch(),
20_000, 20_000,
tip_corner(),
)?; )?;
println!( println!(
" warm base: {it} sweep+respace iterations (last move {last:.1e}) in {:.1} s; {warm_sweeps} per regeneration", " warm base: {it} sweep+respace iterations (last move {last:.1e}) in {:.1} s; {warm_sweeps} per regeneration",
@@ -788,7 +801,7 @@ impl OversetFluid {
/// The patch around the interface `d`. /// The patch around the interface `d`.
pub fn patch_for(&self, d: &[f64]) -> CfdResult<rtx_cfd::mesh::PatchMesh> { pub fn patch_for(&self, d: &[f64]) -> CfdResult<rtx_cfd::mesh::PatchMesh> {
let start = std::time::Instant::now(); let start = std::time::Instant::now();
let (mesh, _) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from( let (mesh, _) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from_tip(
if self.warm_sweeps > 0 { if self.warm_sweeps > 0 {
self.warm_base.as_ref() self.warm_base.as_ref()
} else { } else {
@@ -809,6 +822,7 @@ impl OversetFluid {
} else { } else {
self.sweeps self.sweeps
}, },
tip_corner(),
)?; )?;
self.regen_count.set(self.regen_count.get() + 1); self.regen_count.set(self.regen_count.get() + 1);
self.regen_seconds self.regen_seconds
@@ -100,7 +100,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov
// Every setting the acceptance rule reads, printed once: P5-3 lost a // Every setting the acceptance rule reads, printed once: P5-3 lost a
// day to a floor of 2e-4 against the overnight marches' 1e-6. // day to a floor of 2e-4 against the overnight marches' 1e-6.
println!( println!(
" coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m", " coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m",
cfg.coupler, cfg.coupler,
cfg.reuse, cfg.reuse,
cfg.initial_relaxation, cfg.initial_relaxation,
@@ -117,6 +117,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov
super::overset::bg_convection(), super::overset::bg_convection(),
super::overset::patch_stretch(), super::overset::patch_stretch(),
super::overset::fillet(), super::overset::fillet(),
super::overset::tip_corner(),
); );
// Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the // Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the
@@ -425,13 +425,14 @@ fn fsi2_overset_prescribed_motion() {
let dt = fluid.dt_fluid; let dt = fluid.dt_fluid;
let (d_rigid, l_rigid) = fluid.measure_force(); let (d_rigid, l_rigid) = fluid.measure_force();
println!( println!(
" PRESCRIBED ny = {ny}: dt {dt:.3e}, rigid drag {d_rigid:.2} lift {l_rigid:.2}, replay from t = {t0} (ramp {ramp_w} s) to {t_end}; patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, rounds cap {max_rounds}", " PRESCRIBED ny = {ny}: dt {dt:.3e}, rigid drag {d_rigid:.2} lift {l_rigid:.2}, replay from t = {t0} (ramp {ramp_w} s) to {t_end}; patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m, rounds cap {max_rounds}",
std::env::var("RTX_FSI2O_PATCH_OFFSET").unwrap_or_else(|_| "6".into()), std::env::var("RTX_FSI2O_PATCH_OFFSET").unwrap_or_else(|_| "6".into()),
std::env::var("RTX_FSI2O_PATCH_ROWS").unwrap_or_else(|_| "12".into()), std::env::var("RTX_FSI2O_PATCH_ROWS").unwrap_or_else(|_| "12".into()),
fsi2_harness::overset::patch_convection(), fsi2_harness::overset::patch_convection(),
fsi2_harness::overset::bg_convection(), fsi2_harness::overset::bg_convection(),
fsi2_harness::overset::patch_stretch(), fsi2_harness::overset::patch_stretch(),
fsi2_harness::overset::fillet(), fsi2_harness::overset::fillet(),
fsi2_harness::overset::tip_corner(),
); );
let mut csv = std::env::var("RTX_FSI2O_CSV").ok().map(|p| { let mut csv = std::env::var("RTX_FSI2O_CSV").ok().map(|p| {
use std::io::Write as _; use std::io::Write as _;