add diffusion velocity correction and timing function for dfHighSpeedFoam, correct alpha in rhoEEqn.H and clean the code

applications/solvers/dfHighSpeedFoam/Make/options

@@ -1,4 +1,7 @@
+-include $(GENERAL_RULES)/mplibType
+
 EXE_INC = -std=c++14 \
+    $(PFLAGS) $(PINC) \
     -I$(FOAM_APP)/solvers/compressible/rhoCentralFoam/BCs/lnInclude \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -8,6 +11,7 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(DF_SRC)/CanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel \

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -126,4 +126,45 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(thermo.he());
+fields.add(thermo.he());
+
+const scalar Sct = chemistry.lookupOrDefault("Sct", 1.);
+volScalarField diffAlphaD
+(
+    IOobject
+    (
+        "diffAlphaD",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimEnergy/dimTime/dimVolume, 0)
+);
+volVectorField hDiffCorrFlux
+(
+    IOobject
+    (
+        "hDiffCorrFlux",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector(dimensionSet(1,0,-3,0,0,0,0), Zero)
+);
+volVectorField sumYDiffError
+(
+    IOobject
+    (
+        "sumYDiffError",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero)
+);

applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C

@@ -42,8 +42,7 @@ Description
 #include "directionInterpolate.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
-//add
-//#include "fvOptions.H"
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +56,11 @@ int main(int argc, char *argv[])
     #include "createDynamicFvMesh.H"
     #include "createFields.H"
     #include "createTimeControls.H"
-    //#include "createFvOptions.H"
+
+    double time_monitor_flow;
+    double time_monitor_chem;
+    double time_monitor_Y;
+    clock_t start, end;
 
     turbulence->validate();
 
@@ -200,11 +203,14 @@ int main(int argc, char *argv[])
         volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
 
         // --- Solve density
-        //solve(fvm::ddt(rho) + fvc::div(phi));
         #include "rhoEqn.H"
 
+        start = std::clock();
         // --- Solve momentum
         #include "rhoUEqn.H"
+        end = std::clock();
+        time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC);
+
 
         // --- Solve species
         #include "rhoYEqn.H"
@@ -216,6 +222,10 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
+        Info<< "MonitorTime_chem = " << time_monitor_chem << " s" << nl << endl;
+        Info<< "MonitorTime_Y = " << time_monitor_Y << " s" << nl << endl;
+        Info<< "MonitorTime_flow = " << time_monitor_flow << " s" << nl << endl;
+
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;

applications/solvers/dfHighSpeedFoam/rhoEEqn.H

@@ -10,9 +10,9 @@ surfaceScalarField sigmaDotU
 
 solve
 (
-    fvm::ddt(rhoE)
-  + fvc::div(phiEp)
-  - fvc::div(sigmaDotU)
+      fvm::ddt(rhoE)
+    + fvc::div(phiEp)
+    - fvc::div(sigmaDotU)
 );
 
 e = rhoE/rho - 0.5*magSqr(U);
@@ -27,9 +27,14 @@ if (!inviscid)
 {
     fvScalarMatrix eEqn
     (
-        fvm::ddt(rho, e) - fvc::ddt(rho, e)
-      - fvc::laplacian(turbulence->alphaEff(), e)
+           fvm::ddt(rho, e) - fvc::ddt(rho, e)
+        // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha
+         - fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), e)
+         + diffAlphaD
+      ==
+           fvc::div(hDiffCorrFlux)
     );
+
     eEqn.solve("e");
 
     ha = e + p/rho;

applications/solvers/dfHighSpeedFoam/rhoYEqn.H

@@ -1,3 +1,7 @@
+hDiffCorrFlux = Zero;
+diffAlphaD = Zero;
+sumYDiffError = Zero;
+
 tmp<fv::convectionScheme<scalar>> mvConvection
 (
     fv::convectionScheme<scalar>::New
@@ -9,44 +13,58 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 
+forAll(Y, i)
 {
+    sumYDiffError += chemistry.rhoD(i)*fvc::grad(Y[i]);
+}
+const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
+
+{
+    start = std::clock();
     chemistry.solve(mesh.time().deltaTValue());
+    label flag_mpi_init;
+    MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+    end = std::clock();
+    time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC);
 
     volScalarField Yt(0*Y[0]);
 
+    start = std::clock();
     forAll(Y, i)
     {
+        volScalarField& Yi = Y[i];
+        hDiffCorrFlux += chemistry.hai(i)*(chemistry.rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError);
+        diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry.hai(i), Yi);
+
         if (i != inertIndex)
         {
-            volScalarField& Yi = Y[i];
-
+            tmp<volScalarField> DEff = chemistry.rhoD(i) + turbulence->mut()/Sct;
             fvScalarMatrix YiEqn
             (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
+                 fvm::ddt(rho, Yi)
+               + mvConvection->fvmDiv(phi, Yi)
               ==
-                chemistry.RR(i)
-              //+ fvOptions(rho, Yi)
+                 chemistry.RR(i)
             );
 
             if (!inviscid)
             {
-                YiEqn -= fvm::laplacian(turbulence->muEff(), Yi);
+                YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi);
             }
 
             YiEqn.relax();
 
-            //fvOptions.constrain(YiEqn);
-
             YiEqn.solve("Yi");
 
-            //fvOptions.correct(Yi);
-
             Yi.max(0.0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
+
+    end = std::clock();
+    time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC);
 }

examples/dfHighSpeedFoam/det_1D_H2_01mm/constant/CanteraTorchProperties

@@ -19,6 +19,11 @@ FoamFile
 chemistry           on;
 CanteraMechanismFile "H2_Li.xml";
 transportModel "UnityLewis";
+loadbalancing
+{
+        active  false;
+        //log   true;
+}
 odeCoeffs
 {}
 inertSpecie        "N2";

examples/dfHighSpeedFoam/det_1D_H2_01mm/system/blockMeshDict

@@ -41,7 +41,6 @@ boundary
 (
     Left
     {
-        // type patch;
         type symmetryPlane;
         faces
         (
@@ -50,7 +49,6 @@ boundary
     }
     Right
     {
-        // type patch;
         type symmetryPlane;
         faces
         (

examples/dfHighSpeedFoam/det_1D_H2_01mm/system/controlDict

@@ -38,7 +38,7 @@ writeFormat     ascii;
 
 writePrecision  8;
 
-writeCompression off;
+writeCompression on;
 
 timeFormat      general;
 

examples/dfHighSpeedFoam/det_1D_H2_01mm/system/fvSchemes

@@ -29,33 +29,29 @@ gradSchemes
 
 divSchemes
 {
-    	default         none;
-		div(tauMC)		Gauss linear;
-		div(phi,O2)		Gauss vanLeer01;
-		div(phi,H2)		Gauss vanLeer01;
-		div(phi,H2O)	Gauss vanLeer01;
-		div(phi,AR)		Gauss vanLeer01;
-		div(phi,Yi)		Gauss vanLeer01;
-		div(phi,Yi_h)	Gauss vanLeer;
+    default         		none;
+	div(tauMC)				Gauss linear;
+	div(hDiffCorrFlux)    	Gauss cubic;
+	div(phi,Yi_h)			Gauss vanLeer;
 }
 
 laplacianSchemes
 {
-		default			Gauss linear uncorrected;	
+	default			Gauss linear uncorrected;	
 }
 
 interpolationSchemes
 {
-    	default				linear;
-		reconstruct(rho)	Minmod;
-		reconstruct(U)		MinmodV;
-		reconstruct(T)		Minmod;
-		reconstruct(Yi)		Minmod;
+    default				linear;
+	reconstruct(rho)	Minmod;
+	reconstruct(U)		MinmodV;
+	reconstruct(T)		Minmod;
+	reconstruct(Yi)		Minmod;
 }
 
 snGradSchemes
 {
-		default				uncorrected;
+	default				uncorrected;
 }