From 501b45f9d0a0e36e42fc8de21fc613b24b79a2d3 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Tue, 22 Sep 2026 07:39:54 -0500 Subject: [PATCH] =?UTF-8?q?R4-i=20+=20R5=20design=20pass:=20RTX=5FFSI2O=5F?= =?UTF-8?q?REGEN=5FONCE=20(one=20patch=20regeneration=20per=20coupled=20st?= =?UTF-8?q?ep;=20dies=20at=20step=200=20=E2=80=94=20the=20registered=20tes?= =?UTF-8?q?t=20is=20void),=20the=20step=20CSV's=20subit/dres=5Fy/dres=5Fno?= =?UTF-8?q?rm/fx=5Fnodal/fy=5Fnodal=20columns;=20the=20cut=20predictor's?= =?UTF-8?q?=20per-face=20term=20probe=20(enable=5Fterm=5Fprobe)=20and=20th?= =?UTF-8?q?e=20curved=20instrument's=20exact=20side/wall=20viscous=20integ?= =?UTF-8?q?rals=20=E2=80=94=20the=20static=20convex=20wall's=20flat=20resi?= =?UTF-8?q?dual=20is=20the=20viscous=20closure's=20first-order=20relative?= =?UTF-8?q?=20accuracy=20on=20O(1/h)=20fluxes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../embedded3/step/cut_predictor.rs | 27 +- .../incompressible/embedded3/step/mod.rs | 22 +- .../tests/embedded3_wall_position_curved.rs | 230 ++++++++++++++++-- .../rtx-fsi/tests/fsi2_harness/overset.rs | 32 +++ .../tests/fsi2_harness/overset_march.rs | 36 ++- 5 files changed, 316 insertions(+), 31 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs index c6d3cb8..ed7282c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs @@ -111,7 +111,8 @@ impl Solver { let cv = mask.cv_geometry(c, p); // S2-7b: the convective sides' own apertures (host prototype). let conv_parts = if mask.conv_sides_exact { - mask.cut().and_then(|cut| mask.exact_cv_side_parts(cut, c, p)) + mask.cut() + .and_then(|cut| mask.exact_cv_side_parts(cut, c, p)) } else { None }; @@ -159,7 +160,9 @@ impl Solver { // The two half faces' own apertures times their faces' values. let v = |q: [i64; 3]| val(d, q).unwrap_or(0.0); ( - 0.5 * (parts[d][1].0 * v(add(cell_minus, ed, 1)) + parts[d][1].1 * v(add(cell_plus, ed, 1))) * a_d, + 0.5 * (parts[d][1].0 * v(add(cell_minus, ed, 1)) + + parts[d][1].1 * v(add(cell_plus, ed, 1))) + * a_d, 0.5 * (parts[d][0].0 * v(cell_minus) + parts[d][0].1 * v(cell_plus)) * a_d, ) } @@ -260,7 +263,8 @@ impl Solver { return 0.0; } let cvq = mask.cv_geometry(c, q); - let ubq = mask.surface_velocity_at(body, foot_of(q, lat.face_position(c, q)), c, t_old); + let ubq = + mask.surface_velocity_at(body, foot_of(q, lat.face_position(c, q)), c, t_old); // Explicit, with a coefficient ∝ 1/d_f: take the gradient from // the faces at least TRANSVERSE_DISTANCE_FLOOR h off the wall — // a FULL neighbour too, over its own distance along the cut @@ -426,6 +430,23 @@ impl Solver { + wall_rhs) / (inertia + shear + wall_implicit); let v_alpha = fraction * h[c] * area[c]; + if let Some(tp) = self.term_probe.borrow_mut().as_mut() { + let m = rho * v_eff; + let terms = [ + -conv / m, + diff / m, + (shear * (ub - u_star) - shear_explicit) / m, + pressure / m, + source / m, + (wall_rhs - wall_implicit * u_star) / m, + (u_star - u0) / dt, + ]; + let v = &mut tp[c]; + if v.len() <= idx_f { + v.resize(idx_f + 1, [0.0; 7]); + } + v[idx_f] = terms; + } (u_star, rho * (v_eff - v_alpha) * (u_star - u0) / dt) } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index d73f582..cb755e8 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -230,7 +230,8 @@ impl Default for Parameters { .is_ok_and(|v| v == "off"), // ON by default since S2-7 (`=h` reproduces the records before it). wall_exchange_axis: std::env::var("RTX_E3_WALL_EXCHANGE").map_or(true, |v| v != "h"), - wall_exchange_foot: std::env::var("RTX_E3_WALL_EXCHANGE").is_ok_and(|v| v == "axisfoot"), + wall_exchange_foot: std::env::var("RTX_E3_WALL_EXCHANGE") + .is_ok_and(|v| v == "axisfoot"), conv_sides_exact: std::env::var("RTX_E3_CONV_SIDES").is_ok_and(|v| v == "exact"), wall_flux_true_normal: std::env::var("RTX_E3_WALL_FLUX").is_ok_and(|v| v == "true"), wall_foot_centroid: std::env::var("RTX_E3_WALL_FOOT").is_ok_and(|v| v == "centroid"), @@ -238,7 +239,8 @@ impl Default for Parameters { momentum_volume_tiled: std::env::var("RTX_E3_MOMENTUM_VOLUME") .is_ok_and(|v| v == "tiled"), cv_sides_exact: std::env::var("RTX_E3_CV_SIDES").is_ok_and(|v| v == "exact"), - wall_order2_centroid: std::env::var("RTX_E3_WALL_ORDER2").is_ok_and(|v| v == "centroid"), + wall_order2_centroid: std::env::var("RTX_E3_WALL_ORDER2") + .is_ok_and(|v| v == "centroid"), // ON by default since S2-5 (`=0` reproduces the records before it). diffusion_centroid: std::env::var("RTX_E3_DIFFUSION_CENTROID") .map_or(true, |v| v != "0"), @@ -269,6 +271,11 @@ pub struct Solver { /// cut predictor only): the box route counts it, the operator route /// does not. pub(super) floor_source: std::cell::Cell<[f64; 3]>, + /// R5 design-pass instrument (host cut predictor only): per unknown + /// face of each component, the last step's momentum terms as + /// accelerations `[−conv, diff, wall shear, pressure, source, solid + /// exchange, total] / (ρ V_eff)`; `None` (the default) records nothing. + pub(super) term_probe: std::cell::RefCell; 3]>>, /// The closure lag of the pressure corrections over the last step /// (`Mask::closure_lag` summed over the correctors; host path). pub(super) pressure_lag: std::cell::Cell<[f64; 3]>, @@ -308,6 +315,7 @@ impl Solver { params, momentum_source: None, floor_source: std::cell::Cell::new([0.0; 3]), + term_probe: std::cell::RefCell::new(None), pressure_lag: std::cell::Cell::new([0.0; 3]), vol_old: Vec::new(), apertures_old: None, @@ -402,6 +410,16 @@ impl Solver { .expect("embedded mask") } + /// R5 instrument: record the cut predictor's momentum terms per face + /// from the next step on (host path only; a no-op on the device path). + pub fn enable_term_probe(&self) { + *self.term_probe.borrow_mut() = Some([Vec::new(), Vec::new(), Vec::new()]); + } + /// The recorded terms per component face index (see `term_probe`). + #[must_use] + pub fn term_probe(&self) -> Option<[Vec<[f64; 7]>; 3]> { + self.term_probe.borrow().clone() + } /// The small-cell floor's momentum source over the last step (a force /// on the fluid; host cut predictor only, zero otherwise). #[must_use] 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 83726f8..5db5f1b 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_curved.rs @@ -58,9 +58,17 @@ impl Exact { let (dx, dy) = (x - CENTRE.0, y - CENTRE.1); let r = (dx * dx + dy * dy).sqrt().max(1e-12); let ut = if r < R1 { - if rigid || !outer_drives() { OMEGA * r } else { 0.0 } + if rigid || !outer_drives() { + OMEGA * r + } else { + 0.0 + } } else if r > R2 { - if rigid || outer_drives() { OMEGA * r } else { 0.0 } + if rigid || outer_drives() { + OMEGA * r + } else { + 0.0 + } } else { self.u_theta(r) }; @@ -145,12 +153,14 @@ fn reading(n: usize, rigid: bool) { for k in 0..nz { for j in 0..ny { for i in 0..=nx { - field.u[g.uface(k, j, i)] = ex.velocity(i as f64 * h, (j as f64 + 0.5) * h, rigid).0; + field.u[g.uface(k, j, i)] = + ex.velocity(i as f64 * h, (j as f64 + 0.5) * h, rigid).0; } } for j in 0..=ny { for i in 0..nx { - field.v[g.vface(k, j, i)] = ex.velocity((i as f64 + 0.5) * h, j as f64 * h, rigid).1; + field.v[g.vface(k, j, i)] = + ex.velocity((i as f64 + 0.5) * h, j as f64 * h, rigid).1; } } } @@ -287,7 +297,8 @@ fn reading(n: usize, rigid: bool) { for i in 0..=nx { let f = g.uface(0, j, i); if mask.u_kind(f) == rtx_cfd::solvers::incompressible::embedded3::FaceKind::Ghost { - let e = (field.u[f] - ex.velocity(i as f64 * h, (j as f64 + 0.5) * h, rigid).0) / (OMEGA * R1); + let e = (field.u[f] - ex.velocity(i as f64 * h, (j as f64 + 0.5) * h, rigid).0) + / (OMEGA * R1); gsq += e * e; gn += 1; gmax = gmax.max(e.abs()); @@ -303,7 +314,13 @@ fn reading(n: usize, rigid: bool) { ); println!( " {} n {n}: walls' offsets {off_in:+.4} h (inner) {off_out:+.4} h (outer), positive = inside the fluid; fit A {a:.5} B {b:.5} (exact {:.5} {:.5}, {} points); pressure error of ρ(ΩR1)²: full cells {:.3e} ({cnt}), cut cells {:.3e} mean {:+.3e} ({n_cut}), fraction < 0.5 {:.3e} ({n_small}), ≥ 0.5 {:.3e} ({n_large}), inner wall {:.3e} ({n_in}), outer wall {:.3e} ({n_out}); merged {}; residual {last_res:.1e}", - if rigid { "rigid" } else if outer_drives() { "outer-driven" } else { "couette" }, + if rigid { + "rigid" + } else if outer_drives() { + "outer-driven" + } else { + "couette" + }, ex.a, ex.b, pts.len(), @@ -374,20 +391,29 @@ fn probe(n: usize, rigid: bool) { let g = Grid::cubic(nx, ny, nz, h); let mut field = Field::new(g); solver.initialize(&mut field); - let tables = solver.mask().expect("mask").face_shift_tables().expect("shift tables").clone(); + let tables = solver + .mask() + .expect("mask") + .face_shift_tables() + .expect("shift tables") + .clone(); for k in 0..nz { for j in 0..ny { for i in 0..=nx { let f = g.uface(k, j, i); let t = &tables[0][3 * f..3 * f + 3]; - field.u[f] = ex.velocity(i as f64 * h + t[0], (j as f64 + 0.5) * h + t[1], rigid).0; + field.u[f] = ex + .velocity(i as f64 * h + t[0], (j as f64 + 0.5) * h + t[1], rigid) + .0; } } for j in 0..=ny { for i in 0..nx { let f = g.vface(k, j, i); let t = &tables[1][3 * f..3 * f + 3]; - field.v[f] = ex.velocity((i as f64 + 0.5) * h + t[0], j as f64 * h + t[1], rigid).1; + field.v[f] = ex + .velocity((i as f64 + 0.5) * h + t[0], j as f64 * h + t[1], rigid) + .1; } } } @@ -399,25 +425,119 @@ fn probe(n: usize, rigid: bool) { } { 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()); + mask.impose( + body, + &mut field.u, + &mut field.v, + &mut field.w, + solver.time(), + ); } let dt = 0.5 * h * h / (6.0 * MU); + // R5 design pass: the predictor's terms per face (viscous group = + // diffusion + wall shear + solid exchange, zero on the exact field; + // inertial group = convection + pressure + source, zero on the exact + // steady field): which group carries the cut layer's residual. + solver.enable_term_probe(); solver.advance(&mut field, dt); + let terms = solver.term_probe().expect("term probe"); let mask = solver.mask().expect("mask"); let pp = &field.p_prime; let a_max = OMEGA * OMEGA * R2; let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9; - // [wall][band]: bands α<¼, ¼–½, ½–¾, ¾–1, full next to cut. - let mut acc: [[Vec; 5]; 2] = Default::default(); + // [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(); + // 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 + // inside it (the wall-shear group) — from the analytic gradient of the + // Taylor–Couette field, as accelerations on the same `V_eff` the + // predictor uses. Their sum is the quadrature error (∇²u = 0). + let grad = |c: usize, x: f64, y: f64| -> [f64; 2] { + let (xp, yp) = (x - CENTRE.0, y - CENTRE.1); + let r = (xp * xp + yp * yp).sqrt().max(1e-12); + let (f, fp) = if rigid { + (OMEGA, 0.0) + } else { + (ex.a + ex.b / (r * r), -2.0 * ex.b / (r * r * r)) + }; + if c == 0 { + [-fp * (xp / r) * yp, -f - fp * (yp / r) * yp] + } else { + [f + fp * (xp / r) * xp, fp * (yp / r) * xp] + } + }; + let in_gap = |x: f64, y: f64| { + 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 { + let m = 4096; + let len = ((x1 - x0).powi(2) + (y1 - y0).powi(2)).sqrt(); + let mut s = 0.0; + for q in 0..m { + 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); + 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 mut s = 0.0; + let m = 1 << 20; + for (r, sign) in [(R1, -1.0), (R2, 1.0)] { + 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 { + let g = grad(c, 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; + } + } + } + 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); + let v_eff = a.max(0.1) * h * h * h; + (d * h / (RHO * v_eff), w * h / (RHO * v_eff)) + }; let bin_of = |a: f64, near: bool| -> Option { if a < 1.0 { Some(((a * 4.0).floor() as usize).min(3)) } else if near { Some(4) } else { - None + Some(5) } }; + // [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] { + [ + total / a_max, + (t[1] + t[2] + t[5]) / a_max, + (t[0] + t[3] + t[4]) / a_max, + t[1] / a_max, + t[2] / a_max, + t[5] / a_max, + (t[1] + t[5] - exact.0) / a_max, + (t[2] - exact.1) / a_max, + (exact.0 + exact.1) / a_max, + ] + }; let wall_of = |x: f64, y: f64| usize::from(r_of(x, y) >= 0.5 * (R1 + R2)); for j in 0..ny { for i in 1..nx { @@ -426,9 +546,24 @@ fn probe(n: usize, rigid: bool) { continue; } let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i)); - let Some(b) = bin_of(mask.a_u(f), is_cut(cm) || is_cut(cp)) else { continue }; + let Some(b) = bin_of(mask.a_u(f), is_cut(cm) || is_cut(cp)) else { + continue; + }; let star = field.u[f] + (dt / RHO) * mask.grad_weight(0, f) * (pp[cp] - pp[cm]) / h; - acc[wall_of(i as f64 * h, (j as f64 + 0.5) * h)][b].push((star - field.u_old[f]) / dt / a_max); + let t = terms[0].get(f).copied().unwrap_or([0.0; 7]); + let exact = exact_terms( + 0, + (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), + ); + acc[wall_of(i as f64 * h, (j as f64 + 0.5) * h)][b].push(groups( + &t, + (star - field.u_old[f]) / dt, + exact, + )); } } for j in 1..ny { @@ -438,9 +573,24 @@ fn probe(n: usize, rigid: bool) { continue; } let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i)); - let Some(b) = bin_of(mask.a_v(f), is_cut(cm) || is_cut(cp)) else { continue }; + let Some(b) = bin_of(mask.a_v(f), is_cut(cm) || is_cut(cp)) else { + continue; + }; let star = field.v[f] + (dt / RHO) * mask.grad_weight(1, f) * (pp[cp] - pp[cm]) / h; - acc[wall_of((i as f64 + 0.5) * h, j as f64 * h)][b].push((star - field.v_old[f]) / dt / a_max); + let t = terms[1].get(f).copied().unwrap_or([0.0; 7]); + let exact = exact_terms( + 1, + 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), + ); + acc[wall_of((i as f64 + 0.5) * h, j as f64 * h)][b].push(groups( + &t, + (star - field.v_old[f]) / dt, + exact, + )); } } let scale = RHO * (OMEGA * R1).powi(2); @@ -464,23 +614,55 @@ fn probe(n: usize, rigid: bool) { continue; } let e = (pp[c] - lvl) / scale; - let w = if is_cut(c) { wall_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h) } else { 2 }; + let w = if is_cut(c) { + wall_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h) + } else { + 2 + }; psq[w] += e * e; pn[w] += 1; } } - let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::() / v.len().max(1) as f64).sqrt(); - let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full next to cut"]; - let mode = if rigid { "rigid" } else if outer_drives() { "outer-driven" } else { "couette" }; + let rms = |v: &[[f64; 9]], 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"]; + let mode = if rigid { + "rigid" + } else if outer_drives() { + "outer-driven" + } else { + "couette" + }; for (w, wname) in ["inner wall", "outer wall"].iter().enumerate() { let mut line = format!(" probe {mode} n {n} {wname}:"); for (b, name) in names.iter().enumerate() { - line += &format!(" {name} {} rms {:.3e};", acc[w][b].len(), rms(&acc[w][b])); + line += &format!( + " {name} {} rms {:.3e} (visc {:.3e} inert {:.3e}; diff {:.3e} shear {:.3e} exch {:.3e}; ERR sides {:.3e} wall {:.3e} quad {:.1e});", + acc[w][b].len(), + rms(&acc[w][b], 0), + rms(&acc[w][b], 1), + rms(&acc[w][b], 2), + rms(&acc[w][b], 3), + rms(&acc[w][b], 4), + rms(&acc[w][b], 5), + rms(&acc[w][b], 6), + rms(&acc[w][b], 7), + rms(&acc[w][b], 8) + ); } - line += &format!(" p' at its cut cells {:.3e} ({})", (psq[w] / pn[w].max(1) as f64).sqrt(), pn[w]); + line += &format!( + " p' at its cut cells {:.3e} ({})", + (psq[w] / pn[w].max(1) as f64).sqrt(), + pn[w] + ); println!("{line}"); } - println!(" probe {mode} n {n} interior p' {:.3e} ({})", (psq[2] / pn[2].max(1) as f64).sqrt(), pn[2]); + println!( + " probe {mode} n {n} interior p' {:.3e} ({})", + (psq[2] / pn[2].max(1) as f64).sqrt(), + pn[2] + ); } #[test] diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs index 194e150..b6604d7 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs @@ -240,6 +240,15 @@ pub struct OversetFluid { /// The last fluid step's mass defects (patch acceptor ring, background /// fringe) and the patch's max divergence. pub last_defects: Cell<(f64, f64, f64)>, + /// R4-i `RTX_FSI2O_REGEN_ONCE`: one patch regeneration per coupled + /// step — the first pass's meshes (one per substep) are reused by every + /// later subiteration of the same step (the wall polygon and velocity + /// still follow the candidate). Active between `begin_coupled_step` + /// calls only; `restore` rewinds the cursor to the step's first substep. + pub regen_once: bool, + pub regen_hold: bool, + pub regen_cache: Vec, + pub regen_cursor: usize, } impl OversetFluid { @@ -503,6 +512,10 @@ impl OversetFluid { correctors_total: Cell::new(0), last_d: zero_d, last_defects: Cell::new((0.0, 0.0, 0.0)), + regen_once: std::env::var("RTX_FSI2O_REGEN_ONCE").is_ok(), + regen_hold: false, + regen_cache: Vec::new(), + regen_cursor: 0, warm_base: None, warm_sweeps, }) @@ -918,6 +931,12 @@ impl OversetFluid { self.last_d = d.to_vec(); return Ok(()); } + if self.regen_once && self.regen_hold && self.regen_cursor < self.regen_cache.len() { + let mesh = self.regen_cache[self.regen_cursor].clone(); + self.regen_cursor += 1; + self.last_d = d.to_vec(); + return self.solver.set_patch_mesh(mesh); + } let mesh = match self.patch_for(d) { Ok(m) => m, Err(e) => { @@ -925,10 +944,22 @@ impl OversetFluid { return Err(e); } }; + if self.regen_once && self.regen_hold { + self.regen_cache.push(mesh.clone()); + self.regen_cursor += 1; + } self.last_d = d.to_vec(); self.solver.set_patch_mesh(mesh) } + /// R4-i: a new coupled step begins — the first pass regenerates, the + /// later passes reuse (`RTX_FSI2O_REGEN_ONCE`; a no-op otherwise). + pub fn begin_coupled_step(&mut self) { + self.regen_hold = self.regen_once; + self.regen_cache.clear(); + self.regen_cursor = 0; + } + /// The last geometry, for the offline reproduction /// (`patch_cylinder_flag_deformed.rs`, `RTX_CF_EDGES_FILE`): one edge /// per block, `x y` per line, when `RTX_FSI2O_DUMP_DIR` is set. @@ -1198,6 +1229,7 @@ impl OversetFluid { let t0 = std::time::Instant::now(); self.solver.restore(&saved.0); self.field = saved.1.clone(); + self.regen_cursor = 0; self.t_restore .set(self.t_restore.get() + t0.elapsed().as_secs_f64()); } diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs index 5f8ba2a..458ac41 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs @@ -137,6 +137,11 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov std::env::var("RTX_FSI2O_FICT_MASS").unwrap_or_else(|_| "0".into()), std::env::var("RTX_FSI2O_ROBIN_ALPHA").unwrap_or_else(|_| "0".into()), ); + if std::env::var("RTX_FSI2O_REGEN_ONCE").is_ok() { + println!( + " patch regeneration: ONCE per coupled step (R4-i, RTX_FSI2O_REGEN_ONCE): the first pass's mesh is reused by the later subiterations" + ); + } // Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the // march with the saved state; `RTX_FSI2O_SAVE=dir` saves it). @@ -399,7 +404,11 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov // cell imbalance. let mut step_csv = std::env::var("RTX_FSI2O_STEP_CSV").ok().map(|p| { let mut f = std::fs::File::create(p).expect("step csv"); - writeln!(f, "t,ux,uy,drag,lift,defect_patch,defect_bg,patch_div").unwrap(); + // R4-i columns: the step's subiterations, the last pass's signed + // y-residual Σ_k (d_new − d_candidate)_y and its norm, and the + // nodal load the structure was given (Σ F_x, Σ F_y) — against the + // wall-integral force of the same state (drag, lift). + writeln!(f, "t,ux,uy,drag,lift,defect_patch,defect_bg,patch_div,subit,dres_y,dres_norm,fx_nodal,fy_nodal").unwrap(); f }); let (t_fluid, t_structure) = (std::cell::Cell::new(0.0_f64), std::cell::Cell::new(0.0_f64)); @@ -433,6 +442,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov if robin_alpha > 0.0 { robin_datum.replace(fluid.borrow().inner_tractions()); } + fluid.borrow_mut().begin_coupled_step(); let saved = fluid.borrow().snapshot(); type PassResult = (DynamicState, Vec<(NodeId, Vector3)>, f64, usize); let latest: RefCell> = RefCell::new(None); @@ -498,6 +508,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov let tol_step = cfg.tol_floor.max(cfg.rtol * increment); let retry_at = (5.0 * tol_step).max(0.1 * increment); let acceptable = (cfg.stall_accept * tol_step).max(0.1 * increment); + let mut step_iterations = 0usize; let mut outcome = if let Some(iqn) = iqn.as_mut() { iqn.set_tolerance(tol_step).unwrap(); iqn.solve(&d_predicted, pass) @@ -533,6 +544,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov Ok(c) => { total_subiterations += c.iterations; max_subiterations = max_subiterations.max(c.iterations); + step_iterations = c.iterations; } Err( rtx_fsi::FsiError::CouplingNotConverged { @@ -549,6 +561,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov worst_stall = worst_stall.max(residual); total_subiterations += iterations; max_subiterations = max_subiterations.max(iterations); + step_iterations = iterations; } Err(e) => { println!( @@ -560,6 +573,25 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov } } let (new_state, nodal, conservation, faces) = latest.borrow_mut().take().expect("pass ran"); + // R4-i: the last pass's residual (the structure's answer against the + // candidate the fluid's committed state sits at) and the nodal load. + let (dres_y, dres_norm, fx_nodal, fy_nodal) = { + let d_new = extract(&new_state); + let d_cand = &fluid.borrow().last_d; + let mut sy = 0.0_f64; + let mut n2 = 0.0_f64; + for k in 0..d_new.len() / 2 { + sy += d_new[2 * k + 1] - d_cand[2 * k + 1]; + n2 += (d_new[2 * k] - d_cand[2 * k]).powi(2) + + (d_new[2 * k + 1] - d_cand[2 * k + 1]).powi(2); + } + let (mut fx, mut fy) = (0.0_f64, 0.0_f64); + for (_, f) in &nodal { + fx += f.x; + fy += f.y; + } + (sy, n2.sqrt(), fx, fy) + }; flag_state = new_state; committed_nodal = nodal; prev_area.set(fluid.borrow().shared.read().unwrap().area()); @@ -590,7 +622,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov let (dp, db, pd) = fluid.borrow().last_defects.get(); writeln!( f, - "{t_now:.6},{ux:.6e},{uy:.6e},{drag_now:.6e},{lift_now:.6e},{dp:.3e},{db:.3e},{pd:.3e}" + "{t_now:.6},{ux:.6e},{uy:.6e},{drag_now:.6e},{lift_now:.6e},{dp:.3e},{db:.3e},{pd:.3e},{step_iterations},{dres_y:.3e},{dres_norm:.3e},{fx_nodal:.6e},{fy_nodal:.6e}" ) .unwrap(); }