Take the IB force from the momentum fluxes across the body's faces

docs/documentation/case.md

@@ -859,6 +859,8 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 | 15–17   | angles about x/y/z |
 | 18–20   | centroid, x/y/z |
 
+The force and torque are the momentum the fluid exchanges with the body: the inviscid and viscous fluxes of the latest RK stage, summed over the faces between body cells and fluid cells, with the torque taken about the centroid using each face's centre. Collision and body forces (`bf_x` etc. times `mass`) are added to them.
+
 The file carries no header line, because every record sits at a computed byte offset and a header would shift them all. Each record is exactly 353 bytes including its newline (`I10` followed by nineteen `1X,ES17.9E3` fields), so the whole file loads with `numpy.loadtxt` and a single body or step can be read without scanning it:
 
 ```

examples/2D_ibm_poiseuille_nn/README.md

@@ -2,10 +2,9 @@
 
 Validates the **immersed boundary (IBM) + non-Newtonian viscosity interaction**:
 the channel's no-slip walls are two rectangular IB slabs instead of domain
-boundary conditions, so the flow exercises the per-stencil-sample Herschel-Bulkley
-viscosity `mu_eff` used by the IBM ghost-point and force machinery
-(`s_compute_viscous_stress_tensor` in `m_viscous.fpp`, consumed by
-`s_compute_ib_forces` in `m_ibm.fpp`). Companion to
+boundary conditions, so the flow exercises the Herschel-Bulkley viscous fluxes at
+the IB walls, which `s_accumulate_ib_face_fluxes` in `m_ibm.fpp` also sums into
+the IB force. Companion to
 `examples/2D_poiseuille_thickening_nn`, which validates the same fluid against the
 same analytic profile with BC walls.
 
@@ -20,7 +19,7 @@ face; the slab centroids sit just *inside* the domain (a centroid exactly on the
 boundary is owned by no rank and its `ib_state` force record is never written).
 The domain BCs behind the slabs are no-slip walls (`bc_y = -16`), which keep the
 body-forced dead fluid inside the slabs benign. `patch_ib%mass = 0` so the
-reported IB force is the pure pressure+viscous volume integration (no
+reported IB force is the fluid force alone (no
 `bf_x*mass` bookkeeping term); `ib_state_wrt = T` writes it at every save.
 
 Fluid and forcing match the BC-walled template (single Papanastasiou-regularized
@@ -68,7 +67,7 @@ and `nn1` modes into two scratch copies and compare:
 with `nn = 1, tau0 = 0, K = 0.02` is analytically the same fluid as a Newtonian
 one with `mu = 0.02`. Running both modes with the same fixed `dt = 6e-5` to
 `t = 0.3` gives **bitwise identical** velocity fields (max abs and rel L2
-difference `0.0`) *and* bitwise identical IBM-integrated wall forces — the
+difference `0.0`) *and* bitwise identical IB wall forces — the
 non-Newtonian IBM path reduces exactly to the Newtonian one.
 
 **B — analytic power-law profile (n = 1.5).** Relative L2 error of the steady
@@ -81,13 +80,11 @@ saves 2.0e-3; profile bluntness (mean/peak) 0.640 vs. the `n = 1.5` theory
 The BC-walled template achieves 1.46% on the same fluid; the extra error is the
 diffuse-wall representation, not the viscosity model.
 
-**C — IBM-integrated wall force.** The volume-integrated x-force converges to
-**8.05e-4** per wall (both walls identical by symmetry; plateaued by `t = 0.9`)
-vs. the analytic `tau_w*L_x = 1.0e-3` — a ratio of **0.80**. The deficit is the
-known coarseness of the volume-integration force estimator (second-order
-finite differences of ghost/dead-cell states inside the body), not the
-viscosity model: in Validation A the same integral is bitwise identical between
-the Newtonian and non-Newtonian code paths.
+**C — IB wall force.** The IB force, the momentum flux the solver applies across
+the faces between slab and fluid cells, reaches **9.97e-4** per wall at
+`t = 0.9` (both walls identical by symmetry) vs. the analytic
+`tau_w*L_x = 1.0e-3`, a ratio of **0.997**. The earlier volume-integration
+estimator, which differenced ghost-cell states inside the body, gave 0.971.
 
 ## References
 

examples/2D_ibm_poiseuille_nn/case.py

@@ -1,8 +1,8 @@
 #!/usr/bin/env python3
 """
 2D power-law Poiseuille channel whose no-slip walls are IMMERSED BOUNDARIES,
-validating the IBM + non-Newtonian viscosity interaction (the per-stencil-sample
-Herschel-Bulkley viscosity used by the IBM ghost-point and force machinery).
+validating the IBM + non-Newtonian viscosity interaction (the Herschel-Bulkley
+viscous fluxes at the IB walls, which also make up the IB force).
 
 Two rectangular IB slabs form the channel: the flow gap is y in [0.05, 0.25]
 (half-height H = 0.1, centerline y = 0.15) inside a y-domain [0, 0.3]. Each slab
@@ -93,7 +93,7 @@
     # the domain in x and half outside the domain in y, so the only IB surface
     # the flow sees is the flat face at y = 0.05 (slab 1) / y = 0.25 (slab 2).
     # mass = 0 zeroes the bf_x*mass bookkeeping term in the reported IB force,
-    # so ib_state files hold the pure pressure+viscous surface integration.
+    # so ib_state files hold the fluid force alone.
     "ib": "T",
     "num_ibs": 2,
     "ib_state_wrt": "T",

examples/2D_ibm_poiseuille_nn/check_equivalence.py

@@ -6,7 +6,7 @@
 with mu = K, and for a single fluid the non-Newtonian mixture arithmetic reduces
 to the same 1/K — so a run with IBM_NN_MODE=newtonian (mu = 0.02, NN code path
 OFF) and a run with IBM_NN_MODE=nn1 (K = 0.02, nn = 1, NN code path ON,
-including the per-stencil-sample IBM viscosity) must produce the same fields to
+including the IB force's viscous face fluxes) must produce the same fields to
 near round-off. Both modes use the SAME fixed dt, so the comparison is at
 matched time step with identical dt histories.
 

examples/2D_ibm_poiseuille_nn/compare_analytic.py

@@ -11,9 +11,9 @@
 each wall, since IBM walls are sharp only to ~half a cell).
 
 Validation C: the analytic steady x-force per wall per unit depth is
-tau_w * L_x with tau_w = rho*g*H; compared against the IBM volume-integrated
+tau_w * L_x with tau_w = rho*g*H; compared against the IB face-flux
 force in restart_data/ib_state_<last>.dat (written because ib_state_wrt = T;
-patch_ib%mass = 0 so the record is the pure pressure+viscous integration).
+patch_ib%mass = 0 so the record is the fluid force alone).
 
 Run:  ./mfc.sh run examples/2D_ibm_poiseuille_nn/case.py -n 2
       ./build/venv/bin/python3 examples/2D_ibm_poiseuille_nn/compare_analytic.py

examples/2D_ibm_viscous_drag_over_cylinder/readme.md

@@ -16,6 +16,8 @@ Of particular note in getting the numerics of this case to match closely to thos
 
 To test the required convergence, comparisons to 4 different pressure and viscous drag force integration methods was made: surface integration and volume integration both computed at 2nd and 4th order. These values were compared in a post-processing step. In this case file, we make use of `fd_order: 4`, which enabled 4th order finite differencing techniques to be used in the IB drag force volume integraiton method of MFC. The following plot demonstrates that the MFC calculation matches the post-processing impliementation of 4th-order, and that this technique gives the appropriate convergence.
 
+The plot and the value below come from MFC's earlier IB force, a volume integral of the pressure gradient and viscous-stress divergence over the body's cells. MFC now takes the IB force from the momentum fluxes the solver applies across the faces between body and fluid cells, since the volume integral missed part of the friction drag (#2013).
+
 ![Screenshot](drag_coeffs.png)
 
 The final MFC-obtained value of the drag coefficient is $C_{d}^{MFC} = 1.549$, which differs only slightly from the value from Canuto and Taira of $C_{d}^{Canuto} = 1.540$.

src/simulation/m_derived_variables.fpp

@@ -38,9 +38,7 @@ contains
         ! to be implemented in the subroutine s_compute_finite_difference_coefficients.
 
         ! Allocating centered finite-difference coefficients. The coefficient (second) index is extended by fd_number beyond
-        ! the interior on each side: s_compute_ib_forces evaluates the viscous-stress stencil centered on ghost-adjacent
-        ! cells (i+l for l in -fd_number:fd_number) when an IB sits near a domain boundary, so the coefficient array must
-        ! cover those centers too, not just the interior 0:m.
+        ! the interior on each side, as s_compute_finite_difference_coefficients fills it.
         if (probe_wrt .or. ib) then
             @:ALLOCATE(fd_coeff_x(-fd_number:fd_number,-fd_number:m + fd_number))
             if (n > 0) then