CFD Online Discussion Forums - Implement source term in rhoSimpleFoam

Greetings!

Recently I'm working on adding source term to rhoSimpleFoam.
I read through the source code and did some test, but still I have some problems unsolved. Please if someone could help:

1. Is it correct to add source term, i.e., mySource (volVectorField, constant but updated at every iteration) to the three locations?
a. In UEqn.H:

Code:

solve(UEqn == -fvc::grad(p)+mySource);

b. In pEqn.H, NonOrthogonal correction

Code:

fvScalarMatrix pEqn

        (

            fvc::div(phiHbyA)

          - fvm::laplacian(rhorAUf, p)

          + fvc::div(mySource*rAU)



          ==

            fvOptions(psi, p, rho.name())

        );

c. In pEqn.H, corrector

Code:

U = HbyA - rAU*fvc::grad(p) + mySource;

I saw a thread

HTML Code:

https://www.cfd-online.com/Forums/openfoam/65293-add-source-term-simplefoam.html

related to this. Basically they are discussing adding source term to a. and c., but why not b.?

2. Why in pEqn.H phiHbyA is interpolated like this

Code:

surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));

instead of

Code:

surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));

like the way rhorAUf did, or they are the same.

3. Why rho*HbyA is interpolated to the cell surface in the first place? may be it is more of a algorithm question but I just cant figure it out.