diff --git a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..e4c79b40cc5685ddbeeab18c32e8611d745656b8
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
@@ -0,0 +1,19 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index d8fdff55ae8432751a153864d5b8ad9e7b76dfe0..dea35b76575477e953c521d77f6bb12cf28c3dc0 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& e = thermo.he();
     const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field U\n" << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/eEqn.H
deleted file mode 100644
index a1ea77157311f3f475957a725a237c7baa8fbed3..0000000000000000000000000000000000000000
--- a/applications/solvers/compressible/rhoSimpleFoam/eEqn.H
+++ /dev/null
@@ -1,18 +0,0 @@
-{
-    // Kinetic + pressure energy
-    volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
-
-    fvScalarMatrix eEqn
-    (
-        fvm::div(phi, e)
-      - fvm::Sp(fvc::div(phi), e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-        fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
-    );
-
-    eEqn.relax();
-    eEqn.solve();
-
-    thermo.correct();
-}
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..ff467c0382cfcb444f66d64cbbcb549c0dac1736
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
@@ -0,0 +1,21 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    pZones.addEnergySource(thermo, rho, EEqn);
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
index 96bdb2a122b0b346cb7f0602e90d0dfef6e2612b..4fff74d224cadb920dec4f6c42b936dff55f0452 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
@@ -5,6 +5,7 @@
         rhoThermo::New(mesh)
     );
     rhoThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& e = thermo.he();
 
     Info<< "Reading field U\n" << endl;
     volVectorField U
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H
deleted file mode 100644
index 5d0f174623b04cd20f8d643792d07dcf04ab3a25..0000000000000000000000000000000000000000
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H
+++ /dev/null
@@ -1,20 +0,0 @@
-{
-    // Kinetic + pressure energy
-    volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
-
-    fvScalarMatrix eEqn
-    (
-        fvm::div(phi, e)
-      - fvm::Sp(fvc::div(phi), e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-        fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
-    );
-
-    pZones.addEnergySource(thermo, rho, eEqn);
-
-    eEqn.relax();
-    eEqn.solve();
-
-    thermo.correct();
-}
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 22a58da1cc85f2fb159e7bcd7874ab34b6f6786a..3898295b915061669daa94fa773218a5877ef357 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 0eee9129b568412c8e00c76671748d479c54ca8b..312196583eaae9a90ce43b22350f2ecb5d4ac383 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
index 561c5e2dc95685ef3b061f1dfb8c37aa817fab1d..f707570f60cd0e0a73162bd7ebabb4857b2b325e 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
         // Velocity-pressure-enthalpy SIMPLEC corrector
         {
             #include "UEqn.H"
-            #include "eEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 
diff --git a/applications/solvers/compressible/sonicFoam/EEqn.H b/applications/solvers/compressible/sonicFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..15979381a93aeaad192e1422428a91e8cf268e1a
--- /dev/null
+++ b/applications/solvers/compressible/sonicFoam/EEqn.H
@@ -0,0 +1,10 @@
+{
+    solve
+    (
+        fvm::ddt(rho, e) + fvm::div(phi, e)
+      + fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho))
+      - fvm::laplacian(turbulence->alphaEff(), e)
+    );
+
+    thermo.correct();
+}
diff --git a/applications/solvers/compressible/sonicFoam/createFields.H b/applications/solvers/compressible/sonicFoam/createFields.H
index 3c39d3e7aad309ccefaec4c861d33decc49001a3..c3421369e7c87da7f28a15f25f812eed340e6519 100644
--- a/applications/solvers/compressible/sonicFoam/createFields.H
+++ b/applications/solvers/compressible/sonicFoam/createFields.H
@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "e");
 
     volScalarField& p = thermo.p();
     volScalarField& e = thermo.he();
diff --git a/applications/solvers/compressible/sonicFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/eEqn.H
deleted file mode 100644
index cfd908ded81ed56b906b12b111f94d3f77b5ab1a..0000000000000000000000000000000000000000
--- a/applications/solvers/compressible/sonicFoam/eEqn.H
+++ /dev/null
@@ -1,12 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, e)
-      + fvm::div(phi, e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-      - (fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho)))
-    );
-
-    thermo.correct();
-}
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 99ca26fce35e974a868a4a05a9376ba1cdb7f07c..4d700b3a41ee5a84b9e05663d46fd9a70c8d07de 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -17,7 +17,8 @@ surfaceScalarField phid
 
 volScalarField Dp("Dp", rho*rAU);
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+// Non-orthogonal pressure corrector loop
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
@@ -28,7 +29,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve();
 
-    if (nonOrth == nNonOrthCorr)
+    if (pimple.finalNonOrthogonalIter())
     {
         phi = pEqn.flux();
     }
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..994c2862b8fbffd162e1b94b4ea2676403d2515d
--- /dev/null
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H
@@ -0,0 +1,11 @@
+{
+    solve
+    (
+        fvm::ddt(rho, e) + fvm::div(phi, e)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + fvc::div(phi/fvc::interpolate(rho) + mesh.phi(), p, "div(phiv,p)")
+      - fvm::laplacian(turbulence->alphaEff(), e)
+    );
+
+    thermo.correct();
+}
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H
deleted file mode 100644
index 0c17353ec45ee6026ab67e88b1ff3c882d8d3936..0000000000000000000000000000000000000000
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H
+++ /dev/null
@@ -1,12 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, e)
-      + fvm::div(phi, e)
-      - fvm::laplacian(turbulence->alphaEff(), e)
-     ==
-      - p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi())
-    );
-
-    thermo.correct();
-}
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
index 163baa5b69a537442b503e3485c7448a972088ce..2f7ac9962a0e6c485310268d11d029ebb0608b6d 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "motionSolver.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -62,20 +65,24 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
-
-        #include "eEqn.H"
-
-
-        // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         rho = thermo.rho();
 
         runTime.write();
diff --git a/applications/solvers/compressible/sonicFoam/sonicFoam.C b/applications/solvers/compressible/sonicFoam/sonicFoam.C
index 954404b6a47f225dd398d5935ba1b969931f1048..71d032d4223486a92ecdde1acdcd0a533fb5397a 100644
--- a/applications/solvers/compressible/sonicFoam/sonicFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicFoam.C
@@ -33,6 +33,7 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,21 +55,29 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            #include "eEqn.H"
-            #include "pEqn.H"
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         rho = thermo.rho();
 
         runTime.write();
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index 19a530a6aaa52c9d41198f4d3aa6a42346ea9af6..1ac1ae4ae389cf45612b9364fb783212f6bc7f32 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,57 +55,59 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
 
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(rho, U)
-          + fvm::div(phi, U)
-          - fvm::laplacian(mu, U)
-        );
-
-        solve(UEqn == -fvc::grad(p));
-
-
-        // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            volScalarField rAU(1.0/UEqn.A());
-            U = rAU*UEqn.H();
-
-            surfaceScalarField phid
+            fvVectorMatrix UEqn
             (
-                "phid",
-                psi
-               *(
-                    (fvc::interpolate(U) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, rho, U, phi)
-                )
+                fvm::ddt(rho, U)
+              + fvm::div(phi, U)
+              - fvm::laplacian(mu, U)
             );
 
-            phi = (rhoO/psi)*phid;
-            volScalarField Dp("Dp", rho*rAU);
+            solve(UEqn == -fvc::grad(p));
 
-            fvScalarMatrix pEqn
-            (
-                fvm::ddt(psi, p)
-              + fvc::div(phi)
-              + fvm::div(phid, p)
-              - fvm::laplacian(Dp, p)
-            );
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                volScalarField rAU(1.0/UEqn.A());
+                U = rAU*UEqn.H();
+
+                surfaceScalarField phid
+                (
+                    "phid",
+                    psi
+                   *(
+                       (fvc::interpolate(U) & mesh.Sf())
+                     + fvc::ddtPhiCorr(rAU, rho, U, phi)
+                   )
+                );
+
+                phi = (rhoO/psi)*phid;
+                volScalarField Dp("Dp", rho*rAU);
+
+                fvScalarMatrix pEqn
+                (
+                    fvm::ddt(psi, p)
+                  + fvc::div(phi)
+                  + fvm::div(phid, p)
+                  - fvm::laplacian(Dp, p)
+                );
 
-            pEqn.solve();
+                pEqn.solve();
 
-            phi += pEqn.flux();
+                phi += pEqn.flux();
 
-            #include "compressibleContinuityErrs.H"
+                #include "compressibleContinuityErrs.H"
 
-            U -= rAU*fvc::grad(p);
-            U.correctBoundaryConditions();
+                U -= rAU*fvc::grad(p);
+                U.correctBoundaryConditions();
+            }
         }
 
         rho = rhoO + psi*p;
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3ad3c0f03e3dd9307f2dde9ed5b3689c16ea73b4
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H
@@ -0,0 +1,20 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + fvm::div(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 07ec2300c0b7b8ab4f8a7f19074bc4065bd7150a..dfa88b59a9aebb32486807f80705b1381925946c 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index f2edc4cdba0bf5d47170ef1f06d36c7666dffb3b..f6b183e3fcbf95c113c1a81d7da789c70478ecaa 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -5,6 +5,7 @@
         rhoThermo::New(mesh)
     );
     rhoThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.he();
     const volScalarField& psi = thermo.psi();
 
 
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
deleted file mode 100644
index b1d9e4e8b0237ca0b9f28f6d36c8f5162121c846..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
+++ /dev/null
@@ -1,16 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-    );
-
-    hEqn.relax();
-    hEqn.solve();
-
-    thermo.correct();
-}
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake b/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake
new file mode 100755
index 0000000000000000000000000000000000000000..0fe8e8f4ad18ed67c9f57ed5c360bc03d39f69cf
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake buoyantSimpleRadiationFoam
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..e4c79b40cc5685ddbeeab18c32e8611d745656b8
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H
@@ -0,0 +1,19 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+}
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 2a355a21de09b87745cbacfddfa7d4c546bdd8bf..4a9387c1882b17967d218a95e68553f29d52b64c 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..0d1f41d0c30326a6647e51a0be4cc958921d69b8
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H
@@ -0,0 +1,22 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        radiation->Sh(thermo)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+    radiation->correct();
+}
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files
similarity index 98%
rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files
rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files
index 79ebcd90aa566bb9b87eee17b45aef18df998a97..a8347c552591472f0adfd57571aebeb0d2cb440e 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files
@@ -1,4 +1,3 @@
 buoyantSimpleRadiationFoam.C
 
 EXE = $(FOAM_APPBIN)/buoyantSimpleRadiationFoam
-
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options
similarity index 94%
rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options
index f2ebc095ffc09ece9eca7b3713ef4ecbdc3dd407..f26046adb2a08affb14512022e8082ea6fc73b51 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options
@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I../buoyantSimpleFoam \
+    -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
similarity index 98%
rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index 07718aa41d8cfefbe9ccfcf2006a0d9151c20fbf..d02853591ebfa8ca324256c1a3a66bbb2baa75ab 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -62,7 +62,7 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
             #include "pEqn.H"
         }
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index d2fff48837c72ff56fcf3cc1aa08bb0a41323ad9..bb7a65cb1d28ac920674974f91563447637874e1 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -20,7 +21,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.he();
     const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field U\n" << endl;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
deleted file mode 100644
index 24ed135c08160d6d3d9cfa8bfb8f4396919bb35a..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
+++ /dev/null
@@ -1,15 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::div(phi, h)
-      - fvm::Sp(fvc::div(phi), h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
-    );
-
-    hEqn.relax();
-    hEqn.solve();
-
-    thermo.correct();
-}
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
deleted file mode 100644
index 2c422a5cab526e69b5f38b2450413c24390f7d50..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
+++ /dev/null
@@ -1,19 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::div(phi, h)
-      - fvm::Sp(fvc::div(phi), h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
-      + radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}
diff --git a/etc/controlDict b/etc/controlDict
index 2199d7ec7a917e167476d6e7eafd74d64fbe1539..caf1e22933c01ed70186292923b921efad22ff83 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -884,6 +884,7 @@ DebugSwitches
     wallHeatTransfer    0;
     wallLayerCells      0;
     wallModel           0;
+    warnUnboundedGauss  1;
     waveTransmissive    0;
     wedge               0;
     weighted            0;
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.C b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
index e3754f0a582dbcd70230669093e9dcb5a950f530..8234e573c560002f57968323173e71b11541252d 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
@@ -119,7 +119,7 @@ singleStepCombustion<CombThermoType, ThermoType>::R
     {
         const label fNorm = singleMixturePtr_->specieProd()[specieI];
         const volScalarField fres(singleMixturePtr_->fres(specieI));
-        wSpecie /= max(fNorm*(Y - fres), 1e-2);
+        wSpecie /= max(fNorm*(Y - fres), scalar(1e-2));
 
         return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
     }
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 22988f32169a94d03186bcb7328e24f85c13294e..17abbc6363da0a3e689afa77a991d73def2cc1e3 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -345,6 +345,7 @@ convectionSchemes = finiteVolume/convectionSchemes
 $(convectionSchemes)/convectionScheme/convectionSchemes.C
 $(convectionSchemes)/gaussConvectionScheme/gaussConvectionSchemes.C
 $(convectionSchemes)/multivariateGaussConvectionScheme/multivariateGaussConvectionSchemes.C
+$(convectionSchemes)/boundedConvectionScheme/boundedConvectionSchemes.C
 
 laplacianSchemes = finiteVolume/laplacianSchemes
 $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
new file mode 100644
index 0000000000000000000000000000000000000000..1e5c9117638b6b4b9ac59cf161c9b76631a2de35
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "boundedConvectionScheme.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvMatrices.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+boundedConvectionScheme<Type>::interpolate
+(
+    const surfaceScalarField& phi,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return scheme_().interpolate(phi, vf);
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+boundedConvectionScheme<Type>::flux
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return scheme_().flux(faceFlux, vf);
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type> >
+boundedConvectionScheme<Type>::fvmDiv
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return
+        scheme_().fvmDiv(faceFlux, vf)
+      - fvm::Sp(fvc::surfaceIntegrate(faceFlux), vf);
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+boundedConvectionScheme<Type>::fvcDiv
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return
+        scheme_().fvcDiv(faceFlux, vf)
+      - fvc::surfaceIntegrate(faceFlux)*vf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
new file mode 100644
index 0000000000000000000000000000000000000000..9426104d3398fd922c6bba88bd89ef9f3923fb72
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::boundedConvectionScheme
+
+Description
+    Bounded form of the selected convection scheme.
+
+    Boundedness is achieved by subtracting div(phi)*vf or Sp(div(phi), vf)
+    which is non-conservative if div(phi) != 0 but conservative otherwise.
+
+    Can be used for convection of bounded scalar properties in steady-state
+    solvers to improve stability if insufficient convergence of the pressure
+    equation causes temporary divergence of the flux field.
+
+SourceFiles
+    boundedConvectionScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef boundedConvectionScheme_H
+#define boundedConvectionScheme_H
+
+#include "convectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class boundedConvectionScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class boundedConvectionScheme
+:
+    public fv::convectionScheme<Type>
+{
+    // Private data
+
+        tmp<fv::convectionScheme<Type> > scheme_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        boundedConvectionScheme(const boundedConvectionScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const boundedConvectionScheme&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("bounded");
+
+
+    // Constructors
+
+        //- Construct from flux and Istream
+        boundedConvectionScheme
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            convectionScheme<Type>(mesh, faceFlux),
+            scheme_
+            (
+                fv::convectionScheme<Type>::New(mesh, faceFlux, is)
+            )
+        {}
+
+
+    // Member Functions
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > flux
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<fvMatrix<Type> > fvmDiv
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDiv
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "boundedConvectionScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C
similarity index 59%
rename from src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H
rename to src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C
index e156c58d141b23e834db6602b3ac96564dd07493..33385d4ecfabfceeb3d5c53bd1f813f8134c9ae5 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C
@@ -23,44 +23,17 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "boundedConvectionScheme.H"
+#include "fvMesh.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class CloudType>
-const Foam::volVectorField& Foam::VirtualMassForce<CloudType>::DUcDt() const
+namespace Foam
 {
-    if (DUcDtPtr_)
-    {
-        return *DUcDtPtr_;
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "const volVectorField& VirtualMassForce<CloudType>::DUcDt()"
-            "const"
-        )   << "DUcDt field not allocated" << abort(FatalError);
-
-        return *reinterpret_cast<const volVectorField*>(0);
-    }
-}
-
-
-template<class CloudType>
-inline const Foam::interpolation<Foam::vector>&
-Foam::VirtualMassForce<CloudType>::DUcDtInterp() const
+namespace fv
 {
-    if (!DUcDtInterpPtr_.valid())
-    {
-        FatalErrorIn
-        (
-            "inline const Foam::interpolation<Foam::vector>&"
-            "Foam::VirtualMassForce<CloudType>::DUcDtInterp() const"
-        )   << "Carrier pahase DUcDt interpolation object not set"
-            << abort(FatalError);
-    }
-
-    return DUcDtInterpPtr_();
+    makeFvConvectionScheme(boundedConvectionScheme)
+}
 }
-
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
index 1b27e361a9f87ce0752295d76cd1fdd60f469f60..61d6c22d37aadd055fa8546a2ddbab561760b6ad 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,10 @@ namespace Foam
 namespace fv
 {
 
+//- Temporary debug switch to provide warning about backward-compatibility
+//  issue with setting div schemes for steady-state
+extern int warnUnboundedGauss;
+
 /*---------------------------------------------------------------------------*\
                        Class gaussConvectionScheme Declaration
 \*---------------------------------------------------------------------------*/
@@ -103,7 +107,29 @@ public:
             (
                 surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
             )
-        {}
+        {
+            is.rewind();
+            word bounded(is);
+
+            if
+            (
+                warnUnboundedGauss
+             && word(mesh.ddtScheme("default")) == "steadyState"
+             && bounded != "bounded"
+            )
+            {
+                fileNameList controlDictFiles(findEtcFiles("controlDict"));
+
+                IOWarningIn("gaussConvectionScheme", is)
+                    << "Unbounded 'Gauss' div scheme used in "
+                       "steady-state solver, use 'bounded Gauss' "
+                       "to ensure boundedness.\n"
+                    << "    To remove this warning switch off "
+                    << "'boundedGauss' in "
+                    << controlDictFiles[controlDictFiles.size()-1]
+                    << endl;
+            }
+        }
 
 
     // Member Functions
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C
index 55a5e6049d90eadbf4201717b999ee276e500654..c5ea8dbac7c3890e4a36b904b0e79a26ecf86453 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,6 +32,11 @@ namespace Foam
 {
 namespace fv
 {
+    int warnUnboundedGauss
+    (
+        Foam::debug::debugSwitch("warnUnboundedGauss", true)
+    );
+
     makeFvConvectionScheme(gaussConvectionScheme)
 }
 }
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 321fba33ad63dd0c4dc09975a6ea8f097f7e9bdd..29e696dc81a734dd6b14150203d9dc6ff6425bc1 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -470,8 +470,8 @@ public:
             //- Total rotational kinetic energy in the system
             inline scalar rotationalKineticEnergyOfSystem() const;
 
-            //- Penetration for percentage of the current total mass
-            inline scalar penetration(const scalar& prc) const;
+            //- Penetration for fraction [0-1] of the current total mass
+            inline scalar penetration(const scalar& fraction) const;
 
             //- Mean diameter Dij
             inline scalar Dij(const label i, const label j) const;
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index a6bdf47cc46e750f79091d14201fad830c9e4a2a..2e793ac71e3ff9211c6761bc1176c17c1f2fbcaa 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvmSup.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -336,106 +337,130 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 
     reduce(d, maxOp<scalar>());
 
-    return d;
+    return max(0.0, d);
 }
 
 
 template<class CloudType>
 inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
 (
-    const scalar& prc
+    const scalar& fraction
 ) const
 {
-    scalar distance = 0.0;
-    scalar mTot = 0.0;
+    if ((fraction < 0) || (fraction > 1))
+    {
+        FatalErrorIn
+        (
+            "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
+            "("
+                "const scalar&"
+            ") const"
+        )
+            << "fraction should be in the range 0 < fraction < 1"
+            << exit(FatalError);
+    }
 
-    label np = this->size();
+    scalar distance = 0.0;
 
-    // arrays containing the parcels mass and
-    // distance from injector in ascending order
-    scalarField mass(np);
-    scalarField dist(np);
+    const label nParcel = this->size();
+    globalIndex globalParcels(nParcel);
+    const label nParcelSum = globalParcels.size();
 
-    if (np > 0)
+    if (nParcelSum == 0)
     {
-        label n = 0;
+        return distance;
+    }
 
-        // first arrange the parcels in ascending order
-        // the first parcel is closest to its injection position
-        // and the last one is most far away.
-        forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
-        {
-            const parcelType& p = iter();
-            scalar mi = p.nParticle()*p.mass();
-            scalar di = mag(p.position() - p.position0());
-            mTot += mi;
+    // lists of parcels mass and distance from initial injection point
+    List<scalar> mass(nParcel, 0.0);
+    List<scalar> dist(nParcel, 0.0);
 
-            // insert at the last place
-            mass[n] = mi;
-            dist[n] = di;
+    label i = 0;
+    scalar mSum = 0.0;
+    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+    {
+        const parcelType& p = iter();
+        scalar m = p.nParticle()*p.mass();
+        scalar d = mag(p.position() - p.position0());
+        mSum += m;
 
-            label i = 0;
-            bool found = false;
+        mass[i] = m;
+        dist[i] = d;
 
-            // insert the parcel in the correct place
-            // and move the others
-            while ((i < n) && (!found))
-            {
-                if (di < dist[i])
-                {
-                    found = true;
-                    for (label j=n; j>i; j--)
-                    {
-                        mass[j] = mass[j-1];
-                        dist[j] = dist[j-1];
-                    }
-                    mass[i] = mi;
-                    dist[i] = di;
-                }
-                i++;
-            }
-            n++;
-        }
+        i++;
     }
 
-    reduce(mTot, sumOp<scalar>());
+    // calculate total mass across all processors
+    reduce(mSum, sumOp<scalar>());
+
+    // flatten the mass list
+    List<scalar> allMass(nParcelSum, 0.0);
+    SubList<scalar>
+    (
+        allMass,
+        globalParcels.localSize(Pstream::myProcNo()),
+        globalParcels.offset(Pstream::myProcNo())
+    ).assign(mass);
+
+    // flatten the distance list
+    SortableList<scalar> allDist(nParcelSum, 0.0);
+    SubList<scalar>
+    (
+        allDist,
+        globalParcels.localSize(Pstream::myProcNo()),
+        globalParcels.offset(Pstream::myProcNo())
+    ).assign(dist);
+
+    // sort allDist distances into ascending order
+    // note: allMass masses are left unsorted
+    allDist.sort();
 
-    if (np > 0)
+    if (nParcelSum > 1)
     {
-        scalar mLimit = prc*mTot;
-        scalar mOff = (1.0 - prc)*mTot;
+        const scalar mLimit = fraction*mSum;
+        const labelList& indices = allDist.indices();
 
-        if (np > 1)
+        if (mLimit > (mSum - allMass[indices.last()]))
         {
-            // 'prc' is large enough that the parcel most far
-            // away will be used, no need to loop...
-            if (mLimit > mTot - mass[np-1])
-            {
-                distance = dist[np-1];
-            }
-            else
+            distance = allDist.last();
+        }
+        else
+        {
+            // assuming that 'fraction' is generally closer to 1 than 0, loop
+            // through in reverse distance order
+            const scalar mThreshold = (1.0 - fraction)*mSum;
+            scalar mCurrent = 0.0;
+            label i0 = 0;
+
+            forAllReverse(indices, i)
             {
-                scalar mOffSum = 0.0;
-                label i = np;
+                label indI = indices[i];
+
+                mCurrent += allMass[indI];
 
-                while ((mOffSum < mOff) && (i>0))
+                if (mCurrent > mThreshold)
                 {
-                    i--;
-                    mOffSum += mass[i];
+                    i0 = i;
+                    break;
                 }
-                distance =
-                    dist[i+1]
-                  + (dist[i] - dist[i+1])*(mOffSum - mOff)
-                   /mass[i+1] ;
             }
-        }
-        else
-        {
-            distance = dist[0];
+
+            if (i0 == indices.size() - 1)
+            {
+                distance = allDist.last();
+            }
+            else
+            {
+                // linearly interpolate to determine distance
+                scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
+                distance = allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
+            }
         }
     }
-
-    reduce(distance, maxOp<scalar>());
+    else
+    {
+        distance = allDist.first();
+    }
 
     return distance;
 }
diff --git a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H
index 01c80c16e21ee5a1bc8c8a518b8a7addeb37e6a0..d7a0fcf0d03056df7411dac483549f62240cbca0 100644
--- a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H
@@ -89,7 +89,7 @@ public:
             virtual scalar rotationalKineticEnergyOfSystem() const = 0;
 
             //- Penetration for percentage of the current total mass
-//            virtual scalar penetration(const scalar& prc) const = 0;
+//            virtual scalar penetration(const scalar& fraction) const = 0;
 
             //- Mean diameter Dij
             virtual scalar Dij(const label i, const label j) const = 0;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
index a1d121ee2990e0bbd078f8f882f22de9e1d99dc8..c9d61f6d34d606eb2f72618b435f17cea7bcdaa2 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PressureGradientForce.H"
+#include "fvcDdt.H"
 #include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -33,12 +34,13 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
-    UName_(this->coeffs().lookup("U")),
-    gradUPtr_(NULL)
+    ParticleForce<CloudType>(owner, mesh, dict, forceType, true),
+    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
+    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -50,7 +52,7 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce
 :
     ParticleForce<CloudType>(pgf),
     UName_(pgf.UName_),
-    gradUPtr_(NULL)
+    DUcDtInterpPtr_(NULL)
 {}
 
 
@@ -66,18 +68,48 @@ Foam::PressureGradientForce<CloudType>::~PressureGradientForce()
 template<class CloudType>
 void Foam::PressureGradientForce<CloudType>::cacheFields(const bool store)
 {
+    static word fName("DUcDt");
+
+    bool fieldExists = this->mesh().template foundObject<volVectorField>(fName);
+
     if (store)
     {
-        const volVectorField& U = this->mesh().template
-            lookupObject<volVectorField>(UName_);
-        gradUPtr_ = fvc::grad(U).ptr();
+        if (!fieldExists)
+        {
+            const volVectorField& Uc = this->mesh().template
+                lookupObject<volVectorField>(UName_);
+
+            volVectorField* DUcDtPtr = new volVectorField
+            (
+                fName,
+                fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
+            );
+
+            DUcDtPtr->store();
+        }
+
+        const volVectorField& DUcDt = this->mesh().template
+            lookupObject<volVectorField>(fName);
+
+        DUcDtInterpPtr_.reset
+        (
+            interpolation<vector>::New
+            (
+                this->owner().solution().interpolationSchemes(),
+                DUcDt
+            ).ptr()
+        );
     }
     else
     {
-        if (gradUPtr_)
+        DUcDtInterpPtr_.clear();
+
+        if (fieldExists)
         {
-            delete gradUPtr_;
-            gradUPtr_ = NULL;
+            const volVectorField& DUcDt = this->mesh().template
+                lookupObject<volVectorField>(fName);
+
+            const_cast<volVectorField&>(DUcDt).checkOut();
         }
     }
 }
@@ -95,11 +127,24 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
 {
     forceSuSp value(vector::zero, 0.0);
 
-    const volTensorField& gradU = *gradUPtr_;
-    value.Su() = mass*p.rhoc()/p.rho()*(p.U() & gradU[p.cell()]);
+    vector DUcDt =
+        DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
+
+    value.Su() = mass*p.rhoc()/p.rho()*DUcDt;
 
     return value;
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::PressureGradientForce<CloudType>::massAdd
+(
+    const typename CloudType::parcelType&,
+    const scalar
+) const
+{
+    return 0.0;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
index 663723a82093a8bf9526af9e84a227d9a8c0d2ff..b1986348fba94f22f64ada538c6008935d6032b0 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "ParticleForce.H"
 #include "volFields.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,13 +54,15 @@ class PressureGradientForce
 :
     public ParticleForce<CloudType>
 {
-    // Private data
+protected:
+
+    // Protected data
 
         //- Name of velocity field
         const word UName_;
 
-        //- Velocity gradient field
-        const volTensorField* gradUPtr_;
+        //- Rate of change of carrier phase velocity interpolator
+        autoPtr<interpolation<vector> > DUcDtInterpPtr_;
 
 
 public:
@@ -75,7 +78,8 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& forceType = typeName
         );
 
         //- Construct copy
@@ -99,8 +103,8 @@ public:
 
         // Access
 
-            //- Return const access to the velocity gradient field
-            inline const volTensorField& gradU() const;
+            //- Return the rate of change of carrier phase velocity interpolator
+            inline const interpolation<vector>& DUcDtInterp() const;
 
 
         // Evaluation
@@ -117,6 +121,13 @@ public:
                 const scalar Re,
                 const scalar muc
             ) const;
+
+            //- Return the added mass
+            virtual scalar massAdd
+            (
+                const typename CloudType::parcelType& p,
+                const scalar mass
+            ) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H
index c9bd24c62c66a0df23ebf75ad245fc3e2c571b52..6c085241dedf4289c9cc7b183141a57a0e401e8e 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,23 +26,20 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class CloudType>
-const Foam::volTensorField& Foam::PressureGradientForce<CloudType>::gradU()
-const
+inline const Foam::interpolation<Foam::vector>&
+Foam::PressureGradientForce<CloudType>::DUcDtInterp() const
 {
-    if (gradUPtr_)
-    {
-        return *gradUPtr_;
-    }
-    else
+    if (!DUcDtInterpPtr_.valid())
     {
         FatalErrorIn
         (
-            "const volTensorField& PressureGradientForce<CloudType>::gradU()"
-            "const"
-        )   << "gradU field not allocated" << abort(FatalError);
-
-        return *reinterpret_cast<const volTensorField*>(0);
+            "inline const Foam::interpolation<Foam::vector>&"
+            "Foam::PressureGradientForce<CloudType>::DUcDtInterp() const"
+        )   << "Carrier phase DUcDt interpolation object not set"
+            << abort(FatalError);
     }
+
+    return DUcDtInterpPtr_();
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
index 3dca602dd443f59548e4a671fd8686b827ae17dc..c701e07a00d17d64de5fe69c9b0d7f8bbfb19ba6 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
@@ -24,8 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VirtualMassForce.H"
-#include "fvcDdt.H"
-#include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -34,14 +32,12 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
-    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
-    Cvm_(readScalar(this->coeffs().lookup("Cvm"))),
-    DUcDtPtr_(NULL),
-    DUcDtInterpPtr_(NULL)
+    PressureGradientForce<CloudType>(owner, mesh, dict, forceType),
+    Cvm_(readScalar(this->coeffs().lookup("Cvm")))
 {}
 
 
@@ -51,11 +47,8 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce
     const VirtualMassForce& vmf
 )
 :
-    ParticleForce<CloudType>(vmf),
-    UName_(vmf.UName_),
-    Cvm_(vmf.Cvm_),
-    DUcDtPtr_(NULL),
-    DUcDtInterpPtr_(NULL)
+    PressureGradientForce<CloudType>(vmf),
+    Cvm_(vmf.Cvm_)
 {}
 
 
@@ -71,36 +64,7 @@ Foam::VirtualMassForce<CloudType>::~VirtualMassForce()
 template<class CloudType>
 void Foam::VirtualMassForce<CloudType>::cacheFields(const bool store)
 {
-    if (store && !DUcDtPtr_)
-    {
-        const volVectorField& Uc = this->mesh().template
-            lookupObject<volVectorField>(UName_);
-
-        DUcDtPtr_ = new volVectorField
-        (
-            "DUcDt",
-            fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
-        );
-
-        DUcDtInterpPtr_.reset
-        (
-            interpolation<vector>::New
-            (
-                this->owner().solution().interpolationSchemes(),
-                *DUcDtPtr_
-            ).ptr()
-        );
-    }
-    else
-    {
-        DUcDtInterpPtr_.clear();
-
-        if (DUcDtPtr_)
-        {
-            delete DUcDtPtr_;
-            DUcDtPtr_ = NULL;
-        }
-    }
+    PressureGradientForce<CloudType>::cacheFields(store);
 }
 
 
@@ -114,12 +78,10 @@ Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
     const scalar muc
 ) const
 {
-    forceSuSp value(vector::zero, 0.0);
+    forceSuSp value =
+        PressureGradientForce<CloudType>::calcCoupled(p, dt, mass, Re, muc);
 
-    vector DUcDt =
-        DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
-
-    value.Su() = mass*p.rhoc()/p.rho()*Cvm_*DUcDt;
+    value.Su() *= Cvm_;
 
     return value;
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
index c5277ec6963a7b09c0b20a91ef575d4f1f22ad17..ea0d5809324da2d5bd2024ea6a40e34aafa5c527 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
@@ -28,7 +28,6 @@ Description
     Calculates particle virtual mass force
 
 SourceFiles
-    VirtualMassForceI.H
     VirtualMassForce.C
 
 \*---------------------------------------------------------------------------*/
@@ -36,9 +35,7 @@ SourceFiles
 #ifndef VirtualMassForce_H
 #define VirtualMassForce_H
 
-#include "ParticleForce.H"
-#include "volFields.H"
-#include "interpolation.H"
+#include "PressureGradientForce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,22 +49,13 @@ namespace Foam
 template<class CloudType>
 class VirtualMassForce
 :
-    public ParticleForce<CloudType>
+    public PressureGradientForce<CloudType>
 {
     // Private data
 
-        //- Name of velocity field
-        const word UName_;
-
         //- Virtual mass coefficient - typically 0.5
         scalar Cvm_;
 
-        //- Rate of change of carrier phase velocity
-        volVectorField* DUcDtPtr_;
-
-        //- Rate of change of carrier phase velocity interpolator
-        autoPtr<interpolation<vector> > DUcDtInterpPtr_;
-
 
 public:
 
@@ -82,7 +70,8 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& forceType = typeName
         );
 
         //- Construct copy
@@ -104,15 +93,6 @@ public:
 
     // Member Functions
 
-        // Access
-
-            //- Return the rate of change of carrier phase velocity
-            inline const volVectorField& DUcDt() const;
-
-            //- Return the rate of change of carrier phase velocity interpolator
-            inline const interpolation<vector>& DUcDtInterp() const;
-
-
         // Evaluation
 
             //- Cache fields
@@ -143,8 +123,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "VirtualMassForceI.H"
-
 #ifdef NoRepository
     #include "VirtualMassForce.C"
 #endif
diff --git a/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.C b/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.C
index b675bac3a37b371e765a3147b502007667f3b4df..0e11217066c8e6e42eae1a6250d58961a3b71ddb 100644
--- a/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.C
+++ b/src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.C
@@ -123,7 +123,7 @@ temperatureThermoBaffleFvPatchScalarField
      && !owner_
     )
     {
-        Info << "Creating thermal baffle..." <<  nbrMesh << endl;
+        Info << "Creating thermal baffle" <<  nbrMesh << endl;
         baffle_.reset(baffle::New(thisMesh, dict).ptr());
         owner_ = true;
         dict.lookup("thermoType") >> solidThermoType_;
diff --git a/src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.H b/src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.H
index f0c3181a90d61374cf519c481b864015c67f9330..a8949e89df88e16c167c4b6ef7e0190c094e8123 100644
--- a/src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.H
+++ b/src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,11 +132,9 @@ public:
             const word& modelType,
             const fvMesh& mesh,
             const dictionary& dict
-
         );
 
 
-
     //- Destructor
     virtual ~thermoBaffle2D();
 
@@ -182,6 +180,7 @@ public:
                 //- Return sensible enthalpy/internal energy
                 inline tmp<volScalarField> he() const;
 
+
         // Evolution
 
             //- Pre-evolve  thermal baffle
@@ -191,7 +190,6 @@ public:
             virtual void evolveRegion();
 
 
-
        // I-O
 
             //- Provide some feedback
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index a2ad3c88f7d340257ddb0e1f0202ad19e1efbbce..13d6f0462c763525fc27551cfd574ec434ac1559 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -164,6 +164,89 @@ Foam::basicThermo::~basicThermo()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::basicThermo::validate
+(
+    const word& app,
+    const word& a
+) const
+{
+    if (!(he().name() == a))
+    {
+        FatalErrorIn(app)
+            << "Supported energy type is " << a
+            << ", thermodynamics package provides " << he().name()
+            << exit(FatalError);
+    }
+}
+
+void Foam::basicThermo::validate
+(
+    const word& app,
+    const word& a,
+    const word& b
+) const
+{
+    if (!(he().name() == a || he().name() == b))
+    {
+        FatalErrorIn(app)
+            << "Supported energy types are " << a << " and " << b
+            << ", thermodynamics package provides " << he().name()
+            << exit(FatalError);
+    }
+}
+
+void Foam::basicThermo::validate
+(
+    const word& app,
+    const word& a,
+    const word& b,
+    const word& c
+) const
+{
+    if
+    (
+       !(
+            he().name() == a
+         || he().name() == b
+         || he().name() == c
+        )
+    )
+    {
+        FatalErrorIn(app)
+            << "Supported energy types are " << a << ", " << b << " and " << c
+            << ", thermodynamics package provides " << he().name()
+            << exit(FatalError);
+    }
+}
+
+void Foam::basicThermo::validate
+(
+    const word& app,
+    const word& a,
+    const word& b,
+    const word& c,
+    const word& d
+) const
+{
+    if
+    (
+       !(
+            he().name() == a
+         || he().name() == b
+         || he().name() == c
+         || he().name() == d
+        )
+    )
+    {
+        FatalErrorIn(app)
+            << "Supported energy types are " << a << ", " << b
+            << ", " << c << " and " << d
+            << ", thermodynamics package provides " << he().name()
+            << exit(FatalError);
+    }
+}
+
+
 Foam::volScalarField& Foam::basicThermo::p()
 {
     return p_;
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.H b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
index c9559d4ef0a8feb996c5b8a77b5c8a664d1c4fcd..d25c9c1ef494b691b9a514df0814f90f6b516ff2 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -111,6 +111,45 @@ public:
 
     // Member functions
 
+        //- Check that the thermodynamics package is consistent
+        //  with energy forms supported by the application
+        void validate
+        (
+            const word& app,
+            const word&
+        ) const;
+
+        //- Check that the thermodynamics package is consistent
+        //  with energy forms supported by the application
+        void validate
+        (
+            const word& app,
+            const word&,
+            const word&
+        ) const;
+
+        //- Check that the thermodynamics package is consistent
+        //  with energy forms supported by the application
+        void validate
+        (
+            const word& app,
+            const word&,
+            const word&,
+            const word&
+        ) const;
+
+        //- Check that the thermodynamics package is consistent
+        //  with energy forms supported by the application
+        void validate
+        (
+            const word& app,
+            const word&,
+            const word&,
+            const word&,
+            const word&
+        ) const;
+
+
         //- Update properties
         virtual void correct() = 0;
 
diff --git a/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
index c439348fc045046f9d5e1032b45df61db6b11332..af89aec73a7c433d77ab1ae350304b75c4033833 100644
--- a/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
+++ b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
@@ -142,6 +142,24 @@ makeBasicMixture
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+makeBasicMixture
+(
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    eConstThermo,
+    perfectGas
+);
+
+makeBasicMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    sensibleInternalEnergy,
+    eConstThermo,
+    perfectGas
+);
+
 makeBasicMixture
 (
     pureMixture,
diff --git a/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C b/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
index 56baf3fa0bea99de957b70977791e37d3d17172a..b5a48fc315479656918a4a463f52a721a69e4e08 100644
--- a/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
+++ b/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
@@ -84,27 +84,27 @@ makeThermo
 
 /* * * * * * * * * * * * * * Internal-energy-based * * * * * * * * * * * * * */
 
-// makeThermo
-// (
-//     psiThermo,
-//     hePsiThermo,
-//     pureMixture,
-//     constTransport,
-//     sensibleInternalEnergy,
-//     eConstThermo,
-//     perfectGas
-// );
-
-// makeThermo
-// (
-//     psiThermo,
-//     hePsiThermo,
-//     pureMixture,
-//     sutherlandTransport,
-//     sensibleInternalEnergy,
-//     eConstThermo,
-//     perfectGas
-// );
+makeThermo
+(
+    psiThermo,
+    hePsiThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    eConstThermo,
+    perfectGas
+);
+
+makeThermo
+(
+    psiThermo,
+    hePsiThermo,
+    pureMixture,
+    sutherlandTransport,
+    sensibleInternalEnergy,
+    eConstThermo,
+    perfectGas
+);
 
 makeThermo
 (
diff --git a/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.C b/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.C
index 8d6c67c53836c68e6f8e6dc0ae3964386eca3757..8348899c45907478c5576e8929d2bbd30d9aae10 100644
--- a/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.C
+++ b/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.C
@@ -43,6 +43,21 @@ constSolidRad<thermo>::constSolidRad(const dictionary& dict)
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class thermo>
+void Foam::constSolidRad<thermo>::constSolidRad::write(Ostream& os) const
+{
+    thermo::write(os);
+
+    dictionary dict("radiation");
+    dict.add("kappaRad", kappaRad_);
+    dict.add("sigmaS", sigmaS_);
+    dict.add("emissivity", emissivity_);
+    os  << indent << dict.dictName() << dict;
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class thermo>
diff --git a/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.H b/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.H
index af239b91a7464bda2a29ad02e0c5a10ac233100f..60977f6e83156dbd57db430ca581055b31959cc7 100644
--- a/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.H
+++ b/src/thermophysicalModels/solidSpecie/radiation/const/constSolidRad.H
@@ -58,7 +58,7 @@ Ostream& operator<<
 );
 
 /*---------------------------------------------------------------------------*\
-                      Class constSolidRad Declaration
+                        Class constSolidRad Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class thermo>
@@ -116,6 +116,9 @@ public:
         //- Return  emissivity[]
         inline scalar emissivity(scalar T) const;
 
+        //- Write to Ostream
+        void write(Ostream& os) const;
+
 
     // Member operators
 
diff --git a/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.C b/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.C
index a1a38cc3d0016d7ac2617728e4dfd16315dbe29d..e2c5ae6719a8de95a351eeac612b8f8ad2ff1981 100644
--- a/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.C
+++ b/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.C
@@ -39,6 +39,22 @@ Foam::constAnIsoSolidTransport<thermo>::constAnIsoSolidTransport
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class thermo>
+void Foam::constAnIsoSolidTransport<thermo>::constAnIsoSolidTransport::write
+(
+    Ostream& os
+) const
+{
+    thermo::write(os);
+
+    dictionary dict("transport");
+    dict.add("kappa", kappa_);
+    os  << indent << dict.dictName() << dict;
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class thermo>
diff --git a/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.H b/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.H
index 693afd432aa85eb37772d84966450ec168c7f409..e09d43a1cdacbf04e5de5858928c5c8ec039ac3e 100644
--- a/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.H
+++ b/src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.H
@@ -62,7 +62,7 @@ Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                     Class constAnIsoSolidTransport Declaration
+                   Class constAnIsoSolidTransport Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class thermo>
@@ -79,11 +79,7 @@ class constAnIsoSolidTransport
     // Private Member Functions
 
         //- Construct from components
-        inline constAnIsoSolidTransport
-        (
-            const thermo& t,
-            const vector kappa
-        );
+        inline constAnIsoSolidTransport(const thermo& t, const vector kappa);
 
 
 public:
@@ -110,6 +106,9 @@ public:
         //- Un-isotropic thermal conductivity [W/mK]
         inline vector Kappa(const scalar T) const;
 
+        //- Write to Ostream
+        void write(Ostream& os) const;
+
 
     // Member operators
 
diff --git a/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.C b/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.C
index 5e7f34cf27a4708904b080ea4ae164a452623e1a..e2fc5096bd13f241ce62f705fad2d279891c5315 100644
--- a/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.C
+++ b/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.C
@@ -39,6 +39,22 @@ Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class thermo>
+void Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport::write
+(
+    Ostream& os
+) const
+{
+    thermo::write(os);
+
+    dictionary dict("transport");
+    dict.add("kappa", kappa_);
+    os  << indent << dict.dictName() << dict;
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class thermo>
diff --git a/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.H b/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.H
index 998a60f5ce0b4bc6129aebc8bf1b07ed8fe898cf..b32ae358ee3096c53f8bf637d06598f695e7e3df 100644
--- a/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.H
+++ b/src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.H
@@ -63,7 +63,7 @@ Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                         Class constIsoSolidTransport Declaration
+                   Class constIsoSolidTransport Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class thermo>
@@ -73,18 +73,14 @@ class constIsoSolidTransport
 {
     // Private data
 
-        //- Constant isotropic thermal conductivity.
+        //- Constant isotropic thermal conductivity
         scalar kappa_;
 
 
     // Private Member Functions
 
         //- Construct from components
-        inline constIsoSolidTransport
-        (
-            const thermo& t,
-            const scalar kappa
-        );
+        inline constIsoSolidTransport(const thermo& t, const scalar kappa);
 
 
 public:
@@ -99,8 +95,7 @@ public:
         );
 
         //- Construct from Istream
-        //constIsoSolidTransport(Istream&);
-        constIsoSolidTransport(const dictionary&);
+        constIsoSolidTransport(const dictionary& dict);
 
 
     // Member functions
@@ -111,6 +106,9 @@ public:
         //- Un-isotropic thermal conductivity [W/mK]
         inline vector Kappa(const scalar T) const;
 
+        //- Write to Ostream
+        void write(Ostream& os) const;
+
 
     // Member operators
 
diff --git a/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.C b/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.C
index b479bd4b66db99cd4801eb7bf0fd3e490fd794b1..10070c48ba43c9a200674caa7976591c08b666c7 100644
--- a/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.C
+++ b/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.C
@@ -46,6 +46,24 @@ Foam::exponentialSolidTransport<thermo>::exponentialSolidTransport
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class thermo>
+void Foam::exponentialSolidTransport<thermo>::exponentialSolidTransport::write
+(
+    Ostream& os
+) const
+{
+    thermo::write(os);
+
+    dictionary dict("transport");
+    dict.add("kappa0", kappa0_);
+    dict.add("n0", n0_);
+    dict.add("Tref", Tref_);
+    os  << indent << dict.dictName() << dict;
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class thermo>
diff --git a/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H b/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H
index b89215b42179a67e7de0338c0550f95b5497f6ac..1ace8a3c4518a02f3332c086072890758ee84af5 100644
--- a/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H
+++ b/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H
@@ -60,7 +60,7 @@ Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                        Class exponentialSolidTransport Declaration
+                  Class exponentialSolidTransport Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class thermo>
@@ -116,6 +116,9 @@ public:
         //- Thermal conductivity [W/mK]
         inline vector Kappa(const scalar T) const;
 
+        //- Write to Ostream
+        void write(Ostream& os) const;
+
 
     // Member operators
 
diff --git a/src/thermophysicalModels/solidThermo/heSolidThermo/heSolidThermo.H b/src/thermophysicalModels/solidThermo/heSolidThermo/heSolidThermo.H
index 88a8d8aab6147731c2eab0b75508a6c9238f9f94..3a154f7fee5cf50a207ad78abd63bf823f9a8b16 100644
--- a/src/thermophysicalModels/solidThermo/heSolidThermo/heSolidThermo.H
+++ b/src/thermophysicalModels/solidThermo/heSolidThermo/heSolidThermo.H
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class heSolidThermo Declaration
+                        Class heSolidThermo Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class MixtureType, class BasicSolidThermo>
@@ -62,6 +62,7 @@ class heSolidThermo
         //- Construct as copy (not implemented)
         heSolidThermo(const heSolidThermo<MixtureType, BasicSolidThermo>&);
 
+
 public:
 
     //- Runtime type information
@@ -104,30 +105,17 @@ public:
 
         // Per patch calculation
 
-            //- Anisotropic thermal conductivity [W//m/K]
-            virtual tmp<vectorField> Kappa
-            (
-                const label patchI
-            ) const;
+            //- Anisotropic thermal conductivity [W/m/K]
+            virtual tmp<vectorField> Kappa(const label patchI) const;
 
              //- Absorption coefficient [1/m]
-            virtual tmp<scalarField> kappaRad
-            (
-                const label patchI
-            ) const;
+            virtual tmp<scalarField> kappaRad(const label patchI) const;
 
             //- Scatter coefficient
-            virtual tmp<scalarField> sigmaS
-            (
-                const label patchI
-            ) const;
+            virtual tmp<scalarField> sigmaS(const label patchI) const;
 
             //- Emissivity coefficient [1/m]
-            virtual tmp<scalarField> emissivity
-            (
-                const label patchI
-            ) const;
-
+            virtual tmp<scalarField> emissivity(const label patchI) const;
 };
 
 
@@ -135,7 +123,7 @@ public:
 
 } // End namespace Foam
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
 #   include "heSolidThermo.C"
diff --git a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
index ee80730baa6c7092aa99dfaf538df270f3c8a837..9d186eb039b6eb9d6f9b67437f81817be18e4488 100644
--- a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
@@ -292,7 +292,7 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
+      //***HGW - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_
@@ -314,7 +314,7 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
-      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
+      //***HGW - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.C
index f85d85791a68b52bae276c56d972347cf6ca61d2..e33fa644c235bd6fc63f217e772cbbb62136189f 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.C
@@ -51,7 +51,7 @@ temperatureThermoBaffle1DFvPatchScalarField
     baffleActivated_(true),
     thickness_(p.size()),
     Qs_(p.size()),
-    solid_(NULL)
+    solidPtr_(NULL)
 {}
 
 
@@ -70,7 +70,7 @@ temperatureThermoBaffle1DFvPatchScalarField
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
-    solid_(ptf.solid_)
+    solidPtr_(ptf.solidPtr_)
 {}
 
 
@@ -88,19 +88,19 @@ temperatureThermoBaffle1DFvPatchScalarField
     baffleActivated_(readBool(dict.lookup("baffleActivated"))),
     thickness_(scalarField("thickness", dict, p.size())),
     Qs_(scalarField("Qs", dict, p.size())),
-    solid_(new solidThermoData(dict))
+    solidPtr_(new solidType(dict))
 {
     if (!isA<mappedPatchBase>(this->patch().patch()))
     {
         FatalErrorIn
         (
             "temperatureThermoBaffle1DFvPatchScalarField::"
-            "temperatureThermoBaffle1DFvPatchScalarField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<scalar, volMesh>& iF,\n"
-            "    const dictionary& dict\n"
-            ")\n"
+            "temperatureThermoBaffle1DFvPatchScalarField"
+            "("
+                "const fvPatch&,\n"
+                "const DimensionedField<scalar, volMesh>&, "
+                "const dictionary&"
+            ")"
         )   << "\n    patch type '" << patch().type()
             << "' not type '" << mappedPatchBase::typeName << "'"
             << "\n    for patch " << patch().name()
@@ -141,7 +141,7 @@ temperatureThermoBaffle1DFvPatchScalarField
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
-    solid_(ptf.solid_)
+    solidPtr_(ptf.solidPtr_)
 {}
 
 
@@ -158,7 +158,7 @@ temperatureThermoBaffle1DFvPatchScalarField
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
-    solid_(ptf.solid_)
+    solidPtr_(ptf.solidPtr_)
 {}
 
 
@@ -197,10 +197,8 @@ void temperatureThermoBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
     int oldTag = UPstream::msgType();
     UPstream::msgType() = oldTag+1;
 
-    const mappedPatchBase& mpp = refCast<const mappedPatchBase>
-    (
-        patch().patch()
-    );
+    const mappedPatchBase& mpp =
+        refCast<const mappedPatchBase>(patch().patch());
 
     const label patchI = patch().index();
 
@@ -208,8 +206,7 @@ void temperatureThermoBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
 
     if (baffleActivated_)
     {
-        const fvPatch& nbrPatch =
-            patch().boundaryMesh()[mpp.samplePolyPatch().index()];
+        const fvPatch& nbrPatch = patch().boundaryMesh()[nbrPatchI];
 
         const compressible::turbulenceModel& model =
             db().template lookupObject<compressible::turbulenceModel>
@@ -248,10 +245,7 @@ void temperatureThermoBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
         const temperatureThermoBaffle1DFvPatchScalarField& nbrField =
         refCast<const temperatureThermoBaffle1DFvPatchScalarField>
         (
-            nbrPatch.template lookupPatchField<volScalarField, scalar>
-            (
-                TName_
-            )
+            nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
         );
 
         scalarField nbrTi(nbrField.patchInternalField());
@@ -278,7 +272,8 @@ void temperatureThermoBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
         // Create fields for solid properties
         forAll(KDeltaw, i)
         {
-            KDeltaw[i] = solid_().kappa((Tp[i] + nbrTw[i])/2.0)/thickness_[i];
+            KDeltaw[i] =
+                solidPtr_().kappa((Tp[i] + nbrTw[i])/2.0)/thickness_[i];
         }
 
         const scalarField q
@@ -362,7 +357,7 @@ write(Ostream& os) const
     os.writeKeyword("baffleActivated")
         << baffleActivated_ << token::END_STATEMENT << nl;
     Qs_.writeEntry("Qs", os);
-    solid_().write(os);
+    solidPtr_->write(os);
 }
 
 
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.H b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.H
index 51222848474bafd14cad7e23766c3e4d311f91c9..524fb9c95a6747bee85b9f5f0c7104c322d90e52 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/temperatureThermoBaffle1D/temperatureThermoBaffle1DFvPatchScalarField.H
@@ -48,7 +48,7 @@ namespace compressible
 {
 
 /*---------------------------------------------------------------------------*\
-        Class temperatureThermoBaffle1DFvPatchScalarField Declaration
+         Class temperatureThermoBaffle1DFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class solidType>
@@ -58,74 +58,6 @@ class temperatureThermoBaffle1DFvPatchScalarField
 {
     // Private data
 
-    class solidThermoData
-    {
-            // Solid thermo
-            autoPtr<solidType> solidPtr_;
-
-            // Solid dictionaries
-            dictionary  specieDict_;
-            dictionary  transportDict_;
-            dictionary  radiationDict_;
-            dictionary  thermoDict_;
-            dictionary  densityDict_;
-
-
-    public:
-
-        // Constructor from components
-        solidThermoData(const dictionary& dict)
-        :
-            solidPtr_(new solidType(dict)),
-            specieDict_(dict.subDict("specie")),
-            transportDict_(dict.subDict("transport")),
-            radiationDict_(dict.subDict("radiation")),
-            thermoDict_(dict.subDict("thermodynamics")),
-            densityDict_(dict.subDict("equationOfState"))
-        {}
-
-
-        // Null constructor
-        solidThermoData()
-        :
-            solidPtr_(),
-            specieDict_(),
-            transportDict_(),
-            radiationDict_(),
-            thermoDict_(),
-            densityDict_()
-        {}
-
-
-        // Destructor
-        virtual ~solidThermoData()
-        {}
-
-
-        // Member Functions
-
-            void write(Ostream& os) const
-            {
-                os.writeKeyword("specie");
-                os << specieDict_  << nl;
-                os.writeKeyword("transport");
-                os << transportDict_  << nl;
-                os.writeKeyword("radiation");
-                os << radiationDict_ <<  nl;
-                os.writeKeyword("thermodynamics");
-                os << thermoDict_ << nl;
-                os.writeKeyword("density");
-                os << densityDict_ << nl;
-            }
-
-
-            scalar kappa(const scalar T) const
-            {
-                return solidPtr_().kappa(T);
-            }
-    };
-
-
         //- Name of the temperature field
         word TName_;
 
@@ -138,8 +70,16 @@ class temperatureThermoBaffle1DFvPatchScalarField
         //- Superficial heat source [W/m2]
         scalarField Qs_;
 
-        //- Solid thermo
-        autoPtr<solidThermoData> solid_;
+        // Solid thermo
+        autoPtr<solidType> solidPtr_;
+
+
+    // Private Member Functions
+
+        scalar kappa(const scalar T) const
+        {
+            return solidPtr_().kappa(T);
+        }
 
 
 public:
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
index 29757f1ea6eb381255799639c79c6c4352cc54a0..fbae992c25fd4f773a95087dfa40ca4fe3bef3e7 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
@@ -79,7 +79,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), "undefined", "undefined-K"),
     heatSource_(hsPower),
-    q_(p.size(), 0.0)
+    q_(p.size(), 0.0),
+    QrName_("undefinedQr")
 {}
 
 
@@ -95,7 +96,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
     temperatureCoupledBase(patch(), ptf.KMethod(), ptf.kappaName()),
     heatSource_(ptf.heatSource_),
-    q_(ptf.q_, mapper)
+    q_(ptf.q_, mapper),
+    QrName_(ptf.QrName_)
 {}
 
 
@@ -110,7 +112,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), dict),
     heatSource_(heatSourceTypeNames_.read(dict.lookup("heatSource"))),
-    q_("q", dict, p.size())
+    q_("q", dict, p.size()),
+    QrName_(dict.lookupOrDefault<word>("Qr", "none"))
 {
     fvPatchField<scalar>::operator=(patchInternalField());
     gradient() = 0.0;
@@ -126,7 +129,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(thftpsf),
     temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()),
     heatSource_(thftpsf.heatSource_),
-    q_(thftpsf.q_)
+    q_(thftpsf.q_),
+    QrName_(thftpsf.QrName_)
 {}
 
 
@@ -140,7 +144,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField
     fixedGradientFvPatchScalarField(thftpsf, iF),
     temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()),
     heatSource_(thftpsf.heatSource_),
-    q_(thftpsf.q_)
+    q_(thftpsf.q_),
+    QrName_(thftpsf.QrName_)
 {}
 
 
@@ -183,17 +188,25 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
 
     const scalarField& Tp = *this;
 
+    scalarField qr(this->size(), 0.0);
+
+    //- qr is negative going into the domain
+    if (QrName_ != "none")
+    {
+        qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
+    }
+
     switch (heatSource_)
     {
         case hsPower:
         {
             const scalar Ap = gSum(patch().magSf());
-            gradient() = q_/(Ap*kappa(Tp));
+            gradient() = (q_/Ap + qr)/kappa(Tp);
             break;
         }
         case hsFlux:
         {
-            gradient() = q_/kappa(Tp);
+            gradient() = (q_ + qr)/kappa(Tp);
             break;
         }
         default:
@@ -220,6 +233,7 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::write
         << token::END_STATEMENT << nl;
     temperatureCoupledBase::write(os);
     q_.writeEntry("q", os);
+    os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H
index 254ed52b756f1ed0521f1d50b33514a6af9b6c89..740183e35902541cc29b212efcd879cdf71266ce 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H
@@ -37,6 +37,7 @@ Description
             heatSource      flux;        // power [W]; flux [W/m2]
             q               uniform 10;  // heat power or flux
             kappa           fluidThermo; // calculate kappa=alphaEff*thermo.Cp
+            Qr              none;        // name of the radiative flux
             value           uniform 300; // initial temperature value
         }
 
@@ -93,6 +94,9 @@ private:
         //- Heat power [W] or flux [W/m2]
         scalarField q_;
 
+        //- Name of radiative in flux field
+        word QrName_;
+
 
 public:
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
index c1bf8f9b66ba51ade8ac05f9a358369fd315de4c..f0770e8bd9f126c2b835c19726d70b454ea2a2d4 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
@@ -53,7 +53,7 @@ boundaryField
 
     "(region0_to.*)"
     {
-        type            buoyantPressure;
+        type            fixedFluxPressure;
         value           $internalField;
     }
 }
diff --git a/tutorials/compressible/rhoCentralFoam/shockTube/constant/polyMesh/boundary b/tutorials/compressible/rhoCentralFoam/shockTube/constant/polyMesh/boundary
index a38765097f42199308ac35221e7114385d1eef6f..ef087376e6eafda48b5827e2678d2885571effad 100644
--- a/tutorials/compressible/rhoCentralFoam/shockTube/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoCentralFoam/shockTube/constant/polyMesh/boundary
@@ -26,6 +26,7 @@ FoamFile
     empty
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          400;
         startFace       101;
     }
diff --git a/tutorials/compressible/sonicFoam/laminar/forwardStep/constant/polyMesh/boundary b/tutorials/compressible/sonicFoam/laminar/forwardStep/constant/polyMesh/boundary
index 2abc731c97d64df322dd7036bb0c71065ff7ee86..53cce9211b1a351e8a33abee2738ebd881346c2d 100644
--- a/tutorials/compressible/sonicFoam/laminar/forwardStep/constant/polyMesh/boundary
+++ b/tutorials/compressible/sonicFoam/laminar/forwardStep/constant/polyMesh/boundary
@@ -32,12 +32,14 @@ FoamFile
     bottom
     {
         type            symmetryPlane;
+        inGroups        1(symmetryPlane);
         nFaces          25;
         startFace       10415;
     }
     top
     {
         type            symmetryPlane;
+        inGroups        1(symmetryPlane);
         nFaces          125;
         startFace       10440;
     }
@@ -50,6 +52,7 @@ FoamFile
     defaultFaces
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          10500;
         startFace       10675;
     }
diff --git a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSolution b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSolution
index 45bb949628482ca06e9e8c15979e678f38f2517c..e07aa4bafc6cc262e95eeedb2ae3e2d558b50fc1 100644
--- a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSolution
+++ b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    "p.*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -25,14 +25,14 @@ solvers
         relTol          0;
     }
 
-    "(U|e)"
+    "(U|e).*"
     {
         $p;
         tolerance       1e-05;
         relTol          0;
     }
 
-    rho
+    "rho.*"
     {
         solver          PCG;
         preconditioner  DIC;
@@ -41,9 +41,10 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors     2;
+    nOuterCorrectors 2;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary
index e4b48d35f4f621504a7e32d598b928415b8823d7..7fc4ff3fe5825c9dbf9d6dee48678f01ce9b4baf 100644
--- a/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary
+++ b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary
@@ -26,6 +26,7 @@ FoamFile
     empty
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          4000;
         startFace       1001;
     }
diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSolution b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSolution
index b894099c4c44ae64498d01ccd4a3268bd45e43af..9ed6e1a98226d5a50b3c5764e423fa6461607b28 100644
--- a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSolution
+++ b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    "(p|U|e)"
+    "(p|U|e).*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -25,7 +25,7 @@ solvers
         relTol          0;
     }
 
-    rho
+    "rho.*"
     {
         solver          PCG;
         preconditioner  DIC;
@@ -34,9 +34,10 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors     2;
+    nOuterCorrectors     2;
+    nCorrectors          1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/polyMesh/boundary b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/polyMesh/boundary
index 59a99aa3f008e969e793042560e7658577622325..b6f0d11024aaf164309fa920c645cd84269497b5 100644
--- a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/polyMesh/boundary
+++ b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/polyMesh/boundary
@@ -32,6 +32,7 @@ FoamFile
     SYMP3
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          80000;
         startFace       80170;
     }
diff --git a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSolution b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSolution
index 530a5bbd34bf71661fd8fc91815d685d199a6cc7..afc24d24a2e9106bb45a2a1013a763e503715c94 100644
--- a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSolution
+++ b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSolution
@@ -17,12 +17,12 @@ FoamFile
 
 solvers
 {
-    rho
+    "rho.*"
     {
         solver          diagonal;
     }
 
-    p
+    "p.*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -30,22 +30,23 @@ solvers
         relTol          0;
     }
 
-    "(U|e)"
+    "(U|e).*"
     {
         $p;
         tolerance       1e-9;
     }
 
-    "(k|epsilon)"
+    "(k|epsilon).*"
     {
         $p;
         tolerance       1e-10;
     }
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors     2;
+    nOuterCorrectors 1;
+    nCorrectors      2;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/compressible/sonicFoam/ras/prism/constant/polyMesh/boundary b/tutorials/compressible/sonicFoam/ras/prism/constant/polyMesh/boundary
index 635f185a35776b3a5631fe2c030ec6a810b4eb5e..0b093a1e0f3b727f85b1763c7da717e397dd025f 100644
--- a/tutorials/compressible/sonicFoam/ras/prism/constant/polyMesh/boundary
+++ b/tutorials/compressible/sonicFoam/ras/prism/constant/polyMesh/boundary
@@ -50,6 +50,7 @@ FoamFile
     defaultFaces
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          13272;
         startFace       13463;
     }
diff --git a/tutorials/compressible/sonicFoam/ras/prism/system/fvSolution b/tutorials/compressible/sonicFoam/ras/prism/system/fvSolution
index 3f46bc0a2508c392334e84ad38e3884d205739d6..56886cc6e1f7050ebf816607205bbb01cbefb472 100644
--- a/tutorials/compressible/sonicFoam/ras/prism/system/fvSolution
+++ b/tutorials/compressible/sonicFoam/ras/prism/system/fvSolution
@@ -17,12 +17,12 @@ FoamFile
 
 solvers
 {
-    rho
+    "rho.*"
     {
         solver          diagonal;
     }
 
-    p
+    "p.*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -30,22 +30,23 @@ solvers
         relTol          0;
     }
 
-    "(U|e|R)"
+    "(U|e|R).*"
     {
         $p;
         tolerance       1e-05;
     }
 
-    "(k|epsilon)"
+    "(k|epsilon).*"
     {
         $p;
         tolerance       1e-08;
     }
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors     2;
+    nOuterCorrectors 2;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/compressible/sonicLiquidFoam/decompressionTank/constant/polyMesh/boundary b/tutorials/compressible/sonicLiquidFoam/decompressionTank/constant/polyMesh/boundary
index add55d69d69695034b15d2ce266bbfdc08e97fcd..b1ca7310c67aa271b4b8a5f612db06c7051f6512 100644
--- a/tutorials/compressible/sonicLiquidFoam/decompressionTank/constant/polyMesh/boundary
+++ b/tutorials/compressible/sonicLiquidFoam/decompressionTank/constant/polyMesh/boundary
@@ -26,6 +26,7 @@ FoamFile
     axis
     {
         type            symmetryPlane;
+        inGroups        1(symmetryPlane);
         nFaces          120;
         startFace       7500;
     }
@@ -38,12 +39,14 @@ FoamFile
     back
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3725;
         startFace       7625;
     }
     front
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3725;
         startFace       11350;
     }
diff --git a/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSolution b/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSolution
index f63bb3a65b7c73fac12c2fb369f1d84e8287b896..5f012bcea78e2fd716264667f459ac4b97ba0454 100644
--- a/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSolution
+++ b/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    "p.*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -25,14 +25,14 @@ solvers
         relTol          0;
     }
 
-    U
+    "U.*"
     {
         $p;
         tolerance       1e-05;
         relTol          0;
     }
 
-    rho
+    "rho.*"
     {
         solver          PCG;
         preconditioner  DIC;
@@ -41,9 +41,10 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors     2;
+    nOuterCorrectors 2;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/constant/thermophysicalProperties b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/constant/thermophysicalProperties
index 05bbb78407326dfedc574ed0d2d52d51cc7bbb75..410aa5f42a635b71b1594cbe0a3dbdcca2d07d76 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/constant/thermophysicalProperties
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/constant/thermophysicalProperties
@@ -15,7 +15,8 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      heRhoThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleEnthalpy>>>>;
+thermoType heRhoThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleEnthalpy>>>>;
+//thermoType heRhoThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 pRef            100000;
 
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
index a6897b72d5203bc01d9d3ef325c0606a8b79f239..9ca38f0376fc13b4591e7891231c1637df0b7be5 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
@@ -30,23 +30,19 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phi,h)      Gauss upwind;
+    div(phi,e)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(phi,K)      Gauss linear;
+    div(phi,Ekp)    Gauss linear;
     div(R)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(Dp,p_rgh) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
index 42d96d946fb12aceb56cca571d0928b31fc967fc..7127fdc8c6fa54fac6586cea7bf6ae1bcd6aa718 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
@@ -39,7 +39,7 @@ solvers
         relTol          0;
     }
 
-    "(U|h|k|epsilon|R)"
+    "(U|h|e|k|epsilon|R)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -47,7 +47,7 @@ solvers
         relTol          0.1;
     }
 
-    "(U|h|k|epsilon|R)Final"
+    "(U|h|e|k|epsilon|R)Final"
     {
         $U;
         relTol          0;
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
index 599b1530bc87becb9c2a370a29e1f862763a5744..c52c053c84844fb84137797f741b112943bcd03a 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSchemes
@@ -27,12 +27,12 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinear 0.2;
-    div(phi,K)      Gauss limitedLinear 0.2;
-    div(phi,h)      Gauss limitedLinear 0.2;
-    div(phi,k)      Gauss limitedLinear 0.2;
-    div(phi,epsilon) Gauss limitedLinear 0.2;
-    div(phi,omega) Gauss limitedLinear 0.2;
+    div(phi,U)      bounded Gauss limitedLinear 0.2;
+    div(phi,K)      bounded Gauss limitedLinear 0.2;
+    div(phi,h)      bounded Gauss limitedLinear 0.2;
+    div(phi,k)      bounded Gauss limitedLinear 0.2;
+    div(phi,epsilon) bounded Gauss limitedLinear 0.2;
+    div(phi,omega) bounded Gauss limitedLinear 0.2;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
index e24d69c2022ae1a42c3d805c4c58a86dd383b0cd..1cee83b6cf85af9d1fdbc2e3af652bf2f7e0f607 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
@@ -15,32 +15,32 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 1 -1 0 0 0 0 ];
+dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform ( 0.1 0 0 );
+internalField   uniform (0.1 0 0);
 
 boundaryField
 {
     floor
     {
         type            fixedValue;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
     }
     ceiling
     {
         type            fixedValue;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
     }
     inlet
     {
         type            fixedValue;
-        value           uniform ( 0.1 0 0 );
+        value           uniform (0.1 0 0);
     }
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform ( 0 0 0 );
-        value           uniform ( 0 0 0 );
+        inletValue      uniform (0 0 0);
+        value           uniform (0.1 0 0);
     }
     fixedWalls
     {
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
index 48ad5decf746b4bf479e124086474b778f2ca070..406f4e1f0ac403885f8bae4aa9358a1f9a22923a 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
@@ -44,6 +44,7 @@ FoamFile
     fixedWalls
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          4000;
         startFace       4062;
     }
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict
index 2b5f43bcd4973c693a5b682ab0623928e716e44f..0bc3e79b5c915e69f696e9393d0eb79939dcbd1d 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict
@@ -77,7 +77,7 @@ dictionaryReplacement
         {
             "baffle.*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 0;
             }
         }
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffle b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffle
index dcdadf9491ee2761e752752ef5988cb077beb038..15a97c9a89318c729ce0bf5f6c883471ca38c0fb 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffle
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffle
@@ -77,7 +77,7 @@ dictionaryReplacement
         {
             "baffle1.*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           $internalField;
             }
         }
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffleRegion b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffleRegion
index 5ab7e81015a19176e0ca2b696257eb25351aa47b..5657d6092f1996084b37cc79b2efdcf2dbfc7c12 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffleRegion
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffleRegion
@@ -78,7 +78,7 @@ dictionaryReplacement
         {
             "region0_to_.*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           $internalField;
             }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
index 46f978cc82b4be057ce6ea6e09b616c27cbf12c2..573f34703a9caa6d446cff6e25b6aea3aa420249 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes
@@ -27,12 +27,12 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinear 0.2;
-    div(phi,K)      Gauss limitedLinear 0.2;
-    div(phi,h)      Gauss limitedLinear 0.2;
-    div(phi,k)      Gauss limitedLinear 0.2;
-    div(phi,epsilon) Gauss limitedLinear 0.2;
-    div(phi,omega) Gauss limitedLinear 0.2;
+    div(phi,U)      bounded Gauss limitedLinear 0.2;
+    div(phi,K)      bounded Gauss limitedLinear 0.2;
+    div(phi,h)      bounded Gauss limitedLinear 0.2;
+    div(phi,k)      bounded Gauss limitedLinear 0.2;
+    div(phi,epsilon) bounded Gauss limitedLinear 0.2;
+    div(phi,omega) bounded Gauss limitedLinear 0.2;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/g b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/g
index e0ac2653b5b370ad62f6770588121d30cac51627..037e3b47f3c4f7d7527b12ae4005e2a790edcc56 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/g
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/g
@@ -16,7 +16,6 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 1 -2 0 0 0 0];
-value           ( 0 -9.81 0 );
-
+value           (0 -9.81 0);
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/thermophysicalProperties b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/thermophysicalProperties
index 521cf3b4216fb9f3ef8e04176f1e98449d542777..885c152f4f85ec016999a5e6b84b2c4d1118d26b 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/thermophysicalProperties
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleEnthalpy>>>>;
+thermoType hePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleEnthalpy>>>>;
 
 pRef            100000;
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
index 991a5277f6e9ba1bdaff130ed0beee8c1d75ddfb..62bdbb235bf8d1608607b33ca57e82d86c674792 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
@@ -28,25 +28,17 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      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(R)          Gauss linear;
+    div(phi,U)      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((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear uncorrected;
-    laplacian(Dp,p_rgh) Gauss linear uncorrected;
-    laplacian(alphaEff,h) Gauss linear uncorrected;
-    laplacian(DkEff,k) Gauss linear uncorrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;
-    laplacian(DREff,R) Gauss linear uncorrected;
+    default         Gauss linear uncorrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
index c80f4db413b0c2e6e620077102bec99bcf7a9e1b..f7365e5f6d75e9085d3d5a5d407236b25df92bd5 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
@@ -21,15 +21,15 @@ solvers
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-08;
+        tolerance       1e-8;
         relTol          0.01;
     }
 
-    "(U|h|k|epsilon|R)"
+    "(U|h|e|k|epsilon|R)"
     {
         solver          PBiCG;
         preconditioner  DILU;
-        tolerance       1e-05;
+        tolerance       1e-6;
         relTol          0.1;
     }
 }
@@ -45,6 +45,7 @@ SIMPLE
         p_rgh           1e-2;
         U               1e-3;
         h               1e-3;
+        e               1e-3;
 
         // possibly check turbulence fields
         "(k|epsilon|omega)" 1e-3;
@@ -60,8 +61,8 @@ relaxationFactors
     }
     equations
     {
-        U               0.2;
-        h               0.2;
+        U               0.3;
+        "(h|e)"         0.3;
         "(k|epsilon|R)" 0.5;
     }
 }
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
index 3b6deefb697ca2d9996d39bb948132404efc138c..2ceace5b7c74faf00665a93de5b102db1657d533 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
@@ -28,26 +28,17 @@ 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,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
-    div(R)          Gauss linear;
+    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,epsilon) bounded Gauss upwind;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(Dp,p_rgh) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(gammaRad,G) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
index 817fe7274b45e93858f680d4dd2de0884c4e9190..c36da33bd6a4d2ef1b173782b019c72e458d22d6 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
@@ -28,27 +28,18 @@ 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,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
-    div(R)          Gauss linear;
+    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,epsilon) bounded Gauss upwind;
     div(Ji,Ii_h)    Gauss linearUpwind grad(Ii_h);
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(Dp,p_rgh) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(gammaRad,G) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
index df3988847815c66cdee90c06c2a47ee5b5f4272e..67aee46ff314c7d2040712a8bf1d6292e01ae4e9 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -100,7 +100,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
index c1a4e6bf16904cbea90ad006193bbbac231a07d0..118c410a839fda2b79747de2d2fb0e915ae333ff 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
@@ -135,7 +135,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict
index 1a7f5ecdd54c364ca6f18c8657be819c6485d5de..f403d950dc8ef39a68ca43bd8f938851f5f0f33a 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict
@@ -149,7 +149,7 @@ dictionaryReplacement
 
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 0;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict.save b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict.save
index 797ef03aa283f1f94ce3c2b2a68c293d0e64de03..d6817628c621dbb7ccba120c825220251d5f9571 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict.save
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/bottomWater/changeDictionaryDict.save
@@ -100,7 +100,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/changeDictionaryDict
index 18b85ea2e9ec25e3ba690c48c50a105db78370cb..8bb4a2ed547a7ed5f7548ff789ab4e0148c7f5ad 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater/system/topAir/changeDictionaryDict
@@ -136,7 +136,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
index 7abf71ca46718d2f4c250543ffa8ad993c4e918e..3ccca13960c4e159c7194b420960f347b5671c57 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -102,7 +102,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
index c1a4e6bf16904cbea90ad006193bbbac231a07d0..118c410a839fda2b79747de2d2fb0e915ae333ff 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
@@ -135,7 +135,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
index df3988847815c66cdee90c06c2a47ee5b5f4272e..67aee46ff314c7d2040712a8bf1d6292e01ae4e9 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -100,7 +100,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
index c1a4e6bf16904cbea90ad006193bbbac231a07d0..118c410a839fda2b79747de2d2fb0e915ae333ff 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
@@ -135,7 +135,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/changeDictionaryDict
index c9197f1b246665ed83b4152d21e81e6438ca9252..8bb6b439411fba5a903e2714347fbf5e222e3897 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/bottomAir/changeDictionaryDict
@@ -97,7 +97,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
         }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/changeDictionaryDict
index 704d77ce7005ba08e0141c1bc5610b03ff359588..cff14264b6b5b21ec67a2d4d8414c02b52ba141d 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/system/topAir/changeDictionaryDict
@@ -132,7 +132,7 @@ dictionaryReplacement
         {
             ".*"
             {
-                type            buoyantPressure;
+                type            fixedFluxPressure;
                 value           uniform 1e5;
             }
 
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/p_rgh.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/p_rgh.org
index 79523d2e245da4bd7eccf85f4042b2f899721752..396d380daddd7bc4e4a3c9d427a5226c6c3f44af 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/p_rgh.org
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/p_rgh.org
@@ -22,7 +22,7 @@ boundaryField
 {
     walls
     {
-        type            buoyantPressure;
+        type            fixedFluxPressure;
         value           uniform 1e5;
     }
 
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/p_rgh.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/p_rgh.org
index 15bc97e933c2bfab9a6c45565fc1d8a89735aa94..0394387374d5ef9bb77528165e097f872381d9a4 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/p_rgh.org
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/p_rgh.org
@@ -22,7 +22,7 @@ boundaryField
 {
     walls
     {
-        type            buoyantPressure;
+        type            fixedFluxPressure;
         value           uniform 1e5;
     }
 }
diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/0/p_rgh b/tutorials/multiphase/interFoam/ras/waterChannel/0/p_rgh
index d30fb854876908815e5272d4b9383431b9a65684..d2958fafeba6f4d7d01b1649c3b6e39baea65f4a 100644
--- a/tutorials/multiphase/interFoam/ras/waterChannel/0/p_rgh
+++ b/tutorials/multiphase/interFoam/ras/waterChannel/0/p_rgh
@@ -34,7 +34,7 @@ boundaryField
 
     ".*"
     {
-        type            buoyantPressure;
+        type            fixedFluxPressure;
         value           uniform 0;
     }
 }