Omar SobhandClaude Fable 5.1
415c32ab04
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
831 lines
33 KiB
Rust
831 lines
33 KiB
Rust
//! A-P4 (`docs/overset_metal_campaign.md` §2.2 P4, §5.11): Turek–Hron CFD1
|
||
//! (Re = 20, steady) on the OVERSET — the rigid harness's background
|
||
//! (`turek_hron_cfd.rs`: parabolic inflow, outlet, multigrid, upwind) with
|
||
//! the cylinder–flag O-grid as a static patch (no-slip wall, line-implicit
|
||
//! across). Loads by the patch's wall stress (`surface_force`) and by the
|
||
//! background's control-volume momentum balance (the two-route rule).
|
||
//! 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%).
|
||
|
||
use rtx_cfd::mesh::PatchSide;
|
||
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
|
||
use rtx_cfd::solvers::incompressible::{
|
||
AleBoundaries, CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
|
||
EmbeddedPisoSolver, FlowField, MomentumResidual, NormalDiffusion, OversetField,
|
||
OversetParameters, OversetPisoSolver, PatchConvection, PatchField, PoissonSolverKind,
|
||
SideBoundary,
|
||
};
|
||
use rtx_cfd::{CfdConfig, CfdResult};
|
||
|
||
const L: f64 = 2.5;
|
||
const H: f64 = 0.41;
|
||
const RHO: f64 = 1000.0;
|
||
const NU: f64 = 1e-3;
|
||
const U_MEAN: f64 = 0.2;
|
||
const REF_DRAG: f64 = 14.2929;
|
||
const REF_LIFT: f64 = 1.11905;
|
||
|
||
/// Tag of a settled-field set: resolution, scheme, and the overlap depth
|
||
/// when it is not the default (`RTX_OVERSET_ROWS`).
|
||
fn field_tag(ny: usize) -> String {
|
||
let rows = overlap_rows();
|
||
format!(
|
||
"ny{ny}_{}{}",
|
||
if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
|
||
"tvd"
|
||
} else {
|
||
"upwind"
|
||
},
|
||
if rows == OversetParameters::default().overlap_rows {
|
||
String::new()
|
||
} else {
|
||
format!("_rows{rows}")
|
||
}
|
||
)
|
||
}
|
||
|
||
/// `RTX_OVERSET_ROWS`: patch rows kept non-hole below the acceptor row
|
||
/// (the overlap depth; default 4). The band probe of §5.11.
|
||
fn overlap_rows() -> usize {
|
||
std::env::var("RTX_OVERSET_ROWS")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(OversetParameters::default().overlap_rows)
|
||
}
|
||
|
||
fn inflow(y: f64) -> f64 {
|
||
1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
|
||
}
|
||
|
||
struct Cfd1 {
|
||
drag_surface: f64,
|
||
lift_surface: f64,
|
||
drag_cv: f64,
|
||
lift_cv: f64,
|
||
steps: usize,
|
||
seconds: f64,
|
||
rounds_mean: f64,
|
||
dt: f64,
|
||
residual: MomentumResidual,
|
||
}
|
||
|
||
async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
|
||
let h = H / ny as f64;
|
||
let nx = (L / h).round() as usize;
|
||
let mu = RHO * NU;
|
||
let config = CfdConfig::new()
|
||
.with_density(RHO)
|
||
.with_viscosity(mu)
|
||
.with_reference_velocity(U_MEAN)
|
||
.with_reference_length(0.1);
|
||
let mut background = EmbeddedPisoSolver::new(
|
||
config.clone(),
|
||
EmbeddedParameters {
|
||
corrector_steps: 2,
|
||
tolerance: 1e-7,
|
||
boundaries: AleBoundaries {
|
||
left: SideBoundary::Velocity,
|
||
right: SideBoundary::PressureOutlet,
|
||
bottom: SideBoundary::Velocity,
|
||
top: SideBoundary::Velocity,
|
||
},
|
||
poisson_solver: PoissonSolverKind::Multigrid,
|
||
..EmbeddedParameters::default()
|
||
},
|
||
)?;
|
||
background.set_boundary_velocity(|x, y, _| {
|
||
if x <= 0.0 {
|
||
(inflow(y), 0.0)
|
||
} else {
|
||
(0.0, 0.0)
|
||
}
|
||
});
|
||
|
||
let (mesh, _) = cylinder_flag_patch(
|
||
[0.2, 0.2],
|
||
0.05,
|
||
0.01,
|
||
0.6,
|
||
h,
|
||
0.5 * 0.41 / 41.0,
|
||
6.0 * h,
|
||
12,
|
||
4.0,
|
||
500,
|
||
)?;
|
||
// The patch's explicit along-body diffusion limit (its wall row is
|
||
// line-implicit); the harness's combined criterion for the background.
|
||
let mut hs = f64::INFINITY;
|
||
for c in 0..mesh.cell_count() {
|
||
for (f, _) in mesh.cell_faces(c) {
|
||
if mesh.is_sface(f) {
|
||
let d = mesh.faces()[f].d;
|
||
hs = hs.min((d[0] * d[0] + d[1] * d[1]).sqrt());
|
||
}
|
||
}
|
||
}
|
||
let u_peak = 1.5 * 1.5 * U_MEAN;
|
||
let dt_bg = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
|
||
let dt_patch = 0.4 * (hs * hs / (4.0 * NU)).min(hs / u_peak);
|
||
let dt = dt_bg.min(dt_patch);
|
||
|
||
// P4 step 2: `RTX_OVERSET_CFD1_TVD=1` puts the van Albada deferred
|
||
// correction on the patch (the background stays upwind, as recorded).
|
||
let convection = if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
|
||
PatchConvection::TvdVanAlbada
|
||
} else {
|
||
PatchConvection::Upwind
|
||
};
|
||
let mut patch = CurvilinearPisoSolver::new(
|
||
config,
|
||
CurvilinearParameters {
|
||
tolerance: 1e-5,
|
||
convection,
|
||
normal_diffusion: NormalDiffusion::LineImplicit,
|
||
..CurvilinearParameters::default()
|
||
},
|
||
mesh,
|
||
)?;
|
||
patch.set_side_velocity(PatchSide::Inner, |_, _, _| (0.0, 0.0));
|
||
let mut patch_field = PatchField::new(patch.mesh());
|
||
patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));
|
||
|
||
let mut bg_field = FlowField::new(nx, ny, h, h)?;
|
||
for j in 0..ny {
|
||
let u0 = inflow((j as f64 + 0.5) * h);
|
||
for i in 0..=nx {
|
||
bg_field.u[(j, i)] = u0;
|
||
}
|
||
}
|
||
// A steady march: stop a corrector's rounds when two rounds make no
|
||
// progress (the noise floor); measured at ny = 41 without it: 5.0 rounds
|
||
// per corrector on average (the second corrector 8 every step), 2099 s.
|
||
let params = OversetParameters {
|
||
stall_rounds: std::env::var("RTX_OVERSET_STALL")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(2),
|
||
// Cost question (P4): does the second corrector's ~9 rounds buy a
|
||
// measurable load? `RTX_OVERSET_MAX_ROUNDS=3` caps every corrector.
|
||
max_rounds: std::env::var("RTX_OVERSET_MAX_ROUNDS")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(OversetParameters::default().max_rounds),
|
||
overlap_rows: overlap_rows(),
|
||
..OversetParameters::default()
|
||
};
|
||
let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
|
||
let mut field = OversetField {
|
||
background: bg_field,
|
||
patch: patch_field,
|
||
};
|
||
solver.initialize(&mut field)?;
|
||
|
||
let cv = (
|
||
(0.10 / h).round() as usize,
|
||
(0.75 / h).round() as usize,
|
||
(0.05 / h).round() as usize,
|
||
(0.36 / h).round() as usize,
|
||
);
|
||
let cv_force = |field: &OversetField, solver: &OversetPisoSolver| {
|
||
solver
|
||
.background()
|
||
.mask()
|
||
.expect("mask")
|
||
.control_volume_force(
|
||
&field.background.u,
|
||
&field.background.v,
|
||
&field.background.p,
|
||
&field.background.u_old,
|
||
&field.background.v_old,
|
||
dt,
|
||
RHO,
|
||
mu,
|
||
None,
|
||
cv,
|
||
)
|
||
};
|
||
|
||
// `RTX_OVERSET_CFD1_LOAD=dir`: settled fields saved by a previous run
|
||
// (`RTX_OVERSET_CFD1_SAVE`) replace the march — the diagnostics below
|
||
// run offline in seconds instead of the 20–50 min settle.
|
||
let loaded = match std::env::var("RTX_OVERSET_CFD1_LOAD") {
|
||
Ok(dir) => {
|
||
let tag = field_tag(ny);
|
||
let dir = std::path::Path::new(&dir);
|
||
field.background = FlowField::load(&dir.join(format!("bg_{tag}.bin")))?;
|
||
let read = |name: &str| -> Vec<f64> {
|
||
let bytes = std::fs::read(dir.join(format!("patch_{tag}_{name}.bin")))
|
||
.unwrap_or_else(|e| panic!("load patch {name}: {e}"));
|
||
bytes
|
||
.chunks_exact(8)
|
||
.map(|c| f64::from_le_bytes(c.try_into().expect("8 bytes")))
|
||
.collect()
|
||
};
|
||
field.patch.u = read("u");
|
||
field.patch.v = read("v");
|
||
field.patch.p = read("p");
|
||
field.patch.flux = read("flux");
|
||
assert_eq!(field.patch.u.len(), solver.patch().mesh().cell_count());
|
||
assert_eq!(field.patch.flux.len(), solver.patch().mesh().faces().len());
|
||
println!(" loaded settled fields {tag} from {}", dir.display());
|
||
true
|
||
}
|
||
Err(_) => false,
|
||
};
|
||
let start = std::time::Instant::now();
|
||
let flow_through = L / U_MEAN;
|
||
let min_steps = (flow_through / dt).ceil() as usize;
|
||
let mut history: Vec<f64> = Vec::new();
|
||
let mut steps = 0;
|
||
let mut rounds_total = 0usize;
|
||
let mut correctors_total = 0usize;
|
||
let trace_first = std::env::var("RTX_OVERSET_CFD1_TRACE").is_ok();
|
||
if !loaded {
|
||
loop {
|
||
let r = solver.advance(&mut field, dt).await?;
|
||
steps += 1;
|
||
let every: usize = std::env::var("RTX_OVERSET_CFD1_TRACE_EVERY")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(0);
|
||
if (trace_first && steps <= 6) || (every > 0 && steps % every == 0) {
|
||
let pmax = field
|
||
.background
|
||
.p
|
||
.iter()
|
||
.fold(0.0_f64, |m, v| m.max(v.abs()));
|
||
let ppmax = field.patch.p.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
|
||
let upmax = field.patch.u.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
|
||
println!(
|
||
" step {steps}: rounds {:?} converged {} stalled {} bg res {:.2e} patch div {:.2e} patch iters {} conv {} | max|p| bg {pmax:.3e} patch {ppmax:.3e} max|u| patch {upmax:.3e} | defect bg {:.2e} patch {:.2e}",
|
||
r.rounds,
|
||
r.schwarz_converged,
|
||
r.schwarz_stalled,
|
||
r.background_residual,
|
||
r.patch_max_divergence,
|
||
r.patch_poisson_iterations,
|
||
r.patch_converged,
|
||
r.background_mass_defect,
|
||
r.patch_mass_defect
|
||
);
|
||
}
|
||
if steps >= max_steps {
|
||
break;
|
||
}
|
||
rounds_total += r.rounds.iter().sum::<usize>();
|
||
correctors_total += r.rounds.len();
|
||
if steps % 50 == 0 {
|
||
let (fx, _) = cv_force(&field, &solver);
|
||
history.push(fx);
|
||
let load =
|
||
solver
|
||
.patch()
|
||
.surface_force(&field.patch, PatchSide::Inner, solver.time());
|
||
let umax = field
|
||
.background
|
||
.u
|
||
.iter()
|
||
.fold(0.0_f64, |m, v| m.max(v.abs()));
|
||
if steps % 500 == 0 || !umax.is_finite() {
|
||
println!(
|
||
" ny = {ny}: step {steps} t = {:.2} s drag_cv {fx:.4} drag_wall {:.4} lift_wall {:.4} max|u| {umax:.3} rounds {:?} bg res {:.1e} patch div {:.1e} [{:.0} s]",
|
||
solver.time(),
|
||
load.total()[0],
|
||
load.total()[1],
|
||
r.rounds,
|
||
r.background_residual,
|
||
r.patch_max_divergence,
|
||
start.elapsed().as_secs_f64()
|
||
);
|
||
}
|
||
assert!(
|
||
umax.is_finite(),
|
||
"velocity became non-finite at step {steps}"
|
||
);
|
||
if steps >= min_steps && history.len() > 4 {
|
||
let now = history[history.len() - 1];
|
||
let then = history[history.len() - 5];
|
||
if ((now - then) / now).abs() < 1e-4 {
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
assert!(steps < 2_000_000, "CFD1 at ny = {ny} did not settle");
|
||
}
|
||
}
|
||
let seconds = start.elapsed().as_secs_f64();
|
||
let load = solver
|
||
.patch()
|
||
.surface_force(&field.patch, PatchSide::Inner, solver.time());
|
||
let (drag_cv, lift_cv) = cv_force(&field, &solver);
|
||
// P4 momentum-defect measurement: the force the background transmits
|
||
// into the ring (fringe + hole), into the hole alone, and the patch's
|
||
// wall force — consecutive differences are the active region's
|
||
// residual, the fringe ring's momentum defect, and the patch region's.
|
||
let ring = solver
|
||
.overlap()
|
||
.region_force(&field.background, RHO, mu, |c| c != CellClass::Active);
|
||
let hole = solver
|
||
.overlap()
|
||
.region_force(&field.background, RHO, mu, |c| c == CellClass::Hole);
|
||
let wall = load.total();
|
||
println!(
|
||
" overlap rows {} (patch nn {})",
|
||
overlap_rows(),
|
||
solver.patch().mesh().nn()
|
||
);
|
||
println!(
|
||
" momentum routes ny = {ny}: CV box ({drag_cv:.4}, {lift_cv:.4}) | ring outer ({:.4}, {:.4}) | hole boundary ({:.4}, {:.4}) | wall ({:.4}, {:.4}); defects [% of wall drag]: active {:+.2} fringe ring {:+.2} patch region {:+.2}",
|
||
ring.0,
|
||
ring.1,
|
||
hole.0,
|
||
hole.1,
|
||
wall[0],
|
||
wall[1],
|
||
100.0 * (drag_cv - ring.0) / wall[0],
|
||
100.0 * (ring.0 - hole.0) / wall[0],
|
||
100.0 * (hole.0 - wall[0]) / wall[0],
|
||
);
|
||
// P4 option B: the momentum residual of the solver's OWN staggered
|
||
// upwind stencil on every background face at the settled state. Solved
|
||
// faces read zero by construction (the pin below); the prescribed
|
||
// faces' sum is the momentum the stamping injects, in the solver's
|
||
// metric and without the staircase curves' face-formula error.
|
||
let mr = solver.momentum_residual(&field, dt);
|
||
let pct = |b: &rtx_cfd::solvers::incompressible::ResidualBucket| 100.0 * b.fx / wall[0];
|
||
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}) {}/{} | fringe–fringe ({:+.4}, {:+.4}) {:+.2}% {}/{} | fringe–hole ({:+.4}, {:+.4}) {:+.2}% {}/{} | hole–hole skipped {} (ghosts: {} cells, {} faces) | Σ|r| fringe–fringe ({:.4}, {:.4}) fringe–hole ({:.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_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),
|
||
hole.0 - ring.0,
|
||
100.0 * (hole.0 - ring.0) / wall[0],
|
||
);
|
||
// Where along the ring: the prescribed u faces' x-momentum residual in
|
||
// x-bands (cylinder front, cylinder–flag junction, flag, trailing edge).
|
||
let mut bands = [
|
||
(0.0_f64, 0.20, 0.0_f64, 0usize),
|
||
(0.20, 0.30, 0.0, 0),
|
||
(0.30, 0.55, 0.0, 0),
|
||
(0.55, 1.0, 0.0, 0),
|
||
];
|
||
for f in mr.prescribed.iter().filter(|f| f.is_u && f.r.is_finite()) {
|
||
let x = f.i as f64 * h;
|
||
if let Some(b) = bands.iter_mut().find(|b| x >= b.0 && x < b.1) {
|
||
b.2 += f.r;
|
||
b.3 += 1;
|
||
}
|
||
}
|
||
println!(
|
||
" ring x-momentum residual by x-band ny = {ny} (N/m, u faces): {}; level offset δ = {:.3e} Pa on {} + {} interface faces",
|
||
bands
|
||
.iter()
|
||
.map(|b| format!("x {:.2}–{:.2}: {:+.4} ({} faces)", b.0, b.1, b.2, b.3))
|
||
.collect::<Vec<_>>()
|
||
.join(" | "),
|
||
mr.level_offset(h),
|
||
mr.interface_u,
|
||
mr.interface_v
|
||
);
|
||
// `RTX_OVERSET_CFD1_FRONT=1`: every prescribed u face with x < 0.20
|
||
// (the cylinder front, where the ring's source sits) with the pieces
|
||
// of its residual, and the band's sums by piece.
|
||
if std::env::var("RTX_OVERSET_CFD1_FRONT").is_ok() {
|
||
let mut sums = [0.0_f64; 4];
|
||
let mut consistency = 0.0_f64;
|
||
for f in mr
|
||
.prescribed
|
||
.iter()
|
||
.filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
|
||
{
|
||
let pc = f.pieces;
|
||
let recon = pc[0] + pc[1] - pc[2] + pc[3];
|
||
consistency = consistency.max((recon - f.r).abs());
|
||
for k in 0..4 {
|
||
sums[k] += pc[k];
|
||
}
|
||
let cls = |jj: usize, ii: usize| {
|
||
format!("{:?}", solver.overlap().class(jj, ii))
|
||
.chars()
|
||
.next()
|
||
.unwrap()
|
||
};
|
||
println!(
|
||
" front u face j {:2} i {:2} (x {:.3} y {:.3}) {}|{}: r {:+.4e} = time {:+.2e} + conv {:+.2e} − diff {:+.2e} + pres {:+.2e} u_old {:+.4} u {:+.4}",
|
||
f.j,
|
||
f.i,
|
||
f.i as f64 * h,
|
||
(f.j as f64 + 0.5) * h,
|
||
cls(f.j, f.i - 1),
|
||
cls(f.j, f.i),
|
||
f.r,
|
||
pc[0],
|
||
pc[1],
|
||
pc[2],
|
||
pc[3],
|
||
field.background.u_old[(f.j, f.i)],
|
||
field.background.u[(f.j, f.i)]
|
||
);
|
||
}
|
||
println!(
|
||
" front band sums: time {:+.4} conv {:+.4} −diff {:+.4} pres {:+.4} (pieces reconstruct r to {:.1e})",
|
||
sums[0], sums[1], -sums[2], sums[3], consistency
|
||
);
|
||
}
|
||
// `RTX_OVERSET_CFD1_RAW=1`: the same ring residual on the RAW
|
||
// interpolated stamping — the prescribed faces re-stamped from the
|
||
// patch without `balance_fringe_fluxes` (u and u_old alike) — against
|
||
// the balanced one, and how far the balance moved the faces. If the
|
||
// ring's source is the balance's, the raw residual is small.
|
||
if std::env::var("RTX_OVERSET_CFD1_RAW").is_ok() {
|
||
let mut raw = OversetField {
|
||
background: field.background.clone(),
|
||
patch: field.patch.clone(),
|
||
};
|
||
solver
|
||
.overlap()
|
||
.stamp_fringe_faces(&mut raw.background, &field.patch.u, &field.patch.v);
|
||
let (mut moved, mut moved_max, mut scale) = (0.0_f64, 0.0_f64, 0.0_f64);
|
||
for e in &solver.overlap().fringe_u {
|
||
let d = raw.background.u[(e.j, e.i)] - field.background.u[(e.j, e.i)];
|
||
moved += d.abs();
|
||
moved_max = moved_max.max(d.abs());
|
||
scale = scale.max(field.background.u[(e.j, e.i)].abs());
|
||
}
|
||
for e in &solver.overlap().fringe_v {
|
||
let d = raw.background.v[(e.j, e.i)] - field.background.v[(e.j, e.i)];
|
||
moved += d.abs();
|
||
moved_max = moved_max.max(d.abs());
|
||
}
|
||
raw.background.u_old.copy_from(&raw.background.u);
|
||
raw.background.v_old.copy_from(&raw.background.v);
|
||
// Active faces keep u_old = u^n; only the prescribed faces changed.
|
||
for j in 0..field.background.u.nrows() {
|
||
for i in 0..field.background.u.ncols() {
|
||
if solver.background().mask().expect("mask").u_kind(j, i)
|
||
== rtx_cfd::solvers::incompressible::FaceKind::Fluid
|
||
{
|
||
raw.background.u_old[(j, i)] = field.background.u_old[(j, i)];
|
||
raw.background.u[(j, i)] = field.background.u[(j, i)];
|
||
}
|
||
}
|
||
}
|
||
for j in 0..field.background.v.nrows() {
|
||
for i in 0..field.background.v.ncols() {
|
||
if solver.background().mask().expect("mask").v_kind(j, i)
|
||
== rtx_cfd::solvers::incompressible::FaceKind::Fluid
|
||
{
|
||
raw.background.v_old[(j, i)] = field.background.v_old[(j, i)];
|
||
raw.background.v[(j, i)] = field.background.v[(j, i)];
|
||
}
|
||
}
|
||
}
|
||
let mr_raw = solver.momentum_residual(&raw, dt);
|
||
let mut front_raw = 0.0;
|
||
for f in mr_raw
|
||
.prescribed
|
||
.iter()
|
||
.filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
|
||
{
|
||
front_raw += f.r;
|
||
}
|
||
println!(
|
||
" raw-stamping ring residual ny = {ny}: fringe–fringe ({:+.4}, {:+.4}) fringe–hole ({:+.4}, {:+.4}) ring total x {:+.4} (balanced {:+.4}); front band x {:+.4} (balanced {:+.4}); the balance moved the prescribed faces by Σ|Δu| {:.3e} max {:.3e} (max |u| {:.3}); background mass defect raw {:.3e} balanced {:.3e}",
|
||
mr_raw.fringe_fringe.fx,
|
||
mr_raw.fringe_fringe.fy,
|
||
mr_raw.fringe_hole.fx,
|
||
mr_raw.fringe_hole.fy,
|
||
mr_raw.fringe_fringe.fx + mr_raw.fringe_hole.fx,
|
||
mr.fringe_fringe.fx + mr.fringe_hole.fx,
|
||
front_raw,
|
||
mr.prescribed
|
||
.iter()
|
||
.filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
|
||
.map(|f| f.r)
|
||
.sum::<f64>(),
|
||
moved,
|
||
moved_max,
|
||
scale,
|
||
solver.overlap().background_mass_defect(&raw.background),
|
||
solver.overlap().background_mass_defect(&field.background),
|
||
);
|
||
}
|
||
// The box force in the solver's own flux form on the five boxes (must
|
||
// agree to rounding — the gate), and the bookkeeping it allows: the
|
||
// numerical x-momentum source inside the box on the fluid is
|
||
// wall − box; the ring's share is Σr; the rest is the patch region's
|
||
// (hole boundary → wall) plus the two meshes' metric difference.
|
||
let boxes = [
|
||
(0.10, 0.75, 0.05, 0.36),
|
||
(0.09, 0.70, 0.07, 0.34),
|
||
(0.08, 1.00, 0.03, 0.38),
|
||
(0.10, 1.50, 0.05, 0.36),
|
||
(0.10, 0.75, 0.02, 0.39),
|
||
];
|
||
let solver_forces: Vec<(f64, f64)> = boxes
|
||
.iter()
|
||
.map(|&(x0, x1, y0, y1)| {
|
||
solver.solver_metric_force(
|
||
&field,
|
||
dt,
|
||
(
|
||
(x0 / h).round() as usize,
|
||
(x1 / h).round() as usize,
|
||
(y0 / h).round() as usize,
|
||
(y1 / h).round() as usize,
|
||
),
|
||
)
|
||
})
|
||
.collect();
|
||
let sf = solver_forces[0];
|
||
let spread = solver_forces.iter().fold(0.0_f64, |m, f| {
|
||
m.max((f.0 - sf.0).abs()).max((f.1 - sf.1).abs())
|
||
});
|
||
let ring_sum = mr.fringe_fringe.fx + mr.fringe_hole.fx;
|
||
println!(
|
||
" solver-metric box force ny = {ny}: ({:.4}, {:.4}) [5 boxes spread {:.2e}; CV formula {drag_cv:.4}]; x-momentum sources on the fluid inside the box [N/m, % of wall drag]: total wall − box {:+.4} ({:+.2}%) = ring Σr {:+.4} ({:+.2}%) + rest (patch region + metric) {:+.4} ({:+.2}%)",
|
||
sf.0,
|
||
sf.1,
|
||
spread,
|
||
wall[0] - sf.0,
|
||
100.0 * (wall[0] - sf.0) / wall[0],
|
||
ring_sum,
|
||
100.0 * ring_sum / wall[0],
|
||
wall[0] - sf.0 - ring_sum,
|
||
100.0 * (wall[0] - sf.0 - ring_sum) / wall[0],
|
||
);
|
||
// The patch's own momentum balance on its solved cells (scheme
|
||
// fluxes; the unsteady term is omitted — settled state): the balance
|
||
// residual is the gate; flux_force − wall_force = the least-squares
|
||
// pressure's non-conservation δP; the acceptor band's mismatch is
|
||
// then (box + ring Σr) − flux_force, all in N/m.
|
||
let pb = solver.patch().momentum_balance(&field.patch, solver.time());
|
||
let ff = pb.flux_force();
|
||
let fw = pb.wall_force();
|
||
let dp = pb.pressure_defect();
|
||
let bal = pb.balance();
|
||
let hole_flux = sf.0 + ring_sum;
|
||
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}%)",
|
||
pb.cells,
|
||
pb.acc_faces,
|
||
pb.wall_faces,
|
||
bal[0],
|
||
bal[1],
|
||
ff[0],
|
||
ff[1],
|
||
fw[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
|
||
// diagnostics without the march (background in `FlowField::save`'s
|
||
// format; patch u, v, p as raw little-endian f64 vectors).
|
||
if let Ok(dir) = std::env::var("RTX_OVERSET_CFD1_SAVE") {
|
||
let tag = field_tag(ny);
|
||
let dir = std::path::Path::new(&dir);
|
||
std::fs::create_dir_all(dir).expect("save dir");
|
||
field.background.save(&dir.join(format!("bg_{tag}.bin")))?;
|
||
for (name, vals) in [
|
||
("u", &field.patch.u),
|
||
("v", &field.patch.v),
|
||
("p", &field.patch.p),
|
||
("flux", &field.patch.flux),
|
||
] {
|
||
let bytes: Vec<u8> = vals.iter().flat_map(|x| x.to_le_bytes()).collect();
|
||
std::fs::write(dir.join(format!("patch_{tag}_{name}.bin")), bytes).expect("save patch");
|
||
}
|
||
println!(" saved settled fields to {} as {tag}", dir.display());
|
||
}
|
||
// The wall load split and the fringe ring's extent (the tight box in
|
||
// the sensitivity list must stay outside it).
|
||
let (mut jmin, mut jmax, mut imin, mut imax) = (usize::MAX, 0, usize::MAX, 0);
|
||
for e in &solver.overlap().fringe_cells {
|
||
jmin = jmin.min(e.j);
|
||
jmax = jmax.max(e.j);
|
||
imin = imin.min(e.i);
|
||
imax = imax.max(e.i);
|
||
}
|
||
println!(
|
||
" wall load split: pressure ({:.4}, {:.4}) viscous ({:.4}, {:.4}); fringe ring cells i {imin}–{imax} (x {:.3}–{:.3}) j {jmin}–{jmax} (y {:.3}–{:.3})",
|
||
load.pressure[0],
|
||
load.pressure[1],
|
||
load.viscous[0],
|
||
load.viscous[1],
|
||
imin as f64 * h,
|
||
(imax + 1) as f64 * h,
|
||
jmin as f64 * h,
|
||
(jmax + 1) as f64 * h
|
||
);
|
||
// Box sensitivity of the control-volume route: the same balance on
|
||
// other rectangles, all outside the fringe ring (x 0.110–0.640, y
|
||
// 0.110–0.290 at ny = 41 — a first list had a box at x0 = 0.12 cutting
|
||
// through it and read −21 %). A route that moves with the box by more
|
||
// than its own formula error cannot arbitrate the gap.
|
||
for (x0, x1, y0, y1) in [
|
||
(0.10, 0.75, 0.05, 0.36),
|
||
(0.09, 0.70, 0.07, 0.34),
|
||
(0.08, 1.00, 0.03, 0.38),
|
||
(0.10, 1.50, 0.05, 0.36),
|
||
(0.10, 0.75, 0.02, 0.39),
|
||
] {
|
||
let boxc = (
|
||
(x0 / h).round() as usize,
|
||
(x1 / h).round() as usize,
|
||
(y0 / h).round() as usize,
|
||
(y1 / h).round() as usize,
|
||
);
|
||
let (bx, by) = solver
|
||
.background()
|
||
.mask()
|
||
.expect("mask")
|
||
.control_volume_force(
|
||
&field.background.u,
|
||
&field.background.v,
|
||
&field.background.p,
|
||
&field.background.u_old,
|
||
&field.background.v_old,
|
||
dt,
|
||
RHO,
|
||
mu,
|
||
None,
|
||
boxc,
|
||
);
|
||
println!(
|
||
" CV box x {x0:.2}–{x1:.2} y {y0:.2}–{y1:.2}: drag {bx:.4} ({:+.2}% vs wall) lift {by:.4}",
|
||
100.0 * (bx - wall[0]) / wall[0]
|
||
);
|
||
}
|
||
Ok(Cfd1 {
|
||
drag_surface: load.total()[0],
|
||
lift_surface: load.total()[1],
|
||
drag_cv,
|
||
lift_cv,
|
||
steps,
|
||
seconds,
|
||
rounds_mean: rounds_total as f64 / correctors_total.max(1) as f64,
|
||
dt,
|
||
residual: mr,
|
||
})
|
||
}
|
||
|
||
#[tokio::test]
|
||
async fn cfd1_on_the_overset_against_the_featflow_reference() -> CfdResult<()> {
|
||
let resolutions: Vec<usize> = std::env::var("RTX_OVERSET_CFD1_NY").ok().map_or_else(
|
||
|| vec![41usize],
|
||
|list| {
|
||
list.split(',')
|
||
.map(|t| t.trim().parse().expect("ny list"))
|
||
.collect()
|
||
},
|
||
);
|
||
for &ny in &resolutions {
|
||
let max_steps: usize = std::env::var("RTX_OVERSET_CFD1_MAX_STEPS")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(2_000_000);
|
||
let r = run_cfd1(ny, max_steps).await?;
|
||
let rel = |a: f64, b: f64| 100.0 * (a - b) / b;
|
||
println!(
|
||
" CFD1 overset ny = {ny} (h = {:.4}, dt = {:.2e}, patch {}): wall drag {:.4} ({:+.2}%) lift {:.4} ({:+.2}%); control volume drag {:.4} ({:+.2}%) lift {:.4}; routes differ {:.2}%; [{} steps, {:.0} s, Schwarz rounds mean {:.2}] reference {REF_DRAG} / {REF_LIFT}; embedded staircase at ny=41: 15.71 / 15.62 (+10%)",
|
||
H / ny as f64,
|
||
r.dt,
|
||
if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
|
||
"tvd"
|
||
} else {
|
||
"upwind"
|
||
},
|
||
r.drag_surface,
|
||
rel(r.drag_surface, REF_DRAG),
|
||
r.lift_surface,
|
||
rel(r.lift_surface, REF_LIFT),
|
||
r.drag_cv,
|
||
rel(r.drag_cv, REF_DRAG),
|
||
r.lift_cv,
|
||
100.0 * ((r.drag_surface - r.drag_cv) / r.drag_cv).abs(),
|
||
r.steps,
|
||
r.seconds,
|
||
r.rounds_mean
|
||
);
|
||
assert!(r.drag_surface.is_finite() && r.drag_cv.is_finite());
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// The residual diagnostic is the solver's own operator: on every SOLVED
|
||
/// background face away from the ring the momentum residual (time term
|
||
/// included) is zero to rounding; on the solved faces NEXT to the ring it
|
||
/// is a constant per face that cancels in the sum — the composite's
|
||
/// pressure LEVEL offset between the active cells (whose `p'` had its mean
|
||
/// removed) and the fringe cells (re-stamped from the patch, which never
|
||
/// saw that shift); and the ring buckets are populated. This is what makes
|
||
/// the prescribed faces' sum readable as the stamping's momentum injection
|
||
/// in the solver's metric (§5.11, option B).
|
||
#[tokio::test]
|
||
async fn momentum_residual_vanishes_on_the_solved_faces() -> CfdResult<()> {
|
||
let r = run_cfd1(41, 5).await?;
|
||
let mr = &r.residual;
|
||
let scale = r.drag_surface.abs().max(1.0);
|
||
let far = &mr.solved_far;
|
||
assert_eq!(
|
||
far.evaluated, far.total,
|
||
"every far solved face is evaluable"
|
||
);
|
||
assert!(
|
||
far.abs_x <= 1e-9 * scale && far.abs_y <= 1e-9 * scale,
|
||
"solved far: Σ|r| = ({:.3e}, {:.3e}) is not rounding against {scale:.3}",
|
||
far.abs_x,
|
||
far.abs_y
|
||
);
|
||
let near = &mr.solved_near;
|
||
assert_eq!(
|
||
near.evaluated, near.total,
|
||
"every near solved face is evaluable"
|
||
);
|
||
assert!(
|
||
near.fx.abs() <= 1e-9 * scale && near.fy.abs() <= 1e-9 * scale,
|
||
"solved near: Σr = ({:.3e}, {:.3e}) does not cancel against {scale:.3}",
|
||
near.fx,
|
||
near.fy
|
||
);
|
||
// A pure level offset: every active–fringe INTERFACE face carries the
|
||
// same |r| = δ·h and every other near face (one that only reads a
|
||
// prescribed velocity) reads zero, so Σ|r| = N_interface · max|r| on
|
||
// each lattice.
|
||
assert!(
|
||
(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()
|
||
<= 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 = ({}, {})",
|
||
near.abs_x,
|
||
near.abs_y,
|
||
mr.interface_u as f64 * near.max_abs_x,
|
||
mr.interface_v as f64 * near.max_abs_y,
|
||
mr.interface_u,
|
||
mr.interface_v
|
||
);
|
||
assert!(mr.fringe_fringe.evaluated > 0 && mr.fringe_hole.evaluated > 0);
|
||
assert_eq!(
|
||
mr.fringe_fringe.evaluated, mr.fringe_fringe.total,
|
||
"every fringe–fringe face has a fully valid stencil"
|
||
);
|
||
assert_eq!(
|
||
mr.fringe_hole.evaluated, mr.fringe_hole.total,
|
||
"every fringe–hole face has a fully valid stencil with the ghost band"
|
||
);
|
||
Ok(())
|
||
}
|