{
    tmp<fv::convectionScheme<scalar>> hConvection
    (
        fv::convectionScheme<scalar>::New
        (
            mesh,
            phi,
            mesh.divScheme("div(phi,Yi_h)")
        )
    );

    typedef multiComponentMixture<gasHThermoPhysics> MMix;

    MMix &janafComposition = dynamic_cast<MMix&>(composition);

    volScalarField& he = thermo.he();

    volScalarField hsi(he);

    volVectorField hDiffusionSrc("rho*Sum(hYV)", he * Vc * 0.0);

    forAll(Y, i)
    {
        const MMix::thermoType &tti = janafComposition.getLocalThermo(i);

        forAll(hsi, cellI)
        {
            const scalar pI = p[cellI];
            const scalar TI = T[cellI];

            hsi[cellI] = tti.HE(pI,TI);
        }

        forAll(hsi.boundaryFieldRef(), patchI)
        {
            volScalarField::Patch &hsiP = hsi.boundaryFieldRef()[patchI];
            const volScalarField::Patch &pP = p.boundaryField()[patchI];
            const volScalarField::Patch &TP = T.boundaryField()[patchI];

            forAll(hsiP, faceI)
            {
                const scalar pI = pP[faceI];
                const scalar TI = TP[faceI];

                hsiP[faceI] = tti.HE(pI,TI);
            }
        }

        hDiffusionSrc += hsi * (- diff.D(i) * fvc::grad(Y[i]) + Y[i] * Vc);
    }

    tmp<volVectorField> tTauU(diff.mu()*(U&(Foam::dev2(Foam::T(gradU)) + gradU)));

    hDiffusionSrc *= rho;

    volVectorField hWorkSrc1(rho * U);

    volVectorField hWorkSrc2(rho * Vc);

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + hConvection->fvmDiv(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      - dpdt
      - fvm::laplacian(diff.k()/thermo.Cp(), he)
      + fvc::div(hDiffusionSrc)
     ==
        reaction->Sh()
      + fvc::div(tTauU)
      + fvOptions(rho, he)
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct();

    Info<< "min/max(T) = "
        << min(T).value() << ", " << max(T).value() << endl;
}