Thermodynamics and sub-models: Removed "Sp" boundedness corrections on transport, replaced with "bounded Gauss" scheme

2025-11-28 03:28:01 +00:00 · 2012-09-21 14:34:42 +01:00
parent 1d83996c0f
commit a114345eab
62 changed files with 150 additions and 186 deletions
--- 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)
 );

--- 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();
--- 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)
    );

--- 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)
    );

--- 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)")
--- 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();
        }
    }
--- 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)")
--- 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;
--- 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
--- 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*
            (
--- 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)
    );

--- 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
--- 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)*
            (
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H
--- 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_)
--- 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_)
--- 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_)
--- 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_)
--- 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_))
     ==
--- 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
--- 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_)
--- 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_)
--- 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
--- 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_)
--- 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_)
--- 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_))
     ==
--- 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
--- 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
--- 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_)
--- 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_)
--- 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),
--- 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_)
--- 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
--- 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
--- 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;
 }

--- 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;
 }

--- 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;
 }
--- 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;
 }
--- 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;
 }
--- 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;
--- 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;
--- 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;
 }

--- 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;
 }

--- 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;
 }

--- 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;
 }

--- 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
 {
--- 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
--- 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
--- 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
--- 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;
 }

--- 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;
 }

--- 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;
 }

--- 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;
 }

--- 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;
 }

--- 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
--- 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;
    }
--- 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;
    }
--- 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;
    }
--- 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
        (
--- 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;
    }
--- 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
        (