Omar SobhandClaude Fable 5.1
b80af59ca4
CI / Build (macos-latest) (push) Waiting to run
CI / Test (macos-latest) (push) Blocked by required conditions
CI / Test (ubuntu-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (macos-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (ubuntu-latest) (push) Blocked by required conditions
CI / WASM Build + Size Check (push) Blocked by required conditions
CI / Distributed Training Tests (push) Blocked by required conditions
CI / CI Success (push) Blocked by required conditions
CI / Build CPU-Only (Explicit) (push) Failing after 4s
Documentation / Build API Documentation (push) Failing after 6s
Documentation / Build User Guide (push) Successful in 5s
CI / Format Check (push) Failing after 10s
CI / Build (ubuntu-latest) (push) Failing after 1m30s
CI / Clippy Check (push) Failing after 1m50s
Performance Benchmarks / Run Benchmarks (push) Successful in 2m15s
Co-Authored-By: Claude Fable 5.1 <[email protected]>
790 lines
30 KiB
Rust
790 lines
30 KiB
Rust
//! S2-6 instrument: the cut wall's effective position on an OBLIQUE wall
|
||
//! with the flow IN the plane of the cut. Poiseuille flow (body force `F`
|
||
//! along the tangent) in a channel between two embedded parallel planes
|
||
//! `y = s x + c0` and `y = s x + c0 + w`; the sides carry the exact
|
||
//! solution, z is periodic. The exact field is `U(r) t`, `U = F/(2μ)
|
||
//! (G²/4 − r²)`, constant pressure, and the 5-point Laplacian is exact on
|
||
//! it, so on full faces of a central window `u/t_x + F r²/(2μ)` is the
|
||
//! constant `F G_eff²/(8μ)`: the effective gap from the u faces and from
|
||
//! the v faces separately, each an effective wall offset per wall in
|
||
//! units of h. Unlike the z-directed flat-wall test this one exercises the
|
||
//! own-direction coupling of cut faces, the convective terms' cut-face
|
||
//! values and the projection next to the wall (the spurious pressure is
|
||
//! printed).
|
||
use rtx_cfd::solvers::incompressible::ConvectionScheme;
|
||
use rtx_cfd::solvers::incompressible::embedded3::{
|
||
Body, Boundaries, FaceKind, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
|
||
};
|
||
|
||
const MU: f64 = 0.1;
|
||
const LX: f64 = 2.0;
|
||
const LY: f64 = 2.0;
|
||
const W: f64 = 0.5;
|
||
|
||
/// The cut wall's parameters: the environment's (`None`), or the S2-6
|
||
/// closures forced on / off for the gate.
|
||
fn parameters(s26: Option<bool>) -> Parameters {
|
||
let mut p = Parameters {
|
||
corrector_steps: 2,
|
||
tolerance: 1e-10,
|
||
convection_scheme: ConvectionScheme::Upwind,
|
||
wall_scheme: WallScheme::CutCell,
|
||
boundaries: Boundaries {
|
||
z0: Side::Periodic,
|
||
z1: Side::Periodic,
|
||
..Boundaries::default()
|
||
},
|
||
..Parameters::default()
|
||
};
|
||
if let Some(on) = s26 {
|
||
p.diffusion_centroid = true;
|
||
p.wall_distance_oblique = on;
|
||
p.diffusion_transverse = on;
|
||
p.distance_floor_fine = on;
|
||
}
|
||
p
|
||
}
|
||
|
||
struct Reading {
|
||
/// Effective wall offset per wall from the fitted profile of the u / v faces.
|
||
off_u: f64,
|
||
off_v: f64,
|
||
/// The driving force the profile's curvature implies, over F.
|
||
force_u: f64,
|
||
/// 1 − (the fitted streamwise pressure slope)/F: must equal `force_u`.
|
||
force_p: f64,
|
||
/// RMS of the pressure about its linear fit, over F·G: full cells, cut cells.
|
||
p_full: f64,
|
||
p_cut: f64,
|
||
/// Cut cells (fraction < 1) and virtually merged cells in the mask.
|
||
cut_cells: usize,
|
||
merged: usize,
|
||
}
|
||
|
||
/// Least squares of `y = a − b x`: returns (a, b).
|
||
fn fit(points: &[(f64, f64)]) -> (f64, f64) {
|
||
let n = points.len() as f64;
|
||
let (sx, sy) = points
|
||
.iter()
|
||
.fold((0.0, 0.0), |s, p| (s.0 + p.0, s.1 + p.1));
|
||
let (mx, my) = (sx / n, sy / n);
|
||
let (sxx, sxy) = points.iter().fold((0.0, 0.0), |s, p| {
|
||
(s.0 + (p.0 - mx) * (p.0 - mx), s.1 + (p.0 - mx) * (p.1 - my))
|
||
});
|
||
let slope = sxy / sxx;
|
||
(my - slope * mx, -slope)
|
||
}
|
||
|
||
fn reading(n: usize, slope: f64, c0: f64, along_z: bool, s26: Option<bool>) -> Reading {
|
||
// The driving force (`RTX_E3_OBLIQUE_F`): the problem is linear in it
|
||
// but for the convective terms, so a small value switches them off.
|
||
#[allow(non_snake_case)]
|
||
let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(1.0);
|
||
let h = 1.0 / n as f64;
|
||
let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
|
||
let norm = (1.0 + slope * slope).sqrt();
|
||
let (tx, ty) = (1.0 / norm, slope / norm);
|
||
let gap = W / norm;
|
||
let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm;
|
||
let speed = move |x: f64, y: f64| {
|
||
let r = r_of(x, y);
|
||
if r.abs() < 0.5 * gap {
|
||
F / (2.0 * MU) * (0.25 * gap * gap - r * r)
|
||
} else {
|
||
0.0
|
||
}
|
||
};
|
||
let mut solver = Solver::new(
|
||
Fluid {
|
||
density: 1.0,
|
||
viscosity: MU,
|
||
reference_velocity: 1.0,
|
||
reference_length: 1.0,
|
||
},
|
||
parameters(s26),
|
||
);
|
||
// `along_z`: the same channel with the flow along the periodic z (the
|
||
// cross-direction diffusion of w alone: no pressure, no convection).
|
||
solver.set_boundary_velocity(move |x, y, _z, _t| {
|
||
let s = speed(x, y);
|
||
if along_z {
|
||
(0.0, 0.0, s)
|
||
} else {
|
||
(s * tx, s * ty, 0.0)
|
||
}
|
||
});
|
||
// `RTX_E3_OBLIQUE_DRIVE=pressure`: no body force — the sides' exact
|
||
// profile pins the flow rate and the pressure slope alone drives it.
|
||
let by_pressure = std::env::var("RTX_E3_OBLIQUE_DRIVE").is_ok_and(|v| v == "pressure");
|
||
solver.set_momentum_source(move |_, _, _, _| {
|
||
if by_pressure {
|
||
(0.0, 0.0, 0.0)
|
||
} else if along_z {
|
||
(0.0, 0.0, F)
|
||
} else {
|
||
(F * tx, F * ty, 0.0)
|
||
}
|
||
});
|
||
solver.set_body(Body::from_sdf(move |x, y, _z, _t| {
|
||
0.5 * gap - r_of(x, y).abs()
|
||
}));
|
||
let g = Grid::cubic(nx, ny, nz, h);
|
||
let mut field = Field::new(g);
|
||
for k in 0..=nz {
|
||
for j in 0..ny {
|
||
for i in 0..nx {
|
||
if along_z {
|
||
field.w[g.wface(k, j, i)] = speed((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
for k in 0..nz {
|
||
if along_z {
|
||
break;
|
||
}
|
||
for j in 0..ny {
|
||
for i in 0..=nx {
|
||
field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h);
|
||
}
|
||
}
|
||
for j in 0..=ny {
|
||
for i in 0..nx {
|
||
field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h);
|
||
}
|
||
}
|
||
}
|
||
solver.initialize(&mut field);
|
||
let dt = 0.5 * h * h / (6.0 * MU);
|
||
let steps = (2.0 / dt).ceil() as usize;
|
||
for _ in 0..steps {
|
||
solver.advance(&mut field, dt);
|
||
}
|
||
let mask = solver.mask().expect("mask");
|
||
let cut_cells = (0..g.cells())
|
||
.filter(|&c| mask.cell_active(c) && mask.vol(c) < 1.0)
|
||
.count();
|
||
let merged = mask.merged_cells();
|
||
// The sides pin the flow RATE (the exact profile), so a displaced wall
|
||
// appears as a streamwise pressure slope: F_eff = F − dp/ds, and on the
|
||
// full faces of a central window `u/t_x = F_eff/(2μ) (G_eff²/4 − r²)`
|
||
// exactly (the 5-point Laplacian is exact on it). Fit both constants.
|
||
let window = |x: f64| (x - 0.5 * LX).abs() < 0.3;
|
||
let (mut pu, mut pv) = (Vec::new(), Vec::new());
|
||
let (mut pf, mut pc) = (Vec::new(), Vec::new());
|
||
for j in 0..ny {
|
||
for i in 0..nx {
|
||
let (xu, yu) = (i as f64 * h, (j as f64 + 0.5) * h);
|
||
let fu = g.uface(0, j, i);
|
||
let (xw, yw) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
|
||
let fw = g.wface(0, j, i);
|
||
if along_z {
|
||
if window(xw) && r_of(xw, yw).abs() < 0.4 * gap && mask.a_w(fw) >= 1.0 {
|
||
pu.push((r_of(xw, yw).powi(2), field.w[fw]));
|
||
}
|
||
} else if window(xu) && r_of(xu, yu).abs() < 0.4 * gap && mask.a_u(fu) >= 1.0 {
|
||
pu.push((r_of(xu, yu).powi(2), field.u[fu] / tx));
|
||
}
|
||
let (xv, yv) = ((i as f64 + 0.5) * h, j as f64 * h);
|
||
let fv = g.vface(0, j, i);
|
||
if !along_z
|
||
&& slope > 0.0
|
||
&& window(xv)
|
||
&& r_of(xv, yv).abs() < 0.4 * gap
|
||
&& mask.a_v(fv) >= 1.0
|
||
{
|
||
pv.push((r_of(xv, yv).powi(2), field.v[fv] / ty));
|
||
}
|
||
let c = g.cell(0, j, i);
|
||
let (xc, yc) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
|
||
if window(xc) && mask.cell_active(c) && r_of(xc, yc).abs() < 0.5 * gap + h {
|
||
let s_along = xc * tx + yc * ty;
|
||
let full = r_of(xc, yc).abs() < 0.5 * gap - 1.5 * h;
|
||
if full {
|
||
pf.push((s_along, field.p[c]));
|
||
} else {
|
||
pc.push((s_along, field.p[c]));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let profile = |points: &[(f64, f64)]| {
|
||
if points.is_empty() {
|
||
return (f64::NAN, f64::NAN);
|
||
}
|
||
let (a, b) = fit(points);
|
||
let gap_eff = 2.0 * (a / b).sqrt();
|
||
(0.5 * (gap - gap_eff) / h, 2.0 * MU * b / F)
|
||
};
|
||
let (off_u, force_u) = profile(&pu);
|
||
let (off_v, _) = profile(&pv);
|
||
let (p0, minus_slope) = fit(&pf);
|
||
let rms = |points: &[(f64, f64)]| {
|
||
(points
|
||
.iter()
|
||
.map(|(s, p)| (p - (p0 - minus_slope * s)).powi(2))
|
||
.sum::<f64>()
|
||
/ points.len().max(1) as f64)
|
||
.sqrt()
|
||
/ (F * gap)
|
||
};
|
||
Reading {
|
||
off_u,
|
||
off_v,
|
||
force_u,
|
||
force_p: 1.0 + minus_slope / F,
|
||
p_full: rms(&pf),
|
||
p_cut: rms(&pc),
|
||
cut_cells,
|
||
merged,
|
||
}
|
||
}
|
||
|
||
/// The linear-exactness mode: in-plane Couette flow `u = K dist t` over ONE
|
||
/// embedded oblique wall (no force, constant pressure, the sides carry the
|
||
/// exact field). A scheme exact on linear fields returns the wall position
|
||
/// to round-off; the fitted zero of the profile on full faces is the offset.
|
||
fn couette(n: usize, slope: f64, c0: f64, s26: Option<bool>) -> (f64, f64) {
|
||
const K: f64 = 1.0;
|
||
let h = 1.0 / n as f64;
|
||
let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
|
||
let norm = (1.0 + slope * slope).sqrt();
|
||
let (tx, ty) = (1.0 / norm, slope / norm);
|
||
let dist = move |x: f64, y: f64| (y - slope * x - c0) / norm;
|
||
let speed = move |x: f64, y: f64| K * dist(x, y).max(0.0);
|
||
let mut solver = Solver::new(
|
||
Fluid {
|
||
density: 1.0,
|
||
viscosity: MU,
|
||
reference_velocity: 1.0,
|
||
reference_length: 1.0,
|
||
},
|
||
parameters(s26),
|
||
);
|
||
solver.set_boundary_velocity(move |x, y, _z, _t| {
|
||
let s = speed(x, y);
|
||
(s * tx, s * ty, 0.0)
|
||
});
|
||
solver.set_body(Body::from_sdf(move |x, y, _z, _t| dist(x, y)));
|
||
let g = Grid::cubic(nx, ny, nz, h);
|
||
let mut field = Field::new(g);
|
||
for k in 0..nz {
|
||
for j in 0..ny {
|
||
for i in 0..=nx {
|
||
field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h);
|
||
}
|
||
}
|
||
for j in 0..=ny {
|
||
for i in 0..nx {
|
||
field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h);
|
||
}
|
||
}
|
||
}
|
||
solver.initialize(&mut field);
|
||
let dt = 0.5 * h * h / (6.0 * MU);
|
||
let steps = (2.0 / dt).ceil() as usize;
|
||
for _ in 0..steps {
|
||
solver.advance(&mut field, dt);
|
||
}
|
||
let mask = solver.mask().expect("mask");
|
||
// u/t_x = K (dist − δ): fit on full faces of the central window, two to
|
||
// six cells off the wall.
|
||
let mut points = Vec::new();
|
||
let mut worst: f64 = 0.0;
|
||
for j in 0..ny {
|
||
for i in 0..nx {
|
||
let (x, y) = (i as f64 * h, (j as f64 + 0.5) * h);
|
||
let f = g.uface(0, j, i);
|
||
let dd = dist(x, y);
|
||
if (x - 0.5 * LX).abs() < 0.3 && dd > 2.0 * h && dd < 6.0 * h && mask.a_u(f) >= 1.0 {
|
||
points.push((dd, field.u[f] / tx));
|
||
worst = worst.max((field.u[f] / tx - K * dd).abs() / (K * h));
|
||
}
|
||
}
|
||
}
|
||
let (a, minus_b) = fit(&points);
|
||
// y = a − (−b) x with b the slope: the zero sits at dist = −a / b.
|
||
let b = -minus_b;
|
||
(-a / b / h, worst)
|
||
}
|
||
|
||
#[test]
|
||
#[ignore = "S2-6 instrument: linear exactness of the cut wall on an oblique wall (in-plane Couette; a minute on the host)"]
|
||
fn oblique_wall_linear_exactness() {
|
||
for slope in [0.0, 0.25, 0.5, 1.0] {
|
||
for n in [16usize, 32] {
|
||
let (off, worst) = couette(n, slope, 0.53, None);
|
||
println!(
|
||
" couette slope {slope:.2} n {n}: wall offset {off:+.4} h (negative = inside the body); worst full-face error {worst:.4} of K·h"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
#[ignore = "S2-6 instrument: the oblique cut wall's effective position with in-plane flow (minutes on the host)"]
|
||
fn oblique_wall_effective_position() {
|
||
let list = |name: &str| -> Option<Vec<f64>> {
|
||
std::env::var(name)
|
||
.ok()
|
||
.map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect())
|
||
};
|
||
// S2-7: `RTX_E3_OBLIQUE_SLOPES=0.25,0.5,1` restricts the slopes;
|
||
// `RTX_E3_OBLIQUE_C0_SHIFTS=0,0.25,0.5,0.75` sweeps the registration
|
||
// (c0 = 0.53 + shift · h at each n; a shift of 1 is the same
|
||
// registration); `RTX_E3_OBLIQUE_ZFLOW=0` skips the z-flow rows.
|
||
let slopes = list("RTX_E3_OBLIQUE_SLOPES");
|
||
let shifts = list("RTX_E3_OBLIQUE_C0_SHIFTS").unwrap_or_else(|| vec![0.0]);
|
||
let zflow = std::env::var("RTX_E3_OBLIQUE_ZFLOW").map_or(true, |v| v != "0");
|
||
for (slope, c0_base) in [
|
||
(0.0, 0.53),
|
||
(0.0, 0.77),
|
||
(0.25, 0.53),
|
||
(0.5, 0.53),
|
||
(1.0, 0.53),
|
||
] {
|
||
if slopes
|
||
.as_ref()
|
||
.is_some_and(|l| !l.iter().any(|s| (s - slope).abs() < 1e-9))
|
||
{
|
||
continue;
|
||
}
|
||
// `RTX_E3_OBLIQUE_N=64` adds a finer rung to the in-plane mode.
|
||
let extra: Option<usize> = std::env::var("RTX_E3_OBLIQUE_N")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok());
|
||
let mut runs = vec![(16usize, false), (32, false), (16, true), (32, true)];
|
||
if let Some(n) = extra {
|
||
runs = vec![(n, false)];
|
||
}
|
||
if !zflow {
|
||
runs.retain(|r| !r.1);
|
||
}
|
||
for (n, along_z) in runs {
|
||
for &shift in &shifts {
|
||
let c0 = c0_base + shift / n as f64;
|
||
let r = reading(n, slope, c0, along_z, None);
|
||
println!(
|
||
" slope {slope:.2} c0 {c0:.5} n {n} {}: wall offset {:+.4} h (u faces) {:+.4} h (v faces), positive = inside the fluid; F_eff/F {:.5} (profile) {:.5} (pressure slope); pressure about its fit: {:.2e} full cells, {:.2e} near-wall cells (of F·G); cut cells {} merged {}",
|
||
if along_z {
|
||
"z-flow (w faces)"
|
||
} else {
|
||
"in-plane"
|
||
},
|
||
r.off_u,
|
||
r.off_v,
|
||
r.force_u,
|
||
r.force_p,
|
||
r.p_full,
|
||
r.p_cut,
|
||
r.cut_cells,
|
||
r.merged
|
||
);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/// S2-7 Step 2 probe: the discrete operators applied to the EXACT in-plane
|
||
/// Poiseuille field valued at the faces' open-part centroids — after one
|
||
/// step with one corrector: the predictor's acceleration per near-wall
|
||
/// face in units of the source's `F/ρ` (zero for a consistent momentum
|
||
/// operator on the exact field), by aperture band; the cut cells'
|
||
/// divergence of the exact centroid VALUES, of the exact open-part MEANS
|
||
/// (the exact fluxes: must vanish), of the predicted field `u*` and of the
|
||
/// corrected field, over the cell's largest face flux, by fluid-fraction
|
||
/// class; and the spurious pressure one projection creates at the cut
|
||
/// cells (of F·G). `u*` is recovered from the one correction.
|
||
#[allow(clippy::too_many_lines)]
|
||
fn probe(n: usize, slope: f64, c0: f64) {
|
||
#[allow(non_snake_case)]
|
||
let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(1e-4);
|
||
let h = 1.0 / n as f64;
|
||
let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
|
||
let norm = (1.0 + slope * slope).sqrt();
|
||
let (tx, ty) = (1.0 / norm, slope / norm);
|
||
let gap = W / norm;
|
||
let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm;
|
||
let speed = move |x: f64, y: f64| {
|
||
let r = r_of(x, y);
|
||
if r.abs() < 0.5 * gap {
|
||
F / (2.0 * MU) * (0.25 * gap * gap - r * r)
|
||
} else {
|
||
0.0
|
||
}
|
||
};
|
||
let mut params = parameters(None);
|
||
params.corrector_steps = 1;
|
||
let rho = 1.0;
|
||
let mut solver = Solver::new(
|
||
Fluid {
|
||
density: rho,
|
||
viscosity: MU,
|
||
reference_velocity: 1.0,
|
||
reference_length: 1.0,
|
||
},
|
||
params,
|
||
);
|
||
solver.set_boundary_velocity(move |x, y, _z, _t| {
|
||
let s = speed(x, y);
|
||
(s * tx, s * ty, 0.0)
|
||
});
|
||
solver.set_momentum_source(move |_, _, _, _| (F * tx, F * ty, 0.0));
|
||
solver.set_body(Body::from_sdf(move |x, y, _z, _t| 0.5 * gap - r_of(x, y).abs()));
|
||
let g = Grid::cubic(nx, ny, nz, h);
|
||
let mut field = Field::new(g);
|
||
solver.initialize(&mut field);
|
||
let (su, sv) = {
|
||
let mask = solver.mask().expect("mask");
|
||
let t = mask
|
||
.face_shift_tables()
|
||
.expect("the centroid shift tables (diffusion_centroid on)");
|
||
(t[0].clone(), t[1].clone())
|
||
};
|
||
// The exact field at the open parts' centroids (a full face: its centre).
|
||
for k in 0..nz {
|
||
for j in 0..ny {
|
||
for i in 0..=nx {
|
||
let f = g.uface(k, j, i);
|
||
field.u[f] = tx * speed(i as f64 * h + su[3 * f], (j as f64 + 0.5) * h + su[3 * f + 1]);
|
||
}
|
||
}
|
||
for j in 0..=ny {
|
||
for i in 0..nx {
|
||
let f = g.vface(k, j, i);
|
||
field.v[f] = ty * speed((i as f64 + 0.5) * h + sv[3 * f], j as f64 * h + sv[3 * f + 1]);
|
||
}
|
||
}
|
||
}
|
||
// The ghost faces (inside the body, in the imposition band) carry the
|
||
// solver's own reconstruction from the exact fluid field, as in a march.
|
||
{
|
||
let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask"));
|
||
mask.impose(body, &mut field.u, &mut field.v, &mut field.w, solver.time());
|
||
}
|
||
let dt = 0.5 * h * h / (6.0 * MU);
|
||
solver.advance(&mut field, dt);
|
||
let mask = solver.mask().expect("mask");
|
||
let pp = &field.p_prime;
|
||
// A full cell's fraction is not exactly 1 (the interpolant's volume).
|
||
let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9;
|
||
let (mut vmin, mut vmax) = (f64::INFINITY, 0.0f64);
|
||
let window = |x: f64| (x - 0.5 * LX).abs() < 0.6;
|
||
// The predictor's value of a face: the corrected value plus the one
|
||
// correction (unknown faces), the face's value otherwise.
|
||
let star_u = |j: usize, i: usize| -> f64 {
|
||
let f = g.uface(0, j, i);
|
||
if i == 0 || i == nx || mask.u_kind(f) != FaceKind::Fluid {
|
||
return field.u[f];
|
||
}
|
||
let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i));
|
||
field.u[f] + (dt / rho) * mask.grad_weight(0, f) * (pp[cp] - pp[cm]) / h
|
||
};
|
||
let star_v = |j: usize, i: usize| -> f64 {
|
||
let f = g.vface(0, j, i);
|
||
if j == 0 || j == ny || mask.v_kind(f) != FaceKind::Fluid {
|
||
return field.v[f];
|
||
}
|
||
let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i));
|
||
field.v[f] + (dt / rho) * mask.grad_weight(1, f) * (pp[cp] - pp[cm]) / h
|
||
};
|
||
// (a) the predictor's acceleration on near-wall faces, of F/ρ.
|
||
let bin_of = |a: f64, near: bool| -> Option<usize> {
|
||
if a < 1.0 {
|
||
Some(((a * 4.0).floor() as usize).min(3))
|
||
} else if near {
|
||
Some(4)
|
||
} else {
|
||
None
|
||
}
|
||
};
|
||
let mut acc: [Vec<f64>; 5] = Default::default();
|
||
for j in 0..ny {
|
||
for i in 1..nx {
|
||
let f = g.uface(0, j, i);
|
||
if !window(i as f64 * h) || mask.u_kind(f) != FaceKind::Fluid {
|
||
continue;
|
||
}
|
||
let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i));
|
||
let near = is_cut(cm) || is_cut(cp);
|
||
let Some(b) = bin_of(mask.a_u(f), near) else { continue };
|
||
acc[b].push(rho * (star_u(j, i) - field.u_old[f]) / dt / (F * tx));
|
||
}
|
||
}
|
||
if slope > 0.0 {
|
||
for j in 1..ny {
|
||
for i in 0..nx {
|
||
let f = g.vface(0, j, i);
|
||
if !window((i as f64 + 0.5) * h) || mask.v_kind(f) != FaceKind::Fluid {
|
||
continue;
|
||
}
|
||
let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i));
|
||
let near = is_cut(cm) || is_cut(cp);
|
||
let Some(b) = bin_of(mask.a_v(f), near) else { continue };
|
||
acc[b].push(rho * (star_v(j, i) - field.v_old[f]) / dt / (F * ty));
|
||
}
|
||
}
|
||
}
|
||
// (b) the cut cells' divergence under four valuations, of the cell's
|
||
// largest face flux; (c) the spurious pressure at cut cells.
|
||
let sigma_u = |a: f64| a * a * h * h / (12.0 * norm * norm);
|
||
let sigma_v = |a: f64| slope * slope * a * a * h * h / (12.0 * norm * norm);
|
||
let class_of = |vol: f64| -> usize {
|
||
if vol < 0.1 {
|
||
0
|
||
} else if vol < 0.5 {
|
||
1
|
||
} else if vol < 1.0 {
|
||
2
|
||
} else {
|
||
3
|
||
}
|
||
};
|
||
let mut div: [Vec<[f64; 4]>; 4] = Default::default();
|
||
let (mut p_cut, mut p_full) = (Vec::new(), Vec::new());
|
||
let area = h * h;
|
||
for j in 1..ny - 1 {
|
||
for i in 1..nx - 1 {
|
||
let c = g.cell(0, j, i);
|
||
if !window((i as f64 + 0.5) * h) || !mask.cell_active(c) {
|
||
continue;
|
||
}
|
||
let vol = mask.vol(c);
|
||
let near = is_cut(c)
|
||
|| [g.cell(0, j, i - 1), g.cell(0, j, i + 1), g.cell(0, j - 1, i), g.cell(0, j + 1, i)]
|
||
.iter()
|
||
.any(|&q| is_cut(q));
|
||
if !near {
|
||
if (r_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h)).abs() < 0.5 * gap {
|
||
p_full.push(pp[c]);
|
||
vmin = vmin.min(vol);
|
||
vmax = vmax.max(vol);
|
||
}
|
||
continue;
|
||
}
|
||
let vol = if is_cut(c) { vol } else { 1.0 };
|
||
let mut d = [0.0; 4];
|
||
let mut scale: f64 = 0.0;
|
||
for (comp, fi, fj, sign) in [(0usize, i + 1, j, 1.0), (0, i, j, -1.0), (1, i, j + 1, 1.0), (1, i, j, -1.0)] {
|
||
let (a, exact, cur, star, sig, t) = if comp == 0 {
|
||
let f = g.uface(0, fj, fi);
|
||
(mask.a_u(f), field.u_old[f], field.u[f], star_u(fj, fi), sigma_u(mask.a_u(f)), tx)
|
||
} else {
|
||
let f = g.vface(0, fj, fi);
|
||
(mask.a_v(f), field.v_old[f], field.v[f], star_v(fj, fi), sigma_v(mask.a_v(f)), ty)
|
||
};
|
||
if a <= 0.0 {
|
||
continue;
|
||
}
|
||
let mean = exact - F / (2.0 * MU) * sig * t;
|
||
d[0] += sign * a * area * exact;
|
||
d[1] += sign * a * area * mean;
|
||
d[2] += sign * a * area * star;
|
||
d[3] += sign * a * area * cur;
|
||
scale = scale.max((a * area * exact).abs());
|
||
}
|
||
if scale > 0.0 {
|
||
div[class_of(vol)].push([d[0] / scale, d[1] / scale, d[2] / scale, d[3] / scale]);
|
||
}
|
||
if is_cut(c) {
|
||
p_cut.push(pp[c]);
|
||
}
|
||
}
|
||
}
|
||
let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::<f64>() / v.len().max(1) as f64).sqrt();
|
||
let maxabs = |v: &[f64]| v.iter().fold(0.0f64, |m, x| m.max(x.abs()));
|
||
let p0 = p_full.iter().sum::<f64>() / p_full.len().max(1) as f64;
|
||
let p_cut: Vec<f64> = p_cut.iter().map(|p| (p - p0) / (F * gap)).collect();
|
||
let p_full: Vec<f64> = p_full.iter().map(|p| (p - p0) / (F * gap)).collect();
|
||
println!(
|
||
" probe slope {slope:.2} c0 {c0:.5} n {n}: merged {}; full cells' fraction {vmin:.3e}–{vmax:.3e}",
|
||
mask.merged_cells()
|
||
);
|
||
let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full, next to a cut cell"];
|
||
for (b, name) in names.iter().enumerate() {
|
||
println!(
|
||
" predictor acceleration of F/ρ on faces {name}: {} faces, rms {:.3e}, max {:.3e}",
|
||
acc[b].len(),
|
||
rms(&acc[b]),
|
||
maxabs(&acc[b])
|
||
);
|
||
}
|
||
let classes = ["vol<0.1 (merged class)", "0.1–0.5", "0.5–1", "full, next to a cut cell"];
|
||
for (k, name) in classes.iter().enumerate() {
|
||
let col = |m: usize| div[k].iter().map(|d| d[m]).collect::<Vec<_>>();
|
||
println!(
|
||
" divergence / largest face flux, cells {name}: {} cells; exact centroid values rms {:.3e}, exact means rms {:.3e}, predicted u* rms {:.3e} max {:.3e}, corrected rms {:.3e}",
|
||
div[k].len(),
|
||
rms(&col(0)),
|
||
rms(&col(1)),
|
||
rms(&col(2)),
|
||
maxabs(&col(2)),
|
||
rms(&col(3))
|
||
);
|
||
}
|
||
println!(
|
||
" spurious pressure after one projection (of F·G): cut cells rms {:.3e} max {:.3e}; full cells rms {:.3e}",
|
||
rms(&p_cut),
|
||
maxabs(&p_cut),
|
||
rms(&p_full)
|
||
);
|
||
}
|
||
|
||
/// The z-flow control of the Step 2 probe: the same channel with the exact
|
||
/// Poiseuille `w(x, y)` at the w faces' open-part centroids (no pressure,
|
||
/// no own-direction variation): the predictor's acceleration on the cut
|
||
/// w faces in units of F/ρ, by aperture band.
|
||
fn probe_z(n: usize, slope: f64, c0: f64) {
|
||
#[allow(non_snake_case)]
|
||
let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
|
||
.ok()
|
||
.and_then(|v| v.parse().ok())
|
||
.unwrap_or(1e-4);
|
||
let h = 1.0 / n as f64;
|
||
let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
|
||
let norm = (1.0 + slope * slope).sqrt();
|
||
let gap = W / norm;
|
||
let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm;
|
||
let speed = move |x: f64, y: f64| {
|
||
let r = r_of(x, y);
|
||
if r.abs() < 0.5 * gap {
|
||
F / (2.0 * MU) * (0.25 * gap * gap - r * r)
|
||
} else {
|
||
0.0
|
||
}
|
||
};
|
||
let mut params = parameters(None);
|
||
params.corrector_steps = 1;
|
||
let rho = 1.0;
|
||
let mut solver = Solver::new(
|
||
Fluid {
|
||
density: rho,
|
||
viscosity: MU,
|
||
reference_velocity: 1.0,
|
||
reference_length: 1.0,
|
||
},
|
||
params,
|
||
);
|
||
solver.set_boundary_velocity(move |x, y, _z, _t| (0.0, 0.0, speed(x, y)));
|
||
solver.set_momentum_source(move |_, _, _, _| (0.0, 0.0, F));
|
||
solver.set_body(Body::from_sdf(move |x, y, _z, _t| 0.5 * gap - r_of(x, y).abs()));
|
||
let g = Grid::cubic(nx, ny, nz, h);
|
||
let mut field = Field::new(g);
|
||
solver.initialize(&mut field);
|
||
let sw = solver
|
||
.mask()
|
||
.expect("mask")
|
||
.face_shift_tables()
|
||
.expect("shift tables")[2]
|
||
.clone();
|
||
for k in 0..=nz {
|
||
for j in 0..ny {
|
||
for i in 0..nx {
|
||
let f = g.wface(k, j, i);
|
||
field.w[f] = speed((i as f64 + 0.5) * h + sw[3 * f], (j as f64 + 0.5) * h + sw[3 * f + 1]);
|
||
}
|
||
}
|
||
}
|
||
{
|
||
let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask"));
|
||
mask.impose(body, &mut field.u, &mut field.v, &mut field.w, solver.time());
|
||
}
|
||
let dt = 0.5 * h * h / (6.0 * MU);
|
||
solver.advance(&mut field, dt);
|
||
let mask = solver.mask().expect("mask");
|
||
let window = |x: f64| (x - 0.5 * LX).abs() < 0.6;
|
||
let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9;
|
||
let mut acc: [Vec<f64>; 5] = Default::default();
|
||
for j in 0..ny {
|
||
for i in 0..nx {
|
||
let f = g.wface(0, j, i);
|
||
if !window((i as f64 + 0.5) * h) || mask.w_kind(f) != FaceKind::Fluid {
|
||
continue;
|
||
}
|
||
let a = mask.a_w(f);
|
||
let c = g.cell(0, j, i);
|
||
let b = if a < 1.0 {
|
||
((a * 4.0).floor() as usize).min(3)
|
||
} else if is_cut(c) {
|
||
4
|
||
} else {
|
||
continue;
|
||
};
|
||
// z-uniform: the projection leaves w alone (dp'/dz = 0).
|
||
acc[b].push(rho * (field.w[f] - field.w_old[f]) / dt / F);
|
||
}
|
||
}
|
||
let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::<f64>() / v.len().max(1) as f64).sqrt();
|
||
let maxabs = |v: &[f64]| v.iter().fold(0.0f64, |m, x| m.max(x.abs()));
|
||
println!(" probe-z slope {slope:.2} c0 {c0:.5} n {n}:");
|
||
let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full, in a cut cell"];
|
||
for (b, name) in names.iter().enumerate() {
|
||
println!(
|
||
" predictor acceleration of F/ρ on w faces {name}: {} faces, rms {:.3e}, max {:.3e}",
|
||
acc[b].len(),
|
||
rms(&acc[b]),
|
||
maxabs(&acc[b])
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
#[ignore = "S2-7 Step 2 probe: the discrete operators on the exact field (seconds per rung on the host)"]
|
||
fn oblique_operator_probe() {
|
||
let list = |name: &str| -> Option<Vec<f64>> {
|
||
std::env::var(name)
|
||
.ok()
|
||
.map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect())
|
||
};
|
||
let slopes = list("RTX_E3_OBLIQUE_SLOPES").unwrap_or_else(|| vec![0.25, 0.5, 1.0]);
|
||
let shifts = list("RTX_E3_OBLIQUE_C0_SHIFTS").unwrap_or_else(|| vec![0.0]);
|
||
let ns: Vec<usize> = list("RTX_E3_OBLIQUE_NS")
|
||
.map(|l| l.iter().map(|&x| x as usize).collect())
|
||
.unwrap_or_else(|| vec![16, 32]);
|
||
for &slope in &slopes {
|
||
for &n in &ns {
|
||
for &shift in &shifts {
|
||
probe(n, slope, 0.53 + shift / n as f64);
|
||
if std::env::var("RTX_E3_OBLIQUE_ZFLOW").map_or(true, |v| v != "0") {
|
||
probe_z(n, slope, 0.53 + shift / n as f64);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/// The gate (S2-6): with the oblique distance, the transverse centroid
|
||
/// correction and the fine floor the cut wall is linear-exact to 0.02 h on
|
||
/// an oblique wall (without them it sits 0.05–0.07 h inside the body at
|
||
/// every h), and the z-directed Poiseuille offset halves per rung at
|
||
/// slope ½ (without them: −0.087 → −0.083 h).
|
||
#[test]
|
||
fn oblique_wall_position_is_second_order() {
|
||
for slope in [0.5, 1.0] {
|
||
let (fixed, _) = couette(32, slope, 0.53, Some(true));
|
||
let (before, _) = couette(32, slope, 0.53, Some(false));
|
||
println!(
|
||
" couette slope {slope}: offset {fixed:+.4} h with the S2-6 closures, {before:+.4} h without"
|
||
);
|
||
assert!(fixed.abs() < 0.02, "slope {slope}: {fixed}");
|
||
assert!(
|
||
before.abs() > 2.0 * fixed.abs(),
|
||
"the instrument lost its contrast"
|
||
);
|
||
}
|
||
let coarse = reading(16, 0.5, 0.53, true, Some(true)).off_u;
|
||
let fine = reading(32, 0.5, 0.53, true, Some(true)).off_u;
|
||
println!(" z-flow slope 0.5: offset {coarse:+.4} h at n 16, {fine:+.4} h at n 32");
|
||
assert!(coarse.abs() < 0.04, "n 16 offset {coarse}");
|
||
assert!(
|
||
fine.abs() < 0.65 * coarse.abs(),
|
||
"the offset does not halve: {coarse} → {fine}"
|
||
);
|
||
}
|