diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
index 24316a168545717b0785e83fcfa8a6fc01d9ae81..e916b7a56aac74e5c48380b3f70a3850cd4fcda5 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
@@ -1,30 +1,30 @@
+if (mesh.changing())
 {
-    if (mesh.changing())
+    forAll(U.boundaryField(), patchi)
     {
-        forAll(U.boundaryField(), patchi)
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].initEvaluate();
-            }
+            U.boundaryField()[patchi].initEvaluate();
         }
+    }
 
-        forAll(U.boundaryField(), patchi)
+    forAll(U.boundaryField(), patchi)
+    {
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].evaluate();
+            U.boundaryField()[patchi].evaluate();
 
-                phi.boundaryField()[patchi] =
-                    rho.boundaryField()[patchi]
-                   *(
-                       U.boundaryField()[patchi]
-                     & mesh.Sf().boundaryField()[patchi]
-                    );
-            }
+            phi.boundaryField()[patchi] =
+                rho.boundaryField()[patchi]
+               *(
+                   U.boundaryField()[patchi]
+                 & mesh.Sf().boundaryField()[patchi]
+               );
         }
     }
+}
 
+{
     volScalarField pcorr
     (
         IOobject
@@ -40,13 +40,13 @@
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divrhoU
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divrhoU
         );
 
         pcorrEqn.solve();
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
index ebc62e7e592ae319d8019c409d59b26295d4839c..b693bb78607f7db13f95a818edb33c10a4d4c1d2 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
@@ -37,13 +37,13 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi)
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -54,6 +54,6 @@ if (mesh.changing())
             phi -= pcorrEqn.flux();
         }
     }
-}
 
-#include "continuityErrs.H"
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
index 2dca530c145434bf062f6ff8cca9192ec1617af8..c3fde6b3a1a2558a052e6fe1cb1175253b593bff 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
@@ -47,16 +47,15 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readThermodynamicProperties.H"
-    #include "readControls.H"
-    #include "createFields.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
 
-    surfaceScalarField phivAbs("phivAbs", phiv);
-    fvc::makeAbsolute(phivAbs, U);
-
+    #include "readThermodynamicProperties.H"
+    #include "readControls.H"
+    #include "createFields.H"
+    #include "createUf.H"
+    #include "createPcorrTypes.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -75,18 +74,8 @@ int main(int argc, char *argv[])
 
         scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        {
-            // Calculate the relative velocity used to map relative flux phiv
-            volVectorField Urel("Urel", U);
-
-            if (mesh.moving())
-            {
-                Urel -= fvc::reconstruct(fvc::meshPhi(U));
-            }
-
-            // Do any mesh changes
-            mesh.update();
-        }
+        // Do any mesh changes
+        mesh.update();
 
         if (mesh.changing())
         {
@@ -94,7 +83,16 @@ int main(int argc, char *argv[])
                 << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                 << " s" << endl;
 
-            #include "correctPhi.H"
+            if (correctPhi)
+            {
+                // Calculate absolute flux from the mapped surface velocity
+                phiv = mesh.Sf() & Uf;
+
+                #include "correctPhi.H"
+
+                // Make the flux relative to the mesh motion
+                fvc::makeRelative(phiv, U);
+            }
         }
 
         // --- Pressure-velocity PIMPLE corrector loop
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
index 8115b791e9c876108a87c3556df04ad9ba985437..c05d664838b445d2f1107d5fa9de30b38ea54324 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
@@ -1,18 +1,27 @@
+if (mesh.changing())
 {
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
+    forAll(U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].fixesValue())
+        {
+            U.boundaryField()[patchI].initEvaluate();
+        }
+    }
 
-    forAll (p.boundaryField(), i)
+    forAll(U.boundaryField(), patchI)
     {
-        if (p.boundaryField()[i].fixesValue())
+        if (U.boundaryField()[patchI].fixesValue())
         {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+            U.boundaryField()[patchI].evaluate();
+
+            phiv.boundaryField()[patchI] =
+                U.boundaryField()[patchI]
+              & mesh.Sf().boundaryField()[patchI];
         }
     }
+}
 
+{
     volScalarField pcorr
     (
         IOobject
@@ -29,7 +38,7 @@
     );
 
     surfaceScalarField rhof(fvc::interpolate(rho, "div(phiv,rho)"));
-    dimensionedScalar rAUf("(1|A(U))", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index e7f078d1a6578be726065c128ede11fdacf66a44..0f0d5768414323fb1df135a80b7bd1b5e27d2832 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -18,7 +18,7 @@
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + rhorAUf*fvc::ddtCorr(U, phivAbs);
+         + rhorAUf*fvc::ddtCorr(U, Uf);
     fvc::makeRelative(phiv, U);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
@@ -43,7 +43,6 @@
         if (pimple.finalNonOrthogonalIter())
         {
             phiv += (phiGradp + pEqn.flux())/rhof;
-            phivAbs = fvc::absolute(phiv, U);
         }
     }
 
@@ -82,4 +81,10 @@
     U.correctBoundaryConditions();
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
+
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += mesh.Sf()*(phiv - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..216e9b4247aa70fc76afda34f80f7310118d6035
--- /dev/null
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H
@@ -0,0 +1,9 @@
+#include "readTimeControls.H"
+
+scalar maxAcousticCo
+(
+    readScalar(runTime.controlDict().lookup("maxAcousticCo"))
+);
+
+bool correctPhi =
+    pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index a9c3365bdf73fedfa445979e9cf7e5b2f2ebf937..06444163304bdae83032ea072608fe5459b896ac 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -56,12 +56,13 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "readGravitationalAcceleration.H"
+    #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
 
-    #include "readControls.H"
-    #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createUf.H"
+    #include "readControls.H"
     #include "createPrghCorrTypes.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -74,9 +75,6 @@ int main(int argc, char *argv[])
         #include "readControls.H"
         #include "CourantNo.H"
 
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
         #include "setDeltaT.H"
 
         runTime++;
@@ -85,7 +83,7 @@ int main(int argc, char *argv[])
 
         {
             // Store divU from the previous mesh for the correctPhi
-            volScalarField divU(fvc::div(phi));
+            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
 
             scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
@@ -104,13 +102,16 @@ int main(int argc, char *argv[])
 
             if (mesh.changing() && correctPhi)
             {
+                // Calculate absolute flux from the mapped surface velocity
+                phi = mesh.Sf() & Uf;
+
                 #include "correctPhi.H"
+
+                // Make the fluxes relative to the mesh motion
+                fvc::makeRelative(phi, U);
             }
         }
 
-        // Make the fluxes relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
         if (mesh.changing() && checkMeshCourantNo)
         {
             #include "meshCourantNo.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
index c8709d33476c0e4e8425583ac33d369e8caba7be..97c63e6a3eeee94897f27706be16ea18315cdabe 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
@@ -1,28 +1,26 @@
+if (mesh.changing())
 {
-    if (mesh.changing())
+    forAll(U.boundaryField(), patchi)
     {
-        forAll(U.boundaryField(), patchi)
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].initEvaluate();
-            }
+            U.boundaryField()[patchi].initEvaluate();
         }
+    }
 
-        forAll(U.boundaryField(), patchi)
+    forAll(U.boundaryField(), patchi)
+    {
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].evaluate();
+            U.boundaryField()[patchi].evaluate();
 
-                phi.boundaryField()[patchi] =
-                U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
-            }
+            phi.boundaryField()[patchi] =
+            U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
         }
     }
+}
 
-    #include "continuityErrs.H"
-
+{
     volScalarField pcorr
     (
         IOobject
@@ -38,15 +36,13 @@
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
-
-    adjustPhi(phi, U, pcorr);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divU
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
         );
 
         pcorrEqn.solve();
@@ -56,8 +52,8 @@
             phi -= pcorrEqn.flux();
         }
     }
-}
 
-phi.oldTime() = phi;
+    phi.oldTime() = phi;
 
-#include "continuityErrs.H"
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..c86921371f52ec96f0cff265257110280f460546
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -0,0 +1,133 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
+    );
+
+    surfaceScalarField phig
+    (
+        (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    tmp<fvScalarMatrix> p_rghEqnComp1;
+    tmp<fvScalarMatrix> p_rghEqnComp2;
+
+    if (pimple.transonic())
+    {
+        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
+        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
+
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+          + correction
+            (
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
+        p_rghEqnComp1().relax();
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+          + correction
+            (
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
+        p_rghEqnComp2().relax();
+    }
+    else
+    {
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
+    }
+
+    // Cache p_rgh prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        solve
+        (
+            (
+                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
+              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
+            )
+          + p_rghEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+
+            dgdt =
+            (
+                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+            );
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    }
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    // Update densities from change in p_rgh
+    rho1 += psi1*(p_rgh - p_rgh_0);
+    rho2 += psi2*(p_rgh - p_rgh_0);
+
+    rho = alpha1*rho1 + alpha2*rho2;
+
+    K = 0.5*magSqr(U);
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+}
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index d1e164e3cb4f5c869a95b10fe32d68df6755ec6b..980d1be6306103ff2aea735fdd2cb49eca276ce3 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialrDeltaT.H"
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index 670ce727fda8e3457badeb753564cb33eafd2c97..887ef9e34b5f032d62e730a6fd9faae6a97f0c32 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/correctPhi.H
index a2a61cef21dffc6dd15e0d2538817498c6ecdc54..0d324d1357ddf6c827a878d582863ab1c7276248 100644
--- a/applications/solvers/multiphase/interFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/correctPhi.H
@@ -1,20 +1,6 @@
 {
     #include "continuityErrs.H"
 
-    wordList pcorrTypes
-    (
-        p_rgh.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll (p_rgh.boundaryField(), i)
-    {
-        if (p_rgh.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
     volScalarField pcorr
     (
         IOobject
@@ -30,7 +16,7 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     adjustPhi(phi, U, pcorr);
 
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
index deb4edc3881a0240c048fcf28ad308894dc4d6a9..53d04d827f1975ab033efb99498a06f267875f97 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
@@ -37,13 +37,13 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi)
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -54,6 +54,6 @@ if (mesh.changing())
             phi -= pcorrEqn.flux();
         }
     }
-}
 
-#include "continuityErrs.H"
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index fbcbd7c7e4a34864a96b1931a35ceb1ccfcce9e4..fa9b4fffe44c0407fbb51c4639001c272876c75a 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 1a10e45f3fda491e7fb8495a2fda79e389aee125..93cb8a0c3de15eb5e42884b8356aa17019bd1d5f 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 4e6177ef9ad4680fef0369f685416f7d581af2ea..64938174dc06bec138163a2c160a7a176ed2a102 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index 3dd2b50fed60412d4284d3e542ff3ce2d07023fc..c1d13ca92645a04f75a33c5ce0245e5689461103 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "../interFoam/correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
index 57a79b6091cbf34f3b01156f1a61a311a08aad76..e291b7dbcd263dcaa18d55a7336e4fedc999fb9e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
@@ -1,20 +1,6 @@
 {
     #include "continuityErrs.H"
 
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll (p.boundaryField(), i)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
     volScalarField pcorr
     (
         IOobject
@@ -30,9 +16,7 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    //adjustPhi(phi, U, pcorr);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index 0c5c8a2a86440d8c6afd828fd1bd95b40919adf9..aacc4b6d498c6ed237edb6b656107dd184ab8e57 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
     #include "createMRFZones.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
+    #include "createPcorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index bdcb8105810cca9ca8e9046ef443fb3cf6503264..8c73d798f6e9222a2049b540c0a2bf7bfc7c2afc 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createMRFZones.H"
     #include "readTimeControls.H"
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index 9bfeef4674584a5c3fc6f3a0d30b54a524a31e8e..d2c565b915c6d8e5fa8eee67c274be58694b947b 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.C b/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.C
index b125351314b205acd0ee6e1406954230af58d917..271f334498880fd52b8fd1fd954c73992d3a525c 100644
--- a/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.C
+++ b/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.C
@@ -31,6 +31,118 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::shortEdgeFilter2D::addRegion
+(
+    const label regionI,
+    DynamicList<label>& bPointRegions
+) const
+{
+    if (bPointRegions.empty())
+    {
+        bPointRegions.append(regionI);
+    }
+    else if (findIndex(bPointRegions, regionI) == -1)
+    {
+        bPointRegions.append(regionI);
+    }
+}
+
+
+void Foam::shortEdgeFilter2D::assignBoundaryPointRegions
+(
+    List<DynamicList<label> >& boundaryPointRegions
+) const
+{
+    forAllConstIter(EdgeMap<label>, mapEdgesRegion_, iter)
+    {
+        const edge& e = iter.key();
+        const label& regionI = iter();
+
+        const label startI = e.start();
+        const label endI = e.end();
+
+        addRegion(regionI, boundaryPointRegions[startI]);
+        addRegion(regionI, boundaryPointRegions[endI]);
+    }
+}
+
+
+void Foam::shortEdgeFilter2D::updateEdgeRegionMap
+(
+    const MeshedSurface<face>& surfMesh,
+    const List<DynamicList<label> >& boundaryPtRegions,
+    const labelList& surfPtToBoundaryPt,
+    EdgeMap<label>& mapEdgesRegion,
+    labelList& patchSizes
+) const
+{
+    EdgeMap<label> newMapEdgesRegion(mapEdgesRegion.size());
+
+    const edgeList& edges = surfMesh.edges();
+    const labelList& meshPoints = surfMesh.meshPoints();
+
+    patchSizes.setSize(patchNames_.size(), 0);
+    patchSizes = 0;
+
+    forAll(edges, edgeI)
+    {
+        if (surfMesh.isInternalEdge(edgeI))
+        {
+            continue;
+        }
+
+        const edge& e = edges[edgeI];
+
+        const label startI = meshPoints[e[0]];
+        const label endI = meshPoints[e[1]];
+
+        label region = -1;
+
+        const DynamicList<label> startPtRegions =
+            boundaryPtRegions[surfPtToBoundaryPt[startI]];
+        const DynamicList<label> endPtRegions =
+            boundaryPtRegions[surfPtToBoundaryPt[endI]];
+
+        if (startPtRegions.size() > 1 && endPtRegions.size() > 1)
+        {
+            region = startPtRegions[0];
+
+            WarningIn("shortEdgeFilter2D()")
+                << "Both points in edge are in different regions."
+                << " Assigning edge to region " << region
+                << endl;
+        }
+        else if (startPtRegions.size() > 1 || endPtRegions.size() > 1)
+        {
+            region =
+            (
+                startPtRegions.size() > 1
+              ? endPtRegions[0]
+              : startPtRegions[0]
+            );
+        }
+        else if
+        (
+            startPtRegions[0] == endPtRegions[0]
+         && startPtRegions[0] != -1
+        )
+        {
+            region = startPtRegions[0];
+        }
+
+        if (region != -1)
+        {
+            newMapEdgesRegion.insert(e, region);
+            patchSizes[region]++;
+        }
+    }
+
+    mapEdgesRegion.transfer(newMapEdgesRegion);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::shortEdgeFilter2D::shortEdgeFilter2D
@@ -138,19 +250,12 @@ Foam::shortEdgeFilter2D::filter()
 
     // Check if the point is a boundary point. Flag if it is so that
     // it will not be deleted.
-    boolList boundaryPointFlags(points.size(), false);
-    // This has been removed, otherwise small edges on the boundary are not
-    // removed.
-    /*  forAll(boundaryPointFlags, pointI)
-    {
-        forAll(boundaryPoints, bPoint)
-        {
-            if (meshPoints[boundaryPoints[bPoint]] == pointI)
-            {
-                boundaryPointFlags[pointI] = true;
-            }
-        }
-    }*/
+    List<DynamicList<label> > boundaryPointRegions
+    (
+        points.size(),
+        DynamicList<label>()
+    );
+    assignBoundaryPointRegions(boundaryPointRegions);
 
     // Check if an edge has a boundary point. It it does the edge length
     // will be doubled when working out its length.
@@ -232,7 +337,6 @@ Foam::shortEdgeFilter2D::filter()
             shortEdgeFilterFactor_
            /(psEdges.size() + peEdges.size() - 2);
 
-
         edge lookupInPatchEdge
         (
             meshPoints[startVertex],
@@ -283,11 +387,38 @@ Foam::shortEdgeFilter2D::filter()
             (
                 pointsToRemove[meshPoints[startVertex]] == -1
              && pointsToRemove[meshPoints[endVertex]] == -1
-             && !boundaryPointFlags[meshPoints[startVertex]]
              && !flagDegenerateFace
             )
             {
-                pointsToRemove[meshPoints[startVertex]] = meshPoints[endVertex];
+                const DynamicList<label>& startVertexRegions =
+                    boundaryPointRegions[meshPoints[startVertex]];
+                const DynamicList<label>& endVertexRegions =
+                    boundaryPointRegions[meshPoints[endVertex]];
+
+                if (startVertexRegions.size() && endVertexRegions.size())
+                {
+                    if (startVertexRegions.size() > 1)
+                    {
+                        pointsToRemove[meshPoints[endVertex]] =
+                            meshPoints[startVertex];
+                    }
+                    else
+                    {
+                        pointsToRemove[meshPoints[startVertex]] =
+                            meshPoints[endVertex];
+                    }
+                }
+                else if (startVertexRegions.size())
+                {
+                    pointsToRemove[meshPoints[endVertex]] =
+                        meshPoints[startVertex];
+                }
+                else
+                {
+                    pointsToRemove[meshPoints[startVertex]] =
+                        meshPoints[endVertex];
+                }
+
                 ++nPointsToRemove;
             }
         }
@@ -299,6 +430,9 @@ Foam::shortEdgeFilter2D::filter()
     labelList newPointNumbers(points.size(), -1);
     label numberRemoved = 0;
 
+    // Maintain addressing from new to old point field
+    labelList newPtToOldPt(totalNewPoints, -1);
+
     forAll(points, pointI)
     {
         // If the point is NOT going to be removed.
@@ -306,6 +440,7 @@ Foam::shortEdgeFilter2D::filter()
         {
             newPoints[pointI - numberRemoved] = points[pointI];
             newPointNumbers[pointI] =  pointI - numberRemoved;
+            newPtToOldPt[pointI - numberRemoved] = pointI;
         }
         else
         {
@@ -405,16 +540,14 @@ Foam::shortEdgeFilter2D::filter()
         xferCopy(List<surfZone>())
     );
 
-    const Map<int>& fMeshPointMap = fMesh.meshPointMap();
-
-    // Reset patchSizes_
-    patchSizes_.clear();
-    patchSizes_.setSize(patchNames_.size(), 0);
-
-    label equalEdges = 0;
-    label notFound = 0;
-    label matches = 0;
-    label negativeLabels = 0;
+    updateEdgeRegionMap
+    (
+        fMesh,
+        boundaryPointRegions,
+        newPtToOldPt,
+        mapEdgesRegion_,
+        patchSizes_
+    );
 
     forAll(newPointNumbers, pointI)
     {
@@ -427,78 +560,12 @@ Foam::shortEdgeFilter2D::filter()
         }
     }
 
-    // Create new EdgeMap.
-    Info<< "Creating new EdgeMap." << endl;
-    EdgeMap<label> newMapEdgesRegion(mapEdgesRegion_.size());
-
-    for
-    (
-        label bEdgeI = ms_.nInternalEdges();
-        bEdgeI < edges.size();
-        ++bEdgeI
-    )
-    {
-        label p1 = meshPoints[edges[bEdgeI][0]];
-        label p2 = meshPoints[edges[bEdgeI][1]];
-
-        edge e(p1, p2);
-
-        if (mapEdgesRegion_.found(e))
-        {
-            if
-            (
-                newPointNumbers[p1] != -1
-             && newPointNumbers[p2] != -1
-            )
-            {
-                if (newPointNumbers[p1] != newPointNumbers[p2])
-                {
-                    label region = mapEdgesRegion_.find(e)();
-                    newMapEdgesRegion.insert
-                    (
-                        edge
-                        (
-                            fMeshPointMap[newPointNumbers[p1]],
-                            fMeshPointMap[newPointNumbers[p2]]
-                        ),
-                        region
-                    );
-                    patchSizes_[region]++;
-                    matches++;
-                }
-                else
-                {
-                    equalEdges++;
-                }
-            }
-            else
-            {
-                negativeLabels++;
-            }
-        }
-        else
-        {
-            notFound++;
-        }
-    }
-
-    if (debug)
-    {
-        Info<< "EdgeMapping  :" << nl
-            << "    Matches  : " << matches << nl
-            << "    Equal    : " << equalEdges << nl
-            << "    Negative : " << negativeLabels << nl
-            << "    Not Found: " << notFound << endl;
-    }
-
-    mapEdgesRegion_.transfer(newMapEdgesRegion);
-
     ms_.transfer(fMesh);
 
-    Info<< "    Maximum number of chained collapses = " << maxChain << endl;
-
     if (debug)
     {
+        Info<< "    Maximum number of chained collapses = " << maxChain << endl;
+
         writeInfo(Info);
     }
 }
diff --git a/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.H b/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.H
index 53b248fc3d4bf542676f3a1f64d9cc5db9f3664a..df52e25741b9fdba9f522bb5384e9b0446675a89 100644
--- a/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.H
+++ b/applications/utilities/mesh/generation/foamyQuadMesh/shortEdgeFilter2D.H
@@ -70,6 +70,26 @@ class shortEdgeFilter2D
 
     // Private Member Functions
 
+        void addRegion
+        (
+            const label regionI,
+            DynamicList<label>& bPointRegions
+        ) const;
+
+        void assignBoundaryPointRegions
+        (
+            List<DynamicList<label> >& boundaryPointRegions
+        ) const;
+
+        void updateEdgeRegionMap
+        (
+            const MeshedSurface<face>& surfMesh,
+            const List<DynamicList<label> >& boundaryPtRegions,
+            const labelList& surfPtToBoundaryPt,
+            EdgeMap<label>& mapEdgesRegion,
+            labelList& patchSizes
+        ) const;
+
         //- Disallow default bitwise copy construct
         shortEdgeFilter2D(const shortEdgeFilter2D&);
 
diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
index c7ad881527ae51d82f42eedbe579412fc80abd7e..63d91ac1788608a0e09c3813511dd664e3579d75 100644
--- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
+++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
@@ -186,10 +186,20 @@ int main(int argc, char *argv[])
 
     const word setName = args[1];
 
-    word resultInstance = mesh.pointsInstance();
-    const bool specifiedInstance =
-        args.optionFound("overwrite")
-     || args.optionReadIfPresent("resultTime", resultInstance);
+    word meshInstance = mesh.pointsInstance();
+    word fieldsInstance = runTime.timeName();
+
+    const bool overwrite = args.optionFound("overwrite");
+    const bool specifiedInstance = args.optionReadIfPresent
+    (
+        "resultTime",
+        fieldsInstance
+    );
+    if (specifiedInstance)
+    {
+        // Set both mesh and field to this time
+        meshInstance = fieldsInstance;
+    }
 
 
     Info<< "Reading cell set from " << setName << endl << endl;
@@ -355,14 +365,14 @@ int main(int argc, char *argv[])
     // Write mesh and fields to new time
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (!specifiedInstance)
+    if (overwrite || specifiedInstance)
     {
-        runTime++;
+        runTime.setTime(instant(fieldsInstance), 0);
+        subsetter.subMesh().setInstance(meshInstance);
     }
     else
     {
-        runTime.setTime(instant(resultInstance), 0);
-        subsetter.subMesh().setInstance(runTime.timeName());
+        runTime++;
     }
 
     Info<< "Writing subsetted mesh and fields to time " << runTime.timeName()
@@ -374,31 +384,26 @@ int main(int argc, char *argv[])
     forAll(scalarFlds, i)
     {
         scalarFlds[i].rename(scalarNames[i]);
-
         scalarFlds[i].write();
     }
     forAll(vectorFlds, i)
     {
         vectorFlds[i].rename(vectorNames[i]);
-
         vectorFlds[i].write();
     }
     forAll(sphericalTensorFlds, i)
     {
         sphericalTensorFlds[i].rename(sphericalTensorNames[i]);
-
         sphericalTensorFlds[i].write();
     }
     forAll(symmTensorFlds, i)
     {
         symmTensorFlds[i].rename(symmTensorNames[i]);
-
         symmTensorFlds[i].write();
     }
     forAll(tensorFlds, i)
     {
         tensorFlds[i].rename(tensorNames[i]);
-
         tensorFlds[i].write();
     }
 
@@ -406,31 +411,26 @@ int main(int argc, char *argv[])
     forAll(surfScalarFlds, i)
     {
         surfScalarFlds[i].rename(surfScalarNames[i]);
-
         surfScalarFlds[i].write();
     }
     forAll(surfVectorFlds, i)
     {
         surfVectorFlds[i].rename(surfVectorNames[i]);
-
         surfVectorFlds[i].write();
     }
     forAll(surfSphericalTensorFlds, i)
     {
         surfSphericalTensorFlds[i].rename(surfSphericalTensorNames[i]);
-
         surfSphericalTensorFlds[i].write();
     }
     forAll(surfSymmTensorFlds, i)
     {
         surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
-
         surfSymmTensorFlds[i].write();
     }
     forAll(surfTensorNames, i)
     {
         surfTensorFlds[i].rename(surfTensorNames[i]);
-
         surfTensorFlds[i].write();
     }
 
@@ -438,31 +438,26 @@ int main(int argc, char *argv[])
     forAll(pointScalarFlds, i)
     {
         pointScalarFlds[i].rename(pointScalarNames[i]);
-
         pointScalarFlds[i].write();
     }
     forAll(pointVectorFlds, i)
     {
         pointVectorFlds[i].rename(pointVectorNames[i]);
-
         pointVectorFlds[i].write();
     }
     forAll(pointSphericalTensorFlds, i)
     {
         pointSphericalTensorFlds[i].rename(pointSphericalTensorNames[i]);
-
         pointSphericalTensorFlds[i].write();
     }
     forAll(pointSymmTensorFlds, i)
     {
         pointSymmTensorFlds[i].rename(pointSymmTensorNames[i]);
-
         pointSymmTensorFlds[i].write();
     }
     forAll(pointTensorNames, i)
     {
         pointTensorFlds[i].rename(pointTensorNames[i]);
-
         pointTensorFlds[i].write();
     }
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
index 47cd58d2869a340aa756df50b72985f2cf155c54..ae90704fd3524a6b4a51df1ef99c2326c04f2c90 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
@@ -356,25 +356,29 @@ void ensightField
 
     ensightStream& ensightFile = *ensightFilePtr;
 
+    if (Pstream::master())
+    {
+        if (timeIndex == 0)
+        {
+            ensightCaseFile.setf(ios_base::left);
+
+            ensightCaseFile
+                << ensightPTraits<Type>::typeName
+                << " per element:            1       "
+                << setw(15) << vf.name()
+                << (' ' + prepend + "****." + vf.name()).c_str()
+                << nl;
+        }
+
+        ensightFile.write(ensightPTraits<Type>::typeName);
+    }
+
     if (patchNames.empty())
     {
         eMesh.barrier();
 
         if (Pstream::master())
         {
-            if (timeIndex == 0)
-            {
-                ensightCaseFile.setf(ios_base::left);
-
-                ensightCaseFile
-                    << ensightPTraits<Type>::typeName
-                    << " per element:            1       "
-                    << setw(15) << vf.name()
-                    << (' ' + prepend + "****." + vf.name()).c_str()
-                    << nl;
-            }
-
-            ensightFile.write(ensightPTraits<Type>::typeName);
             ensightFile.writePartHeader(1);
         }
 
@@ -570,25 +574,29 @@ void ensightPointField
 
     ensightStream& ensightFile = *ensightFilePtr;
 
+    if (Pstream::master())
+    {
+        if (timeIndex == 0)
+        {
+            ensightCaseFile.setf(ios_base::left);
+
+            ensightCaseFile
+                << ensightPTraits<Type>::typeName
+                << " per node:            1       "
+                << setw(15) << pf.name()
+                << (' ' + prepend + "****." + pf.name()).c_str()
+                << nl;
+        }
+
+        ensightFile.write(ensightPTraits<Type>::typeName);
+    }
+
     if (eMesh.patchNames().empty())
     {
         eMesh.barrier();
 
         if (Pstream::master())
         {
-            if (timeIndex == 0)
-            {
-                ensightCaseFile.setf(ios_base::left);
-
-                ensightCaseFile
-                    << ensightPTraits<Type>::typeName
-                    << " per node:            1       "
-                    << setw(15) << pf.name()
-                    << (' ' + prepend + "****." + pf.name()).c_str()
-                    << nl;
-            }
-
-            ensightFile.write(ensightPTraits<Type>::typeName);
             ensightFile.writePartHeader(1);
         }
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index 6fab5a8b3b388e9eae62b85f3574751f488a581c..0332572eeb0d16cdb907354ac36fd9a139063803 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -260,7 +260,10 @@ void Foam::ensightMesh::correct()
     // faceZones
     if (faceZones_)
     {
-        const wordList faceZoneNamesAll = mesh_.faceZones().names();
+        wordList faceZoneNamesAll = mesh_.faceZones().names();
+        // Need to sort the list of all face zones since the index may vary
+        // from processor to processor...
+        sort(faceZoneNamesAll);
 
         // Find faceZone names which match that requested at command-line
         forAll(faceZoneNamesAll, nameI)
@@ -300,15 +303,16 @@ void Foam::ensightMesh::correct()
         // Count face types in each faceZone
         forAll(faceZoneNamesAll, zoneI)
         {
-            //const word& zoneName = faceZoneNamesAll[zoneI];
+            const word& zoneName = faceZoneNamesAll[zoneI];
+            const label faceZoneId = mesh_.faceZones().findZoneID(zoneName);
 
-            const faceZone& fz = mesh_.faceZones()[zoneI];
+            const faceZone& fz = mesh_.faceZones()[faceZoneId];
 
             if (fz.size())
             {
-                labelList& tris = faceZoneFaceSets_[zoneI].tris;
-                labelList& quads = faceZoneFaceSets_[zoneI].quads;
-                labelList& polys = faceZoneFaceSets_[zoneI].polys;
+                labelList& tris = faceZoneFaceSets_[faceZoneId].tris;
+                labelList& quads = faceZoneFaceSets_[faceZoneId].quads;
+                labelList& polys = faceZoneFaceSets_[faceZoneId].polys;
 
                 tris.setSize(fz.size());
                 quads.setSize(fz.size());
@@ -356,19 +360,20 @@ void Foam::ensightMesh::correct()
         {
             const word& zoneName = faceZoneNamesAll[zoneI];
             nFacePrimitives nfp;
+            const label faceZoneId = mesh_.faceZones().findZoneID(zoneName);
 
             if (faceZoneNames_.found(zoneName))
             {
                 if
                 (
-                    faceZoneFaceSets_[zoneI].tris.size()
-                 || faceZoneFaceSets_[zoneI].quads.size()
-                 || faceZoneFaceSets_[zoneI].polys.size()
+                    faceZoneFaceSets_[faceZoneId].tris.size()
+                 || faceZoneFaceSets_[faceZoneId].quads.size()
+                 || faceZoneFaceSets_[faceZoneId].polys.size()
                 )
                 {
-                    nfp.nTris   = faceZoneFaceSets_[zoneI].tris.size();
-                    nfp.nQuads  = faceZoneFaceSets_[zoneI].quads.size();
-                    nfp.nPolys  = faceZoneFaceSets_[zoneI].polys.size();
+                    nfp.nTris   = faceZoneFaceSets_[faceZoneId].tris.size();
+                    nfp.nQuads  = faceZoneFaceSets_[faceZoneId].quads.size();
+                    nfp.nPolys  = faceZoneFaceSets_[faceZoneId].polys.size();
                 }
             }
 
@@ -593,6 +598,7 @@ void Foam::ensightMesh::writePolysPoints
     const labelList& polys,
     const cellList& cellFaces,
     const faceList& faces,
+    const labelList& faceOwner,
     ensightStream& ensightGeometryFile
 ) const
 {
@@ -602,12 +608,46 @@ void Foam::ensightMesh::writePolysPoints
 
         forAll(cf, faceI)
         {
-            const face& f = faces[cf[faceI]];
+            const label faceId = cf[faceI];
+            const face& f = faces[faceId];  // points of face (in global points)
+            const label np = f.size();
+            bool reverseOrder = false;
+            if (faceId >= faceOwner.size())
+            {
+                // Boundary face.
+                // Nothing should be done for processor boundary.
+                // The current cell always owns them. Given that we
+                // are reverting the
+                // order when the cell is the neighbour to the face,
+                // the orientation of
+                // all the boundaries, no matter if they are "real"
+                // or processorBoundaries, is consistent.
+            }
+            else
+            {
+                if (faceOwner[faceId] != polys[i])
+                {
+                    reverseOrder = true;
+                }
+            }
 
-            List<int> temp(f.size());
+            // If the face owner is the current cell, write the points
+            // in the standard order.
+            // If the face owner is not the current cell, write the points
+            // in reverse order.
+            // EnSight prefers to have all the faces of an nfaced cell
+            // oriented in the same way.
+            List<int> temp(np);
             forAll(f, pointI)
             {
-                temp[pointI] = f[pointI] + 1;
+                if (reverseOrder)
+                {
+                    temp[np-1-pointI] = f[pointI] + 1;
+                }
+                else
+                {
+                    temp[pointI] = f[pointI] + 1;
+                }
             }
             ensightGeometryFile.write(temp);
         }
@@ -624,6 +664,8 @@ void Foam::ensightMesh::writeAllPolys
     if (meshCellSets_.nPolys)
     {
         const cellList& cellFaces = mesh_.cells();
+        const labelList& faceOwner = mesh_.faceOwner();
+
         // Renumber faces to use global point numbers
         faceList faces(mesh_.faces());
         forAll(faces, i)
@@ -714,6 +756,7 @@ void Foam::ensightMesh::writeAllPolys
                 meshCellSets_.polys,
                 cellFaces,
                 faces,
+                faceOwner,
                 ensightGeometryFile
             );
             // Slaves
@@ -723,12 +766,14 @@ void Foam::ensightMesh::writeAllPolys
                 labelList polys(fromSlave);
                 cellList cellFaces(fromSlave);
                 faceList faces(fromSlave);
+                labelList faceOwner(fromSlave);
 
                 writePolysPoints
                 (
                     polys,
                     cellFaces,
                     faces,
+                    faceOwner,
                     ensightGeometryFile
                 );
             }
@@ -736,7 +781,7 @@ void Foam::ensightMesh::writeAllPolys
         else
         {
             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
-            toMaster<< meshCellSets_.polys << cellFaces << faces;
+            toMaster<< meshCellSets_.polys << cellFaces << faces << faceOwner;
         }
     }
 }
@@ -1197,7 +1242,7 @@ void Foam::ensightMesh::write
             pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
 
             // Find the list of master faces belonging to the faceZone,
-            // in loacal numbering
+            // in local numbering
             faceList faceZoneFaces(fz().localFaces());
 
             // Count how many master faces belong to the faceZone. Is there
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H
index 35bc2e772a8d11ddaae26291c5f8800dfbc4460b..e1f7e8384b9b124ffbfe740c3a177ace95699f3b 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H
@@ -192,6 +192,7 @@ private:
             const labelList& polys,
             const cellList& cellFaces,
             const faceList& faces,
+            const labelList& faceOwner,
             ensightStream& ensightGeometryFile
         ) const;
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index f7983e4b43678cfa4e3ff5f6b3903602668336b3..bc5899154266eb4038934b083fb3ec04b6eecd3a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -135,6 +135,12 @@ int main(int argc, char *argv[])
         "specify faceZones to write - eg '( slice \"mfp-.*\" )'."
     );
     argList::addOption
+    (
+        "fields",
+        "wordReList",
+        "specify fields to export (all by default) - eg '( \"U.*\" )'."
+    );
+    argList::addOption
     (
         "cellZone",
         "word",
@@ -224,6 +230,13 @@ int main(int argc, char *argv[])
         zonePatterns = wordReList(args.optionLookup("faceZones")());
     }
 
+    const bool selectedFields = args.optionFound("fields");
+    wordReList fieldPatterns;
+    if (selectedFields)
+    {
+        fieldPatterns = wordReList(args.optionLookup("fields")());
+    }
+
     word cellZoneName;
     const bool doCellZone = args.optionReadIfPresent("cellZone", cellZoneName);
 
@@ -430,6 +443,15 @@ int main(int argc, char *argv[])
             {
                 const word& fieldName = fieldNames[j];
 
+                // Check if the field has to be exported
+                if (selectedFields)
+                {
+                    if (!findStrings(fieldPatterns, fieldName))
+                    {
+                        continue;
+                    }
+                }
+
                 #include "checkData.H"
 
                 if (!variableGood)
diff --git a/src/finiteVolume/cfdTools/general/include/alphaControls.H b/src/finiteVolume/cfdTools/general/include/alphaControls.H
index c172904f4f0fc6e88bb5dc788ed98c4f9a46210e..e7b3cf7794c14d1dc34d0c5d42ec2873ef73a333 100644
--- a/src/finiteVolume/cfdTools/general/include/alphaControls.H
+++ b/src/finiteVolume/cfdTools/general/include/alphaControls.H
@@ -2,11 +2,3 @@ const dictionary& alphaControls = mesh.solverDict(alpha1.name());
 label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
 label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
 Switch MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
-
-if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
-{
-    FatalErrorIn(args.executable())
-        << "Sub-cycling alpha is only allowed for PISO operation, "
-           "i.e. when the number of outer-correctors = 1"
-        << exit(FatalError);
-}
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
index c6e55ef2a060131e2b591bcc772b352257119c5e..3894a5c3e4533dcd2d689383f2dc59931bc4cb90 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
@@ -53,8 +53,12 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     fixedGradientFvPatchScalarField(p, iF),
     curTimeIndex_(-1)
 {
+    // Map value. Set unmapped values and overwrite with mapped ptf
+    fvPatchField<scalar>::operator=(patchInternalField());
+    map(ptf, mapper);
+    // Map gradient. Set unmapped values and overwrite with mapped ptf
     gradient() = 0.0;
-    map(gradient(), mapper);
+    gradient().map(ptf.gradient(), mapper);
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 60ec434b13278cc7273580e4b8572c68ba78b69f..1d7770f050e8c322797f0166ef2185bc4f085b93 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -194,6 +194,36 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 {}
 
 
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const triSurface& surface,
+    const word& sampleSourceName
+)
+:
+    sampledSurface(name, mesh),
+    surface_
+    (
+        IOobject
+        (
+            name,
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                  // registry
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        surface
+    ),
+    sampleSource_(samplingSourceNames_[sampleSourceName]),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index 79272f2dbb22f1ee3681af10277a0a9b56c04127..1ea10b01e9b5f195ec7bce5e757a392440ad396b 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -161,6 +161,15 @@ public:
             const dictionary& dict
         );
 
+        //- Construct from triSurface
+        sampledTriSurfaceMesh
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const triSurface& surface,
+            const word& sampleSourceName
+        );
+
 
     //- Destructor
     virtual ~sampledTriSurfaceMesh();
diff --git a/src/surfMesh/surfaceFormats/stl/STLpoint.H b/src/surfMesh/surfaceFormats/stl/STLpoint.H
index 834a4cd196d5b9dde1f51f3157ae7a0fe9330d37..157e832ba5a66581799ac4049f1b30ef792929f3 100644
--- a/src/surfMesh/surfaceFormats/stl/STLpoint.H
+++ b/src/surfMesh/surfaceFormats/stl/STLpoint.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,7 +28,6 @@ Description
     A vertex point representation for STL files.
 
 SourceFiles
-    STLpointI.H
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/triSurface/triSurface/interfaces/STL/STLpoint.H b/src/triSurface/triSurface/interfaces/STL/STLpoint.H
deleted file mode 100644
index 307fc3db3fe187b96052cf45df210f4631077248..0000000000000000000000000000000000000000
--- a/src/triSurface/triSurface/interfaces/STL/STLpoint.H
+++ /dev/null
@@ -1,89 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 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::STLpoint
-
-Description
-    A vertex point representation for STL files.
-
-SourceFiles
-    STLpointI.H
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef STLpoint_H
-#define STLpoint_H
-
-#include "point.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class STLpoint Declaration
-\*---------------------------------------------------------------------------*/
-
-class STLpoint
-:
-    public Vector<float>
-{
-
-public:
-
-    // Constructors
-
-        //- Construct null
-        inline STLpoint();
-
-        //- Construct from components
-        inline STLpoint(float x, float y, float z);
-
-        //- Construct from point
-        inline STLpoint(const point&);
-
-        //- Construct from istream
-        inline STLpoint(Istream&);
-
-
-    // Member Operators
-
-        inline operator point() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "STLpointI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/triSurface/triSurface/interfaces/STL/STLpointI.H b/src/triSurface/triSurface/interfaces/STL/STLpointI.H
deleted file mode 100644
index f0fdd2c585324ab147602406a26cb2911ec5d2c9..0000000000000000000000000000000000000000
--- a/src/triSurface/triSurface/interfaces/STL/STLpointI.H
+++ /dev/null
@@ -1,64 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-inline STLpoint::STLpoint()
-{}
-
-inline STLpoint::STLpoint(float x, float y, float z)
-:
-    Vector<float>(x, y, z)
-{}
-
-inline STLpoint::STLpoint(const point& pt)
-:
-    Vector<float>(float(pt.x()), float(pt.y()), float(pt.z()))
-{}
-
-inline STLpoint::STLpoint(Istream& is)
-:
-    Vector<float>(is)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline STLpoint::operator point() const
-{
-    return point(x(), y(), z());
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/tutorials/mesh/foamyQuadMesh/OpenCFD/Allrun b/tutorials/mesh/foamyQuadMesh/OpenCFD/Allrun
index 36e81e496f29362613b2a2051d8e905c05cc19fd..be87dce42024b37377a4ce221b83d678d4cff654 100755
--- a/tutorials/mesh/foamyQuadMesh/OpenCFD/Allrun
+++ b/tutorials/mesh/foamyQuadMesh/OpenCFD/Allrun
@@ -4,6 +4,11 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+
+# Make sure surface is oriented properly. Since the letters
+# contain disconnected regions use
+#   surfaceOrient -usePierceTest
+
 cp system/controlDict.mesher system/controlDict
 
 runApplication surfaceFeatureExtract
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/T
similarity index 79%
rename from tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water
rename to tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/T
index 31402a0144512b8aaa33b8625de36207b50d1088..94f61b78252fd4ac59fea878fd1c3fbdd907f6c9 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/T
@@ -10,28 +10,30 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "constant";
-    object      alpha.water;
+    object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 0 0 0 0 0 0];
+dimensions      [0 0 0 1 0 0 0];
 
-internalField   uniform 0;
+internalField   uniform 300;
 
 boundaryField
 {
-    atmosphere
-    {
-        type            inletOutlet;
-        inletValue      uniform 0;
-        value           uniform 0;
-    }
     walls
     {
         type            zeroGradient;
     }
-}
 
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/U b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/U
similarity index 85%
rename from tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/U
rename to tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/U
index 7ed7234005f66eaff11a126399bb2a8850aef720..e18422b1de70bcafb4d5148125b2c2c329d429a5 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/U
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/U
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "constant";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,17 +20,19 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    atmosphere
+    front
     {
-        type            pressureInletOutletVelocity;
-        value           uniform (0 0 0);
+        type            empty;
+    }
+    back
+    {
+        type            empty;
     }
     walls
     {
-        type            fixedValue;
+        type            movingWallVelocity;
         value           uniform (0 0 0);
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water.org b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/alpha.water.org
similarity index 83%
rename from tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water.org
rename to tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/alpha.water.org
index 33b377faf6b5a5babb5405e44c5746ab997f2779..70366ec29f7aae0cb1cc50f62a3d186553321883 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/alpha.water.org
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/alpha.water.org
@@ -10,8 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "constant";
-    object      alpha.water.org;
+    object      alpha.water;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,11 +20,13 @@ internalField   uniform 0;
 
 boundaryField
 {
-    atmosphere
+    front
     {
-        type            inletOutlet;
-        inletValue      uniform 0;
-        value           uniform 0;
+        type            empty;
+    }
+    back
+    {
+        type            empty;
     }
     walls
     {
@@ -33,5 +34,4 @@ boundaryField
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..ca99ff2f92f2265c5f58bb44c9ce79ec78ead77f
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e6;
+
+boundaryField
+{
+    front
+    {
+        type            empty;
+    }
+    back
+    {
+        type            empty;
+    }
+    walls
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p_rgh
similarity index 75%
rename from tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
rename to tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p_rgh
index 1b9491f31eaf3e549e37e506d6f04eb8eedb2b5c..b4a73cde4217f67330b55bfb964f427ae53c8072 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/0/p_rgh
@@ -10,23 +10,28 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "constant";
     object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 0;
+internalField   uniform 1e6;
 
 boundaryField
 {
-    atmosphere
+    front
     {
-        type            totalPressure;
-        rho             rho;
-        psi             none;
-        gamma           1;
-        p0              uniform 0;
-        value           uniform 0;
+        type            empty;
     }
+    back
+    {
+        type            empty;
+    }
+    walls
+    {
+        type            fixedFluxPressure;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allclean b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..46bf04f6aa941c241cd04b7b4fc9fb5f70348f2b
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+foamCleanTutorials cases
+rm -rf 0/alpha.water 0/alpha.water.gz 0/T.air.gz 0/T.water.gz \
+    probes wallPressure pRefProbe
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allrun b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..7ea0ac7744adb48402cd7061e4ef780e8853c3df
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/Allrun
@@ -0,0 +1,13 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+runApplication blockMesh
+cp 0/alpha.water.org 0/alpha.water
+runApplication setFields
+runApplication `getApplication`
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/RASProperties b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..4773e27a41a100b308cec5d0fd481bb6aab1eda0
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        laminar;
+
+turbulence      off;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/dynamicMeshDict b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..f8a9a9bd8b3dd8cd97772d60b537574b1f46ea6c
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/dynamicMeshDict
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   solidBodyMotionFvMesh;
+
+solidBodyMotionFvMeshCoeffs
+{
+    solidBodyMotionFunction SDA;
+    SDACoeffs
+    {
+        CofG            ( 0 0 0 );
+        lamda           50;
+        rollAmax        0.22654;
+        rollAmin        0.10472;
+        heaveA          3.79;
+        swayA           2.34;
+        Q               2;
+        Tp              13.93;
+        Tpn             11.93;
+        dTi             0.059;
+        dTp             -0.001;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/g b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..897615a50df92736e7b9c64bb2e64fba539496fc
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 0 -9.81 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/blockMeshDict.m4 b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000000000000000000000000000000000..8285c3af966ad1948fe17e4883ee8ed5d357dd50
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,145 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    `format'      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// General m4 macros
+
+changecom(//)changequote([,]) dnl>
+define(calc, [esyscmd(perl -e 'use Math::Trig; print ($1)')]) dnl>
+define(VCOUNT, 0)
+define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])
+
+define(hex2D, hex (b$1 b$2 b$3 b$4 f$1 f$2 f$3 f$4))
+define(quad2D, (b$1 b$2 f$2 f$1))
+define(frontQuad, (f$1 f$2 f$3 f$4))
+define(backQuad, (b$1 b$4 b$3 b$2))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// User-defined parameters
+
+convertToMeters 1;
+
+define(l, 1.0)      // Length of tank (x-direction)
+define(b, 40)       // Breadth of tank (y-direction)
+define(h, 30)       // Depth of tank (z-direction)
+
+define(hlc, 5)      // Depth to the top (height) of lower chamfer
+define(huc, 10)     // Height of upper chamfer
+
+define(thetalc, 45) // Angle of lower chamfer to the horizontal
+define(thetauc, 45) // Angle of upper chamfer to the horizontal
+
+define(CofGy, calc(b/2.0))  // Centre of gravity in y-direction
+define(CofGz, 10.0)         // Centre of gravity in z-direction
+
+define(Nl, 1)       // Number of cells in the length (1 for 2D)
+define(Nb, 40)      // Number of cells in the breadth
+define(Nhlc, 6)     // Number of cells in the height of the lower champfer
+define(Nh, 16)      // Number of cells in the height between the chamfers
+define(Nhuc, 12)    // Number of cells in the height of the upper champfer
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Derived parameters
+
+define(blc, calc(hlc/tan(deg2rad(thetalc)))) // Breadth to the top (height) of lower chamfer
+define(buc, calc(huc/tan(deg2rad(thetauc)))) // Breadth of upper chamfer
+
+define(Yl, -CofGy)
+define(Yllc, calc(Yl + blc))
+define(Yluc, calc(Yl + buc))
+
+define(Yr, calc(Yl + b))
+define(Yrlc, calc(Yr - blc))
+define(Yruc, calc(Yr - buc))
+
+define(Zb, -CofGz)
+define(Zlc, calc(Zb + hlc))
+define(Zt, calc(Zb + h))
+define(Zuc, calc(Zt - huc))
+
+define(Xf, calc(l/2.0))
+define(Xb, calc(Xf - l))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Parametric description
+
+vertices
+(
+    (Xb Yllc Zb) vlabel(bllcb)
+    (Xb Yl Zlc)  vlabel(bllc)
+    (Xb Yl Zuc)  vlabel(bluc)
+    (Xb Yluc Zt) vlabel(bluct)
+    (Xb Yrlc Zb) vlabel(brlcb)
+    (Xb Yr Zlc)  vlabel(brlc)
+    (Xb Yr Zuc)  vlabel(bruc)
+    (Xb Yruc Zt) vlabel(bruct)
+
+    (Xf Yllc Zb) vlabel(fllcb)
+    (Xf Yl Zlc)  vlabel(fllc)
+    (Xf Yl Zuc)  vlabel(fluc)
+    (Xf Yluc Zt) vlabel(fluct)
+    (Xf Yrlc Zb) vlabel(frlcb)
+    (Xf Yr Zlc)  vlabel(frlc)
+    (Xf Yr Zuc)  vlabel(fruc)
+    (Xf Yruc Zt) vlabel(fruct)
+);
+
+blocks
+(
+    // block0
+    hex2D(llcb, rlcb, rlc, llc)
+    (Nb Nhlc Nl)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(llc, rlc, ruc, luc)
+    (Nb Nh Nl)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(luc, ruc, ruct, luct)
+    (Nb Nhuc Nl)
+    simpleGrading (1 1 1)
+);
+
+patches
+(
+    patch walls
+    (
+        quad2D(llcb, rlcb)
+        quad2D(rlcb, rlc)
+        quad2D(rlc, ruc)
+        quad2D(ruc, ruct)
+        quad2D(ruct, luct)
+        quad2D(luct, luc)
+        quad2D(luc, llc)
+        quad2D(llc, llcb)
+    )
+
+    empty front
+    (
+        frontQuad(llcb, rlcb, rlc, llc)
+        frontQuad(llc, rlc, ruc, luc)
+        frontQuad(luc, ruc, ruct, luct)
+    )
+
+    empty back
+    (
+        backQuad(llcb, rlcb, rlc, llc)
+        backQuad(llc, rlc, ruc, luc)
+        backQuad(luc, ruc, ruct, luct)
+    )
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..f1bde7bdce369c1bfc963511de82c352d2772aab
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/polyMesh/boundary
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+3
+(
+    walls
+    {
+        type            patch;
+        nFaces          148;
+        startFace       2646;
+    }
+    front
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          1360;
+        startFace       2794;
+    }
+    back
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          1360;
+        startFace       4154;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..ad67959d02f43287554544a0f4f34f975eba3441
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+pMin            pMin [ 1 -1 -2 0 0 0 0 ] 1000;
+
+sigma           sigma [ 1 0 -2 0 0 0 0 ] 0;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.air b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.air
new file mode 100644
index 0000000000000000000000000000000000000000..e61009c10be927d2af1b6981cc75ad368e85d5c3
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.air
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.water b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.water
new file mode 100644
index 0000000000000000000000000000000000000000..1ffcbddad7b44c5aecf0266ac4d48835674b0344
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/thermophysicalProperties.water
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18.0;
+    }
+    equationOfState
+    {
+        R           3000;
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/transportProperties b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..3729fc7648bb9d76cae6c48139d4b359b92e526c
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+    transportModel  Newtonian;
+    nu              nu [ 0 2 -1 0 0 0 0 ] 1e-06;
+    rho             rho [ 1 -3 0 0 0 0 0 ] 998.2;
+}
+
+air
+{
+    transportModel  Newtonian;
+    nu              nu [ 0 2 -1 0 0 0 0 ] 1.48e-05;
+    rho             rho [ 1 -3 0 0 0 0 0 ] 1;
+}
+
+sigma           sigma [ 1 0 -2 0 0 0 0 ] 0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/turbulenceProperties b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..c2c3b28a1b4e8f4a2cae55f58bd61f9b1a67b488
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/controlDict b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..255b465ca02f1260eeee34d836c9dea027b3991a
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/controlDict
@@ -0,0 +1,98 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     compressibleInterDyMFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         40;
+
+deltaT          0.0001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.05;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+maxAlphaCo      0.5;
+
+maxDeltaT       1;
+
+functions
+{
+    probes
+    {
+        type            probes;
+        functionObjectLibs ("libsampling.so");
+        outputControl   outputTime;
+        probeLocations
+        (
+            ( 0 9.95 19.77 )
+            ( 0 -9.95 19.77 )
+        );
+        fields
+        (
+            p
+        );
+    }
+
+    wallPressure
+    {
+        type            surfaces;
+        functionObjectLibs ("libsampling.so");
+        outputControl   outputTime;
+        surfaceFormat   raw;
+        fields
+        (
+            p
+        );
+        interpolationScheme cellPoint;
+
+        surfaces
+        (
+            walls
+            {
+                type        patch;
+                patches     (walls);
+                triangulate false;
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/decomposeParDict b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..46583aaa1d5d6c360b1cd177b31d54802a80888a
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 16;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 4 2 2 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSchemes b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..3d39f5b409eb703d6ffea93e60538665384fb81b
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSchemes
@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss interfaceCompression 1;
+
+    div(rho*phi,U)  Gauss vanLeerV;
+    div(phi,thermo:rho.water) Gauss linear;
+    div(phi,thermo:rho.air) Gauss linear;
+    div(rho*phi,T)  Gauss vanLeer;
+    div(rho*phi,K)  Gauss linear;
+    div((phi+meshPhi),p)  Gauss linear;
+
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+    pcorr;
+    alpha;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSolution b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..2bec428fa53bd8344a668261028dd009bc1920d1
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/fvSolution
@@ -0,0 +1,134 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    alpha.water
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 3;
+        cAlpha          1.5;
+    }
+
+    ".*(rho|rhoFinal)"
+    {
+        solver          diagonal;
+    }
+
+    pcorr
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       1e-05;
+            relTol          0;
+            smoother        DICGaussSeidel;
+            nPreSweeps      0;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration false;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     1;
+        }
+
+        tolerance       1e-05;
+        relTol          0;
+        maxIter         100;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.01;
+        smoother        DIC;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    p_rghFinal
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       2e-09;
+            relTol          0;
+            nVcycles        2;
+            smoother        DICGaussSeidel;
+            nPreSweeps      2;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration true;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     1;
+        }
+
+        tolerance       2e-09;
+        relTol          0;
+        maxIter         20;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+        nSweeps         1;
+    }
+
+    "(T|k|B|nuTilda).*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   no;
+    transonic           no;
+    nOuterCorrectors    1;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors 0;
+    correctPhi          no;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        "U.*"           1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/setFieldsDict b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..410a0a6f10301fb287d4b429143355dd6ccd526d
--- /dev/null
+++ b/tutorials/multiphase/compressibleInterDyMFoam/ras/sloshingTank2D/system/setFieldsDict
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box ( -100 -100 -100 ) ( 100 100 0 );
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
index 7a04bf3624c8cbb92370b2dca60cf45a3f15530c..60308edee51a2c3dcbd68c32fe3d56f46210acfe 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
@@ -24,6 +24,11 @@ solvers
         cAlpha          1;
     }
 
+    ".*(rho|rhoFinal)"
+    {
+        solver          diagonal;
+    }
+
     pcorr
     {
         solver          PCG;
@@ -46,11 +51,6 @@ solvers
         maxIter         100;
     }
 
-    ".*(rho|rhoFinal)"
-    {
-        solver          diagonal;
-    }
-
     p_rgh
     {
         solver          GAMG;
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0.org/U b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0.org/U
index 6d3ce2531bf967c9224c82b1f55f1c483848fa0b..e2088b5e51788f1e0d90ff42f4c0ec40a72d039c 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0.org/U
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/0.org/U
@@ -28,8 +28,8 @@ boundaryField
     }
     walls
     {
-        type            fixedValue;
-        value           uniform (0 0 0);
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
     }
 }