diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 397e893035..e4fa5cd6b3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -5,6 +5,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
   + turbulence->divDevRhoReff(U)
+ ==
+    sources(rho, U)
 );
 
 UEqn().relax();
@@ -13,6 +15,6 @@ sources.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
+    solve(UEqn() == -fvc::grad(p));
     K = 0.5*magSqr(U);
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 5f28a3456d..d1d1ac5859 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -24,6 +24,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relatoiveFlux(fvc::interpolate(psi), phid);
+
     volScalarField Dp("Dp", rho*rAU);
 
     while (pimple.correctNonOrthogonal())
@@ -59,6 +61,8 @@ else
         )
     );
 
+    sources.relatoiveFlux(fvc::interpolate(rho), phiHbyA);
+
     volScalarField Dp("Dp", rho*rAU);
 
     while (pimple.correctNonOrthogonal())