bug: IBM pressure-force loop reads volume fractions instead of pressure in multi-fluid cases

[sbryngelson]

Summary

In src/simulation/m_ibm.fpp, s_compute_ibm_forces_and_torques computes the pressure-gradient contribution to IBM forces using:

do fluid_idx = 0, num_fluids - 1
    local_force_contribution(1) = local_force_contribution(1) &
        - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i+1,j,k) &
         - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i-1,j,k)) / (2._wp*dx)
    ...
end do

For fluid_idx = 0 this correctly reads q_prim_vf(eqn_idx%E) — the total pressure.

For fluid_idx > 0 this reads q_prim_vf(eqn_idx%E + 1), q_prim_vf(eqn_idx%E + 2), etc., which fall in the advection (volume fraction) block (eqn_idx%adv%beg .. eqn_idx%adv%end) — not per-fluid pressures. Finite differencing volume fractions and adding them as a force contribution is incorrect.

Impact

Only affects multi-fluid IBM cases (num_fluids > 1 with immersed boundaries). Single-fluid IBM (num_fluids = 1) is unaffected since the loop runs once with fluid_idx = 0.

Fix

For multi-component total-pressure force, the loop body should use eqn_idx%E directly (no + fluid_idx offset), or the intent of this loop needs clarification for mixture IBM. If partial-pressure contributions are genuinely needed (e.g., model_eqns == 3), the index arithmetic needs to be matched to the correct EOS model.

Discovered

Found during review of PR #1365. Pre-existing in master.

[sbryngelson]

@danieljvickers is this real