diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index b5a699cc87..f796bc9bc8 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,19 +113,23 @@ int main(int argc, char *argv[])
 
             solve(UEqn() == -fvc::grad(p));
 
-            p.boundaryField().updateCoeffs();
             volScalarField rAU(1.0/UEqn().A());
-            U = rAU*UEqn().H();
+            volVectorField HbyA("HbyA", U);
+            HbyA = rAU*UEqn().H();
             UEqn.clear();
-            phi = fvc::interpolate(U) & mesh.Sf();
-            adjustPhi(phi, U, p);
+            surfaceScalarField phiHbyA
+            (
+                "phiHbyA",
+                fvc::interpolate(HbyA) & mesh.Sf()
+            );
+            adjustPhi(phiHbyA, U, p);
 
             // Non-orthogonal pressure corrector loop
             while (simple.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -133,7 +137,7 @@ int main(int argc, char *argv[])
 
                 if (simple.finalNonOrthogonalIter())
                 {
-                    phi -= pEqn.flux();
+                    phi = phiHbyA - pEqn.flux();
                 }
             }
 
@@ -143,7 +147,7 @@ int main(int argc, char *argv[])
             p.relax();
 
             // Momentum corrector
-            U -= rAU*fvc::grad(p);
+            U = HbyA - rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
@@ -174,19 +178,23 @@ int main(int argc, char *argv[])
 
             solve(UaEqn() == -fvc::grad(pa));
 
-            pa.boundaryField().updateCoeffs();
             volScalarField rAUa(1.0/UaEqn().A());
-            Ua = rAUa*UaEqn().H();
+            volVectorField HbyAa("HbyAa", Ua);
+            HbyAa = rAUa*UaEqn().H();
             UaEqn.clear();
-            phia = fvc::interpolate(Ua) & mesh.Sf();
-            adjustPhi(phia, Ua, pa);
+            surfaceScalarField phiHbyAa
+            (
+                "phiHbyAa",
+                fvc::interpolate(HbyAa) & mesh.Sf()
+            );
+            adjustPhi(phiHbyAa, Ua, pa);
 
             // Non-orthogonal pressure corrector loop
             while (simple.correctNonOrthogonal())
             {
                 fvScalarMatrix paEqn
                 (
-                    fvm::laplacian(rAUa, pa) == fvc::div(phia)
+                    fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
                 );
 
                 paEqn.setReference(paRefCell, paRefValue);
@@ -194,7 +202,7 @@ int main(int argc, char *argv[])
 
                 if (simple.finalNonOrthogonalIter())
                 {
-                    phia -= paEqn.flux();
+                    phia = phiHbyAa - paEqn.flux();
                 }
             }
 
@@ -204,7 +212,7 @@ int main(int argc, char *argv[])
             pa.relax();
 
             // Adjoint momentum corrector
-            Ua -= rAUa*fvc::grad(pa);
+            Ua = HbyAa - rAUa*fvc::grad(pa);
             Ua.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index bfc39e572c..a54af98d8d 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -68,17 +68,22 @@ int main(int argc, char *argv[])
         {
             volScalarField rAU(1.0/UEqn.A());
 
-            U = rAU*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rAU, U, phi);
+            volVectorField HbyA("HbyA", U);
+            HbyA = rAU*UEqn.H();
+            surfaceScalarField phiHbyA
+            (
+                "phiHbyA",
+                (fvc::interpolate(HbyA) & mesh.Sf())
+              + fvc::ddtPhiCorr(rAU, U, phi)
+            );
 
-            adjustPhi(phi, U, p);
+            adjustPhi(phiHbyA, U, p);
 
             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);
@@ -86,13 +91,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/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
index 29964a7d4d..9166ddcd32 100644
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
@@ -71,17 +71,22 @@ int main(int argc, char *argv[])
         {
             volScalarField rAU(1.0/UEqn.A());
 
-            U = rAU*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rAU, U, phi);
+            volVectorField HbyA("HbyA", U);
+            HbyA = rAU*UEqn.H();
+            surfaceScalarField phiHbyA
+            (
+                "phiHbyA",
+                (fvc::interpolate(HbyA) & mesh.Sf())
+              + fvc::ddtPhiCorr(rAU, U, phi)
+            );
 
-            adjustPhi(phi, U, p);
+            adjustPhi(phiHbyA, U, p);
 
             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);
@@ -89,13 +94,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/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index 72c0d66290..29f1c48c07 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -1,15 +1,20 @@
 volScalarField rAUrel(1.0/UrelEqn().A());
-Urel = rAUrel*(UrelEqn() == sources(Urel))().H();
+volVectorField HbyA("HbyA", Urel);
+HbyA = rAUrel*(UrelEqn() == sources(Urel))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
     UrelEqn.clear();
 }
 
-phi = (fvc::interpolate(Urel) & mesh.Sf())
-    + fvc::ddtPhiCorr(rAUrel, Urel, phi);
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::ddtPhiCorr(rAUrel, Urel, phi)
+);
 
-adjustPhi(phi, Urel, p);
+adjustPhi(phiHbyA, Urel, p);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
@@ -17,7 +22,7 @@ while (pimple.correctNonOrthogonal())
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAUrel, p) == fvc::div(phi)
+        fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -26,7 +31,7 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi -= pEqn.flux();
+        phi = phiHbyA - pEqn.flux();
     }
 }
 
@@ -35,6 +40,6 @@ while (pimple.correctNonOrthogonal())
 p.relax();
 
 // Momentum corrector
-Urel -= rAUrel*fvc::grad(p);
+Urel = HbyA - rAUrel*fvc::grad(p);
 Urel.correctBoundaryConditions();
 sources.correct(Urel);
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 55062637c2..444ad44503 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,14 +1,19 @@
-U = rAU*(UEqn() == sources(U))().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
 
-phi = (fvc::interpolate(U) & mesh.Sf())
-    + fvc::ddtPhiCorr(rAU, U, phi);
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::ddtPhiCorr(rAU, U, phi)
+);
 
-adjustPhi(phi, U, p);
+adjustPhi(phiHbyA, U, p);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
@@ -16,7 +21,7 @@ while (pimple.correctNonOrthogonal())
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAU, p) == fvc::div(phi)
+        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -25,7 +30,7 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi -= pEqn.flux();
+        phi = phiHbyA - pEqn.flux();
     }
 }
 
@@ -34,6 +39,6 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 sources.correct(U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index dd90826e7d..d9a6f672a0 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,29 +1,34 @@
-U = rAU*(UEqn() == sources(U))().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
 
-phi = (fvc::interpolate(U) & mesh.Sf());
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+);
 
 if (ddtPhiCorr)
 {
-    phi += fvc::ddtPhiCorr(rAU, U, phi);
+    phiHbyA += fvc::ddtPhiCorr(rAU, U, phi);
 }
 
 if (p.needReference())
 {
-    fvc::makeRelative(phi, U);
-    adjustPhi(phi, U, p);
-    fvc::makeAbsolute(phi, U);
+    fvc::makeRelative(phiHbyA, U);
+    adjustPhi(phiHbyA, U, p);
+    fvc::makeAbsolute(phiHbyA, U);
 }
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAU, p) == fvc::div(phi)
+        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -32,7 +37,7 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi -= pEqn.flux();
+        phi = phiHbyA - pEqn.flux();
     }
 }
 
@@ -44,6 +49,6 @@ p.relax();
 // Make the fluxes relative to the mesh motion
 fvc::makeRelative(phi, U);
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 sources.correct(U);
diff --git a/applications/solvers/incompressible/pisoFoam/pisoFoam.C b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
index 9dec6c12e9..cd2fbb6940 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -81,11 +81,16 @@ int main(int argc, char *argv[])
             {
                 volScalarField rAU(1.0/UEqn.A());
 
-                U = rAU*UEqn.H();
-                phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rAU, U, phi);
+                volVectorField HbyA("HbyA", U);
+                HbyA = rAU*UEqn.H();
+                surfaceScalarField phiHbyA
+                (
+                    "phiHbyA",
+                    (fvc::interpolate(HbyA) & mesh.Sf())
+                  + fvc::ddtPhiCorr(rAU, U, phi)
+                );
 
-                adjustPhi(phi, U, p);
+                adjustPhi(phiHbyA, U, p);
 
                 // Non-orthogonal pressure corrector loop
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -94,7 +99,7 @@ int main(int argc, char *argv[])
 
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rAU, p) == fvc::div(phi)
+                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -114,13 +119,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/incompressible/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
index d7cf6cd5b1..2271e85d9b 100644
--- a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
@@ -1,25 +1,30 @@
 volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
 
-U = rAU*(UEqn() == sources(U))().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
 
-phi = (fvc::interpolate(U) & mesh.Sf())
-    + fvc::ddtPhiCorr(rAU, U, phi);
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::ddtPhiCorr(rAU, U, phi)
+);
 
-adjustPhi(phi, U, p_gh);
+adjustPhi(phiHbyA, U, p_gh);
 
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
-{ 
+{
     fvScalarMatrix p_ghEqn
     (
-        fvm::laplacian(rAUf, p_gh) == fvc::div(phi)
+        fvm::laplacian(rAUf, p_gh) == fvc::div(phiHbyA)
     );
 
     p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
@@ -28,7 +33,7 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi -= p_ghEqn.flux();
+        phi = phiHbyA - p_ghEqn.flux();
     }
 }
 
@@ -39,6 +44,6 @@ p_gh.relax();
 
 p = p_gh + (g & (mesh.C() + zeta - refLevel));
 
-U -= rAU*fvc::grad(p_gh);
+U = HbyA - rAU*fvc::grad(p_gh);
 U.correctBoundaryConditions();
 sources.correct(U);
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index b491552647..645b12c58c 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -97,25 +97,30 @@ int main(int argc, char *argv[])
 
                 surfaceScalarField phih0(ghrAUf*mesh.magSf()*fvc::snGrad(h0));
 
+                volVectorField HbyA("HbyA", hU);
                 if (rotating)
                 {
-                    hU = rAU*(hUEqn.H() - (F ^ hU));
+                    HbyA = rAU*(hUEqn.H() - (F ^ hU));
                 }
                 else
                 {
-                    hU = rAU*hUEqn.H();
+                    HbyA = rAU*hUEqn.H();
                 }
 
-                phi = (fvc::interpolate(hU) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rAU, h, hU, phi)
-                    - phih0;
+                surfaceScalarField phiHbyA
+                (
+                    "phiHbyA",
+                    (fvc::interpolate(HbyA) & mesh.Sf())
+                  + fvc::ddtPhiCorr(rAU, h, hU, phi)
+                  - phih0
+                );
 
                 while (pimple.correctNonOrthogonal())
                 {
                     fvScalarMatrix hEqn
                     (
                         fvm::ddt(h)
-                      + fvc::div(phi)
+                      + fvc::div(phiHbyA)
                       - fvm::laplacian(ghrAUf, h)
                     );
 
@@ -123,11 +128,11 @@ int main(int argc, char *argv[])
 
                     if (pimple.finalNonOrthogonalIter())
                     {
-                        phi += hEqn.flux();
+                        phi = phiHbyA + hEqn.flux();
                     }
                 }
 
-                hU -= rAU*h*magg*fvc::grad(h + h0);
+                hU = HbyA - rAU*h*magg*fvc::grad(h + h0);
 
                 // Constrain the momentum to be in the geometry if 3D geometry
                 if (mesh.nGeometricD() == 3)
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
index 93d6537b6a..b7bf83e8cd 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
@@ -1,20 +1,20 @@
 {
-    p.boundaryField().updateCoeffs();
-
     volScalarField rAU(1.0/UEqn().A());
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
-    mrfZones.relativeFlux(phi);
-    adjustPhi(phi, U, p);
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    adjustPhi(phiHbyA, U, p);
+    mrfZones.relativeFlux(phiHbyA);
+    adjustPhi(phiHbyA, U, p);
 
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAU, p) == fvc::div(phi)
+            fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -22,7 +22,7 @@
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 
@@ -32,7 +32,7 @@
     p.relax();
 
     // Momentum corrector
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index b9a9079b92..e30ef7baaa 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -1,19 +1,18 @@
 {
-    p.boundaryField().updateCoeffs();
-
     volScalarField rAUrel(1.0/UrelEqn().A());
-    Urel = rAUrel*UrelEqn().H();
+    volVectorField HbyA("HbyA", Urel);
+    HbyA = rAUrel*UrelEqn().H();
     UrelEqn.clear();
 
-    phi = fvc::interpolate(Urel, "interpolate(HbyA)") & mesh.Sf();
-    adjustPhi(phi, Urel, p);
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    adjustPhi(phiHbyA, Urel, p);
 
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUrel, p) == fvc::div(phi)
+            fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -22,7 +21,7 @@
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 
@@ -32,7 +31,7 @@
     p.relax();
 
     // Momentum corrector
-    Urel -= rAUrel*fvc::grad(p);
+    Urel = HbyA - rAUrel*fvc::grad(p);
     Urel.correctBoundaryConditions();
     sources.correct(Urel);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 78dd40500b..adc0c0761a 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -1,19 +1,18 @@
 {
-    p.boundaryField().updateCoeffs();
-
     volScalarField rAU(1.0/UEqn().A());
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
-    adjustPhi(phi, U, p);
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    adjustPhi(phiHbyA, U, p);
 
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAU, p) == fvc::div(phi)
+            fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -22,7 +21,7 @@
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 
@@ -32,7 +31,7 @@
     p.relax();
 
     // Momentum corrector
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
index fbe81daee5..8193e072c5 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
@@ -1,15 +1,16 @@
+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();
-phi = fvc::interpolate(U) & mesh.Sf();
-adjustPhi(phi, U, p);
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+adjustPhi(phiHbyA, U, p);
 
 while (simple.correctNonOrthogonal())
 {
@@ -17,11 +18,11 @@ while (simple.correctNonOrthogonal())
 
     if (pressureImplicitPorosity)
     {
-        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
+        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phiHbyA));
     }
     else
     {
-        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
+        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phiHbyA));
     }
 
     tpEqn().setReference(pRefCell, pRefValue);
@@ -30,7 +31,7 @@ while (simple.correctNonOrthogonal())
 
     if (simple.finalNonOrthogonalIter())
     {
-        phi -= tpEqn().flux();
+        phi = phiHbyA - tpEqn().flux();
     }
 }
 
@@ -41,11 +42,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/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index fb8893be2c..6fa1965a24 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -13,13 +13,11 @@
 
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
-    volVectorField HbyA(rAU*UEqn.H());
+    volVectorField HbyA("HbyA", rAU*UEqn.H());
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phiv);
 
-    p.boundaryField().updateCoeffs();
-
     surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
 
     phiv -= phiGradp/rhof;