diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 8eabdce99f..0733802b08 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -86,23 +86,28 @@ int main(int argc, char *argv[])
         for (int corr=1; corr<=1; corr++)
         {
             volScalarField rAU(1.0/UEqn.A());
+            volVectorField HbyA("HbyA", U);
+            HbyA = rAU*UEqn.H();
 
-            U = rAU*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rAU, U, phi);
+            surfaceScalarField phiHbyA
+            (
+                "phiHbyA",
+                (fvc::interpolate(HbyA) & mesh.Sf())
+              + fvc::ddtPhiCorr(rAU, U, phi)
+            );
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(rAU, p) == fvc::div(phi)
+                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
             );
 
             pEqn.solve();
 
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
 
             #include "continuityErrs.H"
 
-            U -= rAU*fvc::grad(p);
+            U = HbyA - rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index cb25e7029f..723df7af98 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -34,19 +35,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -54,7 +58,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -62,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index d1b7135649..00603dc7d9 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -9,7 +10,7 @@ if (pimple.transonic())
     (
         "phid",
         fvc::interpolate(psi)
-       *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
     while (pimple.correctNonOrthogonal())
@@ -31,15 +32,19 @@ if (pimple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho)
-         *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -47,7 +52,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -55,7 +60,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 872a0513a3..f5364ae314 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -2,25 +2,29 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
-surfaceScalarField phiU
+surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
 (
+    "phiHbyA",
     fvc::interpolate(rho)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
+  - phig
 );
 
-phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
         fvc::ddt(psi, rho)*gh
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
@@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi += p_rghEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U.correctBoundaryConditions();
     }
 }
 
@@ -41,8 +47,6 @@ p = p_rgh + rho*gh;
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index cb25e7029f..723df7af98 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -34,19 +35,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -54,7 +58,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -62,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index 4e94de243c..ef11c677ad 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -6,27 +6,25 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
     if (pimple.transonic())
     {
-        surfaceScalarField phiv
+        surfaceScalarField phiHbyA
         (
-            (fvc::interpolate(U) & mesh.Sf())
+            "phiHbyA",
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-        phi = fvc::interpolate(rho)*phiv;
+        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
 
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(thermo.psi())*phiv
-        );
+        phiHbyA *= fvc::interpolate(rho);
 
         fvScalarMatrix pDDtEqn
         (
-            fvc::ddt(rho) + fvc::div(phi)
+            fvc::ddt(rho) + fvc::div(phiHbyA)
           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
         );
 
@@ -42,23 +40,26 @@
 
             if (pimple.finalNonOrthogonalIter())
             {
-                phi += pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }
     else
     {
-        phi =
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
             fvc::interpolate(rho)
            *(
-                (fvc::interpolate(U) & mesh.Sf())
+                (fvc::interpolate(HbyA) & mesh.Sf())
               + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
+            )
+        );
 
         fvScalarMatrix pDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
         );
 
         while (pimple.correctNonOrthogonal())
@@ -73,7 +74,7 @@
 
             if (pimple.finalNonOrthogonalIter())
             {
-                phi += pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }
@@ -84,7 +85,7 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index f7bb932fdc..cfe34f2873 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -3,7 +3,8 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -17,7 +18,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -41,12 +42,15 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
-        fvc::interpolate(rho)*
-        (
-            (fvc::interpolate(U) & mesh.Sf())
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
@@ -54,7 +58,7 @@ else
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -62,7 +66,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -81,7 +85,7 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index 831f139e81..c5c7602a43 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -4,7 +4,8 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -18,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -43,13 +44,17 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
-        fvc::interpolate(rho)*
-        (
-            (fvc::interpolate(U) & mesh.Sf())
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
-    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+        )
+    );
+
+    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -57,7 +62,7 @@ else
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -65,7 +70,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -83,7 +88,7 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4d3f1ac76e..56c444cdab 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -4,7 +4,9 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
 UEqn.clear();
 
 bool closedVolume = false;
@@ -14,7 +16,7 @@ if (simple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
     while (simple.correctNonOrthogonal())
@@ -40,14 +42,19 @@ if (simple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    closedVolume = adjustPhi(phi, U, p);
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    closedVolume = adjustPhi(phiHbyA, U, p);
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rho*rAU, p) == fvc::div(phi)
+            fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -56,7 +63,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 }
@@ -67,7 +74,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 // For closed-volume cases adjust the pressure and density levels
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index b1f83304cd..c3e0a43a15 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -1,10 +1,12 @@
+volVectorField HbyA("HbyA", U);
+
 if (pressureImplicitPorosity)
 {
-    U = trTU()&UEqn().H();
+    HbyA = trTU() & UEqn().H();
 }
 else
 {
-    U = trAU()*UEqn().H();
+    HbyA = trAU()*UEqn().H();
 }
 
 UEqn.clear();
@@ -16,7 +18,7 @@ if (simple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
@@ -45,10 +47,15 @@ if (simple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho*U) & mesh.Sf();
-    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho*HbyA) & mesh.Sf()
+    );
 
-    closedVolume = adjustPhi(phi, U, p);
+    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    closedVolume = adjustPhi(phiHbyA, U, p);
 
     while (simple.correctNonOrthogonal())
     {
@@ -56,11 +63,11 @@ else
 
         if (pressureImplicitPorosity)
         {
-            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi));
+            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA));
         }
         else
         {
-            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi));
+            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA));
         }
 
         tpEqn().setReference(pRefCell, pRefValue);
@@ -69,7 +76,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= tpEqn().flux();
+            phi = phiHbyA - tpEqn().flux();
         }
     }
 }
@@ -81,11 +88,11 @@ p.relax();
 
 if (pressureImplicitPorosity)
 {
-    U -= trTU()&fvc::grad(p);
+    U = HbyA - (trTU() & fvc::grad(p));
 }
 else
 {
-    U -= trAU()*fvc::grad(p);
+    U = HbyA - trAU()*fvc::grad(p);
 }
 
 U.correctBoundaryConditions();
diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
index ccc7a1b21e..a4d9325522 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
@@ -7,7 +7,10 @@ volScalarField p0(p);
 
 volScalarField AU(UEqn().A());
 volScalarField AtU(AU - UEqn().H1());
-U = UEqn().H()/AU;
+
+volVectorField HbyA("HbyA", U);
+HbyA = UEqn().H()/AU;
+
 UEqn.clear();
 
 bool closedVolume = false;
@@ -19,7 +22,7 @@ if (simple.transonic())
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi*U) & mesh.Sf()
+            fvc::interpolate(psi*HbyA) & mesh.Sf()
         );
 
         surfaceScalarField phic
@@ -56,13 +59,18 @@ else
 {
     while (simple.correctNonOrthogonal())
     {
-        phi = fvc::interpolate(rho*U) & mesh.Sf();
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            fvc::interpolate(rho*HbyA) & mesh.Sf()
+        );
+
         closedVolume = adjustPhi(phi, U, p);
         phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
 
         fvScalarMatrix pEqn
         (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
         //- fvm::laplacian(rho/AU, p)
           - fvm::laplacian(rho/AtU, p)
         );
@@ -73,7 +81,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -85,8 +93,8 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
-//U -= fvc::grad(p)/AU;
+U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
+//U = HbyA - fvc::grad(p)/AU;
 
 U.correctBoundaryConditions();
 
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index a4311dfd2e..46cc36216c 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -1,14 +1,15 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
 );
@@ -33,5 +34,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 98cc22d15a..cad5bb7217 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -1,14 +1,15 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = UEqn.H()/UEqn.A();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       - fvc::meshPhi(rho, U)
     )
 );
@@ -29,5 +30,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index a91ee8885e..51f04820eb 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -93,17 +93,21 @@ int main(int argc, char *argv[])
             for (int corr=0; corr<nCorr; corr++)
             {
                 volScalarField rAU(1.0/UEqn.A());
+                volVectorField HbyA("HbyA", U);
+                HbyA = rAU*UEqn.H();
 
-                U = rAU*UEqn.H();
-
-                phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rAU, U, phi);
+                surfaceScalarField phiHbyA
+                (
+                    "phiHbyA",
+                    (fvc::interpolate(HbyA) & mesh.Sf())
+                  + fvc::ddtPhiCorr(rAU, U, phi)
+                );
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                 {
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rAU, p) == fvc::div(phi)
+                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -111,13 +115,13 @@ int main(int argc, char *argv[])
 
                     if (nonOrth == nNonOrthCorr)
                     {
-                        phi -= pEqn.flux();
+                        phi = phiHbyA - pEqn.flux();
                     }
                 }
 
                 #include "continuityErrs.H"
 
-                U -= rAU*fvc::grad(p);
+                U = HbyA - rAU*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
         }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 14e87ea844..9db0897efa 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -2,19 +2,24 @@
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phi = (fvc::interpolate(U) & mesh.Sf())
-        + fvc::ddtPhiCorr(rAU, U, phi);
+    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
-    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
-    phi -= buoyancyPhi;
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, U, phi)
+      - phig
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -24,14 +29,14 @@
         if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 5a9ec66707..61d4358a8a 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -2,20 +2,27 @@
     volScalarField rAU("rAU", 1.0/UEqn().A());
     surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(U) & mesh.Sf();
-    adjustPhi(phi, U, p_rgh);
+    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
-    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
-    phi -= buoyancyPhi;
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA -= phig;
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -25,14 +32,14 @@
         if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 397c1c8bf4..462b96b17d 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -8,21 +8,26 @@
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phi = fvc::interpolate(rho)*
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
-        (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
+      + phig
     );
 
-    surfaceScalarField buoyancyPhi(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi += buoyancyPhi;
-
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
     );
 
     while (pimple.correctNonOrthogonal())
@@ -38,14 +43,14 @@
         if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi += p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
             K = 0.5*magSqr(U);
         }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index d2793d7192..586cddb55b 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -5,20 +5,27 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    bool closedVolume = adjustPhi(phi, U, p_rgh);
+    surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
-    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi -= buoyancyPhi;
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA -= phig
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -27,14 +34,14 @@
         if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index a582d09213..76f6f8e30e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -8,24 +8,27 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
+      + phig
     );
 
-    phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
-
     {
         fvScalarMatrix p_rghDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
         );
 
         // Thermodynamic density needs to be updated by psi*d(p) after the
@@ -57,7 +60,11 @@
 
             if (nonOrth == nNonOrthCorr)
             {
-                phi += p_rghEqn.flux();
+                phi = phiHbyA + p_rghEqn.flux();
+                U = HbyA
+                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+                U.correctBoundaryConditions();
+                K = 0.5*magSqr(U);
             }
         }
 
@@ -65,9 +72,6 @@
         thermo.rho() += psi*p_rgh;
     }
 
-    // Correct velocity field
-    U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-    U.correctBoundaryConditions();
     p = p_rgh + rho*gh;
 
     // Update pressure time derivative
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index 8a77a2656e..e6df2ca5cf 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -6,27 +6,23 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*(UEqn == sources(rho, U))().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
-    if (pZones.size() > 0)
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    if (pZones.size() == 0)
     {
-        // ddtPhiCorr not well defined for cases with porosity
-        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    }
-    else
-    {
-        phi =
-            fvc::interpolate(rho)
-           *(
-                (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
+        // ddtPhiCorr only used without porosity
+        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
     }
 
+    phiHbyA *= fvc::interpolate(rho);
+
+
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
      ==
         parcels.Srho()
       + sources(psi, p, rho.name())
@@ -46,7 +42,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 
@@ -58,7 +54,7 @@
     #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index cb8553bbfb..3a61b28e60 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*(UEqn == sources(rho, U))().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*(UEqn == sources(rho, U))().H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -39,19 +40,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
@@ -64,7 +68,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -72,7 +76,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 sources.correct(U);
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 6ed7affe77..4c9b41bba4 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -6,27 +6,23 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*(UEqn == sources(rho, U))().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
-    if (pZones.size() > 0)
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    if (pZones.size() == 0)
     {
-        // ddtPhiCorr not well defined for cases with porosity
-        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    }
-    else
-    {
-        phi =
-            fvc::interpolate(rho)
-           *(
-                (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
+        // ddtPhiCorr only used without porosity
+        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
     }
 
+    phiHbyA *= fvc::interpolate(rho);
+
+
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
      ==
         parcels.Srho()
       + sources(psi, p, rho.name())
@@ -47,7 +43,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 
@@ -57,7 +53,7 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 872a0513a3..f5364ae314 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -2,25 +2,29 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
-surfaceScalarField phiU
+surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
 (
+    "phiHbyA",
     fvc::interpolate(rho)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
+  - phig
 );
 
-phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
         fvc::ddt(psi, rho)*gh
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
@@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi += p_rghEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U.correctBoundaryConditions();
     }
 }
 
@@ -41,8 +47,6 @@ p = p_rgh + rho*gh;
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 5ee551a0ea..c504fdf2ea 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -36,19 +37,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
@@ -58,7 +62,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -66,7 +70,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index 5109788ea1..5c4fa0adc1 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -36,19 +37,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
@@ -58,7 +62,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -66,7 +70,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 6fa1965a24..34898172ce 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -9,11 +9,13 @@
         )/psi;
     }
 
-    surfaceScalarField rhof(fvc::interpolate(rho, "rhof"));
+    surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
-    volVectorField HbyA("HbyA", rAU*UEqn.H());
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phiv);
@@ -22,8 +24,6 @@
 
     phiv -= phiGradp/rhof;
 
-    #include "resetPhivPatches.H"
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H
deleted file mode 100644
index e7d0c2f93e..0000000000
--- a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H
+++ /dev/null
@@ -1,15 +0,0 @@
-fvsPatchScalarFieldField& phiPatches = phi.boundaryField();
-const fvPatchScalarFieldField& rhoPatches = rho.boundaryField();
-const fvPatchVectorFieldField& Upatches = U.boundaryField();
-const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField();
-
-forAll(phiPatches, patchI)
-{
-    if (phi.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phiPatch =
-            refCast<calculatedFvsPatchScalarField>(phiPatches[patchI]);
-
-        phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]);
-    }
-}
diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H
deleted file mode 100644
index 2f62dcfc37..0000000000
--- a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H
+++ /dev/null
@@ -1,19 +0,0 @@
-surfaceScalarField::GeometricBoundaryField& phivPatches =
-    phiv.boundaryField();
-
-const volVectorField::GeometricBoundaryField& Upatches =
-    U.boundaryField();
-
-const surfaceVectorField::GeometricBoundaryField& SfPatches =
-    mesh.Sf().boundaryField();
-
-forAll(phivPatches, patchI)
-{
-    if (phiv.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phivPatch =
-            refCast<calculatedFvsPatchScalarField>(phivPatches[patchI]);
-
-        phivPatch == (Upatches[patchI] & SfPatches[patchI]);
-    }
-}