I guess if we want to use turbulence model, we'd better stick to the original definition of phi = fvc::interpoplate(U) & mesh.Sf(), because, for example, if we want to use k-epsilon model, in k-equation and epsilon equation, the convection term could use phi as well.
Meanwhile, in p-equation, we need to account for the divergence of the face velocity field by taking out the difference between the interpolated velocity and the flux: fvc::ddtPhiCorr(rUA, U, phi)
I guess this is just a minor bug, if we don't consider to use turbulence model, however, it could make a difference when we use turbulence model such as k-Epsilon model, am I right?