diff --git a/applications/solvers/multiphase/cavitatingFoam/CourantNo.H b/applications/solvers/multiphase/cavitatingFoam/CourantNo.H
index 02a6c86db7..d42f6344ce 100644
--- a/applications/solvers/multiphase/cavitatingFoam/CourantNo.H
+++ b/applications/solvers/multiphase/cavitatingFoam/CourantNo.H
@@ -37,7 +37,7 @@ if (mesh.nInternalFaces())
 {
     scalarField sumPhi
     (
-        fvc::surfaceSum(mag(phiv))().internalField()
+        fvc::surfaceSum(mag(phi))().internalField()
     );
 
     CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
@@ -54,7 +54,7 @@ if (mesh.nInternalFaces())
     )*runTime.deltaTValue();
 }
 
-Info<< "phiv Courant Number mean: " << meanCoNum
+Info<< "phi Courant Number mean: " << meanCoNum
     << " max: " << CoNum
     << " acoustic max: " << acousticCoNum
     << endl;
diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 1ca0b9f0a5..8edf0d3934 100644
--- a/applications/solvers/multiphase/cavitatingFoam/UEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
@@ -1,7 +1,7 @@
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
-      + fvm::div(phi, U)
+      + fvm::div(rhoPhi, U)
       + turbulence->divDevRhoReff(rho, U)
     );
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
index c3fde6b3a1..863054c59f 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createUf.H"
     #include "createPcorrTypes.H"
-    #include "compressibleCourantNo.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -86,12 +86,12 @@ int main(int argc, char *argv[])
             if (correctPhi)
             {
                 // Calculate absolute flux from the mapped surface velocity
-                phiv = mesh.Sf() & Uf;
+                phi = mesh.Sf() & Uf;
 
                 #include "correctPhi.H"
 
                 // Make the flux relative to the mesh motion
-                fvc::makeRelative(phiv, U);
+                fvc::makeRelative(phi, U);
             }
         }
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
index c05d664838..4bb1ce9d0d 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
@@ -14,7 +14,7 @@ if (mesh.changing())
         {
             U.boundaryField()[patchI].evaluate();
 
-            phiv.boundaryField()[patchI] =
+            phi.boundaryField()[patchI] =
                 U.boundaryField()[patchI]
               & mesh.Sf().boundaryField()[patchI];
         }
@@ -37,21 +37,21 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    surfaceScalarField rhof(fvc::interpolate(rho, "div(phiv,rho)"));
+    surfaceScalarField rhof(fvc::interpolate(rho, "div(phi,rho)"));
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::ddt(rho) + fvc::div(phiv*rhof)
+            fvm::laplacian(rAUf, pcorr) == fvc::ddt(rho) + fvc::div(phi*rhof)
         );
 
         pcorrEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv -= pcorrEqn.flux()/rhof;
+            phi -= pcorrEqn.flux()/rhof;
         }
     }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index a0e0e6c74d..b07846b6de 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -17,13 +17,13 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+    phi = (fvc::interpolate(HbyA) & mesh.Sf())
          + rhorAUf*fvc::ddtCorr(U, Uf);
-    fvc::makeRelative(phiv, U);
+    fvc::makeRelative(phi, U);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
-    phiv -= phiGradp/rhof;
+    phi -= phiGradp/rhof;
 
     volScalarField rho0(rho - psi*p);
 
@@ -33,7 +33,7 @@
         (
             fvc::ddt(rho)
           + psi*correction(fvm::ddt(p))
-          + fvc::div(phiv, rho)
+          + fvc::div(phi, rho)
           + fvc::div(phiGradp)
           - fvm::laplacian(rhorAUf, p)
         );
@@ -42,7 +42,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv += (phiGradp + pEqn.flux())/rhof;
+            phi += (phiGradp + pEqn.flux())/rhof;
         }
     }
 
@@ -85,6 +85,6 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += n*(phiv/mesh.magSf() - (n & Uf));
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index c4453df391..4ba2402c01 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
     #include "readControls.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
-    #include "compressibleCourantNo.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     pimpleControl pimple(mesh);
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index 605c0deece..2c834dc511 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -39,12 +39,25 @@
         mesh
     );
 
-    #include "createPhiv.H"
-    #include "compressibleCreatePhi.H"
+    #include "createPhi.H"
+
+    // Mass flux (corrected by rhoEqn.H)
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
 
     Info<< "Reading transportProperties\n" << endl;
 
-    incompressibleTwoPhaseMixture twoPhaseProperties(U, phiv);
+    incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alphav(twoPhaseProperties.alpha1());
     alphav.oldTime();
@@ -72,5 +85,5 @@
     // Create incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::turbulenceModel::New(U, phiv, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 210a30fd14..da8047914c 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -17,12 +17,12 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + rhorAUf*fvc::ddtCorr(U, phiv);
+    phi = (fvc::interpolate(HbyA) & mesh.Sf())
+         + rhorAUf*fvc::ddtCorr(U, phi);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
-    phiv -= phiGradp/rhof;
+    phi -= phiGradp/rhof;
 
     while (pimple.correctNonOrthogonal())
     {
@@ -30,7 +30,7 @@
         (
             fvm::ddt(psi, p)
           - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
-          + fvc::div(phiv, rho)
+          + fvc::div(phi, rho)
           + fvc::div(phiGradp)
           - fvm::laplacian(rhorAUf, p)
         );
@@ -39,7 +39,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv += (phiGradp + pEqn.flux())/rhof;
+            phi += (phiGradp + pEqn.flux())/rhof;
         }
     }
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H b/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
index d0bd6e1dad..a01823ad49 100644
--- a/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
@@ -2,12 +2,12 @@
     fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
-      + fvm::div(phiv, rho)
+      + fvm::div(phi, rho)
     );
 
     rhoEqn.solve();
 
-    phi = rhoEqn.flux();
+    rhoPhi = rhoEqn.flux();
 
     Info<< "max-min rho: " << max(rho).value()
         << " " << min(rho).value() << endl;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/0/p b/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
index f484f1a42e..22afb75d57 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
index 6b7b46498b..17346054f6 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
@@ -28,9 +28,9 @@ interpolationSchemes
 divSchemes
 {
     default         none;
-    div(phiv,rho)   Gauss vanLeer;
-    div(phi,U)      Gauss LUST grad(U);
-    div(phiv,k)     Gauss LUST grad(k);
+    div(phi,rho)    Gauss vanLeer;
+    div(rhoPhi,U)   Gauss LUST grad(U);
+    div(phi,k)      Gauss LUST grad(k);
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
index 1785af6f4e..ca2022df7b 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
@@ -78,5 +78,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p b/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
index 58cc484ff7..a85f4b9d0d 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
index f7e2958bf6..64aa6e54b9 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
@@ -29,9 +29,9 @@ divSchemes
 {
     default         none;
 
-    div(phiv,rho)   Gauss vanLeer;
-    div(phi,U)      Gauss LUST grad(U);
-    div(phiv,k)     Gauss LUST grad(k);
+    div(phi,rho)    Gauss vanLeer;
+    div(rhoPhi,U)   Gauss LUST grad(U);
+    div(phi,k)      Gauss LUST grad(k);
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
index 7133194be0..7e1a541e0a 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
@@ -75,5 +75,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p b/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
index f484f1a42e..22afb75d57 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
index 6f20de56af..d253239d92 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
@@ -28,10 +28,10 @@ interpolationSchemes
 divSchemes
 {
     default         none;
-    div(phiv,rho)   Gauss limitedLinear 1;
-    div(phi,U)      Gauss limitedLinearV 1;
-    div(phiv,omega) Gauss limitedLinear 1;
-    div(phiv,k)     Gauss limitedLinear 1;
+    div(phi,rho)    Gauss limitedLinear 1;
+    div(rhoPhi,U)   Gauss limitedLinearV 1;
+    div(phi,omega)  Gauss limitedLinear 1;
+    div(phi,k)      Gauss limitedLinear 1;
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
index 1785af6f4e..d2ebcca130 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
@@ -33,7 +33,6 @@ solvers
         relTol          0;
     }
 
-
     p
     {
         solver          GAMG;
@@ -78,5 +77,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
index bf1202d6ea..ca3c8a6104 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
@@ -74,9 +74,6 @@ PIMPLE
 
 relaxationFactors
 {
-    fields
-    {
-    }
     equations
     {
         ".*" 1;