diff --git a/applications/solvers/df0DFoam/createFields.H b/applications/solvers/df0DFoam/createFields.H index cfe0c17f..5ac2944f 100644 --- a/applications/solvers/df0DFoam/createFields.H +++ b/applications/solvers/df0DFoam/createFields.H @@ -2,7 +2,7 @@ Info<< "Reading thermophysical properties\n" << endl; -CanteraMixture::setEnergyName("ha"); +CanteraMixture::setEnergyName("hs"); // fluidThermo* pThermo = new hePsiThermo(mesh, word::null); fluidThermo* pThermo = new heRhoThermo(mesh, word::null); fluidThermo& thermo = *pThermo; diff --git a/applications/solvers/dfBuoyancyFoam/createFields.H b/applications/solvers/dfBuoyancyFoam/createFields.H index bb9ba1e0..8f124121 100644 --- a/applications/solvers/dfBuoyancyFoam/createFields.H +++ b/applications/solvers/dfBuoyancyFoam/createFields.H @@ -3,7 +3,7 @@ Info<< "Reading thermophysical properties\n" << endl; //psiReactionThermo& thermo = pThermo(); //thermo.validate(args.executable(), "h", "e"); -CanteraMixture::setEnergyName("ha"); +CanteraMixture::setEnergyName("hs"); // fluidThermo* pThermo = new heRhoThermo(mesh, word::null); fluidThermo* pThermo = new heRhoThermo(mesh, word::null); fluidThermo& thermo = *pThermo; diff --git a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H index 4e79db11..d0170daf 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H @@ -5,7 +5,7 @@ if (!inviscid) hDiffCorrFlux = Zero; diffAlphaD = Zero; sumYDiffError = Zero; - + forAll(Y, i) { sumYDiffError += chemistry->rhoD(i)*fvc::grad(Y[i]); @@ -45,7 +45,7 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); 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); + time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC); volScalarField Yt(0.0*Y[0]); @@ -57,12 +57,12 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); if (!inviscid) { - hDiffCorrFlux += chemistry->hai(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError); - diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hai(i), Yi); + hDiffCorrFlux += chemistry->hei(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError); + diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hei(i), Yi); } // if (i != inertIndex) - { + { if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) { rhoYi_rhs = -fvc::div(rhoPhiYi[i]); @@ -72,12 +72,12 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); rhoYi_rhs.ref() += chemistry->wrate(i); Info <<"max reaction rate "<< Yi.name() << " is " << max(chemistry->wrate(i)).value() << endl; } - + rhoYi[i] = rkcoe1[nrk]*rhoYi_save[i] + rkcoe2[nrk]*rhoYi[i] + rkcoe3[nrk]*rhoYi_rhs*runTime.deltaT(); - + Yi=rhoYi[i]/rho; Yi.max(0.0); @@ -86,10 +86,10 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); tmp DEff = chemistry->rhoD(i) + turbulence->mut()/Sct; - rhoYi_rhs = + rhoYi_rhs = ( turbName == "laminar" - ? (fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phi,Yi_h)")) + ? (fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phi,Yi_h)")) : fvc::laplacian(DEff(), Yi) ); @@ -103,7 +103,7 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); Yi.correctBoundaryConditions(); rhoYi[i] = rho*Yi; - } + } else { // original convection term @@ -115,12 +115,12 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); // combustion->R(Yi) // + parcels.SYi(i, Yi) // ); - + solve ( fvm::ddt(rhoYi[i]) + fvc::div(rhoPhiYi[i]) - == + == chemistry->RR(i) + parcels.SYi(i, rhoYi[i]) ); @@ -134,11 +134,11 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); tmp DEff = chemistry->rhoD(i) + turbulence->mut()/Sct; // original term - // YiEqn -= + // YiEqn -= // ( // turbName == "laminar" // ? (fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi)) - // : fvm::laplacian(DEff(), Yi) + // : fvm::laplacian(DEff(), Yi) // ); @@ -165,7 +165,7 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC); // Yi.max(0.0); rhoYi[i] = rho*Yi; - } + } Yt += Yi; } } diff --git a/applications/solvers/dfLowMachFoam/YEqn.H b/applications/solvers/dfLowMachFoam/YEqn.H index 6203bbef..d5c6322a 100644 --- a/applications/solvers/dfLowMachFoam/YEqn.H +++ b/applications/solvers/dfLowMachFoam/YEqn.H @@ -29,7 +29,7 @@ const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); start1 = std::clock(); time_monitor_YEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); - time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC); + time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC); @@ -41,7 +41,7 @@ if (!splitting) { std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - + #ifdef ODE_GPU_SOLVER scalar dt = runTime.deltaTValue(); @@ -88,8 +88,8 @@ 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); + hDiffCorrFlux += chemistry->hei(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError); + diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hei(i), Yi); if (i != inertIndex) { start1 = std::clock(); @@ -111,7 +111,7 @@ : (fvm::laplacian(DEff(), Yi) + parcels.SYi(i, Yi) + combustion->R(Yi)) ) ); - + end1 = std::clock(); time_monitor_YEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); // YiEqn.relax(); diff --git a/applications/solvers/dfLowMachFoam/YEqn_GPU.H b/applications/solvers/dfLowMachFoam/YEqn_GPU.H index 85d1bb77..b9a5fe99 100644 --- a/applications/solvers/dfLowMachFoam/YEqn_GPU.H +++ b/applications/solvers/dfLowMachFoam/YEqn_GPU.H @@ -42,8 +42,8 @@ if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); 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); + hDiffCorrFlux += chemistry->hei(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError); + diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hei(i), Yi); if (i != inertIndex) { start1 = std::clock(); @@ -84,7 +84,7 @@ if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); ) ); #endif - + end1 = std::clock(); time_monitor_YEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); // YiEqn.relax(); @@ -243,5 +243,5 @@ if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); DEBUG_TRACE; fflush(stderr); - + } \ No newline at end of file diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H index a1229a89..32c8cd2c 100644 --- a/applications/solvers/dfLowMachFoam/createFields.H +++ b/applications/solvers/dfLowMachFoam/createFields.H @@ -3,7 +3,7 @@ Info<< "Reading thermophysical properties\n" << endl; -CanteraMixture::setEnergyName("ha"); +CanteraMixture::setEnergyName("hs"); // fluidThermo* pThermo = new hePsiThermo(mesh, word::null); fluidThermo* pThermo = new heRhoThermo(mesh, word::null); fluidThermo& thermo = *pThermo; diff --git a/applications/solvers/dfSprayFoam/YEqn.H b/applications/solvers/dfSprayFoam/YEqn.H index 2302b0ba..74436b4a 100644 --- a/applications/solvers/dfSprayFoam/YEqn.H +++ b/applications/solvers/dfSprayFoam/YEqn.H @@ -26,13 +26,13 @@ const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); 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); + hDiffCorrFlux += chemistry->hei(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError); + diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hei(i), Yi); if (i != inertIndex) { tmp DEff = chemistry->rhoD(i) + turbulence->mut()/Sct; - + fvScalarMatrix YEqn ( turbName == "laminar" @@ -69,7 +69,7 @@ const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); Yi.max(0.0); Yt += Yi; - + } } diff --git a/applications/solvers/dfSprayFoam/createFields.H b/applications/solvers/dfSprayFoam/createFields.H index 79c9171d..20295f4c 100644 --- a/applications/solvers/dfSprayFoam/createFields.H +++ b/applications/solvers/dfSprayFoam/createFields.H @@ -2,7 +2,7 @@ Info<< "Reading thermophysical properties\n" << endl; -CanteraMixture::setEnergyName("ha"); +CanteraMixture::setEnergyName("hs"); // fluidThermo* pThermo = new hePsiThermo(mesh, word::null); fluidThermo* pThermo = new heRhoThermo(mesh, word::null); fluidThermo& thermo = *pThermo; diff --git a/src/design.UML b/src/design.UML index 4f34f0ce..3b132301 100644 --- a/src/design.UML +++ b/src/design.UML @@ -204,7 +204,7 @@ class dfChemistryModel{ -scalar absTol_ -PtrList& Y_ -PtrList rhoD_; - -PtrList hai_; + -PtrList hei_; -mutable scalarList yTemp_ -mutable scalarList dTemp_; -mutable scalarList hrtTemp_; @@ -238,7 +238,7 @@ class dfChemistryModel{ +PtrList& Y() +const hashedWordList& species() const +const volScalarField& rhoD(const label i) const - +volScalarField& hai(const label i) + +volScalarField& hei(const label i) +void correctThermo() } @enduml \ No newline at end of file diff --git a/src/design.svg b/src/design.svg index 43a370f8..6de7dcdc 100644 --- a/src/design.svg +++ b/src/design.svg @@ -1,31 +1,31 @@ -CanteraSolutionThermoPhaseTransportIOdictionarybasicThermoconst word& phaseName_volScalarField& p_volScalarField T_volScalarField alpha_Switch dpdt_Protected datavolScalarField& lookupOrConstruct(const fvMesh& mesh,const char* name) constwordList heBoundaryTypes()wordList heBoundaryBaseTypes()Protected Member Functionsvoid correct() = 0Switch dpdt() const{return dpdt_}Member functionsvolScalarField& p() // [Pa]const volScalarField& p() const // [Pa]tmp<volScalarField> rho() const = 0 // [kg/m^3]tmp<scalarField> rho(patchi) const = 0 // [kg/m^3]volScalarField& he() = 0 // [J/kg]const volScalarField& he() const = 0 // [J/kg]tmp<volScalarField> he(p,T) const = 0 // [J/kg]tmp<scalarField> he(p,T,const labelList& cells) const = 0 // [J/kg]tmp<scalarField> he(p,T,patchi) const = 0 // [J/kg]tmp<volScalarField> hc() const = 0 // [J/kg]tmp<scalarField> THE(h,p,T0,const labelList& cells) const = 0tmp<scalarField> THE(h,p,T0,patchi) const = 0Access to thermodynamic state variablesconst volScalarField& T() const // [K]volScalarField& T() // [K]tmp<volScalarField> Cp() const = 0 // [J/kg/K]tmp<scalarField> Cp(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> Cv() const = 0 // [J/kg/K]tmp<scalarField> Cv(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> gamma() const = 0 // []tmp<scalarField> gamma(p,T,patchi) const = 0 // []tmp<volScalarField> Cpv() const = 0 // [J/kg/K]tmp<scalarField> Cpv(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> CpByCpv() const = 0 // []tmp<scalarField> CpByCpv(p,T,patchi) const = 0 // []tmp<volScalarField> W() const = 0 // [kg/kmol]tmp<scalarField> W(patchi) const = 0 // [kg/kmol]Access to thermodynamic state variablesconst volScalarField& alpha() const // [kg/m/s]const scalarField& alpha(patchi) const // [kg/m/s]tmp<volScalarField> kappa() const = 0 // [W/m/K]tmp<scalarField> kappa(patchi) const = 0 // [W/m/K]tmp<volScalarField> alphahe() const = 0 // [kg/m/s]tmp<scalarField> alphahe(const label patchi) const = 0 // [kg/m/s]tmp<volScalarField> kappaEff(const volScalarField&) const = 0 // [W/m/K]tmp<scalarField> kappaEff(alphat,patchi) const = 0 // [W/m/K]tmp<volScalarField> alphaEff(alphat) const = 0 // [kg/m/s]tmp<scalarField> alphaEff(alphat,patchi) const = 0 // [kg/m/s]Access to transport state variablesfluidThermovoid correctRho(const volScalarField& deltaRho) = 0const volScalarField& psi() const = 0 // [s^2/m^2]tmp<volScalarField> mu() const = 0 // [kg/m/s]tmp<scalarField> mu(patchi) const = 0 // [kg/m/s]tmp<volScalarField> nu() const // [m^2/s]tmp<scalarField> nu(patchi) const // [m^2/s]Member functionscompressibleTransportModelpsiThermovolScalarField psi_volScalarField mu_virtual tmp<volScalarField> rho() const {return p_*psi_}virtual tmp<scalarField> rho(patchi) constvirtual void correctRho(const volScalarField& deltaRho) {}virtual const volScalarField& psi() constvirtual tmp<volScalarField> mu() constvirtual tmp<scalarField> mu(patchi) constProtected dataCanteraMixtureIOdictionary CanteraTorchProperties_const word CanteraMechanismFile_std::shared_ptr<Cantera::Solution> CanteraSolution_std::shared_ptr<Cantera::ThermoPhase> CanteraGas_word transportModelName_std::shared_ptr<Cantera::Transport> CanteraTransport_hashedWordList species_PtrList<volScalarField> Y_const volScalarField& Tref_const volScalarField& pref_mutable scalarList yTemp_typedef CanteraMixture basicMixtureTypetypedef CanteraMixture thermoTypestatic const bool incompressible = falsestatic const bool isochoric = falsescalarList HaTemp_scalarList CpTemp_scalarList muTemp_CanteraMixture& cellMixture(celli)CanteraMixture& patchFaceMixture(patchi, facei)const scalar THE(scalar p, scalar T)const scalar psi(scalar p, scalar T)const scalar mu(scalar p, scalar T)const scalar alphah(scalar p, scalar T)const scalar HE(scalar p, scalar T)const scalar Hc(scalar p, scalar T)const scalar Cp(scalar p, scalar T)const scalar Cv(scalar p, scalar T)const scalar gamma(scalar p, scalar T)const scalar Cpv(scalar p, scalar T)const scalar CpByCpv(scalar p, scalar T)const scalar W() conststatic inline word heName() {return "ha"}PtrList<volScalarField>& Y()const PtrList<volScalarField>& Y() constvolScalarField& Y(const label i)const volScalarField& Y(const label i) constvolScalarField& Y(const word& specieName)const volScalarField& Y(const word& specieName) constconst hashedWordList& species()scalar nSpecies()std::shared_ptr<Cantera::ThermoPhase> CanteraGas()std::shared_ptr<Cantera::Solution> CanteraSolution()std::shared_ptr<Cantera::Transport> CanteraTransport()const word& transportModelName()void calcCp(const scalar T, const scalar p)void calcMu(const scalar T, const scalar p)void calcH(const scalar T, const scalar p)scalar Cp(label i, scalar p, scalar T) constscalar mu(label i, scalar p, scalar T) constscalar Ha(label i, scalar p, scalar T) constscalar Hc(label i) constscalar Hs(label i, scalar p, scalar T) constscalar kappa(label i, scalar p, scalar T) constscalar Wi(label i) constheThermoBasicThermo, MixtureTypevolScalarField he_Protected datavirtual typename MixtureType::basicMixtureType&composition(){return *this}virtual const typename MixtureType::basicMixtureType&composition() const{return *this}virtual word thermoName() const{return MixtureType::thermoType::typeName()}virtual volScalarField& he(){return he_} // [J/kg]virtual const volScalarField& he() const{return he_} // [J/kg]virtual tmp<volScalarField> he(p,T) const // [J/kg]virtual tmp<scalarField> he(p,T,const labelList& cells) const // [J/kg]virtual tmp<scalarField> he(p,T,patchi) const // [J/kg]virtual tmp<volScalarField> hc() const // [J/kg]virtual tmp<scalarField> THE(he, p, T0, const labelList& cells) constvirtual tmp<scalarField> THE(he, p, T0, patchi) constvirtual tmp<scalarField> Cp(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cp() const // [J/kg/K]virtual tmp<scalarField> Cv(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cv() const // [J/kg/K]virtual tmp<volScalarField> gamma() const // []virtual tmp<scalarField> gamma(p, T, patchi) const // [virtual tmp<scalarField> Cpv(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cpv() const // [J/kg/K]virtual tmp<volScalarField> CpByCpv() const // []virtual tmp<scalarField> CpByCpv(p, T,patchi) const // []virtual tmp<volScalarField> W() const // [kg/kmol]virtual tmp<scalarField> W(const label patchi) const // [kg/kmol]virtual tmp<volScalarField> kappa() const // [W/m/K]virtual tmp<scalarField> kappa(patchi) const // [W/m/K]virtual tmp<volScalarField> alphahe() const // [kg/m/s]virtual tmp<scalarField> alphahe(patchi) const // [kg/m/s]virtual tmp<volScalarField> kappaEff(const volScalarField&) const // [W/m/K]virtual tmp<scalarField> kappaEff(alphat, patchi) const // [W/m/K]virtual tmp<volScalarField> alphaEff(alphat) const // [kg/m/s]virtual tmp<scalarField> alphaEff(alphat, patchi) const // [kg/m/s]Member functionshePsiThermoBasicPsiThermo, MixtureTypecalculate() : voidcorrect() : voiddfChemistryModelThermoTypeThermoType& thermo_;CanteraMixture& mixture_std::shared_ptr<Cantera::ThermoPhase> CanteraGas_const fvMesh& mesh_Switch chemistry_scalar relTol_scalar absTol_PtrList<volScalarField>& Y_PtrList<volScalarField> rhoD_;PtrList<volScalarField> hai_;mutable scalarList yTemp_mutable scalarList dTemp_;mutable scalarList hrtTemp_;mutable scalarField cTemp_;PtrList<volScalarField::Internal> RR_hashedWordList species_volScalarField& alpha_volScalarField& T_const volScalarField& p_const volScalarField& rho_;volScalarField& mu_volScalarField& psi_volScalarField Qdot_;Switch torchSwitch_;word torchModelName_;scalarList Xmu_;scalarList Xstd_;scalarList Ymu_;scalarList Ystd_;scalar Tact_;scalar Qdotact_;void setNumerics(Cantera::ReactorNet &r)scalar canteraSolve(const DeltaTType& deltaT)scalar torchSolve(const DeltaTType& deltaT)virtual scalar solve(const scalar deltaT)virtual scalar solve(const scalarField& deltaT)virtual const volScalarField::Internal& RR(const label i) constvirtual volScalarField::Internal& RR(const label i)virtual tmp<volScalarField> Qdot() constPtrList<volScalarField>& Y()const hashedWordList& species() constconst volScalarField& rhoD(const label i) constvolScalarField& hai(const label i)void correctThermo()MixtureTypeBasicThermoThermoTypeCanteraSolutionThermoPhaseTransportIOdictionarybasicThermoconst word& phaseName_volScalarField& p_volScalarField T_volScalarField alpha_Switch dpdt_Protected datavolScalarField& lookupOrConstruct(const fvMesh& mesh,const char* name) constwordList heBoundaryTypes()wordList heBoundaryBaseTypes()Protected Member Functionsvoid correct() = 0Switch dpdt() const{return dpdt_}Member functionsvolScalarField& p() // [Pa]const volScalarField& p() const // [Pa]tmp<volScalarField> rho() const = 0 // [kg/m^3]tmp<scalarField> rho(patchi) const = 0 // [kg/m^3]volScalarField& he() = 0 // [J/kg]const volScalarField& he() const = 0 // [J/kg]tmp<volScalarField> he(p,T) const = 0 // [J/kg]tmp<scalarField> he(p,T,const labelList& cells) const = 0 // [J/kg]tmp<scalarField> he(p,T,patchi) const = 0 // [J/kg]tmp<volScalarField> hc() const = 0 // [J/kg]tmp<scalarField> THE(h,p,T0,const labelList& cells) const = 0tmp<scalarField> THE(h,p,T0,patchi) const = 0Access to thermodynamic state variablesconst volScalarField& T() const // [K]volScalarField& T() // [K]tmp<volScalarField> Cp() const = 0 // [J/kg/K]tmp<scalarField> Cp(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> Cv() const = 0 // [J/kg/K]tmp<scalarField> Cv(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> gamma() const = 0 // []tmp<scalarField> gamma(p,T,patchi) const = 0 // []tmp<volScalarField> Cpv() const = 0 // [J/kg/K]tmp<scalarField> Cpv(p,T,patchi) const = 0 // [J/kg/K]tmp<volScalarField> CpByCpv() const = 0 // []tmp<scalarField> CpByCpv(p,T,patchi) const = 0 // []tmp<volScalarField> W() const = 0 // [kg/kmol]tmp<scalarField> W(patchi) const = 0 // [kg/kmol]Access to thermodynamic state variablesconst volScalarField& alpha() const // [kg/m/s]const scalarField& alpha(patchi) const // [kg/m/s]tmp<volScalarField> kappa() const = 0 // [W/m/K]tmp<scalarField> kappa(patchi) const = 0 // [W/m/K]tmp<volScalarField> alphahe() const = 0 // [kg/m/s]tmp<scalarField> alphahe(const label patchi) const = 0 // [kg/m/s]tmp<volScalarField> kappaEff(const volScalarField&) const = 0 // [W/m/K]tmp<scalarField> kappaEff(alphat,patchi) const = 0 // [W/m/K]tmp<volScalarField> alphaEff(alphat) const = 0 // [kg/m/s]tmp<scalarField> alphaEff(alphat,patchi) const = 0 // [kg/m/s]Access to transport state variablesfluidThermovoid correctRho(const volScalarField& deltaRho) = 0const volScalarField& psi() const = 0 // [s^2/m^2]tmp<volScalarField> mu() const = 0 // [kg/m/s]tmp<scalarField> mu(patchi) const = 0 // [kg/m/s]tmp<volScalarField> nu() const // [m^2/s]tmp<scalarField> nu(patchi) const // [m^2/s]Member functionscompressibleTransportModelpsiThermovolScalarField psi_volScalarField mu_virtual tmp<volScalarField> rho() const {return p_*psi_}virtual tmp<scalarField> rho(patchi) constvirtual void correctRho(const volScalarField& deltaRho) {}virtual const volScalarField& psi() constvirtual tmp<volScalarField> mu() constvirtual tmp<scalarField> mu(patchi) constProtected dataCanteraMixtureIOdictionary CanteraTorchProperties_const word CanteraMechanismFile_std::shared_ptr<Cantera::Solution> CanteraSolution_std::shared_ptr<Cantera::ThermoPhase> CanteraGas_word transportModelName_std::shared_ptr<Cantera::Transport> CanteraTransport_hashedWordList species_PtrList<volScalarField> Y_const volScalarField& Tref_const volScalarField& pref_mutable scalarList yTemp_typedef CanteraMixture basicMixtureTypetypedef CanteraMixture thermoTypestatic const bool incompressible = falsestatic const bool isochoric = falsescalarList HaTemp_scalarList CpTemp_scalarList muTemp_CanteraMixture& cellMixture(celli)CanteraMixture& patchFaceMixture(patchi, facei)const scalar THE(scalar p, scalar T)const scalar psi(scalar p, scalar T)const scalar mu(scalar p, scalar T)const scalar alphah(scalar p, scalar T)const scalar HE(scalar p, scalar T)const scalar Hc(scalar p, scalar T)const scalar Cp(scalar p, scalar T)const scalar Cv(scalar p, scalar T)const scalar gamma(scalar p, scalar T)const scalar Cpv(scalar p, scalar T)const scalar CpByCpv(scalar p, scalar T)const scalar W() conststatic inline word heName() {return "ha"}PtrList<volScalarField>& Y()const PtrList<volScalarField>& Y() constvolScalarField& Y(const label i)const volScalarField& Y(const label i) constvolScalarField& Y(const word& specieName)const volScalarField& Y(const word& specieName) constconst hashedWordList& species()scalar nSpecies()std::shared_ptr<Cantera::ThermoPhase> CanteraGas()std::shared_ptr<Cantera::Solution> CanteraSolution()std::shared_ptr<Cantera::Transport> CanteraTransport()const word& transportModelName()void calcCp(const scalar T, const scalar p)void calcMu(const scalar T, const scalar p)void calcH(const scalar T, const scalar p)scalar Cp(label i, scalar p, scalar T) constscalar mu(label i, scalar p, scalar T) constscalar Ha(label i, scalar p, scalar T) constscalar Hc(label i) constscalar Hs(label i, scalar p, scalar T) constscalar kappa(label i, scalar p, scalar T) constscalar Wi(label i) constheThermoBasicThermo, MixtureTypevolScalarField he_Protected datavirtual typename MixtureType::basicMixtureType&composition(){return *this}virtual const typename MixtureType::basicMixtureType&composition() const{return *this}virtual word thermoName() const{return MixtureType::thermoType::typeName()}virtual volScalarField& he(){return he_} // [J/kg]virtual const volScalarField& he() const{return he_} // [J/kg]virtual tmp<volScalarField> he(p,T) const // [J/kg]virtual tmp<scalarField> he(p,T,const labelList& cells) const // [J/kg]virtual tmp<scalarField> he(p,T,patchi) const // [J/kg]virtual tmp<volScalarField> hc() const // [J/kg]virtual tmp<scalarField> THE(he, p, T0, const labelList& cells) constvirtual tmp<scalarField> THE(he, p, T0, patchi) constvirtual tmp<scalarField> Cp(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cp() const // [J/kg/K]virtual tmp<scalarField> Cv(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cv() const // [J/kg/K]virtual tmp<volScalarField> gamma() const // []virtual tmp<scalarField> gamma(p, T, patchi) const // [virtual tmp<scalarField> Cpv(p, T, patchi) const // [J/kg/K]virtual tmp<volScalarField> Cpv() const // [J/kg/K]virtual tmp<volScalarField> CpByCpv() const // []virtual tmp<scalarField> CpByCpv(p, T,patchi) const // []virtual tmp<volScalarField> W() const // [kg/kmol]virtual tmp<scalarField> W(const label patchi) const // [kg/kmol]virtual tmp<volScalarField> kappa() const // [W/m/K]virtual tmp<scalarField> kappa(patchi) const // [W/m/K]virtual tmp<volScalarField> alphahe() const // [kg/m/s]virtual tmp<scalarField> alphahe(patchi) const // [kg/m/s]virtual tmp<volScalarField> kappaEff(const volScalarField&) const // [W/m/K]virtual tmp<scalarField> kappaEff(alphat, patchi) const // [W/m/K]virtual tmp<volScalarField> alphaEff(alphat) const // [kg/m/s]virtual tmp<scalarField> alphaEff(alphat, patchi) const // [kg/m/s]Member functionshePsiThermoBasicPsiThermo, MixtureTypecalculate() : voidcorrect() : voiddfChemistryModelThermoTypeThermoType& thermo_;CanteraMixture& mixture_std::shared_ptr<Cantera::ThermoPhase> CanteraGas_const fvMesh& mesh_Switch chemistry_scalar relTol_scalar absTol_PtrList<volScalarField>& Y_PtrList<volScalarField> rhoD_;PtrList<volScalarField> hei_;mutable scalarList yTemp_mutable scalarList dTemp_;mutable scalarList hrtTemp_;mutable scalarField cTemp_;PtrList<volScalarField::Internal> RR_hashedWordList species_volScalarField& alpha_volScalarField& T_const volScalarField& p_const volScalarField& rho_;volScalarField& mu_volScalarField& psi_volScalarField Qdot_;Switch torchSwitch_;word torchModelName_;scalarList Xmu_;scalarList Xstd_;scalarList Ymu_;scalarList Ystd_;scalar Tact_;scalar Qdotact_;void setNumerics(Cantera::ReactorNet &r)scalar canteraSolve(const DeltaTType& deltaT)scalar torchSolve(const DeltaTType& deltaT)virtual scalar solve(const scalar deltaT)virtual scalar solve(const scalarField& deltaT)virtual const volScalarField::Internal& RR(const label i) constvirtual volScalarField::Internal& RR(const label i)virtual tmp<volScalarField> Qdot() constPtrList<volScalarField>& Y()const hashedWordList& species() constconst volScalarField& rhoD(const label i) constvolScalarField& hei(const label i)void correctThermo()MixtureTypeBasicThermoThermoType \ No newline at end of file diff --git a/src/dfChemistryModel/dfChemistryModel.C b/src/dfChemistryModel/dfChemistryModel.C index dd3aee81..3c09f12a 100644 --- a/src/dfChemistryModel/dfChemistryModel.C +++ b/src/dfChemistryModel/dfChemistryModel.C @@ -58,7 +58,7 @@ Foam::dfChemistryModel::dfChemistryModel absTol_(this->subDict("odeCoeffs").lookupOrDefault("absTol",1e-15)), Y_(mixture_.Y()), rhoD_(mixture_.nSpecies()), - hai_(mixture_.nSpecies()), + hei_(mixture_.nSpecies()), hc_(mixture_.nSpecies()), yTemp_(mixture_.nSpecies()), dTemp_(mixture_.nSpecies()), @@ -296,16 +296,16 @@ Foam::dfChemistryModel::dfChemistryModel ); } - forAll(hai_, i) + forAll(hei_, i) { - hai_.set + hei_.set ( i, new volScalarField ( IOobject ( - "hai_" + Y_[i].name(), + "hei_" + Y_[i].name(), mesh_.time().timeName(), mesh_, IOobject::NO_READ, @@ -428,7 +428,7 @@ void Foam::dfChemistryModel::correctThermo() forAll(Y_, i) { yTemp_[i] = Y_[i][celli]; - hc += yTemp_[i]*mixture_.Hc(i); + hc += yTemp_[i]*hc_[i]; } if(useThermoTranNN) { @@ -481,12 +481,12 @@ void Foam::dfChemistryModel::correctThermo() rhoD_[4][celli] = data_ptr4[6*celli+3]; //CO2 rhoD_[5][celli] = data_ptr4[6*celli+5]; //N2 - hai_[0][celli] = data_ptr5[6*celli]; // O2 - hai_[1][celli] = data_ptr5[6*celli+4]; //H2O - hai_[2][celli] = data_ptr5[6*celli+1]; //CH4 - hai_[3][celli] = data_ptr5[6*celli+2]; //CO - hai_[4][celli] = data_ptr5[6*celli+3]; //CO2 - hai_[5][celli] = data_ptr5[6*celli+5]; //N2 + hei_[0][celli] = data_ptr5[6*celli]; // O2 + hei_[1][celli] = data_ptr5[6*celli+4]; //H2O + hei_[2][celli] = data_ptr5[6*celli+1]; //CH4 + hei_[3][celli] = data_ptr5[6*celli+2]; //CO + hei_[4][celli] = data_ptr5[6*celli+3]; //CO2 + hei_[5][celli] = data_ptr5[6*celli+5]; //N2 } #endif @@ -552,7 +552,24 @@ void Foam::dfChemistryModel::correctThermo() rhoD_[i][celli] = rho_[celli]*dTemp_[i]; // CanteraGas_->molecularWeight(i) kg/kmol - hai_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + if(mixture_.heName() == "ha") + { + hei_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + } + else if(mixture_.heName() == "hs") + { + hei_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - hc_[i]; + } + else if(mixture_.heName() == "ea") + { + hei_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - p_[celli]/rho_[celli]; + } + else + { + FatalErrorInFunction + << "Wrong CanteraMixture enregy name: " << CanteraMixture::getEnergyName() + << abort(FatalError); + } } } } @@ -593,7 +610,7 @@ void Foam::dfChemistryModel::correctThermo() forAll(Y_, i) { yTemp_[i] = Y_[i].boundaryField()[patchi][facei]; - hc += yTemp_[i]*mixture_.Hc(i); + hc += yTemp_[i]*hc_[i]; } mixture_.setState_TPY(pT[facei], pp[facei], yTemp_.begin()); @@ -624,7 +641,24 @@ void Foam::dfChemistryModel::correctThermo() { rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i]; - hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + if(mixture_.heName() == "ha") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + } + else if(mixture_.heName() == "hs") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - hc_[i]; + } + else if(mixture_.heName() == "ea") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - pp[facei]/prho[facei]; + } + else + { + FatalErrorInFunction + << "Wrong CanteraMixture enregy name: " << CanteraMixture::getEnergyName() + << abort(FatalError); + } } } } @@ -637,7 +671,7 @@ void Foam::dfChemistryModel::correctThermo() forAll(Y_, i) { yTemp_[i] = Y_[i].boundaryField()[patchi][facei]; - hc += yTemp_[i]*mixture_.Hc(i); + hc += yTemp_[i]*hc_[i]; } if(useThermoTranNN) { @@ -681,12 +715,12 @@ void Foam::dfChemistryModel::correctThermo() rhoD_[4].boundaryFieldRef()[patchi][facei] = pdata_ptr4[6*facei+3]; //CO2 rhoD_[5].boundaryFieldRef()[patchi][facei] = pdata_ptr4[6*facei+5]; //N2 - hai_[0].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei]; // O2 - hai_[1].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+4]; //H2O - hai_[2].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+1]; //CH4 - hai_[3].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+2]; //CO - hai_[4].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+3]; //CO2 - hai_[5].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+5]; //N2 + hei_[0].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei]; // O2 + hei_[1].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+4]; //H2O + hei_[2].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+1]; //CH4 + hei_[3].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+2]; //CO + hei_[4].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+3]; //CO2 + hei_[5].boundaryFieldRef()[patchi][facei] = pdata_ptr5[6*facei+5]; //N2 } #endif @@ -746,7 +780,24 @@ void Foam::dfChemistryModel::correctThermo() { rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i]; - hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + if(mixture_.heName() == "ha") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i); + } + else if(mixture_.heName() == "hs") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - hc_[i]; + } + else if(mixture_.heName() == "ea") + { + hei_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i) - pp[facei]/prho[facei]; + } + else + { + FatalErrorInFunction + << "Wrong CanteraMixture enregy name: " << CanteraMixture::getEnergyName() + << abort(FatalError); + } } } } diff --git a/src/dfChemistryModel/dfChemistryModel.H b/src/dfChemistryModel/dfChemistryModel.H index 4e8850e5..59ea2483 100644 --- a/src/dfChemistryModel/dfChemistryModel.H +++ b/src/dfChemistryModel/dfChemistryModel.H @@ -106,7 +106,7 @@ public IOdictionary // species mass diffusion coefficients, [kg/m/s] PtrList rhoD_; // species absolute enthalpy, [J/kg] - PtrList hai_; + PtrList hei_; // species chemistry enthalpy, [J/kg] scalarList hc_; // temp mass fraction @@ -312,11 +312,11 @@ public: const volScalarField& rhoD(const label i) const {return rhoD_[i];} - const volScalarField& hai(const label i) {return hai_[i];} + const volScalarField& hei(const label i) {return hei_[i];} const scalar & hci(const label i) {return hc_[i];} - // update T, psi, mu, alpha, rhoD, hai (if needed) + // update T, psi, mu, alpha, rhoD, hei (if needed) void correctThermo(); ThermoType& thermo() {return thermo_;} diff --git a/src_gpu/dfYEqn.H b/src_gpu/dfYEqn.H index b3715efc..c0ccd5a2 100644 --- a/src_gpu/dfYEqn.H +++ b/src_gpu/dfYEqn.H @@ -131,22 +131,22 @@ public: void yeqn_compute_thermo_alpha(cudaStream_t stream, int num_cells, const double *rhoD, double *thermo_alpha, int num_boundary_surfaces, const double *boundary_rhoD, double *boundary_thermo_alpha); - void yeqn_compute_DEff_via_lewisNumber(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, + void yeqn_compute_DEff_via_lewisNumber(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, double *lewis_number, const double *alpha, const double *mut_sct, double *DEff, const double *boundary_alpha, const double *boundary_mut_sct, double *boundary_DEff); - void yeqn_compute_RR(dfChemistrySolver& chemistrySolver, cudaStream_t stream, const double *h_T, const double *d_T, + void yeqn_compute_RR(dfChemistrySolver& chemistrySolver, cudaStream_t stream, const double *h_T, const double *d_T, const double *p, const double *y, const double *rho, double *RR); void yeqn_fvc_laplacian_scalar(cudaStream_t stream, ncclComm_t comm, const int *neighbor_peer, int num_species, int num_cells, int num_surfaces, int num_boundary_surfaces, const int *lowerAddr, const int *upperAddr, const double *weight, const double *mag_sf, const double *delta_coeffs, const double *volume, - const double *thermo_alpha, const double *hai, const double *vf, double *output, // end for internal + const double *thermo_alpha, const double *hei, const double *vf, double *output, // end for internal int num_patches, const int *patch_size, const int *patch_type, const int *boundary_cell_face, const double *boundary_weight, const double *boundary_mag_sf, const double *boundary_delta_coeffs, const double *boundary_thermo_alpha, const double *boundary_hai, const double *boundary_vf, const int *cyclicNeighbor, const int *patchSizeOffset, double *boundary_output); void yeqn_compute_sumYDiffError_and_hDiffCorrFlux(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, - const double *rhoD, const double *hai, const double *y, const double *grad_y, + const double *rhoD, const double *hei, const double *y, const double *grad_y, double *sumY_diff_error, double *hDiff_corr_flux, const double *boundary_hai, const double *boundary_y, const double *boundary_grad_y, const double *boundary_rhoD, double *boundary_sumY_diff_error, double *boundary_hDiff_corr_flux); diff --git a/src_gpu/dfYEqn.cu b/src_gpu/dfYEqn.cu index b86cb5a6..dbb06846 100644 --- a/src_gpu/dfYEqn.cu +++ b/src_gpu/dfYEqn.cu @@ -56,14 +56,14 @@ __global__ void yeqn_compute_phiUc_internal(int num_cells, int num_surfaces, double sfy = sf[num_surfaces * 1 + index]; double sfz = sf[num_surfaces * 2 + index]; - double w = weight[index]; + double w = weight[index]; double ssfx = (w * (sumY_diff_error[num_cells * 0 + owner] - sumY_diff_error[num_cells * 0 + neighbor]) + sumY_diff_error[num_cells * 0 + neighbor]); double ssfy = (w * (sumY_diff_error[num_cells * 1 + owner] - sumY_diff_error[num_cells * 1 + neighbor]) + sumY_diff_error[num_cells * 1 + neighbor]); double ssfz = (w * (sumY_diff_error[num_cells * 2 + owner] - sumY_diff_error[num_cells * 2 + neighbor]) + sumY_diff_error[num_cells * 2 + neighbor]); phiUc[index] = sfx * ssfx + sfy * ssfy + sfz * ssfz; } - + __global__ void yeqn_compute_phiUc_boundary(int num_boundary_surfaces, const double *boundary_sf, const double *boundary_sumY_diff_error, double *boundary_phiUc) { @@ -81,9 +81,9 @@ __global__ void yeqn_compute_phiUc_boundary(int num_boundary_surfaces, boundary_phiUc[index] = boundary_sfx * boundary_ssfx + boundary_sfy * boundary_ssfy + boundary_sfz * boundary_ssfz; } - + __global__ void yeqn_sumError_and_compute_hDiffCorrFlux(int num_species, int num, - const double *rhoD, const double *hai, const double *y, const double *grady, + const double *rhoD, const double *hei, const double *y, const double *grady, double *sum_rhoD_grady, double *hDiffCorrFlux) { int index = blockDim.x * blockIdx.x + threadIdx.x; @@ -98,7 +98,7 @@ __global__ void yeqn_sumError_and_compute_hDiffCorrFlux(int num_species, int num double sum_rhoD_grady_z = 0; double sum_hai_y = 0; for (int s = 0; s < num_species; s++) { - double hai_value = hai[num * s + index]; + double hai_value = hei[num * s + index]; double rhoD_value = rhoD[num * s + index]; // le = alpha/D double y_value = y[num * s + index]; double grady_x = grady[num * s * 3 + num * 0 + index]; @@ -123,7 +123,7 @@ __global__ void yeqn_sumError_and_compute_hDiffCorrFlux(int num_species, int num __global__ void yeqn_fvc_laplacian_scalar_internal(int num_species, int num_cells, int num_surfaces, const int *lower_index, const int *upper_index, const double *mag_sf, const double *delta_coeffs, const double *weight, - const double *thermo_alpha, const double *hai, const double *vf, double *output) + const double *thermo_alpha, const double *hei, const double *vf, double *output) { int index = blockDim.x * blockIdx.x + threadIdx.x; if (index >= num_surfaces) @@ -142,14 +142,14 @@ __global__ void yeqn_fvc_laplacian_scalar_internal(int num_species, int num_cell // printf("input index: %d, thermo: %.16lf, %.16lf\n", index, thermo_alpha_owner, thermo_alpha_neighbor); double sum_ssf = 0; for (int s = 0; s < num_species; s++) { - double haii_owner = hai[num_cells * s + owner]; - double haii_neighbor = hai[num_cells * s + neighbor]; + double haii_owner = hei[num_cells * s + owner]; + double haii_neighbor = hei[num_cells * s + neighbor]; double gamma = w * (thermo_alpha_owner * haii_owner) + (1 - w) * (thermo_alpha_neighbor * haii_neighbor); double sngrad = delta_coeff * (vf[num_cells * s + neighbor] - vf[num_cells * s + owner]); double ssf = gamma * sngrad * magsf; sum_ssf += ssf; //if (owner == 21 || neighbor == 21) - // printf("hai: %.16lf, %.16lf, gamma: %.16lf, sngrad: %.16lf, ssf: %.16lf\n", haii_owner, haii_neighbor, gamma, sngrad, ssf); + // printf("hei: %.16lf, %.16lf, gamma: %.16lf, sngrad: %.16lf, ssf: %.16lf\n", haii_owner, haii_neighbor, gamma, sngrad, ssf); } // owner @@ -521,7 +521,7 @@ void dfYEqn::process() { update_boundary_coeffs_scalar(dataBase_.stream, dataBase_.num_patches, dataBase_.patch_size.data(), patch_type.data(), dataBase_.d_boundary_delta_coeffs, dataBase_.d_boundary_y + dataBase_.num_boundary_surfaces * s, - dataBase_.d_boundary_weight, + dataBase_.d_boundary_weight, d_value_internal_coeffs + dataBase_.num_boundary_surfaces * s, d_value_boundary_coeffs + dataBase_.num_boundary_surfaces * s, d_gradient_internal_coeffs + dataBase_.num_boundary_surfaces * s, @@ -530,7 +530,7 @@ void dfYEqn::process() { // compute sumYDiffError and hDiffCorrFlux yeqn_compute_sumYDiffError_and_hDiffCorrFlux(dataBase_.stream, dataBase_.num_species, dataBase_.num_cells, dataBase_.num_boundary_surfaces, - dataBase_.d_thermo_rhoD, d_hai, dataBase_.d_y, d_grad_y, + dataBase_.d_thermo_rhoD, d_hai, dataBase_.d_y, d_grad_y, d_sumY_diff_error, dataBase_.d_hDiff_corr_flux, d_boundary_hai, dataBase_.d_boundary_y, d_boundary_grad_y, dataBase_.d_boundary_thermo_rhoD, d_boundary_sumY_diff_error, dataBase_.d_boundary_hDiff_corr_flux); @@ -584,11 +584,11 @@ void dfYEqn::process() { dataBase_.d_rho, dataBase_.d_rho_old, dataBase_.d_y + dataBase_.num_cells * s, dataBase_.d_volume, d_diag, d_source, 1.); // **calculate div weights with limitedLinear scheme** - // compute_limitedLinear_weight(dataBase_.stream, dataBase_.nccl_comm, dataBase_.neighbProcNo.data(), dataBase_.num_surfaces, + // compute_limitedLinear_weight(dataBase_.stream, dataBase_.nccl_comm, dataBase_.neighbProcNo.data(), dataBase_.num_surfaces, // dataBase_.num_cells, dataBase_.num_boundary_surfaces, dataBase_.d_owner, dataBase_.d_neighbor, dataBase_.d_mesh_dis, // dataBase_.d_weight, dataBase_.d_sf, dataBase_.d_y + dataBase_.num_cells * s, dataBase_.d_phi, dataBase_.d_phi_weight, // dataBase_.num_patches, dataBase_.patch_size.data(), patch_type.data(), dataBase_.d_boundary_weight, dataBase_.d_boundary_face_cell, - // dataBase_.d_boundary_y + dataBase_.num_boundary_surfaces * s, dataBase_.d_boundary_sf, dataBase_.d_volume, + // dataBase_.d_boundary_y + dataBase_.num_boundary_surfaces * s, dataBase_.d_boundary_sf, dataBase_.d_volume, // dataBase_.d_boundary_mag_sf, dataBase_.d_boundary_phi, dataBase_.cyclicNeighbor.data(), dataBase_.patchSizeOffset.data(), // dataBase_.d_boundary_delta_coeffs); @@ -646,7 +646,7 @@ void dfYEqn::process() { correct_boundary_conditions_scalar(dataBase_.stream, dataBase_.nccl_comm, dataBase_.neighbProcNo.data(), dataBase_.num_boundary_surfaces, dataBase_.num_patches, dataBase_.patch_size.data(), patch_type.data(), dataBase_.d_boundary_delta_coeffs, dataBase_.d_boundary_face_cell, - dataBase_.d_y + dataBase_.num_cells * s, dataBase_.d_boundary_y + dataBase_.num_boundary_surfaces * s, + dataBase_.d_y + dataBase_.num_cells * s, dataBase_.d_boundary_y + dataBase_.num_boundary_surfaces * s, dataBase_.cyclicNeighbor.data(), dataBase_.patchSizeOffset.data(), dataBase_.d_boundary_weight); } TICK_END_EVENT(YEqn post process for all species correctBC); @@ -694,7 +694,7 @@ void dfYEqn::process() { sync(); } -void dfYEqn::solve(int speciesIndex) { +void dfYEqn::solve(int speciesIndex) { TICK_INIT_EVENT; TICK_START_EVENT; dataBase_.solve(num_iteration, AMGXSetting::u_setting, d_A, dataBase_.d_y + dataBase_.num_cells * speciesIndex, d_source); @@ -723,7 +723,7 @@ void dfYEqn::yeqn_compute_thermo_alpha(cudaStream_t stream, num_boundary_surfaces, boundary_rhoD, boundary_thermo_alpha); } -void dfYEqn::yeqn_compute_DEff_via_lewisNumber(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, +void dfYEqn::yeqn_compute_DEff_via_lewisNumber(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, double *lewis_number, const double *alpha, const double *mut_sct, double *DEff, const double *boundary_alpha, const double *boundary_mut_sct, double *boundary_DEff) { @@ -746,7 +746,7 @@ void dfYEqn::yeqn_fvc_laplacian_scalar(cudaStream_t stream, ncclComm_t comm, con int num_species, int num_cells, int num_surfaces, int num_boundary_surfaces, const int *lowerAddr, const int *upperAddr, const double *weight, const double *mag_sf, const double *delta_coeffs, const double *volume, - const double *thermo_alpha, const double *hai, const double *vf, double *output, // end for internal + const double *thermo_alpha, const double *hei, const double *vf, double *output, // end for internal int num_patches, const int *patch_size, const int *patch_type, const int *boundary_cell_face, const double *boundary_weight, const double *boundary_mag_sf, const double *boundary_delta_coeffs, const double *boundary_thermo_alpha, const double *boundary_hai, const double *boundary_vf, @@ -755,7 +755,7 @@ void dfYEqn::yeqn_fvc_laplacian_scalar(cudaStream_t stream, ncclComm_t comm, con size_t threads_per_block = 1024; size_t blocks_per_grid = (num_surfaces + threads_per_block - 1) / threads_per_block; yeqn_fvc_laplacian_scalar_internal<<>>(num_species, num_cells, num_surfaces, - lowerAddr, upperAddr, mag_sf, delta_coeffs, weight, thermo_alpha, hai, vf, output); + lowerAddr, upperAddr, mag_sf, delta_coeffs, weight, thermo_alpha, hei, vf, output); int offset = 0; for (int i = 0; i < num_patches; i++) { @@ -775,7 +775,7 @@ void dfYEqn::yeqn_fvc_laplacian_scalar(cudaStream_t stream, ncclComm_t comm, con boundary_thermo_alpha, boundary_hai, vf, boundary_vf, output); } else if (patch_type[i] == boundaryConditions::cyclic) { yeqn_fvc_laplacian_scalar_boundary_cyclic<<>>( - num_species, num_cells, num_boundary_surfaces, patch_size[i], offset, patchSizeOffset[cyclicNeighbor[i]], + num_species, num_cells, num_boundary_surfaces, patch_size[i], offset, patchSizeOffset[cyclicNeighbor[i]], boundary_cell_face, boundary_mag_sf, boundary_delta_coeffs, boundary_thermo_alpha, boundary_hai, vf, boundary_vf, output); } else if (patch_type[i] == boundaryConditions::processor @@ -819,7 +819,7 @@ void dfYEqn::yeqn_fvc_laplacian_scalar(cudaStream_t stream, ncclComm_t comm, con } void dfYEqn::yeqn_compute_sumYDiffError_and_hDiffCorrFlux(cudaStream_t stream, int num_species, int num_cells, int num_boundary_surfaces, - const double *rhoD, const double *hai, const double *y, const double *grad_y, + const double *rhoD, const double *hei, const double *y, const double *grad_y, double *sumY_diff_error, double *hDiff_corr_flux, const double *boundary_hai, const double *boundary_y, const double *boundary_grad_y, const double *boundary_rhoD, double *boundary_sumY_diff_error, double *boundary_hDiff_corr_flux) @@ -827,7 +827,7 @@ void dfYEqn::yeqn_compute_sumYDiffError_and_hDiffCorrFlux(cudaStream_t stream, i size_t threads_per_block = 1024; size_t blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; yeqn_sumError_and_compute_hDiffCorrFlux<<>>(num_species, num_cells, - rhoD, hai, y, grad_y, sumY_diff_error, hDiff_corr_flux); + rhoD, hei, y, grad_y, sumY_diff_error, hDiff_corr_flux); blocks_per_grid = (num_boundary_surfaces + threads_per_block - 1) / threads_per_block; yeqn_sumError_and_compute_hDiffCorrFlux<<>>(num_species, num_boundary_surfaces, boundary_rhoD, boundary_hai, boundary_y, boundary_grad_y, boundary_sumY_diff_error, boundary_hDiff_corr_flux); diff --git a/test/corrtest.cpp b/test/corrtest.cpp index 7c13ea3a..3eab56a1 100644 --- a/test/corrtest.cpp +++ b/test/corrtest.cpp @@ -49,11 +49,11 @@ TEST(corrtest,dfHighSpeedFoam){ } TEST(corrtest,dfLowMachFoam_TGV){ - EXPECT_FLOAT_EQ(TGV500,1532.92); // compare the maximum temperature along y direction in 2D TGV after 500 time steps - EXPECT_FLOAT_EQ(TGV400,1297.64); // ..........400 time steps - EXPECT_FLOAT_EQ(TGV300,871.092); - EXPECT_FLOAT_EQ(TGV200,537.614); - EXPECT_FLOAT_EQ(TGV100,363.504); + EXPECT_FLOAT_EQ(TGV500,1535.03); // compare the maximum temperature along y direction in 2D TGV after 500 time steps + EXPECT_FLOAT_EQ(TGV400,1298.64); // ..........400 time steps + EXPECT_FLOAT_EQ(TGV300,874.13098); + EXPECT_FLOAT_EQ(TGV200,538.34802); + EXPECT_FLOAT_EQ(TGV100,363.06601); } TEST(corrtest,2DSandia){ @@ -71,10 +71,10 @@ TEST(corrtest,2DSandia){ } TEST(corrtest,dfLowMachFoam_2DaachenBomb){ - EXPECT_NEAR(aachenBomb1,809.163,30); - EXPECT_NEAR(aachenBomb2,1793.4,30); - EXPECT_NEAR(aachenBomb3,908.059,30); - EXPECT_NEAR(aachenBomb4,2493.09,30); + EXPECT_NEAR(aachenBomb1,809.39398193359375,0.0001); + EXPECT_NEAR(aachenBomb2,1814.1099853515625,0.0001); + EXPECT_NEAR(aachenBomb3,912.61602783203125,0.0001); + EXPECT_NEAR(aachenBomb4,2492.31005859375,0.0001); } float readmaxTH2(){