rtx-cfd: CurvilinearPisoSolver::momentum_balance — the patch's own momentum balance on its solved cells in the scheme's fluxes (outward ρFu_f with the predictor's face value, Laplacian-form μ∇u·S on the solved/acceptor interface and the wall, the least-squares pressure volume sum vs the face-pressure integrals); flux_force, wall_force, pressure_defect δP; overset_cfd1 prints it and the acceptor band's mismatch, saves the patch flux, and RTX_OVERSET_CFD1_LOAD=dir runs the diagnostics offline on saved fields

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/balance.rs

@@ -0,0 +1,149 @@
+//! P4 option B, patch side (`docs/overset_metal_campaign.md` §5.11): the
+//! patch's own momentum balance on its SOLVED cells (everything but the
+//! acceptor row), in the scheme's own fluxes.
+//!
+//! The marched equation on a solved cell is the face-flux form of
+//! `predictor.rs` (outward `ρ F_f u_f` with the scheme's face value,
+//! outward `μ L_f(u)` in Laplacian form) plus the projection's pressure
+//! force, which is the least-squares CELL gradient `−A_c ∇p_c` — not a
+//! face pressure, so the pressure term does not telescope: summed over
+//! the solved cells it need not equal the boundary integral of any face
+//! pressure, and that difference `δP` is the patch's momentum
+//! non-conservation. Everything here is a plain sum of the scheme's
+//! terms; the steady balance `−conv_acc + visc_acc + visc_wall − p_ls = 0`
+//! is the gate (the unsteady term is not stored by `PatchField`; at the
+//! settled state it is small and reported as omitted).
+
+use super::{CurvilinearPisoSolver, PatchConvection, PatchField, SideBc};
+use crate::mesh::{PatchMesh, PatchSide};
+
+/// The pieces of the patch's momentum balance, x / y, N/m.
+#[derive(Debug, Clone, Copy, Default)]
+pub struct PatchBalance {
+    /// Outward `ρ F_f u_f` through the solved/acceptor interface faces.
+    pub conv_acc: [f64; 2],
+    /// Outward `μ ∇u · S` (Laplacian form) through the interface faces.
+    pub visc_acc: [f64; 2],
+    /// `Σ p_f S_out` on the interface faces, `p_f` linear in the two cells.
+    pub p_face_acc: [f64; 2],
+    /// Outward `μ ∇u · S` through the wall faces (Dirichlet wall value).
+    pub visc_wall: [f64; 2],
+    /// `Σ p_f S_out` on the wall faces, `p_f` extrapolated as
+    /// `surface_force` does (cell value + least-squares gradient).
+    pub p_face_wall: [f64; 2],
+    /// `Σ_solved A_c ∇p_c` — the pressure force the scheme applied.
+    pub p_ls: [f64; 2],
+    /// Solved cells.
+    pub cells: usize,
+    /// Interface faces.
+    pub acc_faces: usize,
+    /// Wall faces.
+    pub wall_faces: usize,
+}
+
+impl PatchBalance {
+    /// Steady balance residual of the marched equation on the solved
+    /// region: `−conv_acc + visc_acc + visc_wall − p_ls` (the gate).
+    pub fn balance(&self) -> [f64; 2] {
+        [0, 1].map(|k| -self.conv_acc[k] + self.visc_acc[k] + self.visc_wall[k] - self.p_ls[k])
+    }
+    /// The body force read through the interface in flux form:
+    /// `∮ (σ·n − ρ u u·n)` with `σ` in the scheme's Laplacian form.
+    pub fn flux_force(&self) -> [f64; 2] {
+        [0, 1].map(|k| -self.p_face_acc[k] + self.visc_acc[k] - self.conv_acc[k])
+    }
+    /// The wall force in the scheme's own wall fluxes (Laplacian form).
+    pub fn wall_force(&self) -> [f64; 2] {
+        [0, 1].map(|k| self.p_face_wall[k] - self.visc_wall[k])
+    }
+    /// The pressure non-conservation `p_ls − p_face_acc − p_face_wall`
+    /// (= `flux_force − wall_force` when the balance holds).
+    pub fn pressure_defect(&self) -> [f64; 2] {
+        [0, 1].map(|k| self.p_ls[k] - self.p_face_acc[k] - self.p_face_wall[k])
+    }
+}
+
+impl CurvilinearPisoSolver {
+    /// The momentum balance of the solved cells at time `t` (see the
+    /// module doc). Stationary mesh only.
+    pub fn momentum_balance(&self, field: &PatchField, t: f64) -> PatchBalance {
+        let mesh: &PatchMesh = &self.mesh;
+        let rho = self.config.density;
+        let mu = self.config.viscosity;
+        let ops = &self.ops;
+        let bvel = |side: PatchSide, xy: [f64; 2]| -> Option<(f64, f64)> {
+            match self.params.boundaries.get(side) {
+                SideBc::Velocity => Some(self.boundary_velocity(side, xy[0], xy[1], t)),
+                SideBc::Outlet => None,
+            }
+        };
+        let un = ops.node_values(mesh, &field.u, &|s, xy| bvel(s, xy).map(|v| v.0));
+        let vn = ops.node_values(mesh, &field.v, &|s, xy| bvel(s, xy).map(|v| v.1));
+        let mut b = PatchBalance::default();
+        for c in 0..mesh.cell_count() {
+            if self.is_acceptor(c) {
+                continue;
+            }
+            b.cells += 1;
+            let g = self.pressure_gradient(&field.p, c);
+            let a = mesh.area(c);
+            b.p_ls[0] += a * g[0];
+            b.p_ls[1] += a * g[1];
+            for (f, sign) in mesh.cell_faces(c) {
+                let face = &mesh.faces()[f];
+                let s_out = [sign * face.s[0], sign * face.s[1]];
+                match (face.owner, face.neigh) {
+                    (Some(p), Some(q)) => {
+                        let other = if p == c { q } else { p };
+                        if !self.is_acceptor(other) {
+                            continue;
+                        }
+                        b.acc_faces += 1;
+                        let out = sign * field.flux[f];
+                        let (uf, vf) = match self.params.convection {
+                            PatchConvection::Upwind => {
+                                let up = if out >= 0.0 { c } else { other };
+                                (field.u[up], field.v[up])
+                            }
+                            PatchConvection::TvdVanAlbada => {
+                                let (up, dn) = if out >= 0.0 { (c, other) } else { (other, c) };
+                                let (du, dv) = self.tvd_correction(field, f, up, dn);
+                                (field.u[up] + du, field.v[up] + dv)
+                            }
+                            PatchConvection::None => (0.0, 0.0),
+                        };
+                        b.conv_acc[0] += rho * out * uf;
+                        b.conv_acc[1] += rho * out * vf;
+                        b.visc_acc[0] +=
+                            mu * sign * ops.face_gradient_flux(mesh, f, &field.u, &un, None);
+                        b.visc_acc[1] +=
+                            mu * sign * ops.face_gradient_flux(mesh, f, &field.v, &vn, None);
+                        let pf = face.w * field.p[p] + (1.0 - face.w) * field.p[q];
+                        b.p_face_acc[0] += pf * s_out[0];
+                        b.p_face_acc[1] += pf * s_out[1];
+                    }
+                    _ => {
+                        let side = mesh.side(f).expect("boundary face has a side");
+                        if side != PatchSide::Inner {
+                            continue;
+                        }
+                        b.wall_faces += 1;
+                        let bv = bvel(side, face.centre);
+                        b.visc_wall[0] += mu
+                            * sign
+                            * ops.face_gradient_flux(mesh, f, &field.u, &un, bv.map(|v| v.0));
+                        b.visc_wall[1] += mu
+                            * sign
+                            * ops.face_gradient_flux(mesh, f, &field.v, &vn, bv.map(|v| v.1));
+                        let xc = mesh.centre(c);
+                        let dxf = [face.centre[0] - xc[0], face.centre[1] - xc[1]];
+                        let pf = field.p[c] + g[0] * dxf[0] + g[1] * dxf[1];
+                        b.p_face_wall[0] += pf * s_out[0];
+                        b.p_face_wall[1] += pf * s_out[1];
+                    }
+                }
+            }
+        }
+        b
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/mod.rs

@@ -30,11 +30,13 @@
 //! static one (no mesh-velocity term). A stationary mesh through this path
 //! is bit-identical to the static path.
 
+mod balance;
 mod motion;
 mod operators;
 mod predictor;
 mod projection;
 
+pub use balance::PatchBalance;
 pub use motion::StepGeometry;
 pub use operators::Operators;
 

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/predictor.rs

@@ -164,7 +164,13 @@ impl CurvilinearPisoSolver {
     /// between the upwind cell `up` and the downwind cell `dn`: `w_up psi(r)
     /// (phi_dn − phi_up)` for `u` and `v`, zero when the far-upwind cell
     /// (across `up`'s opposite face) lies outside the patch.
-    fn tvd_correction(&self, field: &PatchField, f: usize, up: usize, dn: usize) -> (f64, f64) {
+    pub(super) fn tvd_correction(
+        &self,
+        field: &PatchField,
+        f: usize,
+        up: usize,
+        dn: usize,
+    ) -> (f64, f64) {
         let mesh = &self.mesh;
         let faces = mesh.cell_faces(up);
         let Some(pos) = faces.iter().position(|&(g, _)| g == f) else {

crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs

@@ -48,8 +48,8 @@ pub use boundary_conditions::{
 };
 pub use curvilinear::{
     CurvilinearParameters, CurvilinearPisoSolver, CurvilinearResult, CurvilinearSolverState,
-    NormalDiffusion, Operators, PatchBoundaries, PatchConvection, PatchField, PatchLoad, SideBc,
+    NormalDiffusion, Operators, PatchBalance, PatchBoundaries, PatchConvection, PatchField,
-    StepGeometry,
+    PatchLoad, SideBc, StepGeometry,
 };
 pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
 pub use embedded_body::{

crates/specialized/rtx-cfd/tests/overset_cfd1.rs

@@ -177,6 +177,40 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
             )
     };
 
+    // `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 = format!(
+                "ny{ny}_{}",
+                if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
+                    "tvd"
+                } else {
+                    "upwind"
+                }
+            );
+            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;
@@ -185,6 +219,7 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
     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;
@@ -221,7 +256,8 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
             if steps % 50 == 0 {
                 let (fx, _) = cv_force(&field, &solver);
                 history.push(fx);
-            let load = solver
+                let load =
+                    solver
                         .patch()
                         .surface_force(&field.patch, PatchSide::Inner, solver.time());
                 let umax = field
@@ -255,6 +291,7 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
             }
             assert!(steps < 2_000_000, "CFD1 at ny = {ny} did not settle");
         }
+    }
     let seconds = start.elapsed().as_secs_f64();
     let load = solver
         .patch()
@@ -373,6 +410,29 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
         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).
@@ -392,6 +452,7 @@ async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
             ("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");