diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index b7b4b04db5c477cefc3dd2dfeab98d3ee6557a11..0f5ec2487b30854d2982e39196ea194bd9d32c4e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -4,7 +4,6 @@ tmp<fvVectorMatrix> UEqn
 (
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
-  - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U)
   + turbulence->divDevRhoReff(U)
 );
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
index 57d72887e60ba1f6dda559ff3353d1294ff3346d..a33c82b0acc5f7e14bd3ed318fc5b3fe8e56204d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
@@ -3,14 +3,10 @@
     (
         fvm::ddt(rho, h)
       + fvm::div(phi, h)
-      - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
         dpdt
-      - (
-            fvc::ddt(rho, K) + fvc::div(phi, K)
-          - (fvc::ddt(rho) + fvc::div(phi))*K
-        )
+      - fvc::ddt(rho, K) + fvc::div(phi, K)
     );
 
     hEqn.relax();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index f1bed4d071e985d92c9057c777ca9e216b174c54..21ec2646be5bb5702e0b2210335592e1fe124ec0 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
index 180a33d9186d3f450ed4bb06c8b4e79bb69ee41f..76adf004cbf2b99a2babf5116b12305786816b67 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
@@ -7,7 +7,6 @@
     fvScalarMatrix TEqn
     (
         fvm::div(phi, T)
-      - fvm::Sp(fvc::div(phi), T)
       - fvm::laplacian(kappaEff, T)
     );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
index e55b6d141d4c5f2d67d0f401bf91cd989f682532..10385aa91adf6a4836f056a9c180c5a4899e9f45 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
@@ -2,7 +2,6 @@
     fvScalarMatrix hEqn
     (
         fvm::div(phi, h)
-      - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turb.alphaEff(), h)
      ==
       - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index e6fc3d6f9a3b3de7559b97ec727c6d00c15c9ef0..d2ea510ced756fd512808be931430aa46287e3a6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -7,23 +7,31 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("Dp", 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 phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA += phig;
+
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi += phig;
-
     // Solve pressure
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference
@@ -37,14 +45,14 @@
         if (nonOrth == nNonOrthCorr)
         {
             // 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((phig - p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H
index 653389b221312286b60cbe000ab7489a28f4587e..fcc8b054a2e658615090961e7a94ae6628ddebd3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H
@@ -2,7 +2,6 @@
     fvScalarMatrix hPorousEqn
     (
         fvm::div(porousPhi, porousH)
-      - fvm::Sp(fvc::div(porousPhi), porousH)
       - fvm::laplacian(turbPorous.alphaEff(), porousH)
      ==
       - fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 80011d1ab593f367c75f6f0e31b37e9deeea9901..3197a72beeed1203b6248d6394e27dcdcf940269 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -5,6 +5,27 @@
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+    );
+    phi = phiHbyA;
+
+    surfaceScalarField phig
+    (
+        (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
 
@@ -27,27 +48,6 @@
           - fvc::Sp(fvc::div(phid2), p_rgh);
     }
 
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
-
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    );
-    phi = phiHbyA;
-
-    surfaceScalarField phig
-    (
-        (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
-    );
-
-    phiHbyA += phig;
-
     // Thermodynamic density needs to be updated by psi*d(p) after the
     // pressure solution - done in 2 parts. Part 1:
     //thermo.rho() -= psi*p_rgh;
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
index 861590b235bf2577f2a0161350263d1a1f90159c..8f38ca9d91fe910cfad67e53608612d277564b3a 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
@@ -6,7 +6,6 @@
     (
         fvm::ddt(alpha1, T1)
       + fvm::div(alphaPhi1, T1)
-      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1)
       - fvm::laplacian(kByCp1, T1)
      ==
         heatTransferCoeff*T2/Cp1/rho1
@@ -18,7 +17,6 @@
     (
         fvm::ddt(alpha2, T2)
       + fvm::div(alphaPhi2, T2)
-      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2)
       - fvm::laplacian(kByCp2, T2)
      ==
         heatTransferCoeff*T1/Cp2/rho2
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
index 5d61f0eaf786dbc2dcdd54e1e8c3b9831e3abe9e..606b4573e47988c921954a011626b7ce80ecfb4b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
@@ -30,7 +30,6 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         (
             fvm::ddt(alpha1, U1)
           + fvm::div(alphaPhi1, U1)
-          - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
 
           + Cvm*rho2*alpha1*alpha2/rho1*
             (
@@ -61,7 +60,6 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         (
             fvm::ddt(alpha2, U2)
           + fvm::div(alphaPhi2, U2)
-          - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
 
           + Cvm*rho2*alpha1*alpha2/rho2*
             (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index 062c5523c9d4c94a469e49e202b65e25ad1b285d..7cc250a66a2ab15a22680ab352f321e731b17c17 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -2,7 +2,6 @@
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
       + turbulence->divDevRhoReff(rho, U)
     );
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H
index 26015141d27b79c84bd2d5cb4706ccd79101dbc2..427a5f0b4cdcc175070bbdafd58b59db3ebe3b6e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H
@@ -6,7 +6,6 @@
     (
         fvm::ddt(alpha1, T1)
       + fvm::div(alphaPhi1, T1)
-      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1)
       - fvm::laplacian(kByCp1, T1)
      ==
         heatTransferCoeff*T2/Cp1/rho1
@@ -18,7 +17,6 @@
     (
         fvm::ddt(alpha2, T2)
       + fvm::div(alphaPhi2, T2)
-      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2)
       - fvm::laplacian(kByCp2, T2)
      ==
         heatTransferCoeff*T1/Cp2/rho2
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
index d4b371915ca4b61c1aae8baf9fc235bb6a2a0ab1..e6478e4d89014fbf5133f58e1b7307525ed39cfd 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
@@ -17,7 +17,6 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
         (
             fvm::ddt(alpha, U)
           + fvm::div(phase.phiAlpha(), U)
-          - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U)
 
           + (alpha/phase.rho())*fluid.Cvm(phase)*
             (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H
similarity index 100%
rename from applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H
rename to applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H
similarity index 100%
rename from applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H
rename to applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H
diff --git a/src/turbulenceModels/compressible/RAS/LRR/LRR.C b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
index 154df3132f03d42b0fb63ea5356bdae86638ca4f..175cb8b463970bb3c76e5e3634f8b3e1155c259d 100644
--- a/src/turbulenceModels/compressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
@@ -351,7 +351,6 @@ void LRR::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
@@ -394,7 +393,6 @@ void LRR::correct()
     (
         fvm::ddt(rho_, R_)
       + fvm::div(phi_, R_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clrr1_*rho_*epsilon_/k_, R_)
diff --git a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index 5f852ab59149da9f267325fe59eabccd33f879cd..b4e906cfb474399c45e4d8c58250ded7ae42a2ed 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -389,7 +389,6 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
@@ -433,7 +432,6 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(rho_, R_)
       + fvm::div(phi_, R_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clg1_*rho_*epsilon_/k_, R_)
diff --git a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 9a2c72b4c95b00dda89d88f3548af9c3d2b3cd5e..4e246491e499f881d0c7686a1a3872a0bb9b3be3 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -314,7 +314,6 @@ void LaunderSharmaKE::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_
@@ -334,7 +333,6 @@ void LaunderSharmaKE::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
diff --git a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index 8302d28df1c27ef1c64f122e81da10b81a91fc32..2048e4d2d0ccaf86a6f747d38b2f1d8e4a0d3b24 100644
--- a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -321,7 +321,6 @@ void RNGkEpsilon::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
       ==
         (C1_ - R)*G*epsilon_/k_
@@ -343,7 +342,6 @@ void RNGkEpsilon::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
       ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
diff --git a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 5a6bf8e8b90ee168c7a30cc81df4c2a77cc85a30..6c5f3aae0c8d7a2396497b5406b77e4fd948772b 100644
--- a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -416,7 +416,6 @@ void SpalartAllmaras::correct()
     (
         fvm::ddt(rho_, nuTilda_)
       + fvm::div(phi_, nuTilda_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), nuTilda_)
       - fvm::laplacian(DnuTildaEff(), nuTilda_)
       - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
      ==
diff --git a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
index 9d186eb039b6eb9d6f9b67437f81817be18e4488..3aa8ee2ec1e7e78a5ee9284dd2509ec21d35c4b6 100644
--- a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
@@ -292,7 +292,6 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      //***HGW - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_
@@ -314,7 +313,6 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      //***HGW - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index e84be63ad6cd86e5ef1de6f37c36d77d4c778b37..229fb3e291314e419cafcb8c627c4a6c712d3b7a 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
@@ -411,7 +411,6 @@ void kOmegaSST::correct()
     (
         fvm::ddt(rho_, omega_)
       + fvm::div(phi_, omega_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
      ==
         rhoGammaF1*GbyMu
@@ -436,7 +435,6 @@ void kOmegaSST::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(F1), k_)
      ==
         min(G, (c1_*betaStar_)*rho_*k_*omega_)
diff --git a/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
index 6206fda9399a8a458512ef9fff00b4eb621a7d44..f248d98138c5c2675d52df2e6c4992fa5ed81726 100644
--- a/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
@@ -331,7 +331,6 @@ void realizableKE::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1*rho_*magS*epsilon_
@@ -356,7 +355,6 @@ void realizableKE::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index 12a6f1b3d111b051302587ca94d4709c24f0077a..2b54d0ef09223a0942b96f6345f4211e18ce13e8 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -361,7 +361,6 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
         (
             fvm::ddt(k_)
           + fvm::div(phi(), k_)
-          - fvm::Sp(fvc::div(phi()), k_)
           - fvm::laplacian(DkEff(F1), k_)
         ==
             min(G, c1_*betaStar_*k_*omega_)
@@ -385,7 +384,6 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
         (
             fvm::ddt(omega_)
           + fvm::div(phi(), omega_)
-          - fvm::Sp(fvc::div(phi()), omega_)
           - fvm::laplacian(DomegaEff(F1), omega_)
         ==
             gamma(F1)*S2
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
index c2c7e5efbe456224c622afbfb168183dfb430dae..14683732b21010fe700b33266d23eb2275706608 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
@@ -349,7 +349,6 @@ void LRR::correct()
     (
         fvm::ddt(epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*(K/epsilon_)*R, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
       ==
@@ -392,7 +391,6 @@ void LRR::correct()
     (
         fvm::ddt(R_)
       + fvm::div(phi_, R_)
-      - fvm::Sp(fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clrr1_*epsilon_/k_, R_)
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index 49326b1ff50c27f6453223ec8f489c590b13cce3..1231be6d5df6f3e1739aa733b23c1dfcdedc56b8 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -396,7 +396,6 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*(k_/epsilon_)*R_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
@@ -440,7 +439,6 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(R_)
       + fvm::div(phi_, R_)
-      - fvm::Sp(fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clg1_*epsilon_/k_, R_)
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 4105cb1876e1664fd74006e7193ccd895f01b951..76046a80e251195b106f9acbdad79b1ee803ece7 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -408,7 +408,6 @@ void SpalartAllmaras::correct()
     (
         fvm::ddt(nuTilda_)
       + fvm::div(phi_, nuTilda_)
-      - fvm::Sp(fvc::div(phi_), nuTilda_)
       - fvm::laplacian(DnuTildaEff(), nuTilda_)
       - Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_))
      ==
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
index 029b9f9db790f4a44b555f4e01182741e16c3fe4..673a8f719f92fd4fd0a9bf76de366639faf09fbe 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
@@ -245,7 +245,6 @@ void kEpsilon::correct()
     (
         fvm::ddt(epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_
@@ -265,7 +264,6 @@ void kEpsilon::correct()
     (
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
index 51e1fbd491963020544af6770dc5c282d574cee8..052e9629c710a25b7d5d2702caed41f397583cc5 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
@@ -254,7 +254,6 @@ void kOmega::correct()
     (
         fvm::ddt(omega_)
       + fvm::div(phi_, omega_)
-      - fvm::Sp(fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(), omega_)
      ==
         alpha_*G*omega_/k_
@@ -274,7 +273,6 @@ void kOmega::correct()
     (
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index e6c49db1c86724808fed61940c3d9ac16246383b..9344ca0d7e3459a31134b041bdd58dac9ab0334b 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -381,7 +381,6 @@ void kOmegaSST::correct()
     (
         fvm::ddt(omega_)
       + fvm::div(phi_, omega_)
-      - fvm::Sp(fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
      ==
         gamma(F1)*S2
@@ -405,7 +404,6 @@ void kOmegaSST::correct()
     (
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(F1), k_)
      ==
         min(G, c1_*betaStar_*k_*omega_)
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C.new b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C.new
index 6e3d815d7ea9c99d7bad988714fd56418db88fbe..5ca2c2be4fe31f5588759550ba54bd58793c3ee4 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C.new
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C.new
@@ -374,10 +374,8 @@ void kOmegaSST::correct()
     (
         fvm::ddt(omega_)
       + fvm::div(phi_, omega_)
-      - fvm::Sp(fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
       + fvm::div(CDkPhiOmega, omega_)
-      - fvm::Sp(fvc::div(CDkPhiOmega), omega_)
      ==
         gamma(F1)*2*S2
       - fvm::Sp(beta(F1)*omega_, omega_)
@@ -395,7 +393,6 @@ void kOmegaSST::correct()
     (
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(F1), k_)
      ==
         min(G, c1_*betaStar_*k_*omega_)
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
index 45f4f69d1302091cc5bec564af3c8d1f3680b97a..2aa10e84c682852b5138b363132ae9c5b07a85fd 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
@@ -716,7 +716,6 @@ void kkLOmega::correct()
     (
         fvm::ddt(kt_)
       + fvm::div(phi_, kt_)
-      - fvm::Sp(fvc::div(phi_), kt_)
       - fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
      ==
         Pkt
@@ -739,7 +738,6 @@ void kkLOmega::correct()
     (
         fvm::ddt(kl_)
       + fvm::div(phi_, kl_)
-      - fvm::Sp(fvc::div(phi_), kl_)
       - fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
      ==
         Pkl
@@ -761,7 +759,6 @@ void kkLOmega::correct()
     (
         fvm::ddt(omega_)
       + fvm::div(phi_, omega_)
-      - fvm::Sp(fvc::div(phi_), omega_)
       - fvm::laplacian
         (
             DomegaEff(alphaTEff),
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
index 495b8930dfc6b0df2b18c350e3725af06a65fe0d..5ea6831ea212813f2b31738384e13792df20c0b5 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
@@ -308,7 +308,6 @@ void realizableKE::correct()
     (
         fvm::ddt(epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1*magS*epsilon_
@@ -332,7 +331,6 @@ void realizableKE::correct()
     (
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::Sp(epsilon_/k_, k_)
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
index 0a14cc31a562d3cf4ef18e5f212213f57ba396d9..e106336cb2dd0a748f8d4dc499995db3c8067024 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -27,12 +27,12 @@ gradSchemes
 
 divSchemes
 {
-    div(phi,U)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phi,e)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,Ekp)    Gauss upwind;
+    div(phi,e)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,Ekp)    bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
index 0a14cc31a562d3cf4ef18e5f212213f57ba396d9..e106336cb2dd0a748f8d4dc499995db3c8067024 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -27,12 +27,12 @@ gradSchemes
 
 divSchemes
 {
-    div(phi,U)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phi,e)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,Ekp)    Gauss upwind;
+    div(phi,e)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,Ekp)    bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
index e1e47211d72678553a22afbf86cea2a1328511c2..0cab2f77d3ac1d8f4f224427ea6d8162e1f3c398 100644
--- a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
@@ -29,14 +29,14 @@ divSchemes
 {
     default             none;
 
-    div(phi,U)          Gauss upwind;
+    div(phi,U)          bounded Gauss upwind;
     div((muEff*dev2(T(grad(U)))))      Gauss linear;
-    div(phi,e)          Gauss upwind;
-    div(phi,epsilon)    Gauss upwind;
-    div(phi,k)          Gauss upwind;
+    div(phi,e)          bounded Gauss upwind;
+    div(phi,epsilon)    bounded Gauss upwind;
+    div(phi,k)          bounded Gauss upwind;
 
     div(phid,p)         Gauss upwind;
-    div(phi,Ekp)        Gauss upwind;
+    div(phi,Ekp)        bounded Gauss upwind;
     div((phi|interpolate(rho)),p)  Gauss upwind;
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
index 6fc42b1ca74d8432bc3222d4dab79fc56b92d787..4db877838c9762d1f878f419fc5aafd5d4e6c26b 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
@@ -28,10 +28,10 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,T)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,T)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
index b36bd5f8fcbd2378f485eef05e5ff10ec0f5aa8b..7a53eb9600ae3fc8926ba2dccddf2883fc109fc5 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
@@ -28,11 +28,11 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,T)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,T)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
index 09a66a2658baa460e70388dc986fc78c6a2381a7..2a65c2d5f05717f0c05c11ed67319e2bedc35748 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
@@ -27,13 +27,13 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,h)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
index 9b8eac50b508961960c847c3930d658aeaf483df..1f8388eeb23a981e52db1780e53bd3c9c62dfc78 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
@@ -27,13 +27,13 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,h)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/fvSchemes
index ca25b1d27ca7656c874002aa7e936abb3fdba428..02d978723c0fafff9f742bcec1a80b86c32ed6b1 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/fvSchemes
@@ -27,13 +27,13 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,h)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
     div(Ji,Ii_h)    Gauss linearUpwind grad(U);
     div((muEff*dev2(T(grad(U))))) Gauss linear;
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/fvSchemes
index ca25b1d27ca7656c874002aa7e936abb3fdba428..02d978723c0fafff9f742bcec1a80b86c32ed6b1 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/fvSchemes
@@ -27,13 +27,13 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,h)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,K)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
     div(Ji,Ii_h)    Gauss linearUpwind grad(U);
     div((muEff*dev2(T(grad(U))))) Gauss linear;
diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes
index d78114f231040e2fe6b409000a5cad91a599e5bb..d4699567a4007dba4e44ca3d318547cff6deaf25 100644
--- a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes
@@ -30,9 +30,9 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinearV 1;
-    div(phi,k)      Gauss limitedLinear 1;
-    div(phi,epsilon) Gauss limitedLinear 1;
+    div(phi,U)      bounded Gauss limitedLinearV 1;
+    div(phi,k)      bounded Gauss limitedLinear 1;
+    div(phi,epsilon) bounded Gauss limitedLinear 1;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSchemes b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSchemes
index bd069db99aa95bbd36f76f1c6b161830660bad3d..db947bf4a67f0327cc0727aa93622a891128aeab 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSchemes
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSchemes
@@ -30,13 +30,13 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,Urel)   Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,omega)  Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,Urel)   bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,omega)  bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
-    div(phi,nuTilda) Gauss upwind;
+    div(phi,nuTilda) bounded Gauss upwind;
     div((nuEff*dev(T(grad(Urel))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes b/tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes
index ffda45b2c7cd530059ffa7516551c5dd40f9695d..d06e77f4d583d16f16ddd54e4df7e1446fe11988 100644
--- a/tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes
+++ b/tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes
@@ -32,12 +32,12 @@ divSchemes
 {
     default         none;
 
-    div(phi,U)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 
-    div(-phi,Ua)    Gauss upwind;
+    div(-phi,Ua)    bounded Gauss upwind;
     div((nuEff*dev(T(grad(Ua))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes
index ae53be662314c563f218adb483ee869c9f4b5537..14ac340c9bfcaf22bad1e12c9bb13ff9248a88f3 100644
--- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes
+++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes
@@ -29,9 +29,9 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linearUpwind grad(U);
-    div(phi,k)      Gauss upwind;
-    div(phi,omega)  Gauss upwind;
+    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,omega)  bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/decomposeParDict b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/decomposeParDict
index a630750b7dcb0b04f91f757072ef8c91d2985e11..fd59ef09a0868db7cd5e8d515d5f1bf382708ed0 100644
--- a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/decomposeParDict
+++ b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/decomposeParDict
@@ -17,7 +17,7 @@ FoamFile
 
 numberOfSubdomains 8;
 
-method          hierarchical; //ptscotch;
+method          hierarchical;
 
 simpleCoeffs
 {
diff --git a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/fvSchemes b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/fvSchemes
index 9050cc767a2b1d7e49c827109723edea508b8060..b17cff7ff37cbc796c88ecf5d348c4bf21eabf6b 100644
--- a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/fvSchemes
+++ b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/fvSchemes
@@ -30,11 +30,11 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linearUpwindV grad(U);
-    div(phi,k)      Gauss upwind;
-    div(phi,omega)  Gauss upwind;
+    div(phi,U)      bounded Gauss linearUpwindV grad(U);
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,omega)  bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
-    div(phi,nuTilda) Gauss upwind;
+    div(phi,nuTilda) bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
index a3d6a26e3e492dc8ad78446e42471f88cd8defef..fbd760de1fb64e6a56ca0b2b9a1f97f1f6a737a4 100644
--- a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -29,10 +29,10 @@ gradSchemes
 
 divSchemes
 {
-    div(phi,U)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
index a3d6a26e3e492dc8ad78446e42471f88cd8defef..fbd760de1fb64e6a56ca0b2b9a1f97f1f6a737a4 100644
--- a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -29,10 +29,10 @@ gradSchemes
 
 divSchemes
 {
-    div(phi,U)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/incompressible/simpleFoam/airFoil2D/system/fvSchemes b/tutorials/incompressible/simpleFoam/airFoil2D/system/fvSchemes
index 3f602ada18f0e18c9a74c80f5744bb056f55dffc..bd2f1c8952d59687b33c496bccf2f7b385e91bc8 100644
--- a/tutorials/incompressible/simpleFoam/airFoil2D/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/airFoil2D/system/fvSchemes
@@ -30,8 +30,8 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linearUpwind grad(U);
-    div(phi,nuTilda) Gauss linearUpwind grad(nuTilda);
+    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
index d23b2aca63f7287a78ee055f96f75bdc0e7a0a45..36a6021c97a5b210efcb5206792ea644ef62b423 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
@@ -27,9 +27,9 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linearUpwindV grad(U);
-    div(phi,k)      Gauss upwind;
-    div(phi,omega)  Gauss upwind;
+    div(phi,U)      bounded Gauss linearUpwindV grad(U);
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,omega)  bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSchemes b/tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSchemes
index c1af3f96495a78e314042566f18171c71f5523d6..51aa595facf2d12988dfdb4e53522cb150de2ee4 100644
--- a/tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSchemes
@@ -30,12 +30,12 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinearV 1;
-    div(phi,k)      Gauss limitedLinear 1;
-    div(phi,epsilon) Gauss limitedLinear 1;
-    div(phi,R)      Gauss limitedLinear 1;
+    div(phi,U)      bounded Gauss limitedLinearV 1;
+    div(phi,k)      bounded Gauss limitedLinear 1;
+    div(phi,epsilon) bounded Gauss limitedLinear 1;
+    div(phi,R)      bounded Gauss limitedLinear 1;
     div(R)          Gauss linear;
-    div(phi,nuTilda) Gauss limitedLinear 1;
+    div(phi,nuTilda) bounded Gauss limitedLinear 1;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
index bffb2ed43b19db2eaf8d7062202cbdc1d3fa501a..540ac90b64c61b981cea3dc2b7fc10f6a8293ebc 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
@@ -30,12 +30,12 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
-    div(phi,nuTilda) Gauss upwind;
+    div(phi,nuTilda) bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes
index bffb2ed43b19db2eaf8d7062202cbdc1d3fa501a..540ac90b64c61b981cea3dc2b7fc10f6a8293ebc 100644
--- a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes
@@ -30,12 +30,12 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
     div(R)          Gauss linear;
-    div(phi,nuTilda) Gauss upwind;
+    div(phi,nuTilda) bounded Gauss upwind;
     div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
index 4c101301a10b3092838b386bbbf2b4441310c6c3..774f507b640eca0d324b6804a502629ad98583eb 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
@@ -27,10 +27,10 @@ gradSchemes
 divSchemes
 {
     default             none;
-    div(phi,U)          Gauss upwind;
+    div(phi,U)          bounded Gauss upwind;
     div((nuEff*dev(T(grad(U)))))    Gauss linear;
-    div(phi,epsilon)    Gauss upwind;
-    div(phi,k)          Gauss upwind;
+    div(phi,epsilon)    bounded Gauss upwind;
+    div(phi,k)          bounded Gauss upwind;
 }
 
 laplacianSchemes
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
index e8088bf7bdd9fab6608cd4a7c43c293bf0d7907d..1319b623261803836c97f2fc155849db5acaa923 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
@@ -20,12 +20,14 @@ FoamFile
     back
     {
         type            symmetryPlane;
+        inGroups        1(symmetryPlane);
         nFaces          9340;
         startFace       265900;
     }
     front
     {
         type            symmetryPlane;
+        inGroups        1(symmetryPlane);
         nFaces          9340;
         startFace       275240;
     }
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/boundary b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/boundary
index 9eb64e86375de42fb8dda163fccc2cf75d8dd9a5..81845d8cb3f4398b699c70a886266f48558c355e 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/boundary
@@ -26,6 +26,7 @@ FoamFile
     frontAndBack
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          25600;
         startFace       25840;
     }
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary
index 067ff53cc9865662f18fb8c2f2fb21a1e3975db0..f1bde7bdce369c1bfc963511de82c352d2772aab 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary
@@ -26,12 +26,14 @@ FoamFile
     front
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          1360;
         startFace       2794;
     }
     back
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          1360;
         startFace       4154;
     }
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/controlDict b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/controlDict
index c939e692fda66968c4d616fcd9eaa83ca87ac8eb..b9a58d7cc2eca79f587a7cf1619b66e952ade567 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/controlDict
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/controlDict
@@ -80,6 +80,7 @@ functions
         (
             p
         );
+        interpolationScheme cellPoint;
 
         surfaces
         (
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/constant/polyMesh/boundary b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/constant/polyMesh/boundary
index 067ff53cc9865662f18fb8c2f2fb21a1e3975db0..f1bde7bdce369c1bfc963511de82c352d2772aab 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/constant/polyMesh/boundary
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/constant/polyMesh/boundary
@@ -26,12 +26,14 @@ FoamFile
     front
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          1360;
         startFace       2794;
     }
     back
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          1360;
         startFace       4154;
     }
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/controlDict b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/controlDict
index a63975e3fb781f9b1699a6b40c83494f87b974c3..252b941da73002ccba779271ab5b8c3c7faaa68f 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/controlDict
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/controlDict
@@ -79,6 +79,7 @@ functions
         (
             p
         );
+        interpolationScheme cellPoint;
 
         surfaces
         (