diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
index 4d0927a6afa5c5bf1bbe3fdfde772851452462ac..363a8fd2dcfdeec4b0b1ae99082d74645ee3f796 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
@@ -109,6 +109,9 @@ int main(int argc, char *argv[])
         surfaceScalarField phiv_pos(U_pos & mesh.Sf());
         surfaceScalarField phiv_neg(U_neg & mesh.Sf());
 
+        fvc::makeRelative(phiv_pos, U);
+        fvc::makeRelative(phiv_neg, U);
+
         volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
         (
@@ -161,17 +164,9 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         mesh.movePoints(motionPtr->newPoints());
-        phiv_pos = U_pos & mesh.Sf();
-        phiv_neg = U_neg & mesh.Sf();
-        fvc::makeRelative(phiv_pos, U);
-        fvc::makeRelative(phiv_neg, U);
-        phiv_neg -= mesh.phi();
-        phiv_pos *= a_pos;
-        phiv_neg *= a_neg;
-        aphiv_pos = phiv_pos - aSf;
-        aphiv_neg = phiv_neg + aSf;
 
-        surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
+        phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+        Info<< phi.boundaryField()[0] << endl;
 
         surfaceVectorField phiUp
         (
@@ -183,6 +178,7 @@ int main(int argc, char *argv[])
         (
             aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
           + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+          + mesh.phi()*(a_pos*p_pos + a_neg*p_neg)
           + aSf*p_pos - aSf*p_neg
         );
 
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..0c17353ec45ee6026ab67e88b1ff3c882d8d3936
--- /dev/null
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H
@@ -0,0 +1,12 @@
+{
+    solve
+    (
+        fvm::ddt(rho, e)
+      + fvm::div(phi, e)
+      - fvm::laplacian(turbulence->alphaEff(), e)
+     ==
+      - p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi())
+    );
+
+    thermo.correct();
+}
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index bac588c73895f5eac7fff922f101cd7edb72b631..3f659d44ec4cfbe0b2a65cd35545f137870c7f66 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -76,34 +76,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu(0.75*Cds*phase2_.rho()*Ur*bp/d);
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*alpha_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index 99b5cacdf2fdc150df0039ce77247d99c628abb5..0a624c6b0a2a8a7bb8e8f4ba552b83eb13b3b967 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 79259cbe05858bdb9138023972809325a2df5d01..f68e5379d2403e64ead7aadf0ac369d5bfa5ed8d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 4f971b4bb7639a29991a074ff7fc48d75f6a5ecb..541868cf870c581b0a50e8bf44e29e7593e11bc8 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -73,15 +73,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index f1126125bb1be131e9f60f34bf83158510683a92..bf279a3990a0d065d2a680896bb78e2212ef9f0b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354f0329c8b9cbc253c92c0e681e85c8..d458f07b754218cb44de74b2cbe483b68626f188 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
index 1ee149c48ec63924503ec6349d91dbee8b846e9f..bce95f7d5a4ccd1ffbd0e0091eff166d1c1a687c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
@@ -21,16 +21,27 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
               + fvm::div(phase.phiAlpha(), U)
               - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U)
             )
-            - fvm::laplacian(alpha*nuEff, U)
-            - fvc::div
-              (
-                  alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
-                  "div(Rc)"
-              )
-          ==
-            - fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
-          //- (alpha*phase.rho())*fluid.lift(phase)
-            + (alpha/phase.rho())*fluid.Svm(phase)
+          - fvm::laplacian(alpha*nuEff, U)
+          - fvc::div
+            (
+                alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
+                "div(Rc)"
+            )
+         ==
+          - fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
+        //- (alpha*phase.rho())*fluid.lift(phase)
+          + (alpha/phase.rho())*fluid.Svm(phase)
+          - fvm::Sp
+            (
+                slamDampCoeff
+               *max
+                (
+                    mag(U.dimensionedInternalField()) - maxSlamVelocity,
+                    dimensionedScalar("U0", dimVelocity, 0)
+                )
+               /pow(mesh.V(), 1.0/3.0),
+                U
+            )
         )
     );
     mrfZones.addCoriolis(alpha, UEqns[phasei]);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
index 484ceae4ac904d5bcd7dc2555f3448e1630b1d75..79fa72d7587eabc198bc6a4d01b132db1db2e5a0 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
@@ -53,6 +53,18 @@
         phi += fvc::interpolate(alpha)*phase.phi();
     }
 
+    scalar slamDampCoeff
+    (
+        fluid.lookupOrDefault<scalar>("slamDampCoeff", 1)
+    );
+
+    dimensionedScalar maxSlamVelocity
+    (
+        "maxSlamVelocity",
+        dimVelocity,
+        fluid.lookupOrDefault<scalar>("maxSlamVelocity", GREAT)
+    );
+
     // dimensionedScalar pMin
     // (
     //     "pMin",
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H
deleted file mode 100644
index 436b68a2f26e553c343f3ba225f6858709bdfd44..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCoeffs.H
+++ /dev/null
@@ -1,80 +0,0 @@
-volScalarField dragCoeff
-(
-    IOobject
-    (
-        "dragCoeff",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0)
-);
-
-volVectorField liftForce
-(
-    IOobject
-    (
-        "liftForce",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedVector("liftForce", dimensionSet(1, -2, -2, 0, 0), vector::zero)
-);
-
-volScalarField heatTransferCoeff
-(
-    IOobject
-    (
-        "heatTransferCoeff",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("heatTransferCoeff", dimensionSet(1, -1, -3, -1, 0), 0)
-);
-
-{
-    volVectorField Ur = U1 - U2;
-    volScalarField magUr = mag(Ur);
-
-    if (dispersedPhase == "1")
-    {
-        dragCoeff = drag1->K(magUr);
-        heatTransferCoeff = heatTransfer1->K(magUr);
-    }
-    else if (dispersedPhase == "2")
-    {
-        dragCoeff = drag2->K(magUr);
-        heatTransferCoeff = heatTransfer2->K(magUr);
-    }
-    else if (dispersedPhase == "both")
-    {
-        dragCoeff =
-        (
-            alpha2*drag1->K(magUr)
-          + alpha1*drag2->K(magUr)
-        );
-
-        heatTransferCoeff =
-        (
-            alpha2*heatTransfer1->K(magUr)
-          + alpha1*heatTransfer2->K(magUr)
-        );
-    }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "dispersedPhase: " << dispersedPhase << " is incorrect"
-            << exit(FatalError);
-    }
-
-    volScalarField alphaCoeff
-    (
-        (alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha)
-    );
-    dragCoeff *= alphaCoeff;
-    heatTransferCoeff *= alphaCoeff;
-
-    liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U));
-}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index e34f219ebc0cbefa3972d31399b26ff852b52fbf..b52b766de9238653fdc699a75c6baa14e914cbc6 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -75,34 +75,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu = 0.75*Cds*phase2_.rho()*Ur*bp/d;
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*phase1_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index ecd25eb386fe6523e6266993e4df67072d1b34f3..25de6b1dc7001cd66e7094996655a3913dca7bed 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 3259b5d0e08b708a80e1f4b66c42322248abefa5..d4d77a2bb35df65363b95393b57ffd8d52d9c19e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -71,15 +71,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 736be885e0fdaa596b12b592989100c2870d5305..07f85b38a1f8e04809689ea6e1d94e13255552a8 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(phase2_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index 6cfc119c1c4ed8e2e79ae7e7332df7a17eba06dc..6319276408b935a09c5102ed726af4fb54004c9c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354f0329c8b9cbc253c92c0e681e85c8..d458f07b754218cb44de74b2cbe483b68626f188 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index 4ca01f1a3a89f3bc9e75697057a8fd42e0fd8920..c1d31263525c8989ce6f36b7fc35af22f099d25f 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -82,7 +82,6 @@ int main(int argc, char *argv[])
             rho = fluid.rho();
             #include "zonePhaseVolumes.H"
 
-            //#include "interfacialCoeffs.H"
             //#include "TEqns.H"
             #include "UEqns.H"
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 744faab968636541faa6d8c32d4d59398012404a..b93ea3a020b37862b4e9bf7acc4a615b67c0048b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -37,8 +37,10 @@
 
         U0s.set(phasei, new volVectorField(phase.U()));
         phi0s.set(phasei, new surfaceScalarField(phase.phi()));
+        mrfZones.absoluteFlux(phi0s[phasei]);
 
         phasei++;
+
     }
 
     surfaceScalarField phi0
@@ -67,6 +69,8 @@
 
         phase.U() = rAUs[phasei]*UEqns[phasei].H();
 
+        phase.phi().oldTime();
+        mrfZones.absoluteFlux(phase.phi().oldTime());
         mrfZones.absoluteFlux(phase.phi());
         phase.phi() =
         (
@@ -74,7 +78,8 @@
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
         );
         mrfZones.relativeFlux(phase.phi());
-        surfaceScalarField pphi0("pphi0", phase.phi());
+        mrfZones.relativeFlux(phase.phi().oldTime());
+        surfaceScalarField pphi0("pphi0", 1.0*phase.phi());
         pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf());
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
@@ -196,7 +201,7 @@
             p.relax();
             mSfGradp = pEqnIncomp.flux()/Dp;
 
-            U =  dimensionedVector("U", dimVelocity, vector::zero);
+            U = dimensionedVector("U", dimVelocity, vector::zero);
 
             phasei = 0;
             forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
@@ -254,9 +259,17 @@
                             phasej++;
                         }
 
-                        phase.U() +=
-                            (1.0/phase.rho())
-                           *rAUs[phasei]*(*dcIter())*U0s[phasej];
+                        // phase.U() +=
+                        //     (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
+                        //    *U0s[phasej];
+
+                        phase.U() += fvc::reconstruct
+                        (
+                            fvc::interpolate
+                            (
+                                (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
+                            )*phi0s[phasej]
+                        );
                     }
                 }
 
diff --git a/applications/test/tensor/Test-tensor.C b/applications/test/tensor/Test-tensor.C
index efc7824dd6577151973afda3a65e028ae323fb3e..6da37d33fffad0bd9cbcb01a64bc64b1b384ec89 100644
--- a/applications/test/tensor/Test-tensor.C
+++ b/applications/test/tensor/Test-tensor.C
@@ -50,10 +50,14 @@ int main()
         << (t1 & t7 & t1.T()) << " " << transform(t1, t7) << endl;
 
     symmTensor st1(1, 2, 3, 4, 5, 6);
+    symmTensor st2(7, 8, 9, 10, 11, 12);
 
     Info<< "Check symmetric transformation "
         << transform(t1, st1) << endl;
 
+    Info<< "Check for dot product of symmetric tensors "
+        << (st1 & st2) << endl;
+
     vector v1(1, 2, 3);
 
     Info<< sqr(v1) << endl;
diff --git a/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.C b/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.C
index 4d1234dcb9762f572e507b0d37815fa577ec00c7..542493f1821818e11a35d28dda688f2e914d4587 100644
--- a/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.C
+++ b/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.C
@@ -55,6 +55,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.H b/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.H
index da45faa82c10d4e2db4a7f6e795e62df0c52107b..56f5aac80bcfbcf73d3991a2068633e108c27987 100644
--- a/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/symmTensorFieldField/symmTensorFieldField.H
@@ -65,6 +65,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
index 4d9e4b6dbcab805efea5895272141d6bb301cbfc..fe5c4a3b79d251888b474fa33944a2dceb8eb3a9 100644
--- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
+++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
@@ -164,6 +164,9 @@ tmp<Field<symmTensor> > transformFieldMask<symmTensor>
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.H b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.H
index b73e270e15c8c638501343b118f1e41789ad229c..f3a8cabd23569bc140fd8de82cdb609be093ad36 100644
--- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.H
+++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.H
@@ -39,6 +39,7 @@ SourceFiles
 #include "vectorField.H"
 #include "sphericalTensor.H"
 #include "symmTensor.H"
+#include "tensor.H"
 
 #define TEMPLATE
 #include "FieldFunctionsM.H"
@@ -69,6 +70,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv)
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, dot)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.C b/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.C
index d5131675270475b3335dd5712f7ebb13c5e8d3de..23b616aa7916bf11f41c4404863aedb9eece0162 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.C
@@ -53,6 +53,9 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv, inv)
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual, transform)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.H b/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.H
index c0a8b81cde8174c185e63641a1799622c2763f46..67e1c84d0acce2f7b52c972c97e96ee8a55bff6f 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricSymmTensorField/GeometricSymmTensorField.H
@@ -65,6 +65,10 @@ UNARY_FUNCTION(symmTensor, symmTensor, inv, inv)
 
 UNARY_OPERATOR(vector, symmTensor, *, hdual, transform)
 
+BINARY_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
+BINARY_TYPE_OPERATOR(tensor, symmTensor, symmTensor, &, '&', dot)
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
index 5f9fbe9fe9483f799fdfb2f4400dedce0343d102..4566cc7550897d5625a52ca7c12c73c1f7425a13 100644
--- a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
+++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Vector.H"
+#include "Tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -184,18 +185,21 @@ inline Vector<Cmpt> operator*(const SymmTensor<Cmpt>& st)
 
 //- Inner-product between two symmetric tensors
 template <class Cmpt>
-inline SymmTensor<Cmpt>
+inline Tensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
 {
-    return SymmTensor<Cmpt>
+    return Tensor<Cmpt>
     (
         st1.xx()*st2.xx() + st1.xy()*st2.xy() + st1.xz()*st2.xz(),
         st1.xx()*st2.xy() + st1.xy()*st2.yy() + st1.xz()*st2.yz(),
         st1.xx()*st2.xz() + st1.xy()*st2.yz() + st1.xz()*st2.zz(),
 
+        st1.xy()*st2.xx() + st1.yy()*st2.xy() + st1.yz()*st2.xz(),
         st1.xy()*st2.xy() + st1.yy()*st2.yy() + st1.yz()*st2.yz(),
         st1.xy()*st2.xz() + st1.yy()*st2.yz() + st1.yz()*st2.zz(),
 
+        st1.xz()*st2.xx() + st1.yz()*st2.xy() + st1.zz()*st2.xz(),
+        st1.xz()*st2.xy() + st1.yz()*st2.yy() + st1.zz()*st2.yz(),
         st1.xz()*st2.xz() + st1.yz()*st2.yz() + st1.zz()*st2.zz()
     );
 }
diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H
index c3dd5228e021a66dada1fd1c8c5743d57fbc92c1..a1554dc916cfc7ea8225ce5788d09cd67ece0f99 100644
--- a/src/OpenFOAM/primitives/Tensor/Tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/Tensor.H
@@ -40,13 +40,15 @@ SourceFiles
 
 #include "Vector.H"
 #include "SphericalTensor.H"
-#include "SymmTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+template<class Cmpt>
+class SymmTensor;
+
 /*---------------------------------------------------------------------------*\
                            Class Tensor Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H
index cbfe6772e154f68f9d92e226b6167fe93d8a05d0..237bb6ceb393057cfd9657252d55a08f6facd25b 100644
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "SymmTensor.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphaair b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphaair
index 459439ad6aeece48c733362fdbe67ab8bd264cd9..11b988820d587db93bb057de863b7915b1075878 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphaair
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphaair
@@ -970,381 +970,381 @@ internalField   nonuniform List<scalar>
 0
 0
 0
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
 1
 1
 1
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphas b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphas
index fc0a0cff8050245fabdf0c455763851e64ea964d..b32a5dff468eecddcda9456f39387e71ccad0631 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphas
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphas
@@ -970,381 +970,381 @@ internalField   nonuniform List<scalar>
 1
 1
 1
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
 0
 0
 0
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphawater b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphawater
index b54346a56ed6b1aab843c9670607daf8451d73b2..e2a093f9e48b3d1acaad653631178114425d2b71 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphawater
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/alphawater
@@ -970,381 +970,381 @@ internalField   nonuniform List<scalar>
 1
 1
 1
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
 0
 0
 0
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/LESProperties b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/LESProperties
index 67d11839160e197d48177cbbe65a2771af91170d..ec20bc01dc69b12eccdc04c57e2e0be12a752262 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/LESProperties
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/LESProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-LESModel        Smagorinsky;
+LESModel        laminar;
 
 printCoeffs     on;
 
@@ -26,61 +26,4 @@ cubeRootVolCoeffs
     deltaCoeff      1;
 }
 
-PrandtlCoeffs
-{
-    delta           cubeRootVol;
-    cubeRootVolCoeffs
-    {
-        deltaCoeff      1;
-    }
-
-    smoothCoeffs
-    {
-        delta           cubeRootVol;
-        cubeRootVolCoeffs
-        {
-            deltaCoeff      1;
-        }
-
-        maxDeltaRatio   1.1;
-    }
-
-    Cdelta          0.158;
-}
-
-vanDriestCoeffs
-{
-    delta           cubeRootVol;
-    cubeRootVolCoeffs
-    {
-        deltaCoeff      1;
-    }
-
-    smoothCoeffs
-    {
-        delta           cubeRootVol;
-        cubeRootVolCoeffs
-        {
-            deltaCoeff      1;
-        }
-
-        maxDeltaRatio   1.1;
-    }
-
-    Aplus           26;
-    Cdelta          0.158;
-}
-
-smoothCoeffs
-{
-    delta           cubeRootVol;
-    cubeRootVolCoeffs
-    {
-        deltaCoeff      1;
-    }
-
-    maxDeltaRatio   1.1;
-}
-
-
 // ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/setFieldsDict b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/setFieldsDict
index 6c6bd80320706f8f536ce49fbb57c46ce7707a84..755148a96ecba41b315ba372b927a8b724005bf8 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/setFieldsDict
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/setFieldsDict
@@ -25,7 +25,7 @@ regions
 (
     boxToCell
     {
-        box (0 0 -0.1) (0.15 0.501 0.1);
+        box (0 0 -0.1) (0.15 0.701 0.1);
         fieldValues
         (
             volScalarFieldValue alphaair 0
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/U b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/U
new file mode 100644
index 0000000000000000000000000000000000000000..b82b7ce18aeae391934094cf3a54196467dce9dc
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uair b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uair
new file mode 100644
index 0000000000000000000000000000000000000000..e6cb214d17a9c5c0f0dd1ddfb392ba8ffef0ae1a
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uair
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      Uair;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Umercury b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Umercury
new file mode 100644
index 0000000000000000000000000000000000000000..71ee76460bd5cc2df7376fb473ec2db97a0f0a48
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Umercury
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      Umercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uoil b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uoil
new file mode 100644
index 0000000000000000000000000000000000000000..4f5e766638972d3237f292fc3fee6492d8ef4503
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uoil
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      Uoil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uwater b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uwater
new file mode 100644
index 0000000000000000000000000000000000000000..4977c8a6d5638865f5f2381e2ef4666d22ec4319
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/Uwater
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      Uwater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaair b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaair
new file mode 100644
index 0000000000000000000000000000000000000000..ba8f81857933dfd828c54adb49ff5dff7c82ac3b
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaair
@@ -0,0 +1,79 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      alphaair;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    rightWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 1;
+    }
+    lowerWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphamercury b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphamercury
new file mode 100644
index 0000000000000000000000000000000000000000..598b083fef717e84b4bd5dc548522810f689c70c
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphamercury
@@ -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       volScalarField;
+    location    "0";
+    object      alphamercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaoil b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaoil
new file mode 100644
index 0000000000000000000000000000000000000000..21ce07a42c449bd6a5e66b79129198a87415f8bb
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphaoil
@@ -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       volScalarField;
+    location    "0";
+    object      alphaoil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphas b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphas
new file mode 100644
index 0000000000000000000000000000000000000000..10e3ea7a0724cf102fe0ee5f8d3b7b7c7c3d7ade
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphas
@@ -0,0 +1,47 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      alphas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphawater b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphawater
new file mode 100644
index 0000000000000000000000000000000000000000..a91b1a596894ac539419d9e3797883d406ef2e4c
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/alphawater
@@ -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       volScalarField;
+    location    "0";
+    object      alphawater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p
new file mode 100644
index 0000000000000000000000000000000000000000..7d12a339569a8ca8a9ef78453ef2963f94f5fba0
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p
@@ -0,0 +1,59 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            multiphaseFixedFluxPressure;
+        value           uniform 0;
+    }
+
+    rightWall
+    {
+        type            multiphaseFixedFluxPressure;
+        value           uniform 0;
+    }
+
+    lowerWall
+    {
+        type            multiphaseFixedFluxPressure;
+        value           uniform 0;
+    }
+
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 0;
+        U               Uair;
+        phi             phiair;
+        rho             rho;
+        psi             none;
+        gamma           1;
+        value           uniform 0;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allclean b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..cd78b26599c3d4f7cd89a2197d3e52eb5f83087b
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+\rm -rf 0
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allrun b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..300fc394aa8d44bba50536bc57f38e54d93da0f0
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/Allrun
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Set application name
+application=`getApplication`
+
+\rm -rf 0
+cp -r 0.org 0
+
+runApplication blockMesh
+runApplication setFields
+runApplication decomposePar
+runParallel $application 4
+runApplication reconstructPar
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/LESProperties b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/LESProperties
new file mode 100644
index 0000000000000000000000000000000000000000..ec20bc01dc69b12eccdc04c57e2e0be12a752262
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/LESProperties
@@ -0,0 +1,29 @@
+/*--------------------------------*- 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      LESProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+LESModel        laminar;
+
+printCoeffs     on;
+
+delta           cubeRootVol;
+
+cubeRootVolCoeffs
+{
+    deltaCoeff      1;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/MRFZones b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..0034a58c2293992131cec848f2dcb0f8ad50a0d4
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/MRFZones
@@ -0,0 +1,20 @@
+/*--------------------------------*- 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      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0()
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/g b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..e0ac2653b5b370ad62f6770588121d30cac51627
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/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 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/motionProperties b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/motionProperties
new file mode 100644
index 0000000000000000000000000000000000000000..93e07019847e525066b39393f4cb5b25e7c424e5
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/motionProperties
@@ -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      motionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+movingFvMesh    staticFvMesh;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/blockMeshDict b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..fd11bf53fdc733aa71c3199b5151879989d99bdd
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/blockMeshDict
@@ -0,0 +1,108 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.146;
+
+vertices
+(
+    (0 0 0)
+    (2 0 0)
+    (2.16438 0 0)
+    (4 0 0)
+    (0 0.32876 0)
+    (2 0.32876 0)
+    (2.16438 0.32876 0)
+    (4 0.32876 0)
+    (0 4 0)
+    (2 4 0)
+    (2.16438 4 0)
+    (4 4 0)
+    (0 0 0.1)
+    (2 0 0.1)
+    (2.16438 0 0.1)
+    (4 0 0.1)
+    (0 0.32876 0.1)
+    (2 0.32876 0.1)
+    (2.16438 0.32876 0.1)
+    (4 0.32876 0.1)
+    (0 4 0.1)
+    (2 4 0.1)
+    (2.16438 4 0.1)
+    (4 4 0.1)
+);
+
+blocks
+(
+    hex (0 1 5 4 12 13 17 16) (92 15 1) simpleGrading (1 1 1)
+    hex (2 3 7 6 14 15 19 18) (76 15 1) simpleGrading (1 1 1)
+    hex (4 5 9 8 16 17 21 20) (92 180 1) simpleGrading (1 1 1)
+    hex (5 6 10 9 17 18 22 21) (8 180 1) simpleGrading (1 1 1)
+    hex (6 7 11 10 18 19 23 22) (76 180 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    leftWall
+    {
+        type wall;
+        faces
+        (
+            (0 12 16 4)
+            (4 16 20 8)
+        );
+    }
+    rightWall
+    {
+        type wall;
+        faces
+        (
+            (7 19 15 3)
+            (11 23 19 7)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (0 1 13 12)
+            (1 5 17 13)
+            (5 6 18 17)
+            (2 14 18 6)
+            (2 3 15 14)
+        );
+    }
+    atmosphere
+    {
+        type patch;
+        faces
+        (
+            (8 20 21 9)
+            (9 21 22 10)
+            (10 22 23 11)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/boundary b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..04429b16df52a30539b01ed1deef6ba5f46a118b
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/polyMesh/boundary
@@ -0,0 +1,52 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    leftWall
+    {
+        type            wall;
+        nFaces          195;
+        startFace       68014;
+    }
+    rightWall
+    {
+        type            wall;
+        nFaces          195;
+        startFace       68209;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          206;
+        startFace       68404;
+    }
+    atmosphere
+    {
+        type            patch;
+        nFaces          176;
+        startFace       68610;
+    }
+    defaultFaces
+    {
+        type            empty;
+        nFaces          68400;
+        startFace       68786;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/transportProperties b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..fad0a7a482eabf8aabe1be23939748a5941cf13d
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/transportProperties
@@ -0,0 +1,248 @@
+/*--------------------------------*- 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
+     {
+         nu         1e-06;
+         kappa      1e-06;
+         Cp         4195;
+         rho        1000;
+
+         diameterModel constant;
+         constantCoeffs
+         {
+             d               1e-3;
+         }
+     }
+
+     oil
+     {
+         nu         1e-06;
+         kappa      1e-06;
+         Cp         4195;
+         rho        500;
+
+         diameterModel constant;
+         constantCoeffs
+         {
+             d               1e-3;
+         }
+     }
+
+     mercury
+     {
+         nu         1.125e-07;
+         kappa      1e-06;
+         Cp         4195;
+         rho        13529;
+
+         diameterModel constant;
+         constantCoeffs
+         {
+             d               1e-3;
+         }
+     }
+
+     air
+     {
+         nu         1.48e-05;
+         kappa      2.63e-2;
+         Cp         1007;
+         rho        1;
+
+         diameterModel constant;
+         constantCoeffs
+         {
+             d              3e-3;
+         }
+     }
+);
+
+sigmas
+(
+    (air water)     0.07
+    (air oil)       0.07
+    (air mercury)   0.07
+    (water oil)     0
+    (water mercury) 0
+    (oil mercury)   0
+);
+
+interfaceCompression
+(
+    (air water)     1
+    (air oil)       1
+    (air mercury)   1
+    (water oil)     0
+    (water mercury) 0
+    (oil mercury)   0
+);
+
+virtualMass
+(
+    (air water)     0
+    (air oil)       0
+    (air mercury)   0
+    (water oil)     0.5
+    (water mercury) 0.5
+    (oil mercury)   0.5
+);
+
+drag
+(
+    (air water)
+    {
+        type blended;
+
+        air
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        water
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+
+    (air oil)
+    {
+        type blended;
+
+        air
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        oil
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+
+    (air mercury)
+    {
+        type blended;
+
+        air
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        mercury
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+
+    (water oil)
+    {
+        type blended;
+
+        water
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        oil
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+
+    (water mercury)
+    {
+        type blended;
+
+        water
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        mercury
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+
+    (oil mercury)
+    {
+        type blended;
+
+        oil
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        mercury
+        {
+            type SchillerNaumann;
+            residualPhaseFraction 0;
+            residualSlip 0;
+        }
+
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
+    }
+);
+
+
+// This is a dummy to support the Smagorinsky model
+transportModel  Newtonian;
+nu              nu [ 0 2 -1 0 0 0 0 ] 0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/turbulenceProperties b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..c2c3b28a1b4e8f4a2cae55f58bd61f9b1a67b488
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/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/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..902244128790ada7cfd7bb009c8fa1cc6ef243fb
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict
@@ -0,0 +1,56 @@
+/*--------------------------------*- 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     multiphaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         6;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+maxAlphaCo      0.5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/decomposeParDict b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..e7e490bf74b1a3f58a34923b4f98569f1b09e483
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/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 4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..5a4bb79d9628ec2e14f036974ba794d9b805986f
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes
@@ -0,0 +1,61 @@
+/*--------------------------------*- 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\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
+
+    "div\(phiAlpha.*,U.*\)" Gauss limitedLinearV 1;
+    div(Rc)                 Gauss linear;
+    "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+    pcorr;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..296b0cd48839b3ce647e02a066bcaaa341e692f1
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution
@@ -0,0 +1,94 @@
+/*--------------------------------*- 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
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-7;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration on;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    pFinal
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       1e-7;
+            relTol          0;
+            nVcycles        2;
+            smoother        GaussSeidel;
+            nPreSweeps      0;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration on;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     1;
+        }
+        tolerance       1e-7;
+        relTol          0;
+        maxIter         20;
+    }
+
+    pcorr
+    {
+        $pFinal;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-8;
+        relTol          0.1;
+        nSweeps         1;
+    }
+
+    UFinal
+    {
+        $U;
+        tolerance       1e-7;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    nAlphaSubCycles 3;
+}
+
+relaxationFactors
+{
+    "U.*"           1;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/setFieldsDict b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..3abdf3f73af90b343bc431de989f70d159098b39
--- /dev/null
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/setFieldsDict
@@ -0,0 +1,65 @@
+/*--------------------------------*- 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 alphaair 1
+    volScalarFieldValue alphawater 0
+    volScalarFieldValue alphaoil 0
+    volScalarFieldValue alphamercury 0
+    volVectorFieldValue U ( 0 0 0 )
+);
+
+regions
+(
+    boxToCell
+    {
+        box ( 0 0 -1 ) ( 0.1461 0.292 1 );
+        fieldValues
+        (
+            volScalarFieldValue alphawater 1
+            volScalarFieldValue alphaoil 0
+            volScalarFieldValue alphamercury 0
+            volScalarFieldValue alphaair 0
+        );
+    }
+    boxToCell
+    {
+        box ( 0.1461 0 -1 ) ( 0.2922 0.292 1 );
+        fieldValues
+        (
+            volScalarFieldValue alphawater 0
+            volScalarFieldValue alphaoil 1
+            volScalarFieldValue alphamercury 0
+            volScalarFieldValue alphaair 0
+        );
+    }
+    boxToCell
+    {
+        box ( 0 0 -1 ) ( 0.1461 0.1 1 );
+        fieldValues
+        (
+            volScalarFieldValue alphawater 0
+            volScalarFieldValue alphaoil 0
+            volScalarFieldValue alphamercury 1
+            volScalarFieldValue alphaair 0
+        );
+    }
+);
+
+
+// ************************************************************************* //