From e9fe155c4bd674fe8235a1e8222f9165c56a4156 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Tue, 22 Sep 2026 08:05:33 -0500 Subject: [PATCH] =?UTF-8?q?R5=20design=20pass:=20the=20curved=20probe's=20?= =?UTF-8?q?least-squares=20reconstruction=20instrument=20(fit=5Fterms;=20R?= =?UTF-8?q?TX=5FE3=5FFIT=5FORDER=3D3=20cubic)=20=E2=80=94=20the=20candidat?= =?UTF-8?q?e=20cut-face=20closure's=20flux=20truncation=20against=20the=20?= =?UTF-8?q?exact=20box=20integrals?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../tests/embedded3_wall_position_curved.rs | 219 ++++++++++++++++-- 1 file changed, 195 insertions(+), 24 deletions(-) diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs index 5db5f1b..f89f73d 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs @@ -8,10 +8,10 @@ //! exact p(r) (the S2-7b form). `RTX_E3_CURVED_MODE=rigid` is the //! linear-exactness mode: solid-body rotation of both cylinders //! (`u = Ω r e_θ`, `p = ρ Ω² r²/2`). -use rtx_cfd::solvers::incompressible::ConvectionScheme; use rtx_cfd::solvers::incompressible::embedded3::{ Body, Boundaries, FaceKind, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, }; +use rtx_cfd::solvers::incompressible::ConvectionScheme; const MU: f64 = 0.1; const RHO: f64 = 1.0; @@ -447,7 +447,7 @@ fn probe(n: usize, rigid: bool) { let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9; // [wall][band]: bands α<¼, ¼–½, ½–¾, ¾–1, full next to cut, interior; // per face [total, viscous group, inertial group] / a_max. - let mut acc: [[Vec<[f64; 9]>; 6]; 2] = Default::default(); + let mut acc: [[Vec<[f64; 11]>; 6]; 2] = Default::default(); // R5 design pass: the EXACT viscous line integrals over the momentum // control volume of a face — the open parts of its four sides (the // side-diffusion group: diffusion + solid exchange) and the wall arc @@ -472,7 +472,8 @@ fn probe(n: usize, rigid: bool) { let r = r_of(x, y); (R1..=R2).contains(&r) }; - let side = |c: usize, x0: f64, y0: f64, x1: f64, y1: f64, n: [f64; 2]| -> f64 { + type GradFn<'a> = &'a dyn Fn(f64, f64) -> [f64; 2]; + let side = |x0: f64, y0: f64, x1: f64, y1: f64, n: [f64; 2], gf: GradFn| -> f64 { let m = 4096; let len = ((x1 - x0).powi(2) + (y1 - y0).powi(2)).sqrt(); let mut s = 0.0; @@ -480,13 +481,13 @@ fn probe(n: usize, rigid: bool) { let t = (q as f64 + 0.5) / m as f64; let (x, y) = (x0 + t * (x1 - x0), y0 + t * (y1 - y0)); if in_gap(x, y) { - let g = grad(c, x, y); + let g = gf(x, y); s += (g[0] * n[0] + g[1] * n[1]) * len / m as f64; } } MU * s }; - let arc = |c: usize, xa: f64, xb: f64, ya: f64, yb: f64| -> f64 { + let arc = |xa: f64, xb: f64, ya: f64, yb: f64, gf: GradFn| -> f64 { let mut s = 0.0; let m = 1 << 20; for (r, sign) in [(R1, -1.0), (R2, 1.0)] { @@ -494,7 +495,7 @@ fn probe(n: usize, rigid: bool) { let th = std::f64::consts::TAU * (q as f64 + 0.5) / m as f64; let (x, y) = (CENTRE.0 + r * th.cos(), CENTRE.1 + r * th.sin()); if x >= xa && x < xb && y >= ya && y < yb { - let g = grad(c, x, y); + let g = gf(x, y); let n = [sign * th.cos(), sign * th.sin()]; s += (g[0] * n[0] + g[1] * n[1]) * r * std::f64::consts::TAU / m as f64; } @@ -502,17 +503,173 @@ fn probe(n: usize, rigid: bool) { } MU * s }; - // (side-diffusion, wall-shear) exact accelerations for the face of - // component `c` whose CV is the box [xa, xb] × [ya, yb], aperture `a`. - let exact_terms = |c: usize, xa: f64, xb: f64, ya: f64, yb: f64, a: f64| -> (f64, f64) { - let d = side(c, xa, ya, xa, yb, [-1.0, 0.0]) - + side(c, xb, ya, xb, yb, [1.0, 0.0]) - + side(c, xa, ya, xb, ya, [0.0, -1.0]) - + side(c, xa, yb, xb, yb, [0.0, 1.0]); - let w = arc(c, xa, xb, ya, yb); + // (side-diffusion, wall-shear) accelerations for the face of + // component `c` whose CV is the box [xa, xb] × [ya, yb], aperture `a`, + // from a gradient field `gf`. + let box_terms = |xa: f64, xb: f64, ya: f64, yb: f64, a: f64, gf: GradFn| -> (f64, f64) { + let d = side(xa, ya, xa, yb, [-1.0, 0.0], gf) + + side(xb, ya, xb, yb, [1.0, 0.0], gf) + + side(xa, ya, xb, ya, [0.0, -1.0], gf) + + side(xa, yb, xb, yb, [0.0, 1.0], gf); + let w = arc(xa, xb, ya, yb, gf); let v_eff = a.max(0.1) * h * h * h; (d * h / (RHO * v_eff), w * h / (RHO * v_eff)) }; + let exact_terms = |c: usize, xa: f64, xb: f64, ya: f64, yb: f64, a: f64| -> (f64, f64) { + box_terms(xa, xb, ya, yb, a, &|x, y| grad(c, x, y)) + }; + // R5-b candidate at the reconstruction level: a weighted least-squares + // QUADRATIC of the exact face values at the open neighbours' centroids + // (a 5 × 5 face neighbourhood) with the wall's velocity at three points + // of the arc inside the box as strong constraints; its gradient gives + // the same box fluxes. Its error against the exact integrals is the + // truncation a second-order cut-face closure would leave. + let fit_terms = + |c: usize, i: usize, j: usize, xa: f64, xb: f64, ya: f64, yb: f64, a: f64| -> (f64, f64) { + let mut pts: Vec<(f64, f64, f64, f64)> = Vec::new(); + for dj in -2i64..=2 { + for di in -2i64..=2 { + let (ii, jj) = (i as i64 + di, j as i64 + dj); + if ii < 0 || jj < 0 { + continue; + } + let (ii, jj) = (ii as usize, jj as usize); + let (ok, f, x, y) = if c == 0 { + if ii > nx || jj >= ny { + continue; + } + let f = g.uface(0, jj, ii); + ( + mask.u_kind(f) == FaceKind::Fluid && mask.a_u(f) > 0.0, + f, + ii as f64 * h, + (jj as f64 + 0.5) * h, + ) + } else { + if ii >= nx || jj > ny { + continue; + } + let f = g.vface(0, jj, ii); + ( + mask.v_kind(f) == FaceKind::Fluid && mask.a_v(f) > 0.0, + f, + (ii as f64 + 0.5) * h, + jj as f64 * h, + ) + }; + if !ok { + continue; + } + let t = &tables[c][3 * f..3 * f + 3]; + let val = if c == 0 { + field.u_old[f] + } else { + field.v_old[f] + }; + pts.push((x + t[0], y + t[1], val, 1.0)); + } + } + // The wall constraints: three points of each arc inside the box. + for (r, _) in [(R1, -1.0), (R2, 1.0)] { + let m = 1 << 14; + let mut ths: Vec = Vec::new(); + for q in 0..m { + let th = std::f64::consts::TAU * (q as f64 + 0.5) / m as f64; + let (x, y) = (CENTRE.0 + r * th.cos(), CENTRE.1 + r * th.sin()); + if x >= xa && x < xb && y >= ya && y < yb { + ths.push(th); + } + } + if ths.is_empty() { + continue; + } + let (t0, t1) = (ths[0], ths[ths.len() - 1]); + for th in [t0, 0.5 * (t0 + t1), t1] { + let (x, y) = (CENTRE.0 + r * th.cos(), CENTRE.1 + r * th.sin()); + let v = ex.velocity(x, y, rigid); + pts.push((x, y, if c == 0 { v.0 } else { v.1 }, 1.0e4)); + } + } + let (xc, yc) = (0.5 * (xa + xb), 0.5 * (ya + yb)); + // `RTX_E3_FIT_ORDER=3`: a cubic (10 terms) instead of the quadratic (6). + let nb = if std::env::var("RTX_E3_FIT_ORDER").is_ok_and(|v| v == "3") { + 10 + } else { + 6 + }; + let basis = |x: f64, y: f64| -> [f64; 10] { + let (u, v) = ((x - xc) / h, (y - yc) / h); + [ + 1.0, + u, + v, + u * u, + u * v, + v * v, + u * u * u, + u * u * v, + u * v * v, + v * v * v, + ] + }; + let mut m = [[0.0f64; 11]; 10]; + for &(x, y, val, w0) in &pts { + let b = basis(x, y); + let rho2 = ((x - xc) / h).powi(2) + ((y - yc) / h).powi(2); + let w = w0 / (1.0 + rho2); + for r in 0..nb { + for cc in 0..nb { + m[r][cc] += w * b[r] * b[cc]; + } + m[r][10] += w * b[r] * val; + } + } + // Gaussian elimination with partial pivoting. + for k in 0..nb { + let piv = (k..nb) + .max_by(|&p1, &p2| m[p1][k].abs().partial_cmp(&m[p2][k].abs()).unwrap()) + .unwrap(); + m.swap(k, piv); + let d = m[k][k]; + if d.abs() < 1e-300 { + continue; + } + for r in 0..nb { + if r != k { + let fct = m[r][k] / d; + for cc in k..11 { + m[r][cc] -= fct * m[k][cc]; + } + } + } + } + let mut coef = [0.0f64; 10]; + for k in 0..nb { + if m[k][k].abs() >= 1e-300 { + coef[k] = m[k][10] / m[k][k]; + } + } + let gf = move |x: f64, y: f64| -> [f64; 2] { + let (u, v) = ((x - xc) / h, (y - yc) / h); + [ + (coef[1] + + 2.0 * coef[3] * u + + coef[4] * v + + 3.0 * coef[6] * u * u + + 2.0 * coef[7] * u * v + + coef[8] * v * v) + / h, + (coef[2] + + coef[4] * u + + 2.0 * coef[5] * v + + coef[7] * u * u + + 2.0 * coef[8] * u * v + + 3.0 * coef[9] * v * v) + / h, + ] + }; + box_terms(xa, xb, ya, yb, a, &gf) + }; let bin_of = |a: f64, near: bool| -> Option { if a < 1.0 { Some(((a * 4.0).floor() as usize).min(3)) @@ -525,7 +682,7 @@ fn probe(n: usize, rigid: bool) { // [total, viscous group, inertial group, diffusion, wall shear, solid exchange] / a_max. // …, then the side-diffusion error (diff + exch − exact), the wall-shear // error (shear − exact) and the exact sum (the quadrature's own error). - let groups = |t: &[f64; 7], total: f64, exact: (f64, f64)| -> [f64; 9] { + let groups = |t: &[f64; 7], total: f64, exact: (f64, f64), fit: (f64, f64)| -> [f64; 11] { [ total / a_max, (t[1] + t[2] + t[5]) / a_max, @@ -536,6 +693,8 @@ fn probe(n: usize, rigid: bool) { (t[1] + t[5] - exact.0) / a_max, (t[2] - exact.1) / a_max, (exact.0 + exact.1) / a_max, + (fit.0 - exact.0) / a_max, + (fit.1 - exact.1) / a_max, ] }; let wall_of = |x: f64, y: f64| usize::from(r_of(x, y) >= 0.5 * (R1 + R2)); @@ -551,18 +710,23 @@ fn probe(n: usize, rigid: bool) { }; let star = field.u[f] + (dt / RHO) * mask.grad_weight(0, f) * (pp[cp] - pp[cm]) / h; let t = terms[0].get(f).copied().unwrap_or([0.0; 7]); - let exact = exact_terms( - 0, + let (xa, xb, ya, yb) = ( (i as f64 - 0.5) * h, (i as f64 + 0.5) * h, j as f64 * h, (j as f64 + 1.0) * h, - mask.a_u(f), ); + let exact = exact_terms(0, xa, xb, ya, yb, mask.a_u(f)); + let fit = if b < 5 { + fit_terms(0, i, j, xa, xb, ya, yb, mask.a_u(f)) + } else { + exact + }; acc[wall_of(i as f64 * h, (j as f64 + 0.5) * h)][b].push(groups( &t, (star - field.u_old[f]) / dt, exact, + fit, )); } } @@ -578,18 +742,23 @@ fn probe(n: usize, rigid: bool) { }; let star = field.v[f] + (dt / RHO) * mask.grad_weight(1, f) * (pp[cp] - pp[cm]) / h; let t = terms[1].get(f).copied().unwrap_or([0.0; 7]); - let exact = exact_terms( - 1, + let (xa, xb, ya, yb) = ( i as f64 * h, (i as f64 + 1.0) * h, (j as f64 - 0.5) * h, (j as f64 + 0.5) * h, - mask.a_v(f), ); + let exact = exact_terms(1, xa, xb, ya, yb, mask.a_v(f)); + let fit = if b < 5 { + fit_terms(1, i, j, xa, xb, ya, yb, mask.a_v(f)) + } else { + exact + }; acc[wall_of((i as f64 + 0.5) * h, j as f64 * h)][b].push(groups( &t, (star - field.v_old[f]) / dt, exact, + fit, )); } } @@ -623,7 +792,7 @@ fn probe(n: usize, rigid: bool) { pn[w] += 1; } } - let rms = |v: &[[f64; 9]], k: usize| { + let rms = |v: &[[f64; 11]], k: usize| { (v.iter().map(|x| x[k] * x[k]).sum::() / v.len().max(1) as f64).sqrt() }; let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full next to cut", "interior"]; @@ -638,7 +807,7 @@ fn probe(n: usize, rigid: bool) { let mut line = format!(" probe {mode} n {n} {wname}:"); for (b, name) in names.iter().enumerate() { line += &format!( - " {name} {} rms {:.3e} (visc {:.3e} inert {:.3e}; diff {:.3e} shear {:.3e} exch {:.3e}; ERR sides {:.3e} wall {:.3e} quad {:.1e});", + " {name} {} rms {:.3e} (visc {:.3e} inert {:.3e}; diff {:.3e} shear {:.3e} exch {:.3e}; ERR sides {:.3e} wall {:.3e} quad {:.1e}; FIT sides {:.3e} wall {:.3e});", acc[w][b].len(), rms(&acc[w][b], 0), rms(&acc[w][b], 1), @@ -648,7 +817,9 @@ fn probe(n: usize, rigid: bool) { rms(&acc[w][b], 5), rms(&acc[w][b], 6), rms(&acc[w][b], 7), - rms(&acc[w][b], 8) + rms(&acc[w][b], 8), + rms(&acc[w][b], 9), + rms(&acc[w][b], 10) ); } line += &format!(