From e9bdd8c390491a31bcf601fbd0dadceb5c393c9f Mon Sep 17 00:00:00 2001 From: pkuLmq <1900011062@pku.edu.cn> Date: Fri, 30 Jun 2023 00:30:14 +0800 Subject: [PATCH 1/2] realize RK2SSP and RK3SSP in dfHighSpeedFoam and add wrate opption for chemistry --- .../solvers/dfHighSpeedFoam/Make/files | 4 +- .../solvers/dfHighSpeedFoam/calculateR.H | 49 +++ .../solvers/dfHighSpeedFoam/createFields.H | 39 ++ .../dfHighSpeedFoam/createFields_rk2.H | 27 ++ .../solvers/dfHighSpeedFoam/dfHighSpeedFoam.C | 341 ------------------ .../dfHighSpeedFoam/dfHighSpeedFoam_new.C | 255 +++++++++++++ applications/solvers/dfHighSpeedFoam/phiCal.H | 45 +++ applications/solvers/dfHighSpeedFoam/preCal.H | 118 ++++++ .../solvers/dfHighSpeedFoam/rhoEEqn.H | 67 +++- applications/solvers/dfHighSpeedFoam/rhoEqn.H | 25 +- .../solvers/dfHighSpeedFoam/rhoUEqn.H | 57 ++- .../solvers/dfHighSpeedFoam/rhoYEqn.H | 122 +++++-- .../dfHighSpeedFoam/updateFieldsSave.H | 12 + 13 files changed, 744 insertions(+), 417 deletions(-) create mode 100644 applications/solvers/dfHighSpeedFoam/calculateR.H create mode 100644 applications/solvers/dfHighSpeedFoam/createFields_rk2.H delete mode 100644 applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C create mode 100644 applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C create mode 100644 applications/solvers/dfHighSpeedFoam/phiCal.H create mode 100644 applications/solvers/dfHighSpeedFoam/preCal.H create mode 100644 applications/solvers/dfHighSpeedFoam/updateFieldsSave.H diff --git a/applications/solvers/dfHighSpeedFoam/Make/files b/applications/solvers/dfHighSpeedFoam/Make/files index 70642de5e..1d93e0527 100644 --- a/applications/solvers/dfHighSpeedFoam/Make/files +++ b/applications/solvers/dfHighSpeedFoam/Make/files @@ -1,3 +1,3 @@ -dfHighSpeedFoam.C +dfHighSpeedFoam_new.C -EXE = $(DF_APPBIN)/dfHighSpeedFoam +EXE = $(DF_APPBIN)/dfHighSpeedFoam_new diff --git a/applications/solvers/dfHighSpeedFoam/calculateR.H b/applications/solvers/dfHighSpeedFoam/calculateR.H new file mode 100644 index 000000000..651b00ad8 --- /dev/null +++ b/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; + diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H index ddffc838e..e01dab723 100644 --- a/applications/solvers/dfHighSpeedFoam/createFields.H +++ b/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(mesh, word::null); diff --git a/applications/solvers/dfHighSpeedFoam/createFields_rk2.H b/applications/solvers/dfHighSpeedFoam/createFields_rk2.H new file mode 100644 index 000000000..25bdeffc3 --- /dev/null +++ b/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 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] + ) + ); +} + diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C deleted file mode 100644 index 67254ec4d..000000000 --- a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C +++ /dev/null @@ -1,341 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 . - -Application - rhoCentralFoam - -Description - Density-based compressible flow solver based on central-upwind schemes of - Kurganov and Tadmor with support for mesh-motion and topology changes. - -\*---------------------------------------------------------------------------*/ - -#include "dfChemistryModel.H" -#include "CanteraMixture.H" -// #include "hePsiThermo.H" -#include "heRhoThermo.H" - -#ifdef USE_PYTORCH -#include -#include -#include //used to convert -#endif - -#ifdef USE_LIBTORCH -#include -#include "DNNInferencer.H" -#endif - -#include "fvCFD.H" -#include "dynamicFvMesh.H" -// #include "psiThermo.H" -#include "rhoThermo.H" -#include "turbulentFluidThermoModel.H" -#include "fixedRhoFvPatchScalarField.H" -#include "directionInterpolate.H" -#include "localEulerDdtScheme.H" -#include "fvcSmooth.H" -#include "PstreamGlobals.H" -#include "CombustionModel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ -#ifdef USE_PYTORCH - pybind11::scoped_interpreter guard{};//start python interpreter -#endif - #define NO_CONTROL - - #include "postProcess.H" - - // #include "setRootCaseLists.H" - #include "listOptions.H" - #include "setRootCase2.H" - #include "listOutput.H" - #include "createTime.H" - #include "createDynamicFvMesh.H" - #include "createFields.H" - #include "createTimeControls.H" - - double time_monitor_flow=0; - double time_monitor_chem=0; - double time_monitor_Y=0; - double time_monitor_AMR=0; - clock_t start, end; - - turbulence->validate(); - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - #include "readFluxScheme.H" - - dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0); - dimensionedScalar refCri("refCri", dimensionSet(1, -4, 0, 0, 0), 0.0); - - // Courant numbers used to adjust the time-step - scalar CoNum = 0.0; - scalar meanCoNum = 0.0; - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.run()) - { - #include "readTimeControls.H" - - //used for AMR - refCri = max(mag(fvc::grad(rho))); - tmp tmagGradrho = mag(fvc::grad(rho)); - volScalarField normalisedGradrho - ( - "normalisedGradrho", - tmagGradrho()/refCri - ); - normalisedGradrho.writeOpt() = IOobject::AUTO_WRITE; - tmagGradrho.clear(); - - if (!LTS) - { - #include "setDeltaT.H" - runTime++; - - // Do any mesh changes - start = std::clock(); - mesh.update(); - end = std::clock(); - time_monitor_AMR += double(end - start) / double(CLOCKS_PER_SEC); - - } - - // --- Directed interpolation of primitive fields onto faces - - surfaceScalarField rho_pos(interpolate(rho, pos)); - surfaceScalarField rho_neg(interpolate(rho, neg)); - - PtrList rhoYi_pos(nspecies); - PtrList rhoYi_neg(nspecies); - forAll(rhoYi_pos,i) - { - rhoYi_pos.set - ( - i, - new surfaceScalarField - ( - IOobject - ( - "rhoYi_pos" + Y[i].name(), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - interpolate(rhoYi[i], pos,"Yi") - ) - ); - } - - forAll(rhoYi_neg,i) - { - rhoYi_neg.set - ( - i, - new surfaceScalarField - ( - IOobject - ( - "rhoYi_neg" + Y[i].name(), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - interpolate(rhoYi[i], neg,"Yi") - ) - ); - } - - surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name())); - surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name())); - - volScalarField rPsi("rPsi", 1.0/psi); - surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name())); - surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name())); - - surfaceScalarField ea_pos(interpolate(ea, pos, T.name())); - surfaceScalarField ea_neg(interpolate(ea, neg, T.name())); - - surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos); - 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(); - } - - volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); - surfaceScalarField cSf_pos - ( - "cSf_pos", - interpolate(c, pos, T.name())*mesh.magSf() - ); - surfaceScalarField cSf_neg - ( - "cSf_neg", - interpolate(c, neg, T.name())*mesh.magSf() - ); - - 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); - - 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); - - // Reuse amaxSf for the maximum positive and negative fluxes - // estimated by the central scheme - amaxSf = max(mag(aphiv_pos), mag(aphiv_neg)); - - #include "centralCourantNo.H" - - if (LTS) - { - #include "setRDeltaT.H" - runTime++; - } - - Info<< "Time = " << runTime.timeName() << nl << endl; - - phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; - - PtrList phiYi(nspecies); - forAll(phiYi,i) - { - phiYi.set - ( - i, - new surfaceScalarField - ( - IOobject - ( - "phiYi_" + Y[i].name(), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i] - ) - ); - } - - 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)))); - - // --- Solve density - #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" - - // --- Solve energy - #include "rhoEEqn.H" - - turbulence->correct(); - - 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<< "MonitorTime_AMR = " << time_monitor_AMR << " s" << nl << endl; - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return 0; -} - -// ************************************************************************* // diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C new file mode 100644 index 000000000..39c195cdb --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C @@ -0,0 +1,255 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Application + rhoCentralFoam + +Description + Density-based compressible flow solver based on central-upwind schemes of + Kurganov and Tadmor with support for mesh-motion and topology changes. + +\*---------------------------------------------------------------------------*/ + +#include "dfChemistryModel.H" +#include "CanteraMixture.H" +// #include "hePsiThermo.H" +#include "heRhoThermo.H" + +#ifdef USE_PYTORCH +#include +#include +#include //used to convert +#endif + +#ifdef USE_LIBTORCH +#include +#include "DNNInferencer.H" +#endif + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +// #include "psiThermo.H" +#include "rhoThermo.H" +#include "turbulentFluidThermoModel.H" +#include "fixedRhoFvPatchScalarField.H" +#include "directionInterpolate.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" +#include "PstreamGlobals.H" +#include "CombustionModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +#ifdef USE_PYTORCH + pybind11::scoped_interpreter guard{};//start python interpreter +#endif + #define NO_CONTROL + + #include "postProcess.H" + + // #include "setRootCaseLists.H" + #include "listOptions.H" + #include "setRootCase2.H" + #include "listOutput.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "createFields.H" + #include "createFields_rk2.H" + #include "createTimeControls.H" + + double time_monitor_flow=0; + double time_monitor_chem=0; + double time_monitor_Y=0; + double time_monitor_AMR=0; + clock_t start, end; + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + #include "readFluxScheme.H" + + dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0); + dimensionedScalar refCri("refCri", dimensionSet(1, -4, 0, 0, 0), 0.0); + + // Courant numbers used to adjust the time-step + scalar CoNum = 0.0; + scalar meanCoNum = 0.0; + + std::vector rkcoe1(3); + std::vector rkcoe2(3); + std::vector rkcoe3(3); + scalar rk; + label nrk=0; + + if(ddtSchemes == "RK2SSP") + { + rkcoe1[0]=1.0; rkcoe2[0]=0.0; rkcoe3[0]=1.0; + rkcoe1[1]=0.5; rkcoe2[1]=0.5; rkcoe3[1]=0.5; + rkcoe1[2]=0.0; rkcoe2[2]=0.0; rkcoe3[2]=0.0; + rk=2; + } + else if(ddtSchemes == "RK3SSP") + { + rkcoe1[0]=1.0; rkcoe2[0]=0.0; rkcoe3[0]=1.0; + rkcoe1[1]=0.75; rkcoe2[1]=0.25; rkcoe3[1]=0.25; + rkcoe1[2]=0.33; rkcoe2[2]=0.66; rkcoe3[2]=0.66; + rk=3; + } + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + //used for AMR + refCri = max(mag(fvc::grad(rho))); + tmp tmagGradrho = mag(fvc::grad(rho)); + volScalarField normalisedGradrho + ( + "normalisedGradrho", + tmagGradrho()/refCri + ); + normalisedGradrho.writeOpt() = IOobject::AUTO_WRITE; + tmagGradrho.clear(); + + if (!LTS) + { + #include "setDeltaT.H" + runTime++; + + // Do any mesh changes + start = std::clock(); + mesh.update(); + end = std::clock(); + time_monitor_AMR += double(end - start) / double(CLOCKS_PER_SEC); + + } + + volScalarField rho_rhs("rho_rhs",rho_save/runTime.deltaT()); + volVectorField rhoU_rhs("rhoU_rhs",rhoU_save/runTime.deltaT()); + volScalarField rhoYi_rhs("rhoYi_rhs",rhoYi_save[0]/runTime.deltaT()); + volScalarField rhoE_rhs("rhoE_rhs",rhoE_save/runTime.deltaT()); + + if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + for( nrk=0 ; nrkcorrect(); + + 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<< "MonitorTime_AMR = " << time_monitor_AMR << " s" << nl << endl; + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + +// ************************************************************************* // diff --git a/applications/solvers/dfHighSpeedFoam/phiCal.H b/applications/solvers/dfHighSpeedFoam/phiCal.H new file mode 100644 index 000000000..e5f569ed8 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/phiCal.H @@ -0,0 +1,45 @@ +phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; + +PtrList phiYi(nspecies); +forAll(phiYi,i) +{ + phiYi.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "phiYi_" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i] + ) + ); +} + +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 new file mode 100644 index 000000000..879170583 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/preCal.H @@ -0,0 +1,118 @@ +// --- Directed interpolation of primitive fields onto faces +surfaceScalarField rho_pos(interpolate(rho, pos)); +surfaceScalarField rho_neg(interpolate(rho, neg)); + +PtrList rhoYi_pos(nspecies); +PtrList rhoYi_neg(nspecies); +forAll(rhoYi_pos,i) +{ + rhoYi_pos.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "rhoYi_pos" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interpolate(rhoYi[i], pos,"Yi") + ) + ); +} + +forAll(rhoYi_neg,i) +{ + rhoYi_neg.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "rhoYi_neg" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interpolate(rhoYi[i], neg,"Yi") + ) + ); +} + +surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name())); +surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name())); + +volScalarField rPsi("rPsi", 1.0/psi); +surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name())); +surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name())); + +surfaceScalarField ea_pos(interpolate(ea, pos, T.name())); +surfaceScalarField ea_neg(interpolate(ea, neg, T.name())); + +surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos); +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(); +} + +volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); +surfaceScalarField cSf_pos +( + "cSf_pos", + interpolate(c, pos, T.name())*mesh.magSf() +); +surfaceScalarField cSf_neg +( + "cSf_neg", + interpolate(c, neg, T.name())*mesh.magSf() +); + +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); + +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); + +// 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 diff --git a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H index 532d064d1..2b96e8862 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H @@ -8,12 +8,23 @@ surfaceScalarField sigmaDotU & (a_pos*U_pos + a_neg*U_neg) ); -solve -( - fvm::ddt(rhoE) - + fvc::div(phiEp) - - fvc::div(sigmaDotU) -); +if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rhoE_rhs = -fvc::div(phiEp)+fvc::div(sigmaDotU); + + rhoE = rkcoe1[nrk]*rhoE_save + + rkcoe2[nrk]*rhoE + + rkcoe3[nrk]*rhoE_rhs*runTime.deltaT(); +} +else +{ + solve + ( + fvm::ddt(rhoE) + + fvc::div(phiEp) + - fvc::div(sigmaDotU) + ); +} ea = rhoE/rho - 0.5*magSqr(U); ea.correctBoundaryConditions(); @@ -25,21 +36,39 @@ rhoE.boundaryFieldRef() == rho.boundaryField()*(ea.boundaryField() + 0.5*magSqr( if (!inviscid) { - fvScalarMatrix eEqn - ( - fvm::ddt(rho, ea) - fvc::ddt(rho, ea) - // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha - - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) - + diffAlphaD - == - fvc::div(hDiffCorrFlux) - ); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoE_rhs = fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) + - diffAlphaD + + fvc::div(hDiffCorrFlux); + + rhoE += rkcoe3[nrk]*rhoE_rhs*runTime.deltaT(); + + ea = rhoE/rho - 0.5*magSqr(U); + ea.correctBoundaryConditions(); + + ha = ea + p/rho; + chemistry->correctThermo(); // before this, we must update ha = e + p/rho + } + else + { + fvScalarMatrix eEqn + ( + fvm::ddt(rho, ea) - fvc::ddt(rho, ea) + // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha + - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) + + diffAlphaD + == + fvc::div(hDiffCorrFlux) + ); - eEqn.solve("ea"); + eEqn.solve("ea"); - ha = ea + p/rho; - chemistry->correctThermo(); - rhoE = rho*(ea + 0.5*magSqr(U)); + ha = ea + p/rho; + chemistry->correctThermo(); + rhoE = rho*(ea + 0.5*magSqr(U)); + } + } Info<< "min/max(T) = " diff --git a/applications/solvers/dfHighSpeedFoam/rhoEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEqn.H index 21be7d5be..d02a21314 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoEqn.H @@ -1,5 +1,20 @@ -solve -( - fvm::ddt(rho) - + fvc::div(phi) -); +if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rho_rhs = -fvc::div(phi); + + rho = rkcoe1[nrk]*rho_save + + rkcoe2[nrk]*rho + + rkcoe3[nrk]*rho_rhs*runTime.deltaT(); + + Info <<"in rk"<< nrk+1 << " finish calculate rho" << nl << endl; +} +else +{ + solve + ( + fvm::ddt(rho) + + fvc::div(phi) + ); +} + + diff --git a/applications/solvers/dfHighSpeedFoam/rhoUEqn.H b/applications/solvers/dfHighSpeedFoam/rhoUEqn.H index f00a6899c..df6903eb9 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoUEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoUEqn.H @@ -1,23 +1,50 @@ -//same with rhoReactingFoam -solve -( - fvm::ddt(rhoU) - + fvc::div(phiUp) -); +if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rhoU_rhs = -fvc::div(phiUp); -U.ref() = rhoU()/rho(); -U.correctBoundaryConditions(); -rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + rhoU = rkcoe1[nrk]*rhoU_save + + rkcoe2[nrk]*rhoU + + rkcoe3[nrk]*rhoU_rhs*runTime.deltaT(); + + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); -if (!inviscid) + if (!inviscid) + { + rhoU_rhs = fvc::laplacian(turbulence->muEff(),U) + fvc::div(tauMC); + rhoU += rkcoe3[nrk]*rhoU_rhs*runTime.deltaT(); + + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + } + + Info <<"in rk"<muEff(), U) - - fvc::div(tauMC) + fvm::ddt(rhoU) + + fvc::div(phiUp) ); - rhoU = rho*U; + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + + if (!inviscid) + { + solve + ( + fvm::ddt(rho, U) - fvc::ddt(rho, U) + - fvm::laplacian(turbulence->muEff(), U) + - fvc::div(tauMC) + ); + rhoU = rho*U; + } + + Info << "\nmin / max mag(U) ; " << gMin(mag(U)()()) << " / " << gMax(mag(U)()()) << nl << endl; } -Info << "\nmin / max mag(U) ; " << gMin(mag(U)()()) << " / " << gMax(mag(U)()()) << nl << endl; + diff --git a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H index 46b0c9a94..26beb31c8 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H @@ -22,13 +22,33 @@ tmp> mvConvection ); { - 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); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + if(chemScheme == "wrate") + { + start = std::clock(); + + chemistry->calculateW(); + + 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); + } + } + else + { + 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]); @@ -44,23 +64,41 @@ tmp> mvConvection } if (i != inertIndex) - { - // original convection term - // fvScalarMatrix YiEqn - // ( - // fvm::ddt(rho, Yi) - // + mvConvection->fvmDiv(phi, Yi) - // == - // combustion->R(Yi) - // ); - - solve - ( - fvm::ddt(rhoYi[i]) - + fvc::div(phiYi[i]) - == - chemistry->RR(i) - ); + { + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoYi_rhs = -fvc::div(phiYi[i]); + + if(chemScheme == "wrate") + { + 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(); + } + else + { + // original convection term + // fvScalarMatrix YiEqn + // ( + // fvm::ddt(rho, Yi) + // + mvConvection->fvmDiv(phi, Yi) + // == + // combustion->R(Yi) + // ); + + solve + ( + fvm::ddt(rhoYi[i]) + + fvc::div(phiYi[i]) + == + chemistry->RR(i) + ); + } Yi=rhoYi[i]/rho; Yi.max(0.0); @@ -70,22 +108,34 @@ tmp> mvConvection const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); tmp DEff = chemistry->rhoD(i) + turbulence->mut()/Sct; - // original term - // YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoYi_rhs = fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phic,Yi)"); - fvScalarMatrix YiEqn - ( - fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi) - - fvm::laplacian(DEff(), Yi) - + mvConvection->fvmDiv(phiUc, Yi) - ); + rhoYi[i] += rkcoe3[nrk]*rhoYi_rhs*runTime.deltaT(); + + Yi = rhoYi[i]/rho; + Yi.max(0.0); + } + else + { + // original term + // YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi); - YiEqn.relax(); + fvScalarMatrix YiEqn + ( + fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi) + - fvm::laplacian(DEff(), Yi) + + mvConvection->fvmDiv(phiUc, Yi) + ); - YiEqn.solve("Yi"); + YiEqn.relax(); - Yi.max(0.0); + YiEqn.solve("Yi"); + Yi.max(0.0); + } + } // original term @@ -93,6 +143,7 @@ tmp> mvConvection // YiEqn.solve("Yi"); // Yi.max(0.0); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) Yi.correctBoundaryConditions(); rhoYi[i] = rho*Yi; Yt += Yi; } @@ -100,6 +151,7 @@ tmp> mvConvection 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); diff --git a/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H b/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H new file mode 100644 index 000000000..2f9e8c33a --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H @@ -0,0 +1,12 @@ +rho_save = rho; + +rhoU_save = rho*U; + +rhoE_save = rho*(ea + 0.5*magSqr(U)); + +forAll(rhoYi_save,i) +{ + rhoYi_save[i] = rho*Y[i]; +} + +Info <<"finish update saved fields"<< endl; \ No newline at end of file From 8a329ba109973a5bc32754db79b04aac02c23df5 Mon Sep 17 00:00:00 2001 From: pkuLmq <1900011062@pku.edu.cn> Date: Fri, 30 Jun 2023 00:33:02 +0800 Subject: [PATCH 2/2] a little fix --- applications/solvers/dfHighSpeedFoam/Make/files | 4 ++-- .../{dfHighSpeedFoam_new.C => dfHighSpeedFoam.C} | 0 2 files changed, 2 insertions(+), 2 deletions(-) rename applications/solvers/dfHighSpeedFoam/{dfHighSpeedFoam_new.C => dfHighSpeedFoam.C} (100%) diff --git a/applications/solvers/dfHighSpeedFoam/Make/files b/applications/solvers/dfHighSpeedFoam/Make/files index 1d93e0527..70642de5e 100644 --- a/applications/solvers/dfHighSpeedFoam/Make/files +++ b/applications/solvers/dfHighSpeedFoam/Make/files @@ -1,3 +1,3 @@ -dfHighSpeedFoam_new.C +dfHighSpeedFoam.C -EXE = $(DF_APPBIN)/dfHighSpeedFoam_new +EXE = $(DF_APPBIN)/dfHighSpeedFoam diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C similarity index 100% rename from applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam_new.C rename to applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C