realize RK2SSP and RK3SSP in dfHighSpeedFoam and add wrate option for chemistry calculation

applications/solvers/dfHighSpeedFoam/calculateR.H

@@ -0,0 +1,49 @@
+start = std::clock();
+
+combustion->correct();
+
+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.0*Y[0]);
+
+start = std::clock();
+forAll(Y, i)
+{
+    volScalarField& Yi = Y[i];
+
+    if (i != inertIndex)
+    {          
+        rhoYi[i].ref() += chemistry->RR(i)*runTime.deltaT();
+        Info <<"max reaction rate "<< Yi.name() << " is " << max(chemistry->RR(i)).value() << endl;
+
+        Yi = rhoYi[i]/rho;
+        Yi.max(0.0);
+
+        Yi.correctBoundaryConditions();
+        rhoYi[i] = rho*Yi;
+        Yt += Yi;
+    }
+}
+
+Y[inertIndex] = scalar(1) - Yt;
+Y[inertIndex].max(0.0);
+rhoYi[inertIndex] = rho*Y[inertIndex];
+
+end = std::clock();
+time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC);
+
+chemistry->correctThermo();
+Info<< "min/max(T) = "
+        << min(T).value() << ", " << max(T).value() << endl;
+
+p.ref() = rho()/psi();
+p.correctBoundaryConditions();
+rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
+
+
+Info << " finish calculate Reaction" << nl << endl;
+

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -1,5 +1,44 @@
 #include "createRDeltaT.H"
 
+word ddtSchemes("Euler");
+if (mesh.schemesDict().readIfPresent("timeScheme", ddtSchemes))
+{
+    if(ddtSchemes == "RK2SSP" || ddtSchemes == "RK3SSP")
+    {
+        Info<< "ddtSchemes: " << ddtSchemes << endl;
+    }
+}
+
+word chemScheme("RR");
+mesh.schemesDict().readIfPresent("chemScheme", chemScheme);
+
+if ((chemScheme == "wrate") || (chemScheme == "RR"))
+{
+    Info<< "chemScheme: " << chemScheme << endl;
+}
+else
+{
+    FatalErrorInFunction
+        << "chemScheme: " << chemScheme
+        << " is not a valid choice. "
+        << "Options are: wrate, RR"
+        << abort(FatalError);
+}
+
+
+if((ddtSchemes != "RK2SSP") && (ddtSchemes != "RK3SSP"))
+{
+    if(chemScheme == "wrate")
+    {
+        FatalErrorInFunction
+            << "This combination is not a valid choice. "
+            << "If you want to use wrate for chemistry, please use RK2SSP or RK3SSP scheme."
+            << abort(FatalError);
+    }
+}
+
+
+
 Info<< "Reading thermophysical properties\n" << endl;
 
 // psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);

applications/solvers/dfHighSpeedFoam/createFields_rk2.H

@@ -0,0 +1,27 @@
+volScalarField rho_save("rho_save",thermo.rho()); 
+
+volVectorField rhoU_save("rhoU_save",rho*U);
+
+volScalarField rhoE_save("rhoE_save",rho*(ea + 0.5*magSqr(U))); 
+
+PtrList<volScalarField> rhoYi_save(nspecies);
+forAll(rhoYi_save,i)
+{
+    rhoYi_save.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoYi_save" + Y[i].name(),
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            rho*Y[i]
+        )
+    );
+}
+