diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
new file mode 100644
index 0000000000000000000000000000000000000000..c894b36b41c4e48a07a650590f9c6fe6fdc7edc2
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
@@ -0,0 +1,49 @@
+{
+    volScalarField k1
+    (
+        "k1",
+        alpha1*(thermo1.alpha()/rho1 + sqr(Ct)*nut2*thermo1.CpByCpv()/Prt)
+    );
+
+    volScalarField k2
+    (
+        "k2",
+        alpha2*(thermo2.alpha()/rho2 + nut2*thermo2.CpByCpv()/Prt)
+    );
+
+    volScalarField& he1 = thermo1.he();
+    volScalarField& he2 = thermo2.he();
+
+    Info<< max(he1) << min(he1) << endl;
+
+    fvScalarMatrix he1Eqn
+    (
+        fvm::ddt(alpha1, he1)
+      + fvm::div(alphaPhi1, he1)
+      - fvm::laplacian(k1, he1)
+     ==
+        heatTransferCoeff*(he1/thermo1.Cp())/rho1
+      - fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1)
+      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))
+    );
+
+    fvScalarMatrix he2Eqn
+    (
+        fvm::ddt(alpha2, he2)
+      + fvm::div(alphaPhi2, he2)
+      - fvm::laplacian(k2, he2)
+     ==
+        heatTransferCoeff*(he2/thermo2.Cp())/rho2
+      - fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2)
+      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))
+    );
+
+    he1Eqn.relax();
+    he1Eqn.solve();
+
+    he2Eqn.relax();
+    he2Eqn.solve();
+
+    thermo1.correct();
+    thermo2.correct();
+}
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
index 16ce2e3674b971c19291acbbfb292aeb2da3795e..d90ba7bbbdf308e0b43489c0e0efd851b5367e28 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -IturbulenceModel \
@@ -8,6 +9,8 @@ EXE_INC = \
     -Iaveraging
 
 EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
     -lincompressibleTransportModels \
     -lcompressiblePhaseModel \
     -lcompressibleEulerianInterfacialModels \
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
deleted file mode 100644
index 8f38ca9d91fe910cfad67e53608612d277564b3a..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
+++ /dev/null
@@ -1,36 +0,0 @@
-{
-    volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt));
-    volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt));
-
-    fvScalarMatrix T1Eqn
-    (
-        fvm::ddt(alpha1, T1)
-      + fvm::div(alphaPhi1, T1)
-      - fvm::laplacian(kByCp1, T1)
-     ==
-        heatTransferCoeff*T2/Cp1/rho1
-      - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)
-      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1
-    );
-
-    fvScalarMatrix T2Eqn
-    (
-        fvm::ddt(alpha2, T2)
-      + fvm::div(alphaPhi2, T2)
-      - fvm::laplacian(kByCp2, T2)
-     ==
-        heatTransferCoeff*T1/Cp2/rho2
-      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
-      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2
-    );
-
-    T1Eqn.relax();
-    T1Eqn.solve();
-
-    T2Eqn.relax();
-    T2Eqn.solve();
-
-    // Update compressibilities
-    psi1 = 1.0/(R1*T1);
-    psi2 = 1.0/(R2*T2);
-}
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
index 606b4573e47988c921954a011626b7ce80ecfb4b..a243db6300e555bfc8b909c74cd11352d7d66abd 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
@@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         }
         else // If not using kinetic theory is using Ct model
         {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
+            nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1;
         }
 
         volTensorField Rc1
         (
-            "Rc1",
+            "Rc",
             (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
         );
 
@@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         volTensorField gradU2T(fvc::grad(U2)().T());
         volTensorField Rc2
         (
-            "Rc2",
+            "Rc",
             (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
         );
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
index a818ee2e9ea7707a7707b25f40026a2bb08ad591..86d9203dce190cfe520a07149dccbe7ef66eab83 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
@@ -1,9 +1,9 @@
-surfaceScalarField alphaPhi1("alphaPhi1", phi1);
-surfaceScalarField alphaPhi2("alphaPhi2", phi2);
+surfaceScalarField alphaPhi1("alphaPhi", phi1);
+surfaceScalarField alphaPhi2("alphaPhi", phi2);
 
 {
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
+    word scheme("div(phi,alpha)");
+    word schemer("div(phir,alpha)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -68,16 +68,22 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
 
         if (g0.value() > 0.0)
         {
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
             ppMagf =
                 fvc::interpolate((1.0/rho1)*rAU1)
-               *fvc::interpolate
-                (
-                    g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-                );
+               *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
+
+            // ppMagf =
+            //     fvc::interpolate((1.0/rho1)*rAU1)
+            //    *fvc::interpolate
+            //     (
+            //         g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            //     );
 
             alpha1Eqn -= fvm::laplacian
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                alpha1f*ppMagf,
                 alpha1,
                 "laplacian(alphaPpMag,alpha1)"
             );
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
index 855a862b6ba51ac11c1fdadbae38e37f5fa35dc2..65372612b68532228120394c2bcc6e0925e44c2d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "rhoThermo.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvsPatchFields.H"
@@ -86,7 +87,7 @@ int main(int argc, char *argv[])
             #include "alphaEqn.H"
             #include "kEpsilon.H"
             #include "interfacialCoeffs.H"
-            #include "TEqns.H"
+            #include "EEqns.H"
             #include "UEqns.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
index e0cd5328c98a5585040b805c5257cf676ab24996..6095ad6845e3a6b08767950eb9298bc27db2e416 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
@@ -12,55 +12,43 @@
         )
     );
 
+    word phase1Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[0]
+      : "phase1"
+    );
+
+    word phase2Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[1]
+      : "phase2"
+    );
+
     autoPtr<phaseModel> phase1 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "1"
+        phase1Name
     );
 
     autoPtr<phaseModel> phase2 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "2"
+        phase2Name
     );
 
+    volScalarField& alpha1 = phase1();
+    volScalarField& alpha2 = phase2();
+    alpha2 = scalar(1) - alpha1;
+
     volVectorField& U1 = phase1->U();
     surfaceScalarField& phi1 = phase1->phi();
-    const dimensionedScalar& nu1 = phase1->nu();
-    const dimensionedScalar& k1 = phase1->kappa();
-    const dimensionedScalar& Cp1 = phase1->Cp();
-    dimensionedScalar rho10
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase1").lookup("rho0")
-    );
-    dimensionedScalar R1
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase1").lookup("R")
-    );
 
     volVectorField& U2 = phase2->U();
     surfaceScalarField& phi2 = phase2->phi();
-    const dimensionedScalar& nu2 = phase2->nu();
-    const dimensionedScalar& k2 = phase2->kappa();
-    const dimensionedScalar& Cp2 = phase2->Cp();
-    dimensionedScalar rho20
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase2").lookup("rho0")
-    );
-    dimensionedScalar R2
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase2").lookup("R")
-    );
 
     dimensionedScalar pMin
     (
@@ -69,102 +57,16 @@
         transportProperties.lookup("pMin")
     );
 
+    rhoThermo& thermo1 = phase1->thermo();
+    rhoThermo& thermo2 = phase2->thermo();
 
-    Info<< "Reading field T1" << endl;
-    volScalarField T1
-    (
-        IOobject
-        (
-            "T1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField& p = thermo1.p();
 
-    Info<< "Reading field T2" << endl;
-    volScalarField T2
-    (
-        IOobject
-        (
-            "T2",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField psi1
-    (
-        IOobject
-        (
-            "psi1",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R1*T1)
-    );
-
-    volScalarField psi2
-    (
-        IOobject
-        (
-            "psi2",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R2*T2)
-    );
-    psi2.oldTime();
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField rho1("rho1", rho10 + psi1*p);
-    volScalarField rho2("rho2", rho20 + psi2*p);
-
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField rho1("rho" + phase1Name, thermo1.rho());
+    const volScalarField& psi1 = thermo1.psi();
 
-    volScalarField alpha2
-    (
-        IOobject
-        (
-            "alpha2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - alpha1
-    );
+    volScalarField rho2("rho" + phase2Name, thermo2.rho());
+    const volScalarField& psi2 = thermo2.psi();
 
     volVectorField U
     (
@@ -179,27 +81,6 @@
         alpha1*U1 + alpha2*U2
     );
 
-    dimensionedScalar Cvm
-    (
-        "Cvm",
-        dimless,
-        transportProperties.lookup("Cvm")
-    );
-
-    dimensionedScalar Cl
-    (
-        "Cl",
-        dimless,
-        transportProperties.lookup("Cl")
-    );
-
-    dimensionedScalar Ct
-    (
-        "Ct",
-        dimless,
-        transportProperties.lookup("Ct")
-    );
-
     surfaceScalarField phi
     (
         IOobject
@@ -226,8 +107,6 @@
         alpha1*rho1 + alpha2*rho2
     );
 
-    #include "createRASTurbulence.H"
-
     Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
 
     volVectorField DDtU1
@@ -248,6 +127,29 @@
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
 
+    dimensionedScalar Cvm
+    (
+        "Cvm",
+        dimless,
+        transportProperties.lookup("Cvm")
+    );
+
+    dimensionedScalar Cl
+    (
+        "Cl",
+        dimless,
+        transportProperties.lookup("Cl")
+    );
+
+    dimensionedScalar Ct
+    (
+        "Ct",
+        dimless,
+        transportProperties.lookup("Ct")
+    );
+
+    #include "createRASTurbulence.H"
+
     IOdictionary interfacialProperties
     (
         IOobject
@@ -297,8 +199,8 @@
     if
     (
         !(
-            dispersedPhase == "1"
-         || dispersedPhase == "2"
+            dispersedPhase == phase1Name
+         || dispersedPhase == phase2Name
          || dispersedPhase == "both"
         )
     )
@@ -337,7 +239,7 @@
     (
         IOobject
         (
-            "rAU1",
+            "rAU" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -387,5 +289,5 @@
     );
 
     Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K1("K1", 0.5*magSqr(U1));
-    volScalarField K2("K2", 0.5*magSqr(U2));
+    volScalarField K1("K" + phase1Name, 0.5*magSqr(U1));
+    volScalarField K2("K" + phase2Name, 0.5*magSqr(U2));
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
index 8f05a1d1b373ccde772bfba4cee1b575574ae9f3..0a782ef51edf99f32eaf986325fd926db596e410 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
@@ -151,7 +151,7 @@
     (
         IOobject
         (
-            "nut2",
+            "nut" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -165,13 +165,13 @@
     (
         IOobject
         (
-            "nuEff1",
+            "nuEff" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        sqr(Ct)*nut2 + nu1
+        sqr(Ct)*nut2 + thermo1.mu()/rho1
     );
 
     Info<< "Calculating field nuEff2\n" << endl;
@@ -179,11 +179,11 @@
     (
         IOobject
         (
-            "nuEff2",
+            "nuEff" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        nut2 + nu2
+        nut2 + thermo2.mu()/rho2
     );
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
index d6ccf90289ad6891248e665609408a68d4e7c766..d53bec5ea4f5ed2dc8393425bb4c4ddfa816c39d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
@@ -38,12 +38,12 @@ volScalarField heatTransferCoeff
     volVectorField Ur(U1 - U2);
     volScalarField magUr(mag(Ur) + residualSlip);
 
-    if (dispersedPhase == "1")
+    if (dispersedPhase == phase1Name)
     {
         dragCoeff = drag1->K(magUr);
         heatTransferCoeff = heatTransfer1->K(magUr);
     }
-    else if (dispersedPhase == "2")
+    else if (dispersedPhase == phase2Name)
     {
         dragCoeff = drag2->K(magUr);
         heatTransferCoeff = heatTransfer2->K(magUr);
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
index 92fe0b0d74ab657d08fb917834ecd7b4afa9783d..f031e058986a9c4ac3c37c83fafeaf2f592f28df 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I../phaseModel/lnInclude
 
 LIB_LIBS = \
-    -lcompressiblePhaseModel
-
+    -lcompressiblePhaseModel \
+    -lfluidThermophysicalModels \
+    -lspecie
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
index 9bc6cadaf66495c673b2e7c69781aeabdaaf850d..95ade6ba05944fc306d26ce7c7ed1ffff4d97d2d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,8 +72,7 @@ Foam::tmp<Foam::volScalarField> Foam::heatTransferModels::RanzMarshall::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    dimensionedScalar Prb =
-        phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa();
+    volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa());
     volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb));
 
     return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d());
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
index 2fcce9913d9408482098065ad511e280e366c83f..7fdabf573092bd4219e75fa41044236e02b7d50e 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/foam/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I../phaseModel/lnInclude \
     -I../interfacialModels/lnInclude
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index 7d0767914ea7ef3d5ef2fa6f11bd9b1d8ef4abda..18db3a70df7a5265075b96ee74e3bc87a073576d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel
     phi1_(phase1.phi()),
     draga_(draga),
 
-    rho1_(phase1.rho()),
-    nu1_(phase1.nu()),
+    rho1_(phase1.rho()[0]), //***HGW
+    nu1_(phase1.nu()()[0]), //***HGW
 
     kineticTheoryProperties_
     (
@@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
     (
         IOobject
         (
-            "mu1",
+            "mu" + phase1.name(),
             U1_.time().timeName(),
             U1_.mesh(),
             IOobject::NO_READ,
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index f7010b1c459cc8419b658797c7d7e4998ba1ab3a..f9a7807d00d46f51c83c8ee2f21679866b898d68 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
@@ -11,10 +11,10 @@
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    volVectorField HbyA1("HbyA1", U1);
+    volVectorField HbyA1("HbyA" + phase1Name, U1);
     HbyA1 = rAU1*U1Eqn.H();
 
-    volVectorField HbyA2("HbyA2", U2);
+    volVectorField HbyA2("HbyA" + phase2Name, U2);
     HbyA2 = rAU2*U2Eqn.H();
 
     mrfZones.absoluteFlux(phi1.oldTime());
@@ -38,14 +38,14 @@
 
     surfaceScalarField phiHbyA1
     (
-        "phiHbyA1",
+        "phiHbyA" + phase1Name,
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
-        "phiHbyA2",
+        "phiHbyA" + phase2Name,
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
@@ -96,8 +96,16 @@
     //}
     //else
     {
-        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
-        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
+        surfaceScalarField phid1
+        (
+            "phid" + phase1Name,
+            fvc::interpolate(psi1)*phi1
+        );
+        surfaceScalarField phid2
+        (
+            "phid" + phase2Name,
+            fvc::interpolate(psi2)*phi2
+        );
 
         pEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
@@ -174,13 +182,15 @@
 
     p = max(p, pMin);
 
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    thermo1.correct();
+    thermo2.correct();
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
 
-    //***HGW if (thermo.dpdt())
+    if (thermo1.dpdt())
     {
         dpdt = fvc::ddt(p);
     }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
index 0ec11392098bd862e78e2790e25cf1f49ec557fd..e441b0417bb82b021ee7e0ae470e030900e43a71 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude
 
 LIB_LIBS = \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
index d80758ac8c7d6cb754af8fe5aeaf9d2e7e2f7b78..64d59b436d5f2e85d2bd263227878101b3b51c76 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel
     const word& phaseName
 )
 :
-    dict_
+    volScalarField
     (
-        transportProperties.subDict("phase" + phaseName)
+        IOobject
+        (
+            "alpha" + phaseName,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha", dimless, 0)
     ),
     name_(phaseName),
-    nu_
-    (
-        "nu",
-        dimensionSet(0, 2, -1, 0, 0),
-        dict_.lookup("nu")
-    ),
-    kappa_
-    (
-        "kappa",
-        dimensionSet(1, 1, -3, -1, 0),
-        dict_.lookup("kappa")
-    ),
-    Cp_
-    (
-        "Cp",
-        dimensionSet(0, 2, -2, -1, 0),
-        dict_.lookup("Cp")
-    ),
-    rho_
-    (
-        "rho",
-        dimDensity,
-        dict_.lookup("rho")
-    ),
+    phaseDict_(transportProperties.subDict(phaseName)),
+    thermo_(rhoThermo::New(mesh, phaseName)),
     U_
     (
         IOobject
@@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel
         mesh
     )
 {
+    thermo_->validate("phaseModel " + phaseName, "h", "e");
+
     const word phiName = "phi" + phaseName;
 
     IOobject phiHeader
@@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel
 
     dPtr_ = diameterModel::New
     (
-        dict_,
+        phaseDict_,
         *this
     );
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
index 90b2d2184e0d3f971beeb959fe46864805073e89..18f5a47e46950dbec54b8a364b170aa5a3d5bef2 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,25 +51,18 @@ class diameterModel;
 \*---------------------------------------------------------------------------*/
 
 class phaseModel
+:
+    public volScalarField
 {
     // Private data
 
-        dictionary dict_;
-
         //- Name of phase
         word name_;
 
-        //- Kinematic viscosity
-        dimensionedScalar nu_;
-
-        //- Thermal conductivity
-        dimensionedScalar kappa_;
-
-        //- Heat capacity
-        dimensionedScalar Cp_;
+        dictionary phaseDict_;
 
-        //- Density
-        dimensionedScalar rho_;
+        //- Thermophysical properties
+        autoPtr<rhoThermo> thermo_;
 
         //- Velocity
         volVectorField U_;
@@ -116,24 +110,34 @@ public:
 
         tmp<volScalarField> d() const;
 
-        const dimensionedScalar& nu() const
+        tmp<volScalarField> nu() const
+        {
+            return thermo_->mu()/thermo_->rho();
+        }
+
+        tmp<volScalarField> kappa() const
+        {
+            return thermo_->kappa();
+        }
+
+        tmp<volScalarField> Cp() const
         {
-            return nu_;
+            return thermo_->Cp();
         }
 
-        const dimensionedScalar& kappa() const
+        const volScalarField& rho() const
         {
-            return kappa_;
+            return thermo_->rho();
         }
 
-        const dimensionedScalar& Cp() const
+        const rhoThermo& thermo() const
         {
-            return Cp_;
+            return thermo_();
         }
 
-        const dimensionedScalar& rho() const
+        rhoThermo& thermo()
         {
-            return rho_;
+            return thermo_();
         }
 
         const volVectorField& U() const
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
index 7f53411529a31a855f706d13d801f5624d1b7c45..bc9a07b0a81105cdbacc53336c054253844f2d2b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
@@ -61,4 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEff2 = nut2 + nu2;
+nuEff2 = nut2 + thermo2.mu()/rho2;
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
index c91cce9a00e2df5ca912a75a57a02d6f813d72c4..d85181cba2e020e57fb2063d09d4af4f14b01807 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
@@ -4,7 +4,6 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar Cmu75 = ::pow(Cmu.value(), 0.75);
     scalar kappa_ = kappa.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
@@ -30,6 +29,8 @@
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -52,7 +53,8 @@
                     /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (nut2w[facei] + nu2_)*magFaceGradU[facei]
+                    (nut2w[facei] + mu2_[facei]/rho2_[facei])
+                   *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);
             }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
index b153a360146deefed6b9b4a9c2a8a8ecb66075d9..9aa032149ca51f65d51d412c4dc13b7d61b9c46c 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
@@ -2,13 +2,14 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar kappa_ = kappa.value();
     scalar E_ = E.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -20,11 +21,14 @@
 
                 // calculate yPlus
                 scalar yPlus =
-                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
+                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
+                   /(mu2_[facei]/rho2_[facei]);
 
                 if (yPlus > 11.6)
                 {
-                    nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                    nutw[facei] =
+                        (mu2_[facei]/rho2_[facei])
+                       *(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
index 34fdc24afc532d56c0fd0ddf0c08479237d2cdb9..ddb923cf8726096119eac8c6b62e4b763877c642 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
@@ -13,7 +13,7 @@
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -52,18 +52,23 @@
 
         if (g0.value() > 0)
         {
-            ppMagf = rAU1f*fvc::interpolate
-            (
-                (1.0/(rho1*(alpha1 + scalar(0.0001))))
-               *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-            );
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
+            // ppMagf = rAU1f*fvc::interpolate
+            // (
+            //     (1.0/(rho1*(alpha1 + scalar(0.0001))))
+            //    *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            // );
+            ppMagf =
+                rAU1f/(alpha1f + scalar(0.0001))
+               *(g0/rho1)*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
 
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
               - fvm::laplacian
                 (
-                    (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                    alpha1f*ppMagf,
                     alpha1,
                     "laplacian(alpha1PpMag,alpha1)"
                 )
diff --git a/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C
index cfc402e14a943938aa50499e7efdde7f9f470b97..6a8aeb9de098056d949073efdc10f0e0aea09f15 100644
--- a/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C
+++ b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,7 +32,7 @@ Description
 #include "volFields.H"
 #include "PatchEdgeFaceWave.H"
 #include "patchEdgeFaceInfo.H"
-#include "patchDist.H"
+#include "patchPatchDist.H"
 
 using namespace Foam;
 
@@ -136,7 +136,7 @@ int main(int argc, char *argv[])
             << " from edges shared with patches " << otherPatchIDs
             << endl;
 
-        patchDist pwd(patch, otherPatchIDs);
+        patchPatchDist pwd(patch, otherPatchIDs);
 
         // Extract as patchField
         volScalarField vsf
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
index ffa2db8b19ff1779d0ca9e25c2e3c22711e51311..6587c8429874e74f9ad64dea668706f24e5ac848 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict
@@ -70,6 +70,9 @@ extrudeModel    linearNormal;
 //- Extrudes into sphere around (0 0 0)
 // extrudeModel   linearRadial;
 
+//- Extrudes into sphere around (0 0 0) with specified radii
+//extrudeModel        radial;
+
 //- Extrudes into sphere with grading according to pressure (atmospherics)
 // extrudeModel   sigmaRadial;
 
@@ -98,6 +101,14 @@ linearDirectionCoeffs
 linearRadialCoeffs
 {
      R              0.1;
+    // Optional inner radius
+    Rsurface        0.01;
+}
+
+radialCoeffs
+{
+    // Radii specified through interpolation table
+    R               table ((0 0.01)(3 0.03)(10 0.1));
 }
 
 sigmaRadialCoeffs
@@ -107,4 +118,5 @@ sigmaRadialCoeffs
     pStrat          1;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
index 15547e985e980721764cbf29367d741a4ea8fd37..2cd1a68ca56dd338283239af10e3a668448a1014 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -69,8 +69,8 @@ structuredCoeffs
     // Renumber in columns (depthFirst) or in layers
     depthFirst true;
 
-    // Optional: reverse ordering
-    //reverse false;
+    // Reverse ordering
+    reverse false;
 }
 
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
index b383a6efe3d6f07576df236a33415dda2b3a380c..6a27bdfb2c624eb1e5b252f7c12edc3914d011ca 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
@@ -119,7 +119,10 @@ structuredCoeffs
 {
     // Patches to do 2D decomposition on. Structured mesh only; cells have
     // to be in 'columns' on top of patches.
-    patches     (bottomPatch);
+    patches     (movingWall);
+
+    // Method to use on the 2D subset
+    method      scotch;
 }
 
 //// Is the case distributed? Note: command-line argument -roots takes
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 9571be6fb003430090800ba535337674560b494e..bcb6083c6abc87bff0ff12d01261fa1cab512daa 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,7 +184,7 @@ int main(int argc, char *argv[])
     argList::addOption
     (
         "sourceTime",
-        "scalar",
+        "scalar|'latestTime'",
         "specify the source time"
     );
     argList::addOption
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H b/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
index 5aec38af8684fc145149df546d0e029e5691e969..843b386e77fd49bda9cec52ae59a69970c1da600 100644
--- a/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,9 +27,6 @@ Class
 Description
     Simplified diagonal-based incomplete LU smoother for asymmetric matrices.
 
-    To improve efficiency, the residual is evaluated after every nSweeps
-    sweeps.
-
 SourceFiles
     DILUSmoother.C
 
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index 5c17f0547f38ab1ff7d6c74df56452004702b6e3..fc2b484cb0a4a478eb2a1591a8f504c5afd74eea 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -367,7 +367,8 @@ public:
         );
 
         //- Construct from reverse addressing: per data item the send
-        //  processor and the receive processor. All processors get same data.
+        //  processor and the receive processor. (note: data is not stored
+        //  sorted per processor so cannot use printLayout).
         mapDistribute
         (
             const labelList& sendProcs,
@@ -635,7 +636,8 @@ public:
             template<class T>
             void receive(PstreamBuffers&, List<T>&) const;
 
-            //- Debug: print layout
+            //- Debug: print layout. Can only be used on maps with sorted
+            //  storage (local data first, then non-local data)
             void printLayout(Ostream& os) const;
 
             //- Correct for topo change.
diff --git a/src/dynamicMesh/createShellMesh/createShellMesh.H b/src/dynamicMesh/createShellMesh/createShellMesh.H
index df1e000be8c90a46ee3391b4a601d669215ef5dc..65edd65cd84b0d9dff43b0a72aab6ae03f0c32fc 100644
--- a/src/dynamicMesh/createShellMesh/createShellMesh.H
+++ b/src/dynamicMesh/createShellMesh/createShellMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,7 +137,7 @@ public:
             //            be in top patch
             //      < 0 : face in opposite orientation as patch face. face will
             //            be in bottom patch
-            //      = 0 : for all side and internal faces
+            //      = 0 : for all side faces
             const labelList& faceToFaceMap() const
             {
                 return faceToFaceMap_;
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
index ac3762db048d5c5476ed42b99e9608d5176ad65b..0329782315c4d4f4af6ef77d1c4da98602b6a6dd 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,48 @@ reconstruct
 
     const fvMesh& mesh = ssf.mesh();
 
+    surfaceVectorField faceVols
+    (
+        mesh.Sf()/(mesh.magSf()*mesh.nonOrthDeltaCoeffs())
+    );
+
+    faceVols.internalField() *= (1.0 -  mesh.weights().internalField());
+    forAll(faceVols.boundaryField(), patchi)
+    {
+        if (faceVols.boundaryField()[patchi].coupled())
+        {
+            faceVols.boundaryField()[patchi] *=
+                (1.0 -  mesh.weights().boundaryField()[patchi]);
+        }
+    }
+
+    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
+    (
+        new GeometricField<GradType, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "volIntegrate("+ssf.name()+')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            inv(surfaceSum(mesh.Sf()*faceVols))&surfaceSum(faceVols*ssf),
+            zeroGradientFvPatchField<GradType>::typeName
+        )
+    );
+
+    treconField().correctBoundaryConditions();
+
+    return treconField;
+}
+/*
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = ssf.mesh();
+
     tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
     (
         new GeometricField<GradType, fvPatchField, volMesh>
@@ -78,6 +120,7 @@ reconstruct
 
     return treconField;
 }
+*/
 
 
 template<class Type>
diff --git a/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H b/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
index 3510e810a4ab831117bf4fa68010d745eda07850..f9fe6ca9e393deb5a967df81c65d7ed2fb98d542 100644
--- a/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
+++ b/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
@@ -59,63 +59,8 @@ inline bool Foam::pointEdgePoint::update
 
     if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
     {
-        // don't propagate small changes unless origins differ
-        scalar originDiff2 = magSqr(origin_-w2.origin());
-        if
-        (
-            originDiff2 < SMALL
-         || ((distSqr_ > SMALL) && (originDiff2/distSqr_ < tol))
-        )
-        {
-            return false;
-        }
-        else
-        {
-            // Different origin. Choose lexical ordering.
-            bool w2less = false;
-            for (direction cmp = 0; cmp < vector::nComponents; cmp++)
-            {
-                scalar d = w2.origin()[cmp]-origin_[cmp];
-                scalar magD = mag(d);
-
-                if
-                (
-                    magD < SMALL
-                 || ((distSqr_ > SMALL) && (magD/distSqr_ > tol))
-                )
-                {
-                    // Small difference. Test next component.
-                }
-                else if (d < 0)
-                {
-                    w2less = true;
-                    break;
-                }
-                else if (d > 0)
-                {
-                    break;
-                }
-            }
-
-            if (w2less)
-            {
-                // update with new values
-                distSqr_ = dist2;
-                origin_ = w2.origin();
-                return true;
-            }
-            else
-            {
-                Pout<< "** same distance. at:" << pt << nl
-                    << "    old distance:" << Foam::sqr(distSqr_) << nl
-                    << "    new distance:" << Foam::sqr(dist2) << nl
-                    << "    old origin:" << origin_ << nl
-                    << "    new origin:" << w2.origin() << nl
-                    << endl;
-                Pout<< "** w2 looses." << endl;
-                return false;
-            }
-        }
+        // don't propagate small changes
+        return false;
     }
     else
     {
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index 73d64b40f3e8c346f47ef372c3d905d93ef5d86d..304501870181cf88709362df778bd3f7d7c7f714 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ License
 #include "OFstream.H"
 #include "Random.H"
 #include "treeDataFace.H"
+#include "treeDataPoint.H"
 #include "indexedOctree.H"
 #include "polyMesh.H"
 #include "polyPatch.H"
@@ -38,6 +39,7 @@ License
 #include "mapDistribute.H"
 #include "SubField.H"
 #include "triPointRef.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,12 +51,13 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::mappedPatchBase::sampleMode,
-        4
+        5
     >::names[] =
     {
         "nearestCell",
         "nearestPatchFace",
         "nearestPatchFaceAMI",
+        "nearestPatchPoint",
         "nearestFace"
     };
 
@@ -72,7 +75,7 @@ namespace Foam
 }
 
 
-const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 4>
+const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 5>
     Foam::mappedPatchBase::sampleModeNames_;
 
 const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
@@ -315,6 +318,77 @@ void Foam::mappedPatchBase::findSamples
             break;
         }
 
+        case NEARESTPATCHPOINT:
+        {
+            Random rndGen(123456);
+
+            const polyPatch& pp = samplePolyPatch();
+
+            if (pp.empty())
+            {
+                forAll(samples, sampleI)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            else
+            {
+                // patch (local) points
+                treeBoundBox patchBb
+                (
+                    treeBoundBox(pp.points(), pp.meshPoints()).extend
+                    (
+                        rndGen,
+                        1e-4
+                    )
+                );
+                patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+                patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+                indexedOctree<treeDataPoint> boundaryTree
+                (
+                    treeDataPoint   // all information needed to search faces
+                    (
+                        mesh.points(),
+                        pp.meshPoints() // selection of points to search on
+                    ),
+                    patchBb,        // overall search domain
+                    8,              // maxLevel
+                    10,             // leafsize
+                    3.0             // duplicity
+                );
+
+                forAll(samples, sampleI)
+                {
+                    const point& sample = samples[sampleI];
+
+                    pointIndexHit& nearInfo = nearest[sampleI].first();
+                    nearInfo = boundaryTree.findNearest
+                    (
+                        sample,
+                        magSqr(patchBb.span())
+                    );
+
+                    if (!nearInfo.hit())
+                    {
+                        nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                        nearest[sampleI].second().second() =
+                            Pstream::myProcNo();
+                    }
+                    else
+                    {
+                        const point& pt = nearInfo.hitPoint();
+
+                        nearest[sampleI].second().first() = magSqr(pt-sample);
+                        nearest[sampleI].second().second() =
+                            Pstream::myProcNo();
+                    }
+                }
+            }
+            break;
+        }
+
         case NEARESTFACE:
         {
             if (samplePatch_.size() && samplePatch_ != "none")
@@ -373,10 +447,11 @@ void Foam::mappedPatchBase::findSamples
     }
 
 
-    // Find nearest
+    // Find nearest. Combine on master.
     Pstream::listCombineGather(nearest, nearestEqOp());
     Pstream::listCombineScatter(nearest);
 
+
     if (debug)
     {
         Info<< "mappedPatchBase::findSamples on mesh " << sampleRegion_
@@ -388,8 +463,9 @@ void Foam::mappedPatchBase::findSamples
 
             Info<< "    " << sampleI << " coord:"<< samples[sampleI]
                 << " found on processor:" << procI
-                << " in local cell/face:" << localI
-                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
+                << " in local cell/face/point:" << localI
+                << " with location:" << nearest[sampleI].first().rawPoint()
+                << endl;
         }
     }
 
@@ -599,7 +675,8 @@ void Foam::mappedPatchBase::calcMapping() const
           + "_mapped.obj"
         );
         Pout<< "Dumping mapping as lines from patch faceCentres to"
-            << " sampled cell/faceCentres to file " << str.name() << endl;
+            << " sampled cell/faceCentres/points to file " << str.name()
+            << endl;
 
         label vertI = 0;
 
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
index f9afc00cdff0a4794cb837219361188a4d66e22e..f72f407d8102a991e070229c3b1d34cc6c097a31 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,10 @@ Description
                                    - patches need not conform
                                    - uses AMI interpolation
         // - nearestFace         : nearest boundary face on any patch
+        // - nearestPatchPoint   : nearest patch point (for coupled points
+        //                         this might be any of the points so you have
+        //                         to guarantee the point data is synchronised
+        //                         beforehand)
         sampleMode nearestCell;
 
         // If sampleMod is nearestPatchFace : patch to find faces of
@@ -99,8 +103,9 @@ public:
         enum sampleMode
         {
             NEARESTCELL,         // nearest cell
-            NEARESTPATCHFACE,    // faces on selected patch
+            NEARESTPATCHFACE,    // nearest face on selected patch
             NEARESTPATCHFACEAMI, // nearest patch face + AMI interpolation
+            NEARESTPATCHPOINT,   // nearest point on selected patch
             NEARESTFACE          // nearest face
         };
 
@@ -112,7 +117,7 @@ public:
             NORMAL              // use face normal + distance
         };
 
-        static const NamedEnum<sampleMode, 4> sampleModeNames_;
+        static const NamedEnum<sampleMode, 5> sampleModeNames_;
 
         static const NamedEnum<offsetMode, 3> offsetModeNames_;
 
@@ -145,6 +150,27 @@ public:
         }
     };
 
+    class maxProcEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first().hit())
+            {
+                if (!x.first().hit())
+                {
+                    x = y;
+                }
+                else if (y.second().second() > x.second().second())
+                {
+                    x = y;
+                }
+            }
+        }
+    };
+
 
 protected:
 
@@ -159,7 +185,7 @@ protected:
         //- What to sample
         const sampleMode mode_;
 
-        //- Patch (only if NEARESTPATCHFACE)
+        //- Patch (if in sampleMode NEARESTPATCH*)
         const word samplePatch_;
 
         //- How to obtain samples
@@ -187,7 +213,7 @@ protected:
             mutable autoPtr<mapDistribute> mapPtr_;
 
 
-        // AMI interpolator
+        // AMI interpolator (only for NEARESTPATCHFACEAMI)
 
             //- Pointer to AMI interpolator
             mutable autoPtr<AMIPatchToPatchInterpolation> AMIPtr_;
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
index 928fa86ad7ddf26b31e342cdaa16e4d15856525c..927534210f1079d551dde3be37d7ff2bcd99a9a8 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,10 @@ inline Foam::label Foam::mappedPatchBase::sampleSize() const
         {
             return samplePolyPatch().size();
         }
+        case NEARESTPATCHPOINT:
+        {
+            return samplePolyPatch().nPoints();
+        }
         case NEARESTFACE:
         {
             const polyMesh& mesh = sampleMesh();
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceDataI.H b/src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceDataI.H
index 9ce5c90c5f621135f8a6fe96cf780e3ffc929c98..388c923a2ae0d63808b88bf3888c71b2987a3fcb 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceDataI.H
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceDataI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -124,8 +124,7 @@ inline bool Foam::topoDistanceData::updateCell
 {
     if (distance_ == -1)
     {
-        data_ = neighbourInfo.data_;
-        distance_ = neighbourInfo.distance_ + 1;
+        operator=(neighbourInfo);
         return true;
     }
     else
@@ -151,7 +150,8 @@ inline bool Foam::topoDistanceData::updateFace
 
     if (distance_ == -1)
     {
-        operator=(neighbourInfo);
+        data_ = neighbourInfo.data_;
+        distance_ = neighbourInfo.distance_ + 1;
         return true;
     }
     else
diff --git a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
index 0593e2d49f4b0a6adfbecb791dcb5a279ea0c069..b0a13a83997061bddae86c0c2520c82d8071ea01 100644
--- a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
+++ b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,7 +198,7 @@ Foam::labelList Foam::structuredRenumber::renumber
         }
         else
         {
-            label layerI = cellData[cellI].distance()-1;
+            label layerI = cellData[cellI].distance();
             if (depthFirst_)
             {
                 orderedToOld[nLayers*cellData[cellI].data()+layerI] = cellI;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
index b07ba8e4080e6d222a80228aa49f6880ee0f895e..25bc08d671014e05362550e22ecd36921d254728 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      T1;
+    object      Tair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -27,7 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
-        phi                phi1;
+        phi                phiair;
         inletValue         $internalField;
         value              $internalField;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
index cc4307a53689ef43412a13c49df25ccc731623c6..f0f366c7dc27735b261b92c9eb7e2570dcaff136 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      T2;
+    object      Twater;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -27,7 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
-        phi                phi2;
+        phi                phiwater;
         inletValue         uniform 300;
         value              $internalField;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
index 2cf5dec2bc7111419b65e3a82852df94d550df37..ac6020947e16d1537f7ed05c0878183567bac15c 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      binary;
     class       volVectorField;
-    object      U1;
+    object      Uair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type               pressureInletOutletVelocity;
-        phi                phi1;
+        phi                phiair;
         value              $internalField;
     }
     walls
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
index 68e9f354a6e5cab2f825bfd1e3b28d601d791e4d..22ed59a0ef3e78143b36a455a5c2d9f41f09df3a 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      binary;
     class       volVectorField;
-    object      U2;
+    object      Uwater;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type               pressureInletOutletVelocity;
-        phi                phi2;
+        phi                phiwater;
         value              $internalField;
     }
     walls
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
similarity index 99%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
index 9a514b78074968c1df0cd3b7baa6ae041996ef8a..a2224f95c3c32f38ce0e3e0dd361b3fabc5bf60f 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
@@ -1908,7 +1908,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
-        phi             phi1;
+        phi             phiair;
         inletValue      uniform 1;
         value           uniform 1;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
index 4b871e21b9da4e1ddd6159564e28dbf7f303b297..62088eca2323cb4b66162a3131679877f60bb945 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alphaair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
-        phi             phi1;
+        phi             phiair;
         inletValue      uniform 1;
         value           uniform 1;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
index 03a3a667a3d8c1ee3e013766f338b7928053cabc..5964adcedc85ce4cfb8757b6dd2f652ba6cd9b88 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dragModel1          SchillerNaumann;
-dragModel2          SchillerNaumann;
+dragModelair            SchillerNaumann;
+dragModelwater          SchillerNaumann;
 
-heatTransferModel1  RanzMarshall;
-heatTransferModel2  RanzMarshall;
+heatTransferModelair    RanzMarshall;
+heatTransferModelwater  RanzMarshall;
 
-dispersedPhase      both;
+dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
 residualSlip            1e-2;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
index 56e0a545c15839416fcd3beca616ce0f609c5879..bf47f69643c9925d3a1ef19c6b4ddc67cf604e0a 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
@@ -38,6 +38,7 @@ FoamFile
     defaultFaces
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3750;
         startFace       3850;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair
new file mode 100644
index 0000000000000000000000000000000000000000..11c033af59bbfbf08d2b11244e50f0fd00c1542f
--- /dev/null
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater
new file mode 100644
index 0000000000000000000000000000000000000000..672b24a98b365f170f3f2567723294befd79a06e
--- /dev/null
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    equationOfState
+    {
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
index d469dc706f932ef77a08bd0aa00f6cd9ec867c91..408d5f97d66a726605f89ae7c8b5a3e18711f2a6 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
@@ -15,16 +15,10 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phase1
-{
-    rho             1;
-    rho0            0;
-    R               287;
-    Cp              1007;
-
-    nu              1.589e-05;
-    kappa           2.63e-2;
+phases (air water);
 
+air
+{
     diameterModel   isothermal;
     isothermalCoeffs
     {
@@ -33,15 +27,9 @@ phase1
     }
 }
 
-phase2
+water
 {
-    rho             1027;
-    rho0            1027;
-    R               1e10;
-    Cp              4195;
-
-    nu              3.55e-7;
-    kappa           0.668;
+    //R               1e10;
 
     diameterModel constant;
     constantCoeffs
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
index f28ade7e84c4a5d839fed48c74b21ffbdca6b371..24d66f91c7d44dfeb7b24b508c2593511edb6ed9 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
@@ -60,21 +60,21 @@ functions
         outputControl   outputTime;
         fields
         (
-            U1
+            Uair
             {
                  mean        on;
                  prime2Mean  off;
                  base        time;
             }
 
-            U2
+            Uwater
             {
                  mean        on;
                  prime2Mean  off;
                  base        time;
             }
 
-            alpha1
+            alphaair
             {
                  mean        on;
                  prime2Mean  off;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
index 34b14edfcce7676bc1e185e87f785fc2ab086c56..d53ef42583a230181a0d5bc0a16d73acdc017ce2 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -29,23 +29,23 @@ divSchemes
 {
     default             none;
 
-    div(phi,alpha1)     Gauss limitedLinear01 1;
-    div(phir,alpha1)    Gauss limitedLinear01 1;
-    div(alphaPhi1,U1)   Gauss limitedLinearV 1;
-    div(alphaPhi2,U2)   Gauss limitedLinearV 1;
-    div(phi1,U1)        Gauss limitedLinearV 1;
-    div(phi2,U2)        Gauss limitedLinearV 1;
-    div((alpha1*Rc1))   Gauss linear;
-    div((alpha2*Rc2))   Gauss linear;
-    div(alphaPhi1,T1)   Gauss limitedLinear 1;
-    div(alphaPhi2,T2)   Gauss limitedLinear 1;
-    div(alphaPhi2,k)    Gauss limitedLinear 1;
-    div(alphaPhi2,epsilon) Gauss limitedLinear 1;
+    div(phi,alpha)     Gauss limitedLinear01 1;
+    div(phir,alpha)    Gauss limitedLinear01 1;
+    div(alphaPhi,Uair)   Gauss limitedLinearV 1;
+    div(alphaPhi,Uwater)   Gauss limitedLinearV 1;
+    div(phiair,Uair)        Gauss limitedLinearV 1;
+    div(phiwater,Uwater)        Gauss limitedLinearV 1;
+    div((alphaair*Rc))   Gauss linear;
+    div((alphawater*Rc))   Gauss linear;
+    div(alphaPhi,hair)   Gauss limitedLinear 1;
+    div(alphaPhi,hwater)   Gauss limitedLinear 1;
+    div(alphaPhiwater,k)    Gauss limitedLinear 1;
+    div(alphaPhiwater,epsilon) Gauss limitedLinear 1;
     div(phi,Theta)      Gauss limitedLinear 1;
-    div(phid1,p)        Gauss upwind;
-    div(phid2,p)        Gauss upwind;
-    div(phi1,K1)        Gauss limitedLinear 1;
-    div(phi2,K2)        Gauss limitedLinear 1;
+    div(phidair,p)        Gauss upwind;
+    div(phidwater,p)        Gauss upwind;
+    div(phiair,Kair)        Gauss limitedLinear 1;
+    div(phiwater,Kwater)        Gauss limitedLinear 1;
 }
 
 laplacianSchemes
@@ -67,7 +67,7 @@ fluxRequired
 {
     default         no;
     p               ;
-    alpha1          ;
+    alphaair        ;
 }
 
 
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
index d09194f954dde0ac71cc6201ac1cb1f475455c51..96e43b3c1d253411d40e8e71b59a2c071ad2f124 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
@@ -47,10 +47,10 @@ solvers
         relTol          0;
     }
 
-    "T.*"
+    "h.*"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        solver          PCG; //PBiCG;
+        preconditioner  DIC; //DILU;
         tolerance       1e-8;
         relTol          0;
     }
@@ -105,7 +105,7 @@ relaxationFactors
     equations
     {
         "U.*"           1;
-        "T.*"           1;
+        "h.*"           1;
         "alpha.*"       1;
         "Theta.*"       1;
         "k.*"           1;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
index 5b6d65bbecb7870ba8626ebd335c4a8f8d75d8d9..85996cf966762c6ca3a7c37a1eaa8ae462ecdb19 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
@@ -27,7 +27,7 @@ regions
         box (0 0 -0.1) (0.15 0.701 0.1);
         fieldValues
         (
-            volScalarFieldValue alpha1 0
+            volScalarFieldValue alphaair 0
         );
     }
 );