diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H
index 15dc7df9..55753a0b 100644
--- a/applications/solvers/dfHighSpeedFoam/createFields.H
+++ b/applications/solvers/dfHighSpeedFoam/createFields.H
@@ -103,6 +103,19 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField pMax
+(
+    IOobject
+    (
+        "pMax",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.p()
+);
+
 volVectorField rhoU
 (
     IOobject
@@ -190,7 +203,7 @@ const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>("combus
 Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
 
 
-if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
+if (combModelName != "flareFGM")
 {
     chemistry->correctThermo();
     Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C
index 8ba02a46..e4487df6 100644
--- a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C
+++ b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C
@@ -58,6 +58,9 @@ Description
 #include "PstreamGlobals.H"
 #include "CombustionModel.H"
 
+#include "fluxSchemes/AUSMDVFlux.H"
+#include "fluxSchemes/KNPFlux.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -104,6 +107,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+        pMax = max(p, pMax);
+
         #include "readTimeControls.H"
 
         //used for AMR
diff --git a/applications/solvers/dfHighSpeedFoam/fluxSchemes/AUSMDVFlux.H b/applications/solvers/dfHighSpeedFoam/fluxSchemes/AUSMDVFlux.H
new file mode 100644
index 00000000..18deac0f
--- /dev/null
+++ b/applications/solvers/dfHighSpeedFoam/fluxSchemes/AUSMDVFlux.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+Function
+    AUSMDVFlux
+Description
+    A function for solving convective flux based on AUSMDV, which mixes the AUSMD and AUSMV scheme
+Author
+    Daoping Zhang (zdp19980509@163.com)
+\*---------------------------------------------------------------------------*/
+
+void AUSMDVFlux(
+    const surfaceScalarField& meshPhi,
+    const surfaceScalarField& magSf,
+    const surfaceVectorField& normal,
+    const surfaceScalarField& rhoOwn,
+    const surfaceScalarField& rhoNei,
+    const surfaceVectorField& UOwn,
+    const surfaceVectorField& UNei,
+    const surfaceScalarField& HOwn,
+    const surfaceScalarField& HNei,
+    const surfaceScalarField& pOwn,
+    const surfaceScalarField& pNei,
+    const surfaceScalarField& gammaOwn,
+    const surfaceScalarField& gammaNei,
+    const surfaceScalarField& aOwn,
+    const surfaceScalarField& aNei,
+    surfaceVectorField& U12,
+    surfaceScalarField& massFlux,
+    surfaceVectorField& momentumFlux,
+    surfaceScalarField& energyFlux
+)
+{
+
+    dimensionedScalar norP("norP", dimless, 1.0);
+    //dimensionedScalar minU("minU", dimVelocity, SMALL);
+
+    surfaceScalarField UvOwn((UOwn & normal) - meshPhi/magSf);
+
+    surfaceScalarField UvNei((UNei & normal) - meshPhi/magSf);
+
+    // Compute split velocity
+    surfaceScalarField alphaOwn(2*(pOwn/rhoOwn)/(pOwn/rhoOwn + pNei/rhoNei));
+    surfaceScalarField alphaNei(2 - alphaOwn);
+
+    surfaceScalarField cm(max(aOwn , aNei));
+
+    surfaceScalarField uPlus(
+      neg0(mag(UvOwn/cm) - 1)*(alphaOwn*(sqr(UvOwn + cm)/(4*cm) - 0.5*(UvOwn + mag(UvOwn))))
+       + 0.5*(UvOwn + mag(UvOwn))
+    );
+    surfaceScalarField uMinus(
+      neg0(mag(UvNei/cm) - 1)*(alphaNei*(-sqr(UvNei - cm)/(4*cm) - 0.5*(UvNei - mag(UvNei))))
+       + 0.5*(UvNei - mag(UvNei))
+    );
+
+    U12 = (alphaOwn*UOwn + alphaNei*UNei)/2.0;
+
+    surfaceScalarField pPlus(
+       neg0(mag(UvOwn/cm) - 1)*pOwn*sqr(UvOwn/cm + 1.0)*(2.0 - UvOwn/cm)/4.0
+       + Foam::pos(mag(UvOwn/cm) - 1)*pOwn*0.5*(1 + Foam::sign(UvOwn))//pos(mag(UvNei))
+    );
+    surfaceScalarField pMinus(
+       neg0(mag(UvNei/cm) - 1)*pNei*sqr(UvNei/cm - 1.0)*(2.0 + UvNei/cm)/4.0
+       + Foam::pos(mag(UvNei/cm) - 1)*pNei*0.5*(1 - Foam::sign(UvNei))//max(UvNei,minU)
+    );
+
+    surfaceScalarField P12(pPlus + pMinus);
+    surfaceScalarField s(0.5*min(norP , 10.0*mag(pNei - pOwn)/min(pOwn,pNei)));
+
+    surfaceScalarField caseA(Foam::neg(UvOwn - aOwn)*Foam::pos(UvNei - aNei));
+    surfaceScalarField caseB(Foam::neg(UvOwn + aOwn)*Foam::pos(UvNei + aNei));
+
+    massFlux = (uPlus*rhoOwn + uMinus*rhoNei)*magSf
+               -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei - rhoOwn)*magSf
+               + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei - rhoOwn)*magSf);
+
+    surfaceVectorField AUSMV((uPlus*rhoOwn*UOwn + uMinus*rhoNei*UNei)*magSf);
+    surfaceVectorField AUSMD(0.5*(massFlux*(UOwn+UNei) - mag(massFlux)*(UNei-UOwn)));
+
+    momentumFlux = (0.5+s)*AUSMV + (0.5-s)*AUSMD + P12*normal*magSf
+                  -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei*UNei - rhoOwn*UOwn)*magSf
+                  + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei*UNei - rhoOwn*UOwn)*magSf);
+
+    energyFlux = 0.5*(massFlux*(HOwn+HNei) - mag(massFlux)*(HNei-HOwn)) + meshPhi*P12
+                  -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei*HNei - rhoOwn*HOwn)*magSf
+                  + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei*HNei - rhoOwn*HOwn)*magSf);
+
+
+}
diff --git a/applications/solvers/dfHighSpeedFoam/fluxSchemes/KNPFlux.H b/applications/solvers/dfHighSpeedFoam/fluxSchemes/KNPFlux.H
new file mode 100644
index 00000000..d96b7739
--- /dev/null
+++ b/applications/solvers/dfHighSpeedFoam/fluxSchemes/KNPFlux.H
@@ -0,0 +1,70 @@
+void KNPFlux(
+    const surfaceScalarField& meshPhi,
+    const surfaceScalarField& magSf,
+    const surfaceVectorField& normal,
+    const surfaceScalarField& rhoOwn,
+    const surfaceScalarField& rhoNei,
+    const surfaceVectorField& UOwn,
+    const surfaceVectorField& UNei,
+    const surfaceScalarField& HOwn,
+    const surfaceScalarField& HNei,
+    const surfaceScalarField& pOwn,
+    const surfaceScalarField& pNei,
+    const surfaceScalarField& gammaOwn,
+    const surfaceScalarField& gammaNei,
+    const surfaceScalarField& aOwn,
+    const surfaceScalarField& aNei,
+    surfaceVectorField& U12,
+    surfaceScalarField& massFlux,
+    surfaceVectorField& momentumFlux,
+    surfaceScalarField& energyFlux
+)
+{
+    surfaceScalarField UvOwn((UOwn & normal) - meshPhi/magSf);
+
+    surfaceScalarField UvNei((UNei & normal) - meshPhi/magSf);
+
+    surfaceScalarField EOwn("EOwn", HOwn - pOwn/rhoOwn);
+    surfaceScalarField ENei("ENei", HNei - pNei/rhoNei);
+
+    surfaceScalarField phiv_pos("phiv_pos", UvOwn*magSf);  //
+    surfaceScalarField phiv_neg("phiv_neg", UvNei*magSf);
+
+
+    surfaceScalarField cSf_pos("cSf_pos",aOwn*magSf);
+    surfaceScalarField cSf_neg("cSf_neg",aNei*magSf);
+    dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
+    surfaceScalarField ap("ap",max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero));
+    surfaceScalarField am("am",min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero));
+
+    surfaceScalarField a_pos("a_pos", ap/(ap - am));
+
+    surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
+
+    surfaceScalarField aSf("aSf", am*a_pos);
+
+    surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
+
+    phiv_pos *= a_pos;
+    phiv_neg *= a_neg;
+
+    U12 = a_pos*UOwn + a_neg*UNei;
+
+    surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
+    surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
+
+    // Reuse amaxSf for the maximum positive and negative fluxes
+    // estimated by the central scheme
+    amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
+
+    massFlux = aphiv_pos*rhoOwn + aphiv_neg*rhoNei;
+
+    momentumFlux = (aphiv_pos*rhoOwn*UOwn + aphiv_neg*rhoNei*UNei)
+        + (a_pos*pOwn + a_neg*pNei)*magSf*normal;
+
+    energyFlux = aphiv_pos*(rhoOwn*EOwn + pOwn)
+        + aphiv_neg*(rhoNei*ENei + pNei)
+        + aSf*pOwn - aSf*pNei
+        + meshPhi*(a_pos*pOwn + a_neg*pNei);
+
+}
diff --git a/applications/solvers/dfHighSpeedFoam/phiCal.H b/applications/solvers/dfHighSpeedFoam/phiCal.H
index e5f569ed..d85b4f94 100644
--- a/applications/solvers/dfHighSpeedFoam/phiCal.H
+++ b/applications/solvers/dfHighSpeedFoam/phiCal.H
@@ -1,4 +1,82 @@
-phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+//phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+
+surfaceVectorField  phiUp
+(
+    IOobject
+    (
+        "phiUp",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", sqr(dimVelocity)*dimDensity*dimArea, Zero)
+);
+
+surfaceVectorField  u12
+(
+    IOobject
+    (
+        "u12",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", dimVelocity, Zero)
+);
+
+surfaceScalarField  phiEp
+(
+    IOobject
+    (
+        "phiEp",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", pow3(dimVelocity)*dimDensity*dimArea, 0)
+);
+
+
+surfaceScalarField  meshPhi
+(
+    IOobject
+    (
+        "meshPhi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", (dimVelocity)*dimArea, 0)
+);
+if(mesh.moving())
+{
+  meshPhi = meshPhi + mesh.phi();
+}
+
+//-- Computing flux
+fluxMap[fluxScheme](
+    meshPhi,
+    mesh.magSf(),
+    normal,
+    rho_pos,rho_neg,
+    U_pos,U_neg,
+    H_pos,H_neg,
+    p_pos,p_neg,
+    gamma_pos,gamma_neg,
+    a_pos,a_neg,
+    u12,
+    phi,
+    phiUp,
+    phiEp
+);
 
 PtrList<surfaceScalarField> phiYi(nspecies);
 forAll(phiYi,i)
@@ -16,30 +94,12 @@ forAll(phiYi,i)
         IOobject::NO_READ,
         IOobject::NO_WRITE
       ),
-      aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i]
+      phi*fvc::interpolate(Y[i], "Yi").ref() 
     )
   );
 }
 
-surfaceVectorField phiUp
-(
-  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
-    + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
-);
-
-surfaceScalarField phiEp
-(
-    "phiEp",
-    aphiv_pos*(rho_pos*(ea_pos + 0.5*magSqr(U_pos)) + p_pos)
-  + aphiv_neg*(rho_neg*(ea_neg + 0.5*magSqr(U_neg)) + p_neg)
-  + aSf*p_pos - aSf*p_neg
-);
 
-// Make flux for pressure-work absolute
-if (mesh.moving())
-{
-  phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
-}
 
 volScalarField muEff("muEff", turbulence->muEff());
 volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
diff --git a/applications/solvers/dfHighSpeedFoam/preCal.H b/applications/solvers/dfHighSpeedFoam/preCal.H
index 87917058..2d959e1b 100644
--- a/applications/solvers/dfHighSpeedFoam/preCal.H
+++ b/applications/solvers/dfHighSpeedFoam/preCal.H
@@ -1,9 +1,11 @@
+surfaceVectorField normal(mesh.Sf()/mesh.magSf());
 // --- Directed interpolation of primitive fields onto faces
 surfaceScalarField rho_pos(interpolate(rho, pos));
 surfaceScalarField rho_neg(interpolate(rho, neg));
 
 PtrList<surfaceScalarField> rhoYi_pos(nspecies);
 PtrList<surfaceScalarField> rhoYi_neg(nspecies);
+
 forAll(rhoYi_pos,i)
 {
   rhoYi_pos.set
@@ -47,6 +49,9 @@ forAll(rhoYi_neg,i)
 surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
 surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
 
+surfaceScalarField rhoE_pos(interpolate(rhoE, pos, T.name()));
+surfaceScalarField rhoE_neg(interpolate(rhoE, neg, T.name()));
+
 volScalarField rPsi("rPsi", 1.0/psi);
 surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
 surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
@@ -60,59 +65,42 @@ surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
 surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
 surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
 
-surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
-surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
 
-// Make fluxes relative to mesh-motion
-if (mesh.moving())
-{
-  phiv_pos -= mesh.phi();
-  phiv_neg -= mesh.phi();
-}
+surfaceScalarField E_pos("E_pos", ea_pos + 0.5*magSqr(U_pos));//rhoE_pos/rho_pos);
+surfaceScalarField E_neg("E_neg", ea_neg + 0.5*magSqr(U_neg));//rhoE_neg/rho_neg);
 
-volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
-surfaceScalarField cSf_pos
+surfaceScalarField H_pos("H_pos", (rhoE_pos + p_pos)/rho_pos);
+surfaceScalarField H_neg("H_neg", (rhoE_neg + p_neg)/rho_neg);
+
+volScalarField gamma_("gamma_", thermo.Cp()/thermo.Cv());
+surfaceScalarField gamma_pos
 (
-  "cSf_pos",
-  interpolate(c, pos, T.name())*mesh.magSf()
+  "gamma_pos",
+   interpolate(gamma_, pos, T.name())
 );
-surfaceScalarField cSf_neg
+surfaceScalarField gamma_neg
 (
-  "cSf_neg",
-  interpolate(c, neg, T.name())*mesh.magSf()
+  "gamma_neg",
+   interpolate(gamma_, neg, T.name())
 );
 
-surfaceScalarField ap
+volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
+surfaceScalarField a_pos
 (
-  "ap",
-  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
+  "a_pos",
+  interpolate(c, pos, T.name())
 );
-surfaceScalarField am
+surfaceScalarField a_neg
 (
-  "am",
-  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
+  "a_neg",
+  interpolate(c, neg, T.name())
 );
 
-surfaceScalarField a_pos("a_pos", ap/(ap - am));
-
-surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
-
-surfaceScalarField aSf("aSf", am*a_pos);
-
-if (fluxScheme == "Tadmor")
-{
-  aSf = -0.5*amaxSf;
-  a_pos = 0.5;
-}
-
-surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
-
-phiv_pos *= a_pos;
-phiv_neg *= a_neg;
-
-surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
-surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
+surfaceScalarField amaxSf("amaxSf",
+    0.5 * mag((U_pos + U_neg) & mesh.Sf()) 
+    + 0.5 * (a_pos + a_neg) * mesh.magSf()
+);
 
-// Reuse amaxSf for the maximum positive and negative fluxes
-// estimated by the central scheme
-amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
\ No newline at end of file
+if(mesh.moving()){
+  amaxSf -= mesh.phi();
+}
\ No newline at end of file
diff --git a/applications/solvers/dfHighSpeedFoam/readFluxScheme.H b/applications/solvers/dfHighSpeedFoam/readFluxScheme.H
index e8c28d65..ac469e94 100644
--- a/applications/solvers/dfHighSpeedFoam/readFluxScheme.H
+++ b/applications/solvers/dfHighSpeedFoam/readFluxScheme.H
@@ -1,16 +1,44 @@
+
+typedef void (*FnPtr)
+(
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceVectorField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceVectorField&,
+    const surfaceVectorField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    const surfaceScalarField&,
+    surfaceVectorField&,
+    surfaceScalarField&,
+    surfaceVectorField&,
+    surfaceScalarField&
+);
+
+std::map<word, FnPtr> fluxMap;
+fluxMap["AUSMDV"] = AUSMDVFlux;
+fluxMap["Kurganov"] = KNPFlux;
+//fluxMap["RoeC"] = RoeCFlux;
+
 word fluxScheme("Kurganov");
-if (mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme))
+mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme);
+if ((fluxScheme == "AUSMDV") || (fluxScheme == "Kurganov"))
 {
-    if ((fluxScheme == "Tadmor") || (fluxScheme == "Kurganov"))
-    {
-        Info<< "fluxScheme: " << fluxScheme << endl;
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "fluxScheme: " << fluxScheme
-            << " is not a valid choice. "
-            << "Options are: Tadmor, Kurganov"
-            << abort(FatalError);
-    }
+    Info<< "fluxScheme: " << fluxScheme << endl;
 }
+else
+{
+    FatalErrorInFunction
+        << "fluxScheme: " << fluxScheme
+        << " is not a valid choice. "
+        << "Options are: Tadmor, Kurganov"
+        << abort(FatalError);
+}
+
diff --git a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
index 3a76804d..144daa0e 100644
--- a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
+++ b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
@@ -5,7 +5,7 @@ surfaceScalarField sigmaDotU
         fvc::interpolate(turbulence->muEff())*mesh.magSf()*fvc::snGrad(U)
       + (mesh.Sf() & fvc::interpolate(tauMC))
     )
-    & (a_pos*U_pos + a_neg*U_neg)
+    & u12
 );
 
 if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP"))
@@ -30,6 +30,7 @@ ea = rhoE/rho - 0.5*magSqr(U);
 ea.correctBoundaryConditions();
 
 ha = ea + p/rho;
+Info << "begin correct thermo" << endl;
 chemistry->correctThermo(); // before this, we must update ha = e + p/rho
 
 rhoE.boundaryFieldRef() == rho.boundaryField()*(ea.boundaryField() + 0.5*magSqr(U.boundaryField()));