rtx-cfd: OversetPisoSolver::solver_metric_force — the box force in the predictor's own flux form (upwind convective, diffusive, cell pressure on the momentum control volumes' faces, minus the unsteady term); box-independent to rounding in the active region (five boxes spread 2e-12 at step 5 vs the CV formula's ±0.5 %); overset_cfd1 prints wall − box = ring Σr + rest and saves the settled fields under RTX_OVERSET_CFD1_SAVE

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

crates/specialized/rtx-cfd/src/solvers/incompressible/overset/residual.rs

@@ -298,3 +298,91 @@ impl OversetPisoSolver {
         out
     }
 }
+
+impl OversetPisoSolver {
+    /// The force on everything inside the box `(i0, i1, j0, j1)` (cell
+    /// index bounds, as `EmbeddedMask::control_volume_force`) in the
+    /// SOLVER'S OWN flux form: the predictor's upwind convective flux,
+    /// its diffusive flux and the cell pressure, on the momentum control
+    /// volumes' faces that make up the box boundary (u volumes `i0 + 1
+    /// ..= i1` × `j0 .. j1`, v volumes `i0 .. i1` × `j0 + 1 ..= j1`),
+    /// minus the unsteady term over the box's evaluable volumes. On the
+    /// solved faces the residual is rounding, so this is box-INDEPENDENT
+    /// to rounding as long as the box stays in the active region — the
+    /// gate — whereas the control-volume formula moves by ±0.5 % between
+    /// boxes. Returns `(fx, fy)`, positive = drag / lift on the body.
+    pub fn solver_metric_force(
+        &self,
+        field: &OversetField,
+        dt: f64,
+        (i0, i1, j0, j1): (usize, usize, usize, usize),
+    ) -> (f64, f64) {
+        let (nx, ny, dx, dy) = self.grid;
+        assert!(
+            i0 >= 1 && i1 + 1 < nx && j0 >= 1 && j1 + 1 < ny,
+            "box must be interior"
+        );
+        let rho = self.background.config().density;
+        let mu = self.background.config().viscosity;
+        let bg = &field.background;
+        let (u, v, p) = (&bg.u_old, &bg.v_old, &bg.p);
+        let upwind = |f: f64, up: f64, down: f64| if f >= 0.0 { up } else { down };
+        // Outward x-momentum flux through the u-volume face at cell i's
+        // centre (n = +x), per unit length.
+        let phi_u_x = |j: usize, i: usize| {
+            let ue = 0.5 * (u[(j, i)] + u[(j, i + 1)]);
+            rho * ue * upwind(ue, u[(j, i)], u[(j, i + 1)]) - mu * (u[(j, i + 1)] - u[(j, i)]) / dx
+                + p[(j, i)]
+        };
+        // Through the u-volume face at v-face row j (n = +y), for u face i.
+        let phi_u_y = |j: usize, i: usize| {
+            let vn = 0.5 * (v[(j, i - 1)] + v[(j, i)]);
+            rho * vn * upwind(vn, u[(j - 1, i)], u[(j, i)]) - mu * (u[(j, i)] - u[(j - 1, i)]) / dy
+        };
+        // y-momentum: through the v-volume face at cell j's centre (n = +y).
+        let phi_v_y = |j: usize, i: usize| {
+            let vn = 0.5 * (v[(j, i)] + v[(j + 1, i)]);
+            rho * vn * upwind(vn, v[(j, i)], v[(j + 1, i)]) - mu * (v[(j + 1, i)] - v[(j, i)]) / dy
+                + p[(j, i)]
+        };
+        // Through the v-volume face at u-face column i (n = +x), for v face j.
+        let phi_v_x = |j: usize, i: usize| {
+            let ue = 0.5 * (u[(j - 1, i)] + u[(j, i)]);
+            rho * ue * upwind(ue, v[(j, i - 1)], v[(j, i)]) - mu * (v[(j, i)] - v[(j, i - 1)]) / dx
+        };
+        let vol = dx * dy;
+        let (mut out_x, mut out_y) = (0.0, 0.0);
+        let (mut dt_x, mut dt_y) = (0.0, 0.0);
+        for j in j0..j1 {
+            out_x += (phi_u_x(j, i1) - phi_u_x(j, i0)) * dy;
+        }
+        for i in i0 + 1..=i1 {
+            out_x += (phi_u_y(j1, i) - phi_u_y(j0, i)) * dx;
+            for j in j0..j1 {
+                let d = bg.u[(j, i)] - bg.u_old[(j, i)];
+                if d.is_finite()
+                    && (self.overlap.class(j, i) != CellClass::Hole
+                        || self.overlap.class(j, i - 1) != CellClass::Hole)
+                {
+                    dt_x += rho * d / dt * vol;
+                }
+            }
+        }
+        for i in i0..i1 {
+            out_y += (phi_v_y(j1, i) - phi_v_y(j0, i)) * dx;
+        }
+        for j in j0 + 1..=j1 {
+            out_y += (phi_v_x(j, i1) - phi_v_x(j, i0)) * dy;
+            for i in i0..i1 {
+                let d = bg.v[(j, i)] - bg.v_old[(j, i)];
+                if d.is_finite()
+                    && (self.overlap.class(j, i) != CellClass::Hole
+                        || self.overlap.class(j - 1, i) != CellClass::Hole)
+                {
+                    dt_y += rho * d / dt * vol;
+                }
+            }
+        }
+        (-out_x - dt_x, -out_y - dt_y)
+    }
+}