diff --git a/applications/solvers/dfHighSpeedFoam/Make/options b/applications/solvers/dfHighSpeedFoam/Make/options index f5f42f89..e8d6d66c 100644 --- a/applications/solvers/dfHighSpeedFoam/Make/options +++ b/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 \ diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H index 563eef8f..5dded90b 100644 --- a/applications/solvers/dfHighSpeedFoam/createFields.H +++ b/applications/solvers/dfHighSpeedFoam/createFields.H @@ -126,4 +126,45 @@ forAll(Y, i) { fields.add(Y[i]); } -fields.add(thermo.he()); \ No newline at end of file +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) +); diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C index 8654cb95..042d1a6a 100644 --- a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C +++ b/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; diff --git a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H index 23eaf01b..c1ba18d6 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H +++ b/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; diff --git a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H index 84221586..f0a50aec 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H @@ -1,3 +1,7 @@ +hDiffCorrFlux = Zero; +diffAlphaD = Zero; +sumYDiffError = Zero; + tmp> mvConvection ( fv::convectionScheme::New @@ -9,39 +13,50 @@ tmp> 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 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; } @@ -49,4 +64,7 @@ tmp> mvConvection Y[inertIndex] = scalar(1) - Yt; Y[inertIndex].max(0.0); + + end = std::clock(); + time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC); } diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/T.orig.gz b/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/T.gz similarity index 100% rename from examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/T.orig.gz rename to examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/T.gz diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/U.orig.gz b/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/U.gz similarity index 100% rename from examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/U.orig.gz rename to examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/U.gz diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/p.orig.gz b/examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/p.gz similarity index 100% rename from examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/p.orig.gz rename to examples/dfHighSpeedFoam/det_1D_H2_01mm/0_orig/p.gz diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/constant/CanteraTorchProperties b/examples/dfHighSpeedFoam/det_1D_H2_01mm/constant/CanteraTorchProperties index d4b25919..78517df2 100755 --- a/examples/dfHighSpeedFoam/det_1D_H2_01mm/constant/CanteraTorchProperties +++ b/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"; diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/blockMeshDict b/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/blockMeshDict index 8199aa9c..1a73e764 100755 --- a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/blockMeshDict +++ b/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 ( diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/controlDict b/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/controlDict index 7bf316e4..57ed5645 100755 --- a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/controlDict +++ b/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/controlDict @@ -38,7 +38,7 @@ writeFormat ascii; writePrecision 8; -writeCompression off; +writeCompression on; timeFormat general; diff --git a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/fvSchemes b/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/fvSchemes index f4f06e9a..3f654d84 100755 --- a/examples/dfHighSpeedFoam/det_1D_H2_01mm/system/fvSchemes +++ b/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; }