diff --git a/applications/solvers/multiphase/MPPICInterFoam/Make/options b/applications/solvers/multiphase/MPPICInterFoam/Make/options
index abc0537451c6337c71984756ddd167241eb2a886..0ad31d6975c01583c8f4eba520a8518038bf5453 100644
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@@ -2,9 +2,9 @@ interFoamPath = $(FOAM_SOLVERS)/multiphase/interFoam
 
 EXE_INC =  \
     -I. \
+    -I../VoF \
     -I./IncompressibleTwoPhaseMixtureTurbulenceModels/lnInclude \
     -I$(interFoamPath) \
-    -I../VoF \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
diff --git a/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3ce881bd06b3c6f7b917d00babd3df5d080816ef
--- /dev/null
+++ b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H
@@ -0,0 +1,258 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Set the off-centering coefficient according to ddt scheme
+    scalar ocCoeff = 0;
+    {
+        tmp<fv::ddtScheme<scalar>> tddtAlpha
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh,
+                mesh.ddtScheme("ddt(alpha)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
+        {
+            if (nAlphaSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                alphaRestart
+             || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
+                   .ocCoeff();
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Only Euler and CrankNicolson ddt schemes are supported"
+                << exit(FatalError);
+        }
+    }
+
+    // Set the time blending factor, 1 for Euler
+    scalar cnCoeff = 1.0/(1.0 + ocCoeff);
+
+    // Standard face-flux compression coefficient
+    surfaceScalarField phic(mixture.cAlpha()*mag(alphaPhic/mesh.magSf()));
+
+    // Add the optional isotropic compression contribution
+    if (icAlpha > 0)
+    {
+        phic *= (1.0 - icAlpha);
+        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+    }
+
+    surfaceScalarField::Boundary& phicBf =
+        phic.boundaryFieldRef();
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phicBf[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalarField> phiCN(alphaPhic);
+
+    // Calculate the Crank-Nicolson off-centred volumetric flux
+    if (ocCoeff > 0)
+    {
+        phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime();
+    }
+
+    if (MULESCorr)
+    {
+        #include "alphaSuSp.H"
+
+        fvScalarMatrix alpha1Eqn
+        (
+            (
+                LTS
+              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
+              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+            )
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phiCN,
+                upwind<scalar>(mesh, phiCN)
+            ).fvmDiv(phiCN, alpha1)
+         - fvm::Sp(fvc::ddt(alphac) + fvc::div(phiCN), alpha1)
+         ==
+            Su + fvm::Sp(Sp + divU, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+            << endl;
+
+        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        alphaPhi = talphaPhiUD();
+
+        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                alphaPhi,
+                talphaPhiCorr0.ref(),
+                zeroField(), zeroField(),
+                1, 0
+            );
+
+            alphaPhi += talphaPhiCorr0();
+        }
+
+        // Cache the upwind-flux
+        talphaPhiCorr0 = talphaPhiUD;
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        #include "alphaSuSp.H"
+
+        surfaceScalarField phir(phic*mixture.nHatf());
+
+        tmp<surfaceScalarField> talphaPhiUn
+        (
+            fvc::flux
+            (
+                phiCN(),
+                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            volScalarField alpha10("alpha10", alpha1);
+
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                talphaPhiUn(),
+                talphaPhiCorr.ref(),
+                Sp,
+                (-Sp*alpha1)(),
+                1,
+                0
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                alphaPhi += talphaPhiCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                alphaPhi += 0.5*talphaPhiCorr();
+            }
+        }
+        else
+        {
+            alphaPhi = talphaPhiUn;
+
+            MULES::explicitSolve
+            (
+                alphac,
+                alpha1,
+                phiCN,
+                alphaPhi,
+                Sp,
+                (Su + divU*min(alpha1(), scalar(1)))(),
+                1,
+                0
+            );
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+    }
+    else
+    {
+        talphaPhiCorr0.clear();
+    }
+
+    if
+    (
+        word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::EulerDdtScheme<vector>::typeName
+    )
+    {
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
+    }
+    else
+    {
+        if (ocCoeff > 0)
+        {
+            // Calculate the end-of-time-step alpha flux
+            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+        }
+
+        // Calculate the end-of-time-step mass flux
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + alphaPhic*rho2f;
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+        << endl;
+}
diff --git a/applications/solvers/multiphase/MPPICInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/MPPICInterFoam/alphaEqnSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..6820b2e596f38d7335013da8d4a1d3559f23541e
--- /dev/null
+++ b/applications/solvers/multiphase/MPPICInterFoam/alphaEqnSubCycle.H
@@ -0,0 +1,43 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;
+mu = mixture.mu();
+
diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
index 60dd2aecc812fef9ebbfd9edd802c326950f84b3..af74f0362a4ec21ce61cf78d7aaa88a08ab5b256 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
@@ -5,7 +5,7 @@
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
-            fvc::div(fvc::absolute(phi, U), p)
+            divU*p
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 34e8762709abfd0644471c5dcd4cb16cef4238f4..e3d43a8944ff9289cba853549b81ff9c1f19d1f5 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -167,6 +167,12 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = alpha1*rho1 + alpha2*rho2;
+
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
+
         runTime.write();
 
         Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index a657caf01ee388fa595b0bc4364dd643e724ab6b..dee1a578f25cea83ed57808459b85a58c5749819 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) OpenCFD Ltd. 2017
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,6 +113,7 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index bcf8358a105009fcb5e900b15ff123b95696af6a..7c530d79e2f82fa80abf6901c0e1acfc8048397e 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -350,9 +350,6 @@ int main(int argc, char *argv[])
 
         triSurface surf = surfPtr();
 
-        Info<< "NB: Feature line extraction is only valid on closed manifold "
-            << "surfaces." << nl;
-
         Info<< nl
             << "Statistics:" << nl;
         surf.writeStats(Info);
@@ -572,7 +569,9 @@ int main(int argc, char *argv[])
             triSurfaceSearch query(surf);
             surfaceIntersection intersect(query, surfaceDict);
 
-            // Remove rounding noise - could make adjustable
+            // Remove rounding noise. For consistency could use 1e-6,
+            // as per extractFromFile implementation
+
             intersect.mergePoints(10*SMALL);
 
             labelPair sizeInfo
diff --git a/src/OpenFOAM/global/constants/mathematical/mathematicalConstants.H b/src/OpenFOAM/global/constants/mathematical/mathematicalConstants.H
index fa4850fb5a98ec9728bfcb70087bb22272ed1f5f..c7faae1bd273a4ab30a5703f0c17699e646e7611 100644
--- a/src/OpenFOAM/global/constants/mathematical/mathematicalConstants.H
+++ b/src/OpenFOAM/global/constants/mathematical/mathematicalConstants.H
@@ -47,13 +47,13 @@ namespace mathematical
 
     static const char* const group = "mathematical";
 
-    const scalar e(M_E);
-    const scalar pi(M_PI);
-    const scalar twoPi(2*pi);
-    const scalar piByTwo(0.5*pi);
+    constexpr scalar e(M_E);
+    constexpr scalar pi(M_PI);
+    constexpr scalar twoPi(2*M_PI);
+    constexpr scalar piByTwo(0.5*M_PI);
 
     //- Euler's constant
-    const scalar Eu(0.57721566490153286060651209);
+    constexpr scalar Eu(0.57721566490153286060651209);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/global/unitConversion/unitConversion.H b/src/OpenFOAM/global/unitConversion/unitConversion.H
index 10263f2a61d90403a8f2f354eb11ca01169bd863..499ecd6aab8eb29825c866a1ffaad44ec2cdb69d 100644
--- a/src/OpenFOAM/global/unitConversion/unitConversion.H
+++ b/src/OpenFOAM/global/unitConversion/unitConversion.H
@@ -44,13 +44,13 @@ namespace Foam
 //- Conversion from degrees to radians
 inline constexpr scalar degToRad(const scalar deg) noexcept
 {
-    return (deg*Foam::constant::mathematical::pi/180.0);
+    return (deg*M_PI/180.0);
 }
 
 //- Conversion from radians to degrees
 inline constexpr scalar radToDeg(const scalar rad) noexcept
 {
-    return (rad*180.0/Foam::constant::mathematical::pi);
+    return (rad*180.0/M_PI);
 }
 
 //- Conversion from atm to Pa
@@ -69,13 +69,13 @@ inline constexpr scalar paToAtm(const scalar pa) noexcept
 //- User literal for degrees to radians conversion (integers)
 inline constexpr scalar operator "" _deg(unsigned long long int deg) noexcept
 {
-    return (deg*Foam::constant::mathematical::pi/180.0);
+    return (deg*M_PI/180.0);
 }
 
 //- User literal for degrees to radians conversion (floats)
 inline constexpr scalar operator "" _deg(long double deg) noexcept
 {
-    return (deg*Foam::constant::mathematical::pi/180.0);
+    return (deg*M_PI/180.0);
 }
 
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
index e8fad80ee703f74bd6b8e278da42bd3e51ee2e1f..e232b8198f4f90a9969a5c983f0364e0b5151770 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
@@ -344,9 +344,11 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::write
 {
     mixedFvPatchScalarField::write(os);
     os.writeKeyword("Tnbr")<< TnbrName_ << token::END_STATEMENT << nl;
+
     os.writeKeyword("qrNbr")<< qrNbrName_ << token::END_STATEMENT << nl;
     os.writeKeyword("qr")<< qrName_ << token::END_STATEMENT << nl;
-    os.writeEntry("thermalInertia", thermalInertia_);
+    os.writeKeyword("thermalInertia")<< thermalInertia_
+        << token::END_STATEMENT << nl;
 
     thicknessLayers_.writeEntry("thicknessLayers", os);
     kappaLayers_.writeEntry("kappaLayers", os);
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
index bffd5f8d2f2ed5737f37063e81d72951755c5ac5..18d2a9367c88b7156f00416475d4142de400cacb 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,7 +67,7 @@ Foam::solidReaction<ReactionThermo>::solidReaction
     const dictionary& dict
 )
 :
-    Reaction<ReactionThermo>(species, thermoDatabase, dict),
+    Reaction<ReactionThermo>(species, thermoDatabase, dict, false),
     pyrolisisGases_(dict.parent().parent().lookup("gaseousSpecies")),
     glhs_(),
     grhs_()
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index e994ff1dce0825a8625875b4e5ffec744fd1559b..645cc4fc806b2b11734b0690ab481863da23bf31 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -111,19 +111,20 @@ void Foam::Reaction<ReactionThermo>::setThermo
     const HashPtrTable<ReactionThermo>& thermoDatabase
 )
 {
+
     typename ReactionThermo::thermoType rhsThermo
     (
         rhs_[0].stoichCoeff
-       *(*thermoDatabase[species_[rhs_[0].index]]).W()
-       *(*thermoDatabase[species_[rhs_[0].index]])
+        *(*thermoDatabase[species_[rhs_[0].index]]).W()
+        *(*thermoDatabase[species_[rhs_[0].index]])
     );
 
     for (label i=1; i<rhs_.size(); ++i)
     {
         rhsThermo +=
             rhs_[i].stoichCoeff
-           *(*thermoDatabase[species_[rhs_[i].index]]).W()
-           *(*thermoDatabase[species_[rhs_[i].index]]);
+        *(*thermoDatabase[species_[rhs_[i].index]]).W()
+        *(*thermoDatabase[species_[rhs_[i].index]]);
     }
 
     typename ReactionThermo::thermoType lhsThermo
@@ -154,7 +155,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     const speciesTable& species,
     const List<specieCoeffs>& lhs,
     const List<specieCoeffs>& rhs,
-    const HashPtrTable<ReactionThermo>& thermoDatabase
+    const HashPtrTable<ReactionThermo>& thermoDatabase,
+    bool initReactionThermo
 )
 :
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
@@ -163,7 +165,10 @@ Foam::Reaction<ReactionThermo>::Reaction
     lhs_(lhs),
     rhs_(rhs)
 {
-    setThermo(thermoDatabase);
+    if (initReactionThermo)
+    {
+        setThermo(thermoDatabase);
+    }
 }
 
 
@@ -325,7 +330,8 @@ Foam::Reaction<ReactionThermo>::Reaction
 (
     const speciesTable& species,
     const HashPtrTable<ReactionThermo>& thermoDatabase,
-    const dictionary& dict
+    const dictionary& dict,
+    bool initReactionThermo
 )
 :
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
@@ -339,7 +345,11 @@ Foam::Reaction<ReactionThermo>::Reaction
         lhs_,
         rhs_
     );
-    setThermo(thermoDatabase);
+
+    if (initReactionThermo)
+    {
+        setThermo(thermoDatabase);
+    }
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index d96134aac9523d59d4b0fdb77db2663123026eab..7b3b8ea571e823e22ecb26e69770ef9ae0c5694f 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -186,18 +186,27 @@ public:
             const speciesTable& species,
             const List<specieCoeffs>& lhs,
             const List<specieCoeffs>& rhs,
-            const HashPtrTable<ReactionThermo>& thermoDatabase
+            const HashPtrTable<ReactionThermo>& thermoDatabase,
+            bool initReactionThermo = true
         );
 
         //- Construct as copy given new speciesTable
-        Reaction(const Reaction<ReactionThermo>&, const speciesTable& species);
+        Reaction
+        (
+            const Reaction<ReactionThermo>&,
+            const speciesTable& species
+        );
 
         //- Construct from dictionary
+        // NOTE: initReactionThermo is used by solidReaction where there is no
+        // need of setting a lhs - rhs thermo type for each reaction. This is
+        // needed for mechanism with reversible reactions
         Reaction
         (
             const speciesTable& species,
             const HashPtrTable<ReactionThermo>& thermoDatabase,
-            const dictionary& dict
+            const dictionary& dict,
+            bool initReactionThermo = true
         );
 
         //- Construct and return a clone
diff --git a/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/Qr b/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/qr
similarity index 98%
rename from tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/Qr
rename to tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/qr
index 654cb3ae748ca6c766c59281d32618d2c859c3c7..ba9e5916fb1ec647ae13cf9a0881323c73dad3b2 100644
--- a/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/Qr
+++ b/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/qr
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      Qr;
+    object      qr;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones b/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
index a099495670d023da0a3f793e0de0533b8203f3a4..f5fd1efacc503a9dee9d6980ecc49005c85c105d 100644
--- a/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
+++ b/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
@@ -29,7 +29,7 @@ FoamFile
 
             filmCoupled     false;
 
-            radFluxName     Qr;
+            qrHSource       on;
 
             minimumDelta    1e-6;
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
index 22ec9ceabe4afb249989673cd7072899f7f3fc78..c22fa78d78a3c02875c27e37baa127e08b3853a1 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 
 odeCoeffs
 {
-    solver          Rosenbrock43;
+    solver          Rosenbrock34;
     absTol          1e-12;
     relTol          0.01;
 }
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/reactions b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/reactions
index e548bae90ffdb1e3d0822d8af80a19f9e2552636..842baa512bed455cba2fafa959164501aa546b47 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/reactions
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/reactions
@@ -1,3 +1,12 @@
+elements
+(
+O
+C
+H
+N
+);
+
+
 species
 (
     O2
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
index 371a3bb6c87f57d6bde8d7e4b6d585c89893168b..4bbb98aa34305e8df8340d9820a4d8ebe928602e 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
@@ -21,6 +21,10 @@ O2
     {
         molWeight       31.9988;
     }
+    elements
+    {
+        O       2;
+    }
     thermodynamics
     {
         Tlow            200;
@@ -42,6 +46,11 @@ H2O
     {
         molWeight       18.0153;
     }
+    elements
+    {
+        O       1;
+        H       2;
+    }
     thermodynamics
     {
         Tlow            200;
@@ -63,6 +72,11 @@ CH4
     {
         molWeight       16.0428;
     }
+    elements
+    {
+        C       1;
+        H       4;
+    }
     thermodynamics
     {
         Tlow            200;
@@ -84,6 +98,11 @@ CO2
     {
         molWeight       44.01;
     }
+    elements
+    {
+        C       1;
+        O       2;
+    }
     thermodynamics
     {
         Tlow            200;
@@ -105,6 +124,10 @@ N2
     {
         molWeight       28.0134;
     }
+    elements
+    {
+        N       2;
+    }
     thermodynamics
     {
         Tlow            200;