rtx-cfd: overset P3b — the reclassification impulse located (the fringe ring is a staircase of the interpolated velocities' mass defect) and removed by a converged fringe flux balance (default on): falsifier max spike 594 → 5.50 N/m at the FSI2 step (staircase 6490), 10.95 / 16.79 at dt/2 / dt/4 (12600 / 25600), rms spike 0.07% of the force, far probe 6 (7900), KE per event 4.9e-3 J/m falling with Δt (2.6 fixed)

OverlapMap::balance_fringe_fluxes: Gauss–Seidel through the prescribed faces of every fringe cell to 1e-12 of the prescribed flux scale (≤ 50 sweeps), after every fringe stamping (3 fixed sweeps 101 N/m, 10 sweeps 5.5 — converged is the rule). OversetParameters: fringe_flux_balance (default on, RTX_OVERSET_NO_BALANCE off), fringe_balance_tolerance, refill_turned_active (measured no effect: 593.7 → 593.8; kept as the record), stall_rounds opt-in. P3b locating trace RTX_OVERSET_TRACE_SP (continuity source by class change in cell volumes/step, stored-pressure jump of turned-active cells): the flipped cells' mass source ≤ 6e-3 cell volumes/step, their stored pressure 5–10% of the range off their neighbours (4.4% on the static MMS — the meshes' discretization disagreement). Knock-outs refuted (RTX_OVERSET_H1 keep own face velocities, H4 no warm start, pressure refill): 593–597 N/m each. S4 MMS with the balance: velocity errors within 0.1% of the pinned values, the background's overlap mass defect 1e-13 by construction, pressure errors unchanged. overset_mms prints pressure diagnostics; overset_falsifier records the balanced ladder (regression guard 20 N/m at dt; RTX_OVERSET_FALSIFIER_STRICT asserts the registered gates — (ii) holds at dt, misses at dt/2, dt/4; (iv) fails: residual ∝ 1/Δt^0.8). Co-Authored-By: Claude Fable 5.1 <[email protected]> Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2
2026-09-05 19:16:16 -07:00
parent 6b8837301f
commit 62df6bd628
5 changed files with 372 additions and 15 deletions
@@ -60,6 +60,20 @@ pub struct OversetParameters {
    /// transient — on the falsifier plate it stopped rounds that were still
    /// converging and raised the max force spike from 594 to 4044 N/m.
    pub stall_rounds: usize,
+    /// When a background cell turns active, refill its pressure from its
+    /// neighbours that stayed active (default on; `RTX_OVERSET_NO_REFILL`
+    /// keeps the patch's interpolated value it held as a fringe cell — the
+    /// P3 falsifier's reclassification impulse, §5.10).
+    pub refill_turned_active: bool,
+    /// Fringe flux balance (default on): after every stamping of the
+    /// prescribed background faces, sweep the fringe cells divergence-free
+    /// through their prescribed faces, to `fringe_balance_tolerance` × the
+    /// prescribed flux scale or 50 sweeps. Off (`RTX_OVERSET_NO_BALANCE`)
+    /// reproduces the A-P3 reclassification impulse (max spike 594 N/m);
+    /// 3 fixed sweeps gave 101, 10 gave 5.5 — converged is the rule.
+    pub fringe_flux_balance: bool,
+    /// Relative stop of the balance sweeps.
+    pub fringe_balance_tolerance: f64,
    /// Patch rows below the acceptor row kept non-hole
    /// (`overlap::DEFAULT_OVERLAP_ROWS`).
    pub overlap_rows: usize,
@@ -73,6 +87,9 @@ impl Default for OversetParameters {
            max_rounds: 20,
            anderson_depth: 3,
            stall_rounds: 0,
+            refill_turned_active: std::env::var("RTX_OVERSET_NO_REFILL").is_err(),
+            fringe_flux_balance: std::env::var("RTX_OVERSET_NO_BALANCE").is_err(),
+            fringe_balance_tolerance: 1e-12,
            overlap_rows: overlap::DEFAULT_OVERLAP_ROWS,
        }
    }
@@ -308,6 +325,106 @@ impl OversetPisoSolver {
        self.background.initialize(&mut field.background)
    }

+    /// P3b locating trace: after the predictor of a reclassification step,
+    /// the continuity source (the divergence of `u*`) on the background
+    /// cells by class change, in units of one cell volume per step
+    /// (`h² / dt`); and the stored pressure of the cells that just turned
+    /// active against the mean of their active neighbours, relative to the
+    /// pressure range.
+    fn trace_reclassification_source(&self, field: &OversetField, old: &[CellClass], dt: f64) {
+        let (nx, ny, dx, dy) = self.grid;
+        let unit = dx * dy / dt;
+        let bg = &field.background;
+        let div = |j: usize, i: usize| {
+            ((bg.u_star[(j, i + 1)] - bg.u_star[(j, i)]) * dy
+                + (bg.v_star[(j + 1, i)] - bg.v_star[(j, i)]) * dx)
+                / unit
+        };
+        let class_now = |j: usize, i: usize| self.overlap.class(j, i);
+        let changed = |j: usize, i: usize| old[j * nx + i] != class_now(j, i);
+        let (mut n_fa, mut max_fa, mut sum_fa) = (0usize, 0.0_f64, 0.0_f64);
+        let (mut n_nb, mut max_nb, mut sum_nb) = (0usize, 0.0_f64, 0.0_f64);
+        let (mut n_ot, mut max_ot) = (0usize, 0.0_f64);
+        let (mut p_lo, mut p_hi) = (f64::INFINITY, f64::NEG_INFINITY);
+        let mut p_jump_max = 0.0_f64;
+        for j in 0..ny {
+            for i in 0..nx {
+                if class_now(j, i) != CellClass::Active {
+                    continue;
+                }
+                p_lo = p_lo.min(bg.p[(j, i)]);
+                p_hi = p_hi.max(bg.p[(j, i)]);
+                let d = div(j, i);
+                if changed(j, i) {
+                    // fringe → active (hole → active is counted here too)
+                    n_fa += 1;
+                    max_fa = max_fa.max(d.abs());
+                    sum_fa += d;
+                    let (mut ps, mut pc) = (0.0, 0usize);
+                    for (jj, ii) in [
+                        (j, i + 1),
+                        (j, i.wrapping_sub(1)),
+                        (j + 1, i),
+                        (j.wrapping_sub(1), i),
+                    ] {
+                        if jj < ny
+                            && ii < nx
+                            && class_now(jj, ii) == CellClass::Active
+                            && !changed(jj, ii)
+                        {
+                            ps += bg.p[(jj, ii)];
+                            pc += 1;
+                        }
+                    }
+                    if pc > 0 {
+                        p_jump_max = p_jump_max.max((bg.p[(j, i)] - ps / pc as f64).abs());
+                    }
+                } else {
+                    let near = [
+                        (j, i + 1),
+                        (j, i.wrapping_sub(1)),
+                        (j + 1, i),
+                        (j.wrapping_sub(1), i),
+                    ]
+                    .iter()
+                    .any(|&(jj, ii)| jj < ny && ii < nx && changed(jj, ii));
+                    if near {
+                        n_nb += 1;
+                        max_nb = max_nb.max(d.abs());
+                        sum_nb += d;
+                    } else {
+                        n_ot += 1;
+                        max_ot = max_ot.max(d.abs());
+                    }
+                }
+            }
+        }
+        println!(
+            "  SP-TRACE t = {:.5}: source (cell volumes/step) — turned active: {n_fa} cells, max |div| {max_fa:.3e}, sum {sum_fa:+.3e}; \
+             their active neighbours: {n_nb} cells, max {max_nb:.3e}, sum {sum_nb:+.3e}; other active: {n_ot} cells, max {max_ot:.3e}; \
+             stored p of turned-active cells vs neighbours: max jump {p_jump_max:.3e} (active p range {:.3e})",
+            self.patch.time() + dt,
+            p_hi - p_lo
+        );
+    }
+
+    /// The fringe flux balance (see `OversetParameters::fringe_flux_balance`).
+    fn balance_fringe(&self, field: &mut FlowField) {
+        if !self.params.fringe_flux_balance {
+            return;
+        }
+        let (_, _, dx, dy) = self.grid;
+        let mut scale = 0.0_f64;
+        for e in &self.overlap.fringe_u {
+            scale = scale.max((field.u[(e.j, e.i)] * dy).abs());
+        }
+        for e in &self.overlap.fringe_v {
+            scale = scale.max((field.v[(e.j, e.i)] * dx).abs());
+        }
+        let tol = self.params.fringe_balance_tolerance * scale.max(1e-300);
+        self.overlap.balance_fringe_fluxes(field, tol, 50);
+    }
+
    /// Subtract the mean of `p'` over the active background cells (the
    /// composite level pin; the fringe cells' Dirichlet values shift with
    /// it so the face corrections across active–fringe faces are unchanged).
@@ -341,6 +458,7 @@ impl OversetPisoSolver {
    fn exchange(&self, field: &mut OversetField) {
        self.overlap
            .stamp_fringe_faces(&mut field.background, &field.patch.u, &field.patch.v);
+        self.balance_fringe(&mut field.background);
        let p = self.overlap.fringe_cell_values(&field.patch.p);
        self.overlap.stamp_fringe_cells(&mut field.background.p, &p);
        let acc = self
@@ -384,6 +502,11 @@ impl OversetPisoSolver {
        //    path to the bit.
        let mut reclassified = 0usize;
        let mut fresh = 0usize;
+        // P3b locating trace (print-only, `RTX_OVERSET_TRACE_SP`): the class
+        // of every background cell before this step, kept only when a
+        // reclassification happens.
+        let trace_sp = std::env::var("RTX_OVERSET_TRACE_SP").is_ok();
+        let mut old_class: Option<Vec<CellClass>> = None;
        if self.pending.take().is_some() {
            let next = self
                .patch
@@ -397,8 +520,20 @@ impl OversetPisoSolver {
                        reclassified += 1;
                        if was == CellClass::Hole && now == CellClass::Active {
                            fresh += 1;
-                            // Fill from the neighbours that already carry a
-                            // pressure (the embedded solver's fresh-cell rule).
+                        }
+                        // A cell that turns ACTIVE takes its pressure from the
+                        // background's own neighbours that stay active, not the
+                        // patch's interpolated value it held as a fringe cell:
+                        // the two meshes' discrete pressures disagree at their
+                        // interface by their discretization error (4% of the
+                        // range on the static MMS, 5–10% next to the falsifier
+                        // plate), and the predictor's ∇p across the flipped
+                        // faces delivered that mismatch as the reclassification
+                        // impulse (P3b locating trace; the mass source of those
+                        // cells was ≤ 6e-3 cell volumes per step, not the cause).
+                        if now == CellClass::Active
+                            && (self.params.refill_turned_active || was == CellClass::Hole)
+                        {
                            let (mut sum, mut count) = (0.0, 0usize);
                            for (jj, ii) in [
                                (j, i + 1),
@@ -408,7 +543,8 @@ impl OversetPisoSolver {
                            ] {
                                if jj < ny
                                    && ii < nx
-                                    && self.overlap.class(jj, ii) != CellClass::Hole
+                                    && self.overlap.class(jj, ii) == CellClass::Active
+                                    && new.class(jj, ii) == CellClass::Active
                                {
                                    sum += field.background.p[(jj, ii)];
                                    count += 1;
@@ -421,11 +557,51 @@ impl OversetPisoSolver {
                    }
                }
            }
-            self.overlap = new;
+            if trace_sp && reclassified > 0 {
+                old_class = Some(
+                    (0..ny * nx)
+                        .map(|k| self.overlap.class(k / nx, k % nx))
+                        .collect(),
+                );
+            }
+            let old_map = std::mem::replace(&mut self.overlap, new);
            self.background
                .set_overlap(self.overlap.background_mask(), self.overlap.fringe_flags());
+            // Knock-out H1 (`RTX_OVERSET_H1`): faces that were the
+            // background's own (fluid in the old map) keep their values
+            // instead of taking the patch's interpolation when they turn
+            // prescribed.
+            let h1 = std::env::var("RTX_OVERSET_H1").is_ok();
+            let (u_keep, v_keep) = if h1 {
+                (
+                    Some(field.background.u.clone()),
+                    Some(field.background.v.clone()),
+                )
+            } else {
+                (None, None)
+            };
            self.overlap
                .stamp_fringe_faces(&mut field.background, &field.patch.u, &field.patch.v);
+            if let (Some(uk), Some(vk)) = (u_keep, v_keep) {
+                let old_active = |jj: usize, ii: usize| old_map.class(jj, ii) == CellClass::Active;
+                for e in &self.overlap.fringe_u {
+                    // u face (j, i) between cells (j, i-1) and (j, i)
+                    if (e.i > 0 && old_active(e.j, e.i - 1)) || (e.i < nx && old_active(e.j, e.i)) {
+                        field.background.u[(e.j, e.i)] = uk[(e.j, e.i)];
+                    }
+                }
+                for e in &self.overlap.fringe_v {
+                    if (e.j > 0 && old_active(e.j - 1, e.i)) || (e.j < ny && old_active(e.j, e.i)) {
+                        field.background.v[(e.j, e.i)] = vk[(e.j, e.i)];
+                    }
+                }
+            }
+            self.balance_fringe(&mut field.background);
+            // Knock-out H4 (`RTX_OVERSET_H4`): no temporal warm start on a
+            // reclassification step.
+            if std::env::var("RTX_OVERSET_H4").is_ok() {
+                self.acceptor_warm.clear();
+            }
            let p = self.overlap.fringe_cell_values(&field.patch.p);
            self.overlap.stamp_fringe_cells(&mut field.background.p, &p);
        }
@@ -435,6 +611,9 @@ impl OversetPisoSolver {

        // 2. Background predictor.
        let bg_start = self.background.begin_step(&mut field.background, dt)?;
+        if let Some(old) = &old_class {
+            self.trace_reclassification_source(field, old, dt);
+        }

        // 3. Correctors: alternating Schwarz on the acceptor p' vector `a`
        //    (patch solve with Dirichlet a → fringe p' → background solve →
@@ -431,6 +431,77 @@ impl OverlapMap {
            .collect()
    }

+    /// Flux balance at the fringe (Chesshire–Henshaw in spirit): make every
+    /// fringe cell divergence-free by adjusting only its PRESCRIBED faces
+    /// (never a face shared with an active cell), spreading each cell's
+    /// imbalance over them by face length, in Gauss–Seidel sweeps (a face
+    /// shared by two fringe cells is corrected by both) until the largest
+    /// fringe-cell imbalance is below `tol` (volume flux) or `max_sweeps`
+    /// is reached. Returns `(sweeps, worst imbalance)`. This is what removes
+    /// the reclassification impulse of A-P3: with the fringe ring a
+    /// staircase of the interpolated velocities' mass defect, every
+    /// row flip injected that defect in one step (§5.10).
+    pub fn balance_fringe_fluxes(
+        &self,
+        field: &mut FlowField,
+        tol: f64,
+        max_sweeps: usize,
+    ) -> (usize, f64) {
+        let (nx, ny, dx, dy) = (self.nx, self.ny, self.dx, self.dy);
+        let is_prescribed_u: std::collections::HashSet<(usize, usize)> =
+            self.fringe_u.iter().map(|e| (e.j, e.i)).collect();
+        let is_prescribed_v: std::collections::HashSet<(usize, usize)> =
+            self.fringe_v.iter().map(|e| (e.j, e.i)).collect();
+        let _ = (nx, ny);
+        let mut worst = f64::INFINITY;
+        let mut sweeps = 0usize;
+        while sweeps < max_sweeps && worst > tol {
+            sweeps += 1;
+            worst = 0.0;
+            for e in &self.fringe_cells {
+                let (j, i) = (e.j, e.i);
+                let div = (field.u[(j, i + 1)] - field.u[(j, i)]) * dy
+                    + (field.v[(j + 1, i)] - field.v[(j, i)]) * dx;
+                // Prescribed faces of this cell with their outward sign and length.
+                let mut faces: Vec<(bool, usize, usize, f64, f64)> = Vec::with_capacity(4);
+                if is_prescribed_u.contains(&(j, i + 1)) {
+                    faces.push((true, j, i + 1, 1.0, dy));
+                }
+                if is_prescribed_u.contains(&(j, i)) {
+                    faces.push((true, j, i, -1.0, dy));
+                }
+                if is_prescribed_v.contains(&(j + 1, i)) {
+                    faces.push((false, j + 1, i, 1.0, dx));
+                }
+                if is_prescribed_v.contains(&(j, i)) {
+                    faces.push((false, j, i, -1.0, dx));
+                }
+                let total_len: f64 = faces.iter().map(|f| f.4).sum();
+                if total_len == 0.0 {
+                    worst = worst.max(div.abs());
+                    continue;
+                }
+                for (is_u, jj, ii, sign, len) in faces {
+                    // outward flux change on this face = −div · len / total_len
+                    let dvel = -sign * div / total_len;
+                    if is_u {
+                        field.u[(jj, ii)] += dvel;
+                    } else {
+                        field.v[(jj, ii)] += dvel;
+                    }
+                    let _ = len;
+                }
+            }
+            for e in &self.fringe_cells {
+                let (j, i) = (e.j, e.i);
+                let div = (field.u[(j, i + 1)] - field.u[(j, i)]) * dy
+                    + (field.v[(j + 1, i)] - field.v[(j, i)]) * dx;
+                worst = worst.max(div.abs());
+            }
+        }
+        (sweeps, worst)
+    }
+
    /// A cell-centred background scalar (e.g. `p'`) at every acceptor.
    pub fn acceptor_scalar(&self, m: &nalgebra::DMatrix<f64>) -> Vec<f64> {
        self.acceptors.iter().map(|a| a.p.value(m)).collect()