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

[sbryngelson]

Fixes #2013

s_compute_ib_forces summed −∇p + ∇·τ over the cells inside each body. Its stencils read ghost cells, whose velocity is the wall value, so the viscous part missed 16–69 % of the friction the flow carries, and D/ib_forces.dat (and the force on moving_ibm = 2 bodies) was 5–36 % low.

The IB force is now the momentum the fluid exchanges with each body: the fluxes the RHS differences (Riemann flux_n plus the viscous/capillary flux_src_n), summed over every face between a body cell and a fluid cell. s_accumulate_ib_face_fluxes does this as each direction's fluxes are formed in s_compute_rhs (one device pass over the faces per direction and stage, only when forces are needed), and s_compute_ib_forces turns the sums into patch_ib%force/torque with the existing collision forces, neighborhood MPI reduction and body forces, all unchanged.

• Torque comes from the same face fluxes, with the lever arm from the centroid of the body's periodic image to the face centre. In 2D the z torque is now kept; the volume sum dropped it (it summed torque components 1:num_dims only).
• A face counts on the rank owning its fluid cell, so no face is counted twice. Moving bodies get the force of the current RK stage, as before.
• The volume sum and s_compute_viscous_stress_tensor, now unused, are removed.
• igr and hypoelastic HLLD (dual pass) do not form one set of face fluxes, so the validator rejects IB forces (ib_state_wrt or moving IBs) with them.
• Docs: the force definition is added to the D/ib_forces.dat section of case.md. The 2D_ibm_poiseuille_nn README and comments are updated, and the 2D_ibm_viscous_drag_over_cylinder readme notes that its plot came from the old estimator.

DRY note. The face loop follows s_accumulate_ib_face_fluxes in #1957 (lumped IB temperature), which sums energy and mass over the same faces. That routine is not on master and #1957 is stacked on #1956, so this PR does not stack on it. It puts the face loop on master instead. #1957 should then add its energy and mass sums to this loop rather than keep its own.

Verification

Cd from the IB force vs Cd_CV, the drag from a momentum balance on a 4D × 4D box from the saved fields (analyze.py, independent of the IB force), M = 0.2, inert moist air with mixture-averaged transport (campaign build: this commit on a scratch copy of our campaign branch, MI300A):

case	cells/D	Cd_IB before	Cd_IB after	Cd_CV	after vs CV	reference
2D cylinder, Re 40	20	1.311	1.505	1.506	−0.1 %	1.52 (Dennis & Chang)
2D cylinder, Re 40	40	1.385	1.512	1.511	+0.1 %	1.52
2D cylinder, Re 40	80	1.437	1.522	1.520	+0.1 %	1.52
2D cylinder, Re 100 (mean)	20	1.137	1.358	1.359	−0.1 %	1.33–1.35
sphere, Re 100 (Tw/T∞ = 1.1)	20	0.815	1.097	1.096	+0.1 %	1.092 (Schiller–Naumann)
2D cylinder, Re 40, issue repro (from t = 0)	20	1.291	1.485	1.488	−0.2 %	1.52

The first five runs restart the original runs from the save at 80 % of t_end (75 % for the sphere) and cover the same averaging window, the last 20 %. Cd_CV of the new runs matches the original runs to 4 digits: the flow is unchanged and only the force output differs. Cd_IB now converges with resolution together with Cd_CV (1.505 / 1.512 / 1.522 at Re 40). Before, it converged at order ≈ 0.5 from 13 % below.

On master alone (CPU, 48–96 ranks, ideal gas with fluid_pp%Re, no chemistry): a 2D cylinder at Re 40, M 0.2, domain [−5D, 10D] × [−6D, 6D], t* = 30:

cells/D	Cd_IB before	Cd_IB after	Cd_CV
20	1.289	1.481	1.471
40	1.365	1.491	1.479

examples/2D_ibm_poiseuille_nn (IB-walled power-law channel, CPU): the wall force is 0.971 of the analytic τ_w L_x before and 0.997 after.

GPU: the kernel ran on MI300A (OpenMP offload, CCE 19) in all the campaign runs above. A master GPU build ran the master cylinder at 20 cells/D: Cd_IB 1.481, the same as on CPU to 4 digits.

Regenerated goldens. Only these change, because only they record the IB force or move a body with it:

• 09FDDC81 2D → 1 Fluid → IBM → Adaptive dt (moving_ibm = 2): the motion changes at the 1e-10 level.
• E085CC5A, 5A22B45F, 4BED9896, C8AD6271, 49893269, 135F548B (2D/3D particle clouds, Box and Hemisphere Shell, with ib_force_wrt): D/ib_forces.dat changes. These runs start impulsively, and the first-step pressure pulse on each particle now gives a force of ~1e-5. The old volume sum gave ~1e-8, because it read the ghost-cell pressure.

./mfc.sh test --no-mpi --no-gpu -j 24 --only IBM: 61 passed, 0 failed after regeneration. ./mfc.sh precheck passes.

Not covered: axisymmetric (cyl_coord) bodies, where the force is Cartesian as before; Euler–Euler bubbles, whose non-conservative momentum sources are not face fluxes; and multi-rank moving-IB runs beyond the regenerated tests.

This PR was made with Claude Code (AI) on LLNL Tuolumne.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

[danieljvickers]

@sbryngelson this is a significant rewrite of a major component of the IBM code. Can you please provide citations and more documentation about where this comes from and why it is valid? I tried to do something similar a year ago and you explicitly turned me away from it in favor of copying what was done in jCODE.

Have we also shown anything about the convergence with a moving case here?

[sbryngelson]

@danieljvickers I dug into this further with Codex. Your moving-body question is important: the stationary tests support a discrete-consistency improvement, but they do not establish moving-body convergence. We also found a pre-existing error in the body RK update that should be addressed before using two-way trajectories to validate the new force.

Standalone reproducers, checked results, the CPU diagnostic patch, and channel inputs are here: reproducers and numerical evidence. Everything below refers to base a6add01 and this PR's head e85a5b0. The small Fortran tests were run locally with GNU Fortran 15.2, runtime/bounds checks, floating-point traps, and both -O0 and -O2.

What the force change does, and the derivation

There was nothing fundamentally wrong with the continuum volume identity: integrating -grad(p) + div(tau) can give the surface traction with an appropriate stress extension. The problem is our discrete implementation. It differentiates the IBM velocity extension to construct stress, then differentiates that stress again with a centered FD operator. The flow update transfers momentum through a different numerical viscous-gradient/reconstruction/flux operator. Increasing fd_order does not make the two operators identical.

For a fixed Cartesian body mask, write the numerical momentum RHS as

$$ V_{i}\dot{\boldsymbol{m}}_{i}=-\sum_{f\in\partial i} A_{f},\widehat{\boldsymbol{\Pi}}_{f}\boldsymbol{n}_{if}+V_{i}\boldsymbol{S}_{i}. $$

Here $\widehat{\boldsymbol{\Pi}}$ denotes the momentum flux actually differenced by the RHS. In this implementation that is the Riemann momentum flux plus the applicable momentum source flux (flux_n + flux_src_n), with the code's signs. Sum over fluid cells. Fluid-fluid face contributions cancel, leaving outer-domain fluxes and body-fluid face contributions:

$$ \dot{\boldsymbol{P}}_{\mathrm{fluid}}= \boldsymbol{I}_{\mathrm{outer}}+\boldsymbol{S}_{\mathrm{fluid}}- \sum_{b}\boldsymbol{F}_{b}^{\mathrm{flux}}. $$

The proposed force is precisely the opposite of the body-face contribution to this fixed-mask numerical momentum update. Torque uses the same face forces and their face-centre lever arms. A compatible volume sum of the same unmasked numerical flux divergence would give the same force; it is not necessary to reject volume integration in general. The important change is reusing the operator that advances the fluid instead of reconstructing a different one from ghost-cell fields. This also explains why the routine must accumulate during each directional sweep: the flux buffers are reused between directions.

That proves a discrete momentum-accounting identity. It does not, by itself, prove physical traction convergence, accuracy of the boundary state, global angular-momentum conservation of the entire fluid scheme, or conservation through a moving-mask update. In the fixed-body impermeable-wall limit, physical convergence still requires the numerical fluxes/boundary treatment to converge to the correct wall momentum exchange.

Local verification and validation of the force

Actual accumulator, synthetic fluxes. I extracted and compiled the PR's actual s_accumulate_ib_face_fluxes, with the GPU wrappers removed and small type/lookup stubs. There are 64 fixtures and 256 binary executions covering 2D/3D, uniform/stretched face areas, two bodies, periodic-image lever arms, source-flux inclusion, pressure-offset invariance, and one domain versus two x partitions with halos and a body crossing the partition boundary. Independent body-cell boundary enumeration checks force and torque; a sum of unmasked flux divergence checks force again. Maximum discrepancy across the reference, partition, and optimization comparisons was 1.32e-13. This is a serial ownership test, not an MPI-runtime or GPU test.

Full MFC channel runs on identical evolving fields. I kept the base's old force and added the PR accumulator as a diagnostic alongside it. The bodies are fixed and have zero mass, so the added sums do not change the solution. In Newtonian mode of 2D_ibm_poiseuille_nn, rho=1, mu=0.02, g=0.05, channel gap 0.2, streamwise length 0.2, 26 x cells, WENO5/mapped/MPWENO, fd_order=4, RK3, the force ratios at nominal t=0.96 were:

y cells	viscous reconstruction	old / transient reference	face / transient reference
48	WENO	1.036065	0.995612
96	WENO	1.024741	0.996559
48	direct (weno_Re_flux=F)	0.840601	0.997345
96	direct	0.843979	0.998752

The reference is the constant-density transient Poiseuille wall force, evaluated at the last RHS's RK3 half-step time. It is an independent analytic reference for this weakly compressible, very-low-Mach case, not an exact solution for every MFC variable. Only the streamwise wall force is assessed: these wall slabs extend through the y-domain boundaries, so normal-force/torque accounting would need those external faces separately. Halving dt changes the old/reference and face/reference ratios by about 1.2e-6 and 0.9e-6. A clean rebuild and repeat of the 48-cell WENO run gave bitwise-identical final conservative fields, force diagnostics, and sampled profile.

A predictive explanation of the old error. For a planar wall with zero interior velocity, the old fourth-order force reduces to

$$ F_{\mathrm{old}}/L_{x}=(\mu/h)(65u_{0}/144+49u_{1}/144-5u_{2}/48+u_{3}/144). $$

For the steady discrete direct-gradient channel solution with $N$ fluid rows, substituting the numerical profile gives

$$ F_{\mathrm{old}}/F_{\mathrm{face}}=61/72-5/(36N). $$

For N=32, this predicts 0.842881944; the full solver at t=1.92 gives 0.842881919 (difference 2.6e-8). Old/steady analytic force is 0.84278504, and face/steady is 0.99988507. The predicted refinement limit is 61/72, not one. This establishes a genuine operator mismatch on a numerical solution, not just an exact-field counterexample.

The error is not universally an underestimate. With WENO reconstruction, the old channel force is slightly high. Putting exact steady parabolic samples into its functional gives 0.4902 of the steady wall force, whereas the computed near-wall samples give 1.0287. Thus near-wall numerical errors can compensate for an incompatible force functional. Flat and circular manufactured-field controls also show strong sensitivity to the interior extension; those controls alone do not prove the new force correct on exact sampled fields.

A completed short cylinder check at Re=40, Ma=0.1, 20 cells/D, low_Mach=2, MPWENO off, domain [-3,7] x [-4,4], and t approximately 2 gave old Cd=1.611360413 and face Cd=1.631049296. Two enclosing boxes give the face force to 3.3e-15 force units when their outer influx and fluid momentum RHS are accounted for separately. That checks the discrete accounting; it shares the solver's fluxes/RHS and is not independent physical validation or a converged cylinder benchmark. Reconstructing the old force from stored fields agrees with its logged Cd to 4.4e-14. Its viscous contribution is 0.583670410, versus the new viscous-source contribution 0.555605418, another reason to avoid a universal claim that the old viscous term is always low. The new Riemann part contains numerical/advective momentum transport as well as pressure.

The larger cylinder/sphere, GPU, and multi-rank results already listed in the PR are additional reported evidence; this local audit has not rerun those campaigns.

What the earlier cylinder result established

The #1614 result (Cd=1.549 versus 1.540) remains a good reported total-drag prediction. Its figure also shows a fourth-order surface estimate around 1.44, while the fourth-order volume estimate is around 1.55 on the same solution. Those are visual readings, not recovered raw data. The horizontal axis is time step, so that figure shows temporal settling at one grid rather than spatial convergence. The earlier 900x900 comparison was also at a different Mach number.

I have not recovered the original fields/analysis script or reproduced that 3000x3000 run, so I cannot identify which pressure/shear errors produced the particular 1.549. The channel demonstrates a possible compensation mechanism, not a quantitative explanation of that historical result. Matching a Python implementation of the same volume functional verified the implementation, but did not test agreement with the fluid's momentum-transfer operator. I have not audited jCODE itself, so this is a statement about the MFC discretization, not a claim that jCODE has the same problem.

Moving bodies: confirmed RK defect and remaining force questions

The current body integrator has two problems that predate this PR:

1. It reverses the current-stage and beginning-of-step RK weights. step_vel is saved at stage 1, but the body uses (v_start + 3*v_current + dt*a)/4 at stage 2 and (2*v_start + v_current + 2*dt*a)/3 at stage 3. The fluid uses the opposite state weighting.
2. Position and angle use the just-updated velocity/angular velocity rather than the input-stage velocities. Correcting the weights alone leaves the position update first order. This kinematic issue also affects RK2.

Relevant source: coefficients, body update.

Even if stage hydrodynamic forces are equal and opposite, the original RK3 integrates the body impulse as

$$ \Delta P_{b}=\Delta t(F_{0}/4+F_{1}/12+2F_{2}/3), $$

while the fluid integrates the opposite exchange as

$$ \Delta P_{f}=-\Delta t(F_{0}/6+F_{1}/6+2F_{2}/3). $$

The residual is dt*(F0-F1)/12, generally nonzero. Constant force hides this velocity/impulse error, but position is still wrong: from rest the first step gives x-x0=a*dt^2 rather than a*dt^2/2.

I compiled the actual generated body routine with controlled force laws and no-op geometry, comparing the original, weights-only, and complete corrections. The complete correction saves input-stage velocities, puts the current/start states under coefficients 1/2 respectively, and uses the saved velocities in the position/angle RHS. The expanded checks cover 432 configurations: RK1/2/3, constant acceleration, linear/nonlinear drag, harmonic oscillation, time-dependent forcing, unequal mass/inertia, nonzero initial states, uniform/nonuniform dt, and equal/opposite exchange with a fluid surrogate. Corrected trajectories match an independent Butcher-tableau implementation within 1.31e-14, recover the expected temporal orders, and have maximum exchange momentum drift 3.55e-15. At 20 steps the original RK3 exchange drifts by 0.004016637, versus 3.33e-16 after correction.

For the oscillator x'=v, v'=-x, initial (x,v)=(1,0), at t=1:

steps	original position error	weights-only error	complete correction
20	1.498e-2	1.418e-2	2.987e-6
40	7.593e-3	7.144e-3	3.626e-7
80	3.822e-3	3.585e-3	4.465e-8
160	1.917e-3	1.796e-3	5.539e-9

The source-extracted correction is verified for translation and 2D constant-inertia rotation. It has not been installed/tested as a production GPU patch; new per-body temporaries need the appropriate GPU-private declarations. General 3D rotation/inertia and full coupled moving-case convergence remain outside these tests. kin_model=0 means the special prescribed models are off, not that the body is stationary; the two-way moving_ibm=2 branch uses it. kin_model=1/2 with one-way motion bypasses this recurrence. Eighteen stubbed branch/stage checks preserve that dispatch, not the physical flapping/pitch model.

Separately, after the flux update, MFC advances the body, rebuilds its mask, and corrects its interior/ghost states. For a split-step fluid indicator $\chi$ and momentum density $m$, the identity is

$$ \Delta P_{f}=\sum_{i} V_{i}\chi_{i}^{\mathrm{old}}(m_{i}^_-m_{i}^{\mathrm{old}}) +\sum_{i} V_{i}(\chi_{i}^{\mathrm{new}}-\chi_{i}^{\mathrm{old}})m_{i}^_ +\sum_{i} V_{i}\chi_{i}^{\mathrm{new}}(m_{i}^{\mathrm{new}}-m_{i}^*). $$

The face accumulator covers the first term's body exchange. Covered/uncovered-cell accounting and any fluid-state correction need to be checked separately, with the actual RK stage bookkeeping. I found no explicit old/new-mask impulse accumulator in this path, but have not measured a resulting physical error. It would also be wrong to transfer every interior ghost reset blindly to the body, since those are artificial states. For two-way motion the force additionally enters the acceleration-dependent ghost-pressure correction, so changing it changes the subsequent flow, not merely the body output.

Citations and what remains before a moving-body claim

Nangia, Johansen, Patankar & Bhalla (2017), A moving control volume approach to computing hydrodynamic forces and torques on immersed bodies provides the continuum force/torque balance, unsteady/moving-control-volume terms, and useful moving/oscillating-body validation cases. Its Peskin-type force equivalence is not a plug-in proof for MFC's ghost-cell implementation. Lee, Kim, Choi & Yang (2011), Sources of spurious force oscillations from an immersed boundary method for moving-body problems is directly relevant to discontinuities when cells change between solid and fluid.

So, to answer the convergence question: we have not yet demonstrated meaningful moving-body convergence for this force change. The regenerated test's 1e-10 motion change does not establish it. I would first correct the body RK bookkeeping in a separate, verified change, then require (i) a prescribed oscillating cylinder crossing multiple cells, with force amplitude/phase/crossing noise compared to a remote unsteady momentum balance under h/dt refinement, and (ii) a freely moving body in a closed/periodic domain, with stage-resolved total momentum and mask/correction accounting, including a lighter-body case. Moving MPI/offload and general 3D rotation need their own checks.

The defensible claim for this PR now is that it makes the reported force compatible with the fixed-mask fluid momentum-flux update, supported by a derivation, accumulator tests and stationary analytic comparisons. That is a substantive improvement; complete moving-body consistency and physical convergence are additional requirements, not consequences of that identity alone.