embedded3: cut predictor convection carries ρ (host + e3_cut.cu; density-scaling pin); operator load route includes the wall exchange (exchange.rs); reconstructed_parts, probe aperture floor knob; dfg_split diagnostic test

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu

@@ -176,7 +176,8 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         if (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d];
         else if (sides[d][0] == SIDE_VELOCITY) diff -= mu * g_minus * a_d * (u0 - beyond_m) / (0.5 * h[d]);
     }
-    conv -= mass_out * u0;
+    /* The mass fluxes above are volume fluxes: the momentum flux carries rho. */
+    conv = rho * (conv - mass_out * u0);
     int cpi = cut_cell(g, cp[0], cp[1], cp[2]);
     int cmi = cut_cell(g, cm[0], cm[1], cm[2]);
     double p_plus = cpi >= 0 ? f.p[cpi] : 0.0;

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs

@@ -207,6 +207,8 @@ impl Mask {
             step_apertures: None,
             step_open: None,
             merge_master: Vec::new(),
+            scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
+            density: 1.0,
         };
         mask.compute_merging(None);
         Ok(mask)
@@ -564,7 +566,8 @@ impl Mask {
     /// without a cut geometry.
     pub fn cut_wall_force(&self, body: &Body, f: &Field, mu: f64, t: f64) -> Option<[f64; 3]> {
         let (p, s) = self.cut_wall_force_parts(body, f, mu, t)?;
-        Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]])
+        let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, None)?;
+        Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]])
     }
 
     /// The reconstructed wall route (S2-1 remedy): on every wall polygon
@@ -585,12 +588,26 @@ impl Mask {
         t: f64,
         planes: Option<(usize, usize)>,
     ) -> Option<[f64; 3]> {
+        let (p, s) = self.cut_wall_force_reconstructed_parts(body, f, mu, t, planes)?;
+        Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]])
+    }
+
+    /// The reconstructed route split into its pressure and shear parts.
+    pub fn cut_wall_force_reconstructed_parts(
+        &self,
+        body: &Body,
+        f: &Field,
+        mu: f64,
+        t: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<([f64; 3], [f64; 3])> {
         let cut = self.cut.as_ref()?;
         let g = self.grid;
         let (k0, k1) = planes.unwrap_or((0, g.nz));
         let h = g.dx.min(g.dy).min(g.dz);
         let (d1, d2) = (h, 2.0 * h);
         let mut force = [0.0; 3];
+        let mut shear = [0.0; 3];
         for (idx, w) in cut.wall.iter().enumerate() {
             let area = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
             if area == 0.0 || !self.cell_fluid[idx] {
@@ -645,10 +662,10 @@ impl Mask {
                 let dn = wall_gradient(t1[c] - ts[c], t2[c] - ts[c]);
                 // Traction on the body = −(fluid stress on the fluid side):
                 // the shear the fluid exerts on the wall along +t.
-                force[c] += mu * dn * area;
+                shear[c] += mu * dn * area;
             }
         }
-        Some(force)
+        Some((force, shear))
     }
 
     /// The cut-cell load route restricted to the cells (and faces) of the
@@ -663,8 +680,13 @@ impl Mask {
         (k0, k1): (usize, usize),
     ) -> Option<[f64; 3]> {
         let (p, s) = self.cut_wall_force_parts_in(body, f, mu, t, Some((k0, k1)))?;
+        let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, Some((k0, k1)))?;
         let lz = (k1 - k0) as f64 * self.grid.dz;
-        Some([(p[0] + s[0]) / lz, (p[1] + s[1]) / lz, (p[2] + s[2]) / lz])
+        Some([
+            (p[0] + s[0] + x[0]) / lz,
+            (p[1] + s[1] + x[1]) / lz,
+            (p[2] + s[2] + x[2]) / lz,
+        ])
     }
 
     /// The cut-cell load route split into its pressure and shear parts.

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs

@@ -0,0 +1,143 @@
+//! The wall-exchange part of the operator load route (S2-1 remedy, second
+//! form). A fluid face control volume next to a prescribed (ghost or
+//! solid) face still exchanges momentum with it: the diffusive flux
+//! `μ g A (U_b − u_f)/h` through the half-aperture face at the cell
+//! centre and the convective flux `ρ m (u_face − u_f)` on the same face.
+//! Both are forces the wall exerts on the fluid that the closure polygon's
+//! shear `μ A_w (u_f − U_b)/d_f` does not carry, so the operator route
+//! read short of the box route by exactly this exchange (a conservation
+//! gap that did not shrink with h). Summed here with the predictor's own
+//! flux formulas, the operator route closes the discrete momentum balance.
+use super::body::Body;
+use super::field::Field;
+use super::wall::{FaceKind, Mask};
+use crate::solvers::incompressible::ConvectionScheme;
+
+impl Mask {
+    /// The momentum the fluid's face control volumes exchange with the
+    /// prescribed faces beside them, as a force on the body (the negative
+    /// of the force on the fluid), over the z planes `planes` (all when
+    /// `None`). `None` without a cut geometry.
+    pub fn cut_wall_exchange_force(
+        &self,
+        body: &Body,
+        f: &Field,
+        mu: f64,
+        rho: f64,
+        t: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<[f64; 3]> {
+        let _ = body;
+        let _ = t;
+        self.cut.as_ref()?;
+        let g = self.grid;
+        let lat = self.lattice();
+        let h = [g.dx, g.dy, g.dz];
+        let area = [g.dy * g.dz, g.dx * g.dz, g.dx * g.dy];
+        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
+        let (k0, k1) = planes.unwrap_or((0, nz));
+        let w_range = if self.periodic_z { 0..nz } else { 1..nz };
+        let vals: [&[f64]; 3] = [&f.u, &f.v, &f.w];
+        let scheme = self.scheme;
+        let kind = |cc: usize, idx: usize| match cc {
+            0 => self.u_kind[idx],
+            1 => self.v_kind[idx],
+            _ => self.w_kind[idx],
+        };
+        let val = |cc: usize, q: [i64; 3]| lat.face(cc, q).map(|i| vals[cc][i]);
+        let ap = |cc: usize, q: [i64; 3]| self.aperture(cc, q);
+        let e = |d: usize| {
+            let mut v = [0i64; 3];
+            v[d] = 1;
+            v
+        };
+        let add =
+            |a: [i64; 3], b: [i64; 3], s: i64| [a[0] + s * b[0], a[1] + s * b[1], a[2] + s * b[2]];
+        let upwind = |m: f64, up: f64, dn: f64| if m >= 0.0 { up } else { dn };
+        let mut force = [0.0; 3];
+        for c in 0..3 {
+            let (ir, jr, kr) = match c {
+                0 => (1..nx, 0..ny, k0..k1),
+                1 => (0..nx, 1..ny, k0..k1),
+                _ => (0..nx, 0..ny, w_range.start.max(k0)..w_range.end.min(k1)),
+            };
+            let ec = e(c);
+            for k in kr {
+                for j in jr.clone() {
+                    for i in ir.clone() {
+                        let p = [i as i64, j as i64, k as i64];
+                        let idx = lat.face(c, p).expect("face");
+                        if kind(c, idx) != FaceKind::Fluid {
+                            continue;
+                        }
+                        let cv = self.cv_geometry(c, p);
+                        let u0 = vals[c][idx];
+                        let cell_minus = add(p, ec, -1);
+                        let cell_plus = p;
+                        for d in 0..3 {
+                            let ed = e(d);
+                            let a_d = area[d];
+                            let up1 = val(c, add(p, ed, 1));
+                            let up2 = val(c, add(p, ed, 2));
+                            let dn1 = val(c, add(p, ed, -1));
+                            let dn2 = val(c, add(p, ed, -2));
+                            let (m_plus, m_minus) = if d == c {
+                                let f_up =
+                                    ap(c, add(p, ec, 1)).unwrap_or(cv.alpha) * up1.unwrap_or(u0);
+                                let f_dn =
+                                    ap(c, add(p, ec, -1)).unwrap_or(cv.alpha) * dn1.unwrap_or(u0);
+                                let f0 = cv.alpha * u0;
+                                (0.5 * (f0 + f_up) * a_d, 0.5 * (f_dn + f0) * a_d)
+                            } else {
+                                let flux = |q: [i64; 3]| {
+                                    ap(d, q).unwrap_or(1.0) * val(d, q).unwrap_or(0.0)
+                                };
+                                (
+                                    0.5 * (flux(add(cell_minus, ed, 1))
+                                        + flux(add(cell_plus, ed, 1)))
+                                        * a_d,
+                                    0.5 * (flux(cell_minus) + flux(cell_plus)) * a_d,
+                                )
+                            };
+                            // Plus side: a prescribed neighbour face.
+                            if let Some(fp) = lat.face(c, add(p, ed, 1)) {
+                                if kind(c, fp) != FaceKind::Fluid {
+                                    let un = vals[c][fp];
+                                    let delta = if scheme == ConvectionScheme::Upwind {
+                                        0.0
+                                    } else if m_plus >= 0.0 {
+                                        scheme.face_correction(dn1, u0, un)
+                                    } else {
+                                        scheme.face_correction(up2, un, u0)
+                                    };
+                                    let u_face = upwind(m_plus, u0, un) + delta;
+                                    let on_fluid = -rho * m_plus * (u_face - u0)
+                                        + mu * cv.ap[d][1] * a_d * (un - u0) / h[d];
+                                    force[c] -= on_fluid;
+                                }
+                            }
+                            // Minus side.
+                            if let Some(fm) = lat.face(c, add(p, ed, -1)) {
+                                if kind(c, fm) != FaceKind::Fluid {
+                                    let ud = vals[c][fm];
+                                    let delta = if scheme == ConvectionScheme::Upwind {
+                                        0.0
+                                    } else if m_minus >= 0.0 {
+                                        scheme.face_correction(dn2, ud, u0)
+                                    } else {
+                                        scheme.face_correction(up1, u0, ud)
+                                    };
+                                    let u_face = upwind(m_minus, ud, u0) + delta;
+                                    let on_fluid = rho * m_minus * (u_face - u0)
+                                        + mu * cv.ap[d][0] * a_d * (ud - u0) / h[d];
+                                    force[c] -= on_fluid;
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        Some(force)
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs

@@ -13,6 +13,16 @@ use super::body::Body;
 use super::field::Field;
 use super::wall::{FaceKind, Mask, linear_fit, stencil_nodes, z_planes};
 
+/// Minimum face aperture for a velocity node to enter a probe's
+/// interpolation (`RTX_E3_PROBE_MIN_APERTURE`, default 0.5; 0 keeps every
+/// fluid face, the reading before S2-1's fix).
+fn probe_min_aperture() -> f64 {
+    std::env::var("RTX_E3_PROBE_MIN_APERTURE")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(0.5)
+}
+
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub struct SurfaceForce {
     pub f: [f64; 3],
@@ -154,10 +164,16 @@ impl Mask {
             let foot_c = (foot.0, foot.1, if zq == z { foot.2 } else { zq });
             let nodes = stencil_nodes((x, y, zq), c, g, self.periodic_z(), |_| None);
             let values: &[f64] = [&f.u, &f.v, &f.w][c];
+            // A cut face's velocity lives on the open part of the face, not
+            // at its centre: a node whose centre lies in the body would put
+            // the wall profile's value h/2 too deep and bias the wall
+            // gradient at O(1). Nodes below the minimum aperture are dropped
+            // and the point's own wall intercept takes their place.
+            let min_aperture = probe_min_aperture();
             let fluid = |idx: usize| match c {
-                0 => self.u_kind(idx) == FaceKind::Fluid,
+                0 => self.u_kind(idx) == FaceKind::Fluid && self.a_u(idx) >= min_aperture,
-                1 => self.v_kind(idx) == FaceKind::Fluid,
+                1 => self.v_kind(idx) == FaceKind::Fluid && self.a_v(idx) >= min_aperture,
-                _ => self.w_kind(idx) == FaceKind::Fluid,
+                _ => self.w_kind(idx) == FaceKind::Fluid && self.a_w(idx) >= min_aperture,
             };
             if nodes.iter().all(|n| fluid(n.idx)) {
                 out[c] = nodes.iter().map(|n| n.weight * values[n.idx]).sum();

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

@@ -9,6 +9,7 @@
 pub mod body;
 pub mod cut;
 pub mod cutwall;
+pub mod exchange;
 pub mod export_vtk;
 pub mod field;
 pub mod grid;

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs

@@ -201,7 +201,9 @@ impl Solver {
         // control volume (zero for a body at rest, the swept rate
         // otherwise) multiplies the face's own value, so a uniform field
         // stays uniform on any wall motion.
-        conv -= mass_out * u0;
+        // The mass fluxes above are volume fluxes: the momentum flux carries ρ
+        // (inertia, diffusion and the pressure are dynamic).
+        let conv = rho * (conv - mass_out * u0);
         let p_plus = lat.cell(cell_plus).map_or(0.0, |ci| field.p[ci]);
         let p_minus = lat.cell(cell_minus).map_or(0.0, |ci| field.p[ci]);
         let pressure = -(p_plus - p_minus) * cv.alpha * area[c];

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

@@ -213,6 +213,11 @@ impl Solver {
                 Mask::build_cut_from(body, g, t, self.params.boundaries, prev)
             }
         }
+        .map(|mut m| {
+            m.scheme = self.params.convection_scheme;
+            m.density = self.fluid.density;
+            m
+        })
         .expect("embedded mask")
     }
 

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs

@@ -92,6 +92,11 @@ pub struct Mask {
     /// (`usize::MAX` = its own row) — a small cell shares its pressure
     /// unknown with its largest active face neighbour in the projection.
     pub(super) merge_master: Vec<usize>,
+    /// The predictor's convection scheme (the exchange route replicates
+    /// its limited fluxes on the faces next to prescribed ones).
+    pub(super) scheme: crate::solvers::incompressible::ConvectionScheme,
+    /// The fluid's density (the exchange route's convective flux).
+    pub(super) density: f64,
 }
 
 /// The z lattice position of a query: the lower plane index, the upper
@@ -502,6 +507,8 @@ impl Mask {
             step_apertures: None,
             step_open: None,
             merge_master: Vec::new(),
+            scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
+            density: 1.0,
         })
     }
 

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

@@ -0,0 +1,101 @@
+//! Density-scaling pin for the cut-cell predictor (host): the same flow at
+//! `ρ` and `1000 ρ` with `μ` scaled alike is the same velocity field and a
+//! pressure scaled by 1000 — every term of the momentum equation carries
+//! `ρ` (the convection term used to be a bare volume flux times velocity,
+//! which starved every ρ = 1000 cut-cell run of convection).
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
+};
+
+fn run(rho: f64, moving: bool) -> Field {
+    let n = 16;
+    let h = 1.0 / n as f64;
+    let g = Grid::cubic(2 * n, n, n, h);
+    let nu = 1e-2;
+    let mut solver = Solver::new(
+        Fluid {
+            density: rho,
+            viscosity: rho * nu,
+            reference_velocity: 1.0,
+            reference_length: 0.3,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-11,
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+            wall_scheme: WallScheme::CutCell,
+            boundaries: Boundaries {
+                x1: Side::PressureOutlet,
+                ..Boundaries::default()
+            },
+            max_surface_speed: if moving { Some(0.5) } else { None },
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(|x, _, _, _| {
+        if x <= 0.0 {
+            (1.0, 0.0, 0.0)
+        } else {
+            (0.0, 0.0, 0.0)
+        }
+    });
+    let xc = move |t: f64| 0.7 + if moving { 0.1 * (3.0 * t).sin() } else { 0.0 };
+    let body = Body::from_sdf(move |x, y, z, t| {
+        ((x - xc(t)).powi(2) + (y - 0.5_f64).powi(2) + (z - 0.5_f64).powi(2)).sqrt() - 0.15
+    })
+    .with_surface_velocity(move |_, _, _, t| {
+        (if moving { 0.3 * (3.0 * t).cos() } else { 0.0 }, 0.0, 0.0)
+    });
+    if moving {
+        solver.set_moving_body(body);
+    } else {
+        solver.set_body(body);
+    }
+    let mut field = Field::new(g);
+    for k in 0..n {
+        for j in 0..n {
+            for i in 0..=2 * n {
+                field.u[g.uface(k, j, i)] = 1.0;
+            }
+        }
+    }
+    solver.initialize(&mut field);
+    let dt = 0.2 * h;
+    for _ in 0..40 {
+        solver.advance(&mut field, dt);
+    }
+    field
+}
+
+fn compare(moving: bool) {
+    let a = run(1.0, moving);
+    let b = run(1000.0, moving);
+    let max = |x: &[f64], y: &[f64], s: f64| {
+        x.iter()
+            .zip(y)
+            .map(|(p, q)| (p - q / s).abs())
+            .fold(0.0, f64::max)
+    };
+    let du = max(&a.u, &b.u, 1.0)
+        .max(max(&a.v, &b.v, 1.0))
+        .max(max(&a.w, &b.w, 1.0));
+    let dp = max(&a.p, &b.p, 1000.0);
+    let pscale = a.p.iter().fold(0.0f64, |m, p| m.max(p.abs()));
+    println!("  moving {moving}: max |Δu| {du:.3e}, max |Δp/1000| {dp:.3e} (p scale {pscale:.3e})");
+    assert!(du < 1e-9, "velocity is not density-invariant: {du:.3e}");
+    assert!(
+        dp < 1e-9 * pscale.max(1.0),
+        "pressure does not scale with density: {dp:.3e}"
+    );
+}
+
+#[test]
+fn cut_cell_flow_is_density_invariant_at_rest() {
+    compare(false);
+}
+
+#[test]
+fn cut_cell_flow_is_density_invariant_moving() {
+    compare(true);
+}

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

@@ -201,6 +201,23 @@ fn dfg_3d_2z_on_the_device() {
         println!("  instant written to {}", path.display());
     }
     let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples");
+    {
+        // S2-1 diagnosis: each wall route split into its pressure and shear parts.
+        let body = solver.body().expect("body");
+        let (po, so) = mask
+            .cut_wall_force_parts(body, &field, RHO * NU, solver.time())
+            .expect("parts");
+        let (pr, sr) = mask
+            .cut_wall_force_reconstructed_parts(body, &field, RHO * NU, solver.time(), None)
+            .expect("parts");
+        println!(
+            "  SPLIT ny {ny}: operator c_D pressure {:.4} + shear {:.4}; reconstructed pressure {:.4} + shear {:.4}",
+            coef * po[0],
+            coef * so[0],
+            coef * pr[0],
+            coef * sr[0]
+        );
+    }
     println!(
         "  FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; reconstructed {:.4}, {:.2e} from CV)  c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {:.5})  Δp {dp:.4}  — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s",
         ((cd - cd_cv) / cd).abs(),

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

@@ -0,0 +1,130 @@
+//! S2-1 diagnosis (host): DFG 3D-2Z at a coarse rung with every load
+//! route split into pressure and shear parts — which part of the wall
+//! routes departs from the box route. `RTX_E3_DFG_NY` (default 31).
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
+};
+
+const H: f64 = 0.41;
+const L: f64 = 2.5;
+const D: f64 = 0.1;
+const CX: f64 = 0.5;
+const CY: f64 = 0.2;
+const U_M: f64 = 0.45;
+const U_BAR: f64 = 4.0 / 9.0 * U_M;
+const RHO: f64 = 1.0;
+const NU: f64 = 1e-3;
+
+fn inflow(y: f64, z: f64) -> f64 {
+    16.0 * U_M * y * z * (H - y) * (H - z) / (H * H * H * H)
+}
+
+#[test]
+#[ignore = "host DFG at ny 31 with the routes split (about half an hour)"]
+fn dfg_routes_split_on_the_host() {
+    let ny: usize = std::env::var("RTX_E3_DFG_NY")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(31);
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let nz = ny;
+    let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU));
+    let mut solver = Solver::new(
+        Fluid {
+            density: RHO,
+            viscosity: RHO * NU,
+            reference_velocity: U_BAR,
+            reference_length: D,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+            wall_scheme: WallScheme::CutCell,
+            boundaries: Boundaries {
+                x1: Side::PressureOutlet,
+                ..Boundaries::default()
+            },
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(|x, y, z, _t| {
+        if x <= 0.0 {
+            (inflow(y, z), 0.0, 0.0)
+        } else {
+            (0.0, 0.0, 0.0)
+        }
+    });
+    solver.set_body(Body::extruded(
+        rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D),
+        H,
+    ));
+    let g = Grid::cubic(nx, ny, nz, h);
+    let mut field = Field::new(g);
+    for k in 0..nz {
+        for j in 0..ny {
+            let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
+            for i in 0..=nx {
+                field.u[g.uface(k, j, i)] = u0;
+            }
+        }
+    }
+    solver.initialize(&mut field);
+    let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H);
+    let steps = (8.0 / dt).ceil() as usize;
+    let start = std::time::Instant::now();
+    let mut last = (0.0, 0.0);
+    for step in 0..steps {
+        let r = solver.advance(&mut field, dt);
+        if (step + 1) % (steps / 20).max(1) == 0 || step + 1 == steps {
+            let t = solver.time();
+            let mask = solver.mask().unwrap();
+            let body = solver.body().unwrap();
+            let mu = RHO * NU;
+            let (po, so) = mask.cut_wall_force_parts(body, &field, mu, t).unwrap();
+            let ex = mask
+                .cut_wall_exchange_force(body, &field, mu, RHO, t, None)
+                .unwrap();
+            let (pr, sr) = mask
+                .cut_wall_force_reconstructed_parts(body, &field, mu, t, None)
+                .unwrap();
+            let margin = 3.0 * D;
+            let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2);
+            let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2);
+            let bx = (
+                ci(CX - margin),
+                ci(CX + margin),
+                cj(CY - 0.15),
+                cj(CY + 0.15),
+                0,
+                nz,
+            );
+            let fcv = mask.control_volume_force_with_walls(&field, dt, RHO, mu, None, bx, true);
+            let cd = |f: [f64; 3]| coef * f[0];
+            println!(
+                "  t {t:7.3}: c_D operator {:.4} (p {:.4} + s {:.4} + exchange {:.4}) | reconstructed {:.4} (p {:.4} + s {:.4}) | box {:.4}; residual {:.1e} [{:.0} s]",
+                cd(po) + cd(so) + cd(ex),
+                cd(po),
+                cd(so),
+                cd(ex),
+                cd(pr) + cd(sr),
+                cd(pr),
+                cd(sr),
+                cd(fcv),
+                r.final_residual,
+                start.elapsed().as_secs_f64()
+            );
+            let now = (cd(po) + cd(so) + cd(ex), cd(fcv));
+            if (now.0 - last.0).abs() < 1e-4 * now.0.abs()
+                && (now.1 - last.1).abs() < 1e-4 * now.1.abs()
+                && t > 2.0
+            {
+                println!("  settled");
+                break;
+            }
+            last = now;
+        }
+    }
+}