diff --git a/CMakeLists.txt b/CMakeLists.txt index a9faa7801..d3e7a48c4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4) include(config/petscCompilers.cmake) # Set the project details -project(ablateLibrary VERSION 0.13.01) +project(ablateLibrary VERSION 0.13.02) # Load the Required 3rd Party Libaries pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc) diff --git a/CMakePresets.json b/CMakePresets.json index 6fe2a08e9..723987490 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -34,36 +34,6 @@ "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}" } }, - { - "name": "OptMainPetsc", - "displayName": "CLion Opt Config", - "description": "Default build for ABLATE in CLion", - "binaryDir": "${sourceDir}/cmake-build-MainPetsc", - "cacheVariables": { - "CMAKE_BUILD_TYPE": "Release", - "COMPILE_MPI_COMMAND": "$env{PETSC_DIR}/arch-ablate-opt/bin/mpirun" - }, - "environment": { - "PETSC_DIR": "/home/klbud/Software/PetscUpdating/petsc", - "PETSC_ARCH": "arch-ablate-opt", - "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}" - } - }, - { - "name": "OptFluentPetsc", - "displayName": "CLion Opt Config", - "description": "Default build for ABLATE in CLion", - "binaryDir": "${sourceDir}/cmake-build-FluentPetsc", - "cacheVariables": { - "CMAKE_BUILD_TYPE": "Release", - "COMPILE_MPI_COMMAND": "$env{PETSC_DIR}/arch-ablate-opt/bin/mpirun" - }, - "environment": { - "PETSC_DIR": "/home/klbud/Software/PetscUpdating/petscW", - "PETSC_ARCH": "arch-ablate-opt", - "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}" - } - }, { "name": "local-ablate-opt-info", "displayName": "CLion RelWithDebugInfo Config", diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index e3f814ec7..9f9702f5c 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -99,7 +99,6 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV VecGetArrayRead(locGradVecs[field.subId], &locGradArrays[field.subId]) >> utilities::PetscUtilities::checkError; } } - ComputeFluxSourceTerms(dm, ds, totDim, @@ -499,8 +498,8 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD PetscScalar *auxL = nullptr, *auxR = nullptr; // Precompute the offsets to pass into the rhsFluxFunctionDescriptions - std::vector fluxComponentSize(rhsFunctions.size()); - std::vector fluxId(rhsFunctions.size()); + std::vector> fluxComponentSize(rhsFunctions.size()); + std::vector> fluxId(rhsFunctions.size()); std::vector> uOff(rhsFunctions.size()); std::vector> aOff(rhsFunctions.size()); @@ -509,9 +508,11 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD PetscDSGetComponentOffsets(ds, &uOffTotal) >> utilities::PetscUtilities::checkError; for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { - const auto& field = subDomain->GetField(rhsFunctions[fun].field); - fluxComponentSize[fun] = field.numberComponents; - fluxId[fun] = field.id; + for (std::size_t f = 0; f < rhsFunctions[fun].updateFields.size(); f++) { + const auto& field = subDomain->GetField(rhsFunctions[fun].updateFields[f]); + fluxComponentSize[fun].push_back(field.numberComponents); + fluxId[fun].push_back(field.id); + } for (std::size_t f = 0; f < rhsFunctions[fun].inputFields.size(); f++) { uOff[fun].push_back(uOffTotal[rhsFunctions[fun].inputFields[f]]); } @@ -534,7 +535,6 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD DMLabel regionLabel = nullptr; PetscInt regionValue = 0; domain::Region::GetLabel(solverRegion, subDomain->GetDM(), regionLabel, regionValue); - // March over each face in this region for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { const PetscInt face = faceRange.points ? faceRange.points[f] : f; @@ -574,33 +574,36 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD // March over each source function for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { + PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) PetscArrayzero(flux, totDim) >> utilities::PetscUtilities::checkError; const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; rhsFluxFunctionDescription.function(dim, fg, uOff[fun].data(), uL, uR, aOff[fun].data(), auxL, auxR, flux, rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; + // add the fluxes back to the cell + for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { + PetscInt cellLabelValue = regionValue; + PetscScalar *fL = nullptr, *fR = nullptr; + DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[0], fluxId[fun][updateFieldIdx], locFArray, &fL) >> utilities::PetscUtilities::checkError; + } - // add the flux back to the cell - PetscScalar *fL = nullptr, *fR = nullptr; - PetscInt cellLabelValue = regionValue; - DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[0], fluxId[fun], locFArray, &fL) >> utilities::PetscUtilities::checkError; - } - - cellLabelValue = regionValue; - DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[1], fluxId[fun], locFArray, &fR) >> utilities::PetscUtilities::checkError; - } + cellLabelValue = regionValue; + DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[1], fluxId[fun][updateFieldIdx], locFArray, &fR) >> utilities::PetscUtilities::checkError; + } - for (PetscInt d = 0; d < fluxComponentSize[fun]; ++d) { - if (fL) fL[d] -= flux[d] / cgL->volume; - if (fR) fR[d] += flux[d] / cgR->volume; + for (PetscInt d = 0; d < (fluxComponentSize[fun][updateFieldIdx]); ++d) { + if (fL) fL[d] -= flux[fluxOffset + d] / cgL->volume; + if (fR) fR[d] += flux[fluxOffset + d] / cgR->volume; + } + fluxOffset += fluxComponentSize[fun][updateFieldIdx]; } } } diff --git a/src/finiteVolume/cellInterpolant.hpp b/src/finiteVolume/cellInterpolant.hpp index 97fe51c64..6bfbdc7c9 100644 --- a/src/finiteVolume/cellInterpolant.hpp +++ b/src/finiteVolume/cellInterpolant.hpp @@ -26,7 +26,7 @@ class CellInterpolant { DiscontinuousFluxFunction function; void* context; - PetscInt field; + std::vector updateFields; std::vector inputFields; std::vector auxFields; }; diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index 45125ab65..f5607d6a4 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -1,37 +1,55 @@ #include "compressibleFlowSolver.hpp" #include #include "compressibleFlowFields.hpp" +#include "finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.hpp" +#include "finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp" #include "finiteVolume/processes/evTransport.hpp" #include "finiteVolume/processes/navierStokesTransport.hpp" #include "finiteVolume/processes/speciesTransport.hpp" #include "utilities/vectorUtilities.hpp" +/** + * The compact argument in this constructor is used to define whether we use seperate processes for the Euler (NS) transport, species transport, and EV transport. + * a value of 0, is all seperate + * a value of 1, is combined Euler and species transport, with seperate EV transport + * a value of 2, is combined Euler and species and The single field Progress EV transport (soot number density in this case) + */ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, const std::shared_ptr& eosIn, const std::shared_ptr& parameters, const std::shared_ptr& transport, const std::shared_ptr& fluxCalculatorIn, std::vector> additionalProcesses, std::vector> boundaryConditions, - const std::shared_ptr& evTransport) - : FiniteVolumeSolver(std::move(solverId), std::move(region), std::move(options), - utilities::VectorUtilities::Merge( - { - // create assumed processes for compressible flow - std::make_shared( - parameters, eosIn, fluxCalculatorIn, transport, utilities::VectorUtilities::Find(additionalProcesses)), - std::make_shared(eosIn, fluxCalculatorIn, transport, parameters), - std::make_shared(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport), - }, - additionalProcesses), - std::move(boundaryConditions)) {} + const std::shared_ptr& evTransport, int compact) + : FiniteVolumeSolver( + std::move(solverId), std::move(region), std::move(options), + compact ? (compact == 1 ? utilities::VectorUtilities::Merge({std::make_shared( + parameters, eosIn, fluxCalculatorIn, transport, + utilities::VectorUtilities::Find(additionalProcesses)), + std::make_shared(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport)}, + additionalProcesses) + : utilities::VectorUtilities::Merge({std::make_shared( + parameters, eosIn, fluxCalculatorIn, transport, evTransport ? evTransport : transport, + utilities::VectorUtilities::Find(additionalProcesses))}, + additionalProcesses)) + : utilities::VectorUtilities::Merge( + { + // create assumed processes for compressible flow + std::make_shared( + parameters, eosIn, fluxCalculatorIn, transport, utilities::VectorUtilities::Find(additionalProcesses)), + std::make_shared(eosIn, fluxCalculatorIn, transport, parameters), + std::make_shared(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport), + }, + additionalProcesses), + std::move(boundaryConditions)) {} ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, const std::shared_ptr& eosIn, const std::shared_ptr& parameters, const std::shared_ptr& transport, const std::shared_ptr& fluxCalculatorIn, std::vector> boundaryConditions, - const std::shared_ptr& evTransport) - : CompressibleFlowSolver(std::move(solverId), std::move(region), std::move(options), eosIn, parameters, transport, fluxCalculatorIn, {}, std::move(boundaryConditions), evTransport) {} + const std::shared_ptr& evTransport, int compact) + : CompressibleFlowSolver(std::move(solverId), std::move(region), std::move(options), eosIn, parameters, transport, fluxCalculatorIn, {}, std::move(boundaryConditions), evTransport, compact) {} #include "registrar.hpp" REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, "compressible finite volume flow", ARG(std::string, "id", "the name of the flow field"), @@ -41,4 +59,5 @@ REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, " OPT(ablate::finiteVolume::fluxCalculator::FluxCalculator, "fluxCalculator", "the flux calculators (defaults to none)"), OPT(std::vector, "additionalProcesses", "any additional processes besides euler/yi/ev transport"), OPT(std::vector, "boundaryConditions", "the boundary conditions for the flow field"), - OPT(ablate::eos::transport::TransportModel, "evTransport", "when provided, this model will be used for ev transport instead of default")); \ No newline at end of file + OPT(ablate::eos::transport::TransportModel, "evTransport", "when provided, this model will be used for ev transport instead of default"), + OPT(int, "compact", "Integer value describing whether to treat all the transport seperately, partially combined, or fully combined (see commented code above constructor for values)")); \ No newline at end of file diff --git a/src/finiteVolume/compressibleFlowSolver.hpp b/src/finiteVolume/compressibleFlowSolver.hpp index 9c378637e..696e70765 100644 --- a/src/finiteVolume/compressibleFlowSolver.hpp +++ b/src/finiteVolume/compressibleFlowSolver.hpp @@ -31,7 +31,8 @@ class CompressibleFlowSolver : public FiniteVolumeSolver { CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, const std::shared_ptr& eos, const std::shared_ptr& parameters, const std::shared_ptr& transport, const std::shared_ptr& = {}, std::vector> additionalProcesses = {}, - std::vector> boundaryConditions = {}, const std::shared_ptr& evTransport = {}); + std::vector> boundaryConditions = {}, const std::shared_ptr& evTransport = {}, + int compact = 0); /** * Constructor without ev or additional processes @@ -48,7 +49,7 @@ class CompressibleFlowSolver : public FiniteVolumeSolver { CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, const std::shared_ptr& eos, const std::shared_ptr& parameters, const std::shared_ptr& transport, const std::shared_ptr& = {}, std::vector> boundaryConditions = {}, - const std::shared_ptr& evTransport = {}); + const std::shared_ptr& evTransport = {}, int compact = 0); ~CompressibleFlowSolver() override = default; }; } // namespace ablate::finiteVolume diff --git a/src/finiteVolume/faceInterpolant.cpp b/src/finiteVolume/faceInterpolant.cpp index 7be985a2b..b9680e25c 100644 --- a/src/finiteVolume/faceInterpolant.cpp +++ b/src/finiteVolume/faceInterpolant.cpp @@ -314,8 +314,8 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV ablate::domain::Region::GetLabel(solverRegion, subDomain->GetDM(), regionLabel, regionValue); // Precompute the offsets to pass into the rhsFluxFunctionDescriptions - std::vector fluxComponentSize(rhsFunctions.size()); - std::vector fluxId(rhsFunctions.size()); + std::vector> fluxComponentSize(rhsFunctions.size()); + std::vector> fluxId(rhsFunctions.size()); std::vector> uOff(rhsFunctions.size()); std::vector> aOff(rhsFunctions.size()); std::vector> uOff_x(rhsFunctions.size()); @@ -328,9 +328,11 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV PetscDSGetComponentDerivativeOffsets(subDomain->GetDiscreteSystem(), &uGradOffTotal) >> utilities::PetscUtilities::checkError; for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { - const auto& field = subDomain->GetField(rhsFunctions[fun].field); - fluxComponentSize[fun] = field.numberComponents; - fluxId[fun] = field.id; + for (std::size_t f = 0; f < rhsFunctions[fun].updateFields.size(); f++) { + const auto& field = subDomain->GetField(rhsFunctions[fun].updateFields[f]); + fluxComponentSize[fun].push_back(field.numberComponents); + fluxId[fun].push_back(field.id); + } for (std::size_t f = 0; f < rhsFunctions[fun].inputFields.size(); f++) { uOff[fun].push_back(uOffTotal[rhsFunctions[fun].inputFields[f]]); uOff_x[fun].push_back(uGradOffTotal[rhsFunctions[fun].inputFields[f]]); @@ -387,7 +389,7 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV // March over each source function for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { PetscArrayzero(flux.data(), totDim) >> utilities::PetscUtilities::checkError; - + PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; rhsFluxFunctionDescription.function(dim, fg, @@ -402,30 +404,32 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV flux.data(), rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; + for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { + // add the flux back to the cell + PetscScalar *fL = nullptr, *fR = nullptr; + PetscInt cellLabelValue = regionValue; + DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[0], rhsFunctions[fun].updateFields[updateFieldIdx], locFArray, &fL) >> utilities::PetscUtilities::checkError; + } - // add the flux back to the cell - PetscScalar *fL = nullptr, *fR = nullptr; - PetscInt cellLabelValue = regionValue; - DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[0], rhsFunctions[fun].field, locFArray, &fL) >> utilities::PetscUtilities::checkError; - } - - cellLabelValue = regionValue; - DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[1], rhsFunctions[fun].field, locFArray, &fR) >> utilities::PetscUtilities::checkError; - } + cellLabelValue = regionValue; + DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[1], rhsFunctions[fun].updateFields[updateFieldIdx], locFArray, &fR) >> utilities::PetscUtilities::checkError; + } - for (PetscInt d = 0; d < fluxComponentSize[fun]; ++d) { - if (fL) fL[d] -= flux[d] / cgL->volume; - if (fR) fR[d] += flux[d] / cgR->volume; + for (PetscInt d = 0; d < fluxComponentSize[fun][updateFieldIdx]; ++d) { + if (fL) fL[d] -= flux[d + fluxOffset] / cgL->volume; + if (fR) fR[d] += flux[d + fluxOffset] / cgR->volume; + } + fluxOffset += fluxComponentSize[fun][updateFieldIdx]; } } } diff --git a/src/finiteVolume/faceInterpolant.hpp b/src/finiteVolume/faceInterpolant.hpp index e2637bdc3..0a0cabd9c 100644 --- a/src/finiteVolume/faceInterpolant.hpp +++ b/src/finiteVolume/faceInterpolant.hpp @@ -85,7 +85,7 @@ class FaceInterpolant { ContinuousFluxFunction function; void* context; - PetscInt field; + std::vector updateFields; std::vector inputFields; std::vector auxFields; }; diff --git a/src/finiteVolume/finiteVolumeSolver.cpp b/src/finiteVolume/finiteVolumeSolver.cpp index 80d7d74e4..9883e0e42 100644 --- a/src/finiteVolume/finiteVolumeSolver.cpp +++ b/src/finiteVolume/finiteVolumeSolver.cpp @@ -290,13 +290,15 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets PetscFunctionReturn(0); } -void ablate::finiteVolume::FiniteVolumeSolver::RegisterRHSFunction(CellInterpolant::DiscontinuousFluxFunction function, void* context, const std::string& field, +void ablate::finiteVolume::FiniteVolumeSolver::RegisterRHSFunction(CellInterpolant::DiscontinuousFluxFunction function, void* context, const std::vector& fields, const std::vector& inputFields, const std::vector& auxFields) { - // map the field, inputFields, and auxFields to locations - auto& fieldId = subDomain->GetField(field); + CellInterpolant::DiscontinuousFluxFunctionDescription functionDescription{.function = function, .context = context}; - // Create the FVMRHS Function - CellInterpolant::DiscontinuousFluxFunctionDescription functionDescription{.function = function, .context = context, .field = fieldId.id}; + // map the field, inputFields, and auxFields to locations + for (auto& field : fields) { + auto& fieldId = subDomain->GetField(field); + functionDescription.updateFields.push_back(fieldId.id); + } for (auto& inputField : inputFields) { auto& inputFieldId = subDomain->GetField(inputField); @@ -311,13 +313,16 @@ void ablate::finiteVolume::FiniteVolumeSolver::RegisterRHSFunction(CellInterpola discontinuousFluxFunctionDescriptions.push_back(functionDescription); } -void ablate::finiteVolume::FiniteVolumeSolver::RegisterRHSFunction(ablate::finiteVolume::FaceInterpolant::ContinuousFluxFunction function, void* context, const std::string& field, +void ablate::finiteVolume::FiniteVolumeSolver::RegisterRHSFunction(ablate::finiteVolume::FaceInterpolant::ContinuousFluxFunction function, void* context, const std::vector& updateFields, const std::vector& inputFields, const std::vector& auxFields) { // map the field, inputFields, and auxFields to locations - auto& fieldId = subDomain->GetField(field); + FaceInterpolant::ContinuousFluxFunctionDescription functionDescription{.function = function, .context = context}; + for (auto& field : updateFields) { + auto& fieldId = subDomain->GetField(field); + functionDescription.updateFields.push_back(fieldId.id); + } // Create the FVMRHS Function - FaceInterpolant::ContinuousFluxFunctionDescription functionDescription{.function = function, .context = context, .field = fieldId.id}; for (auto& inputField : inputFields) { auto& inputFieldId = subDomain->GetField(inputField); diff --git a/src/finiteVolume/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index bfa46ab3d..fd40dc2cf 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -111,14 +111,14 @@ class FiniteVolumeSolver : public solver::CellSolver, PetscErrorCode ComputeBoundary(PetscReal time, Vec locX, Vec locX_t) override; /** - * Register a FVM rhs discontinuous flux function + * Register a FVM rhs discontinuous flux function for multiple fields * @param function * @param context * @param field * @param inputFields * @param auxFields */ - void RegisterRHSFunction(CellInterpolant::DiscontinuousFluxFunction function, void* context, const std::string& field, const std::vector& inputFields, + void RegisterRHSFunction(CellInterpolant::DiscontinuousFluxFunction function, void* context, const std::vector& field, const std::vector& inputFields, const std::vector& auxFields); /** @@ -129,7 +129,7 @@ class FiniteVolumeSolver : public solver::CellSolver, * @param inputFields * @param auxFields */ - void RegisterRHSFunction(FaceInterpolant::ContinuousFluxFunction function, void* context, const std::string& field, const std::vector& inputFields, + void RegisterRHSFunction(FaceInterpolant::ContinuousFluxFunction function, void* context, const std::vector& fields, const std::vector& inputFields, const std::vector& auxFields); /** diff --git a/src/finiteVolume/processes/CMakeLists.txt b/src/finiteVolume/processes/CMakeLists.txt index b84779afd..f2a6239c4 100644 --- a/src/finiteVolume/processes/CMakeLists.txt +++ b/src/finiteVolume/processes/CMakeLists.txt @@ -16,6 +16,8 @@ target_sources(ablateLibrary thermophoreticDiffusion.cpp surfaceForce.cpp soot.cpp + compactCompressibleNSSpeciesTransport.cpp + compactCompressibleNSSpeciesSingleProgressTransport.cpp PUBLIC process.hpp @@ -35,4 +37,6 @@ target_sources(ablateLibrary thermophoreticDiffusion.hpp surfaceForce.hpp soot.hpp + compactCompressibleNSSpeciesTransport.hpp + compactCompressibleNSSpeciesSingleProgressTransport.hpp ) diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.cpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.cpp new file mode 100644 index 000000000..4a9fb3d3f --- /dev/null +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.cpp @@ -0,0 +1,929 @@ +#include "compactCompressibleNSSpeciesSingleProgressTransport.hpp" +#include +#include "finiteVolume/compressibleFlowFields.hpp" +#include "finiteVolume/fluxCalculator/ausm.hpp" +#include "finiteVolume/processes/evTransport.hpp" +#include "finiteVolume/processes/navierStokesTransport.hpp" +#include "finiteVolume/processes/speciesTransport.hpp" +#include "parameters/emptyParameters.hpp" +#include "utilities/constants.hpp" +#include "utilities/mathUtilities.hpp" +#include "utilities/petscUtilities.hpp" + +ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::CompactCompressibleNSSpeciesSingleProgressTransport( + const std::shared_ptr ¶metersIn, std::shared_ptr eosIn, std::shared_ptr fluxCalcIn, + std::shared_ptr baseTransport, std::shared_ptr evTransport, + std::shared_ptr pgs) + : advectionData(), fluxCalculator(fluxCalcIn), eos(std::move(eosIn)), transportModel(std::move(baseTransport)), SingleProgressTransportModel(std::move(evTransport)) { + auto parameters = ablate::parameters::EmptyParameters::Check(parametersIn); + + if (fluxCalculator) { + // cfl + advectionData.cfl = parameters->Get("cfl", 0.5); + + // extract the difference function from fluxDifferencer object + advectionData.fluxCalculatorFunction = fluxCalculator->GetFluxCalculatorFunction(); + advectionData.fluxCalculatorCtx = fluxCalculator->GetFluxCalculatorContext(); + + advectionData.numberSpecies = (PetscInt)eos->GetSpeciesVariables().size(); + } + + timeStepData.advectionData = &advectionData; + timeStepData.pgs = std::move(pgs); + + if (transportModel) { + // Add in the time stepping + diffusionTimeStepData.conductionStabilityFactor = parameters->Get("conductionStabilityFactor", 0.0); + diffusionTimeStepData.viscousStabilityFactor = parameters->Get("viscousStabilityFactor", 0.0); + diffusionTimeStepData.diffusiveStabilityFactor = parameters->Get("speciesStabilityFactor", 0.0); + diffusionTimeStepData.speciesDiffusionCoefficient.resize(eos->GetSpeciesVariables().size()); + + diffusionData.numberSpecies = (PetscInt)eos->GetSpeciesVariables().size(); + diffusionData.numberEV = 1; + diffusionData.speciesSpeciesSensibleEnthalpy.resize(eos->GetSpeciesVariables().size()); + diffusionData.speciesDiffusionCoefficient.resize(eos->GetSpeciesVariables().size()); + } +} + +void ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::Setup(ablate::finiteVolume::FiniteVolumeSolver &flow) { + // loop over convserved EV's and pull out SingleProgress, if there's someething there other than it, error out. + const auto &evConservedFields = flow.GetSubDomain().GetFields(domain::FieldLocation::SOL, CompressibleFlowFields::EV_TAG); + for (auto &evConservedField : evConservedFields) { + if (evConservedField.name != CompressibleFlowFields::DENSITY_PROGRESS_FIELD) { + throw std::invalid_argument("The ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport process expects only the progress extra variable."); + } + // set up the progress Extra Variables and the euler/species transport + auto nonConservedName = evConservedField.name.substr(CompressibleFlowFields::CONSERVED.length()); // non Conserved form just removes density prefix + if (!flow.GetSubDomain().ContainsField(nonConservedName)) { + throw std::invalid_argument("The ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport process expects the conserved (" + evConservedField.name + + ") and non-conserved (" + nonConservedName + ") extra variables to be in the flow."); + } + // Register the euler,eulerYi, and species source terms + if (fluxCalculator) { + // I don't know why we wouldn't push through the old temperature fields, maybe slower for perfect gas/idealized gas's but when there is a temperature iterative method this should be better + // If it is worse for perfect gas's, going to need to add in an option switch -klb + flow.RegisterRHSFunction(AdvectionFlux, + &advectionData, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD, evConservedField.name}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD, evConservedField.name}, + {CompressibleFlowFields::TEMPERATURE_FIELD}); + + // Set the ComputeCFLTimestepFrom flow Process through + flow.RegisterComputeTimeStepFunction(ComputeCflTimeStep, &timeStepData, "cfl"); + + advectionData.numberEV = evConservedField.numberComponents; + advectionData.computeTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + advectionData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, flow.GetSubDomain().GetFields()); + advectionData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); + advectionData.computePressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); + } + + // if there are any coefficients for diffusion, compute diffusion, for here we will follow suit in allowing multiple diffusion functions + if (transportModel) { + // Store the required data for the low level c functions + diffusionData.muFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Viscosity, flow.GetSubDomain().GetFields()); + diffusionData.kFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Conductivity, flow.GetSubDomain().GetFields()); + diffusionData.diffFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Diffusivity, flow.GetSubDomain().GetFields()); + diffusionData.computeSpeciesSensibleEnthalpyFunction = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeciesSensibleEnthalpy, flow.GetSubDomain().GetFields()); + if (SingleProgressTransportModel) { + diffusionData.numberEV = evConservedField.numberComponents; + diffusionData.evDiffusionCoefficient.resize(diffusionData.numberEV); + diffusionData.evDiffFunction = SingleProgressTransportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Diffusivity, flow.GetSubDomain().GetFields()); + // ExtraVariable + if (diffusionData.evDiffFunction.function) { + if (diffusionData.evDiffFunction.propertySize == 1) { + flow.RegisterRHSFunction( + DiffusionEVFlux, &diffusionData, {evConservedField.name}, {CompressibleFlowFields::EULER_FIELD}, {nonConservedName, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else if (diffusionData.evDiffFunction.propertySize == diffusionData.numberEV) { + flow.RegisterRHSFunction(DiffusionEVFluxVariableDiffusionCoefficient, + &diffusionData, + {evConservedField.name}, + {CompressibleFlowFields::EULER_FIELD}, + {nonConservedName, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else + throw std::invalid_argument( + "The ev diffusion property size must be 1 or number of ev in ablate::finiteVolume::processes::CompatcCompressibleNSSpeciesSingleProgressTransport!"); + } + } + + /* For now we will do 3 different diffusive flux calls just since it doesn't seem like their splitting costs too much time in redundant calculations */ + if (diffusionData.muFunction.function || diffusionData.kFunction.function) { + // Register the Diffusion Source term + flow.RegisterRHSFunction(DiffusionFlux, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::TEMPERATURE_FIELD, CompressibleFlowFields::VELOCITY_FIELD}); + } + // Species + if (diffusionData.diffFunction.function) { + // Specify a different rhs function depending on if the diffusion flux is constant + if (diffusionData.diffFunction.propertySize == 1) { + flow.RegisterRHSFunction(DiffusionEnergyFlux, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + flow.RegisterRHSFunction(DiffusionSpeciesFlux, + &diffusionData, + {CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else if (diffusionData.diffFunction.propertySize == diffusionData.numberSpecies) { + flow.RegisterRHSFunction(DiffusionEnergyFluxVariableDiffusionCoefficient, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + flow.RegisterRHSFunction(DiffusionSpeciesFluxVariableDiffusionCoefficient, + &diffusionData, + {CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else { + throw std::invalid_argument("The diffusion property size must be 1 or number of species in ablate::finiteVolume::processes::SpeciesTransport."); + } + } + + // Check to see if time step calculations should be added for viscosity or conduction + if (diffusionTimeStepData.conductionStabilityFactor > 0 || diffusionTimeStepData.viscousStabilityFactor > 0) { + diffusionTimeStepData.kFunction = diffusionData.kFunction; + diffusionTimeStepData.muFunction = diffusionData.muFunction; + diffusionTimeStepData.specificHeat = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpecificHeatConstantVolume, flow.GetSubDomain().GetFields()); + diffusionTimeStepData.density = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Density, flow.GetSubDomain().GetFields()); + + diffusionTimeStepData.numberSpecies = eos->GetSpeciesVariables().size(); + + if (diffusionTimeStepData.conductionStabilityFactor > 0) { + flow.RegisterComputeTimeStepFunction(ComputeConductionTimeStep, &diffusionTimeStepData, "cond"); + } + if (diffusionTimeStepData.viscousStabilityFactor > 0) { + flow.RegisterComputeTimeStepFunction(ComputeViscousDiffusionTimeStep, &diffusionTimeStepData, "visc"); + } + } + if (diffusionTimeStepData.diffusiveStabilityFactor > 0) { + diffusionTimeStepData.numberSpecies = diffusionData.numberSpecies; + diffusionTimeStepData.diffFunction = diffusionData.diffFunction; + flow.RegisterComputeTimeStepFunction(ComputeViscousSpeciesDiffusionTimeStep, &diffusionTimeStepData, "spec"); + } + } + + // Setup up aux updates and normalizations + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::VELOCITY_FIELD)) { + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxVelocityField, + nullptr, + std::vector{CompressibleFlowFields::VELOCITY_FIELD}, + {CompressibleFlowFields::EULER_FIELD}); + } + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::TEMPERATURE_FIELD)) { + computeTemperatureFunction = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxTemperatureField, + &computeTemperatureFunction, + std::vector{CompressibleFlowFields::TEMPERATURE_FIELD}, + {}); + } + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::PRESSURE_FIELD)) { + computePressureFunction = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); + flow.RegisterAuxFieldUpdate( + ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxPressureField, &computePressureFunction, std::vector{CompressibleFlowFields::PRESSURE_FIELD}, {}); + } + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::SpeciesTransport::UpdateAuxMassFractionField, + &advectionData.numberSpecies, + std::vector{CompressibleFlowFields::YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}); + // clean up the species + flow.RegisterPostEvaluate(ablate::finiteVolume::processes::SpeciesTransport::NormalizeSpecies); + // limit the EV's + if (evConservedField.Tagged(CompressibleFlowFields::PositiveRange)) { + const auto &conservedFieldName = evConservedField.name; + flow.RegisterPostEvaluate( + [conservedFieldName](TS ts, ablate::solver::Solver &solver) { ablate::finiteVolume::processes::EVTransport::PositiveExtraVariables(ts, solver, conservedFieldName); }); + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::EVTransport::UpdatePositiveEVField, + &advectionData.numberEV, + std::vector{nonConservedName}, + {CompressibleFlowFields::EULER_FIELD, evConservedField.name}); + } else if (evConservedField.Tagged(CompressibleFlowFields::BoundRange)) { + const auto &conservedFieldName = evConservedField.name; + flow.RegisterPostEvaluate( + [conservedFieldName](TS ts, ablate::solver::Solver &solver) { ablate::finiteVolume::processes::EVTransport::BoundExtraVariables(ts, solver, conservedFieldName); }); + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::EVTransport::UpdateBoundEVField, + &advectionData.numberEV, + std::vector{nonConservedName}, + {CompressibleFlowFields::EULER_FIELD, evConservedField.name}); + } else if (evConservedField.Tagged(CompressibleFlowFields::MinusOneToOneRange)) { + const auto &conservedFieldName = evConservedField.name; + flow.RegisterPostEvaluate( + [conservedFieldName](TS ts, ablate::solver::Solver &solver) { ablate::finiteVolume::processes::EVTransport::BoundExtraVariablesMinusOneToOne(ts, solver, conservedFieldName); }); + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::EVTransport::UpdateMinusOneToOneBoundEVField, + &advectionData.numberEV, + std::vector{nonConservedName}, + {CompressibleFlowFields::EULER_FIELD, evConservedField.name}); + } else { + // Allow for the entire range + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::EVTransport::UpdateEVField, + &advectionData.numberEV, + std::vector{nonConservedName}, + {CompressibleFlowFields::EULER_FIELD, evConservedField.name}); + } + } // End EV loop +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::AdvectionFlux(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt *uOff, + const PetscScalar *fieldL, const PetscScalar *fieldR, const PetscInt *aOff, + const PetscScalar *auxL, const PetscScalar *auxR, PetscScalar *flux, void *ctx) { + PetscFunctionBeginUser; + + auto advectionData = (AdvectionData *)ctx; + + const int EULER_FIELD = 0; + const int RHOYI_FIELD = 1; + const int RHOEV_FIELD = 2; + + // Compute the norm + PetscReal norm[3]; + utilities::MathUtilities::NormVector(dim, fg->normal, norm); + const PetscReal areaMag = utilities::MathUtilities::MagVector(dim, fg->normal); + + // Decode the left and right states + PetscReal densityL; + PetscReal normalVelocityL; + PetscReal velocityL[3]; + PetscReal internalEnergyL; + PetscReal aL; + PetscReal pL; + + // decode the left side + { + densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + PetscReal temperatureL; + + PetscCall(advectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + .33 * auxR[aOff[0]], &temperatureL, advectionData->computeTemperature.context.get())); + + // Get the velocity in this direction + normalVelocityL = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velocityL[d] = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityL; + normalVelocityL += velocityL[d] * norm[d]; + } + PetscCall(advectionData->computeInternalEnergy.function(fieldL, temperatureL, &internalEnergyL, advectionData->computeInternalEnergy.context.get())); + PetscCall(advectionData->computeSpeedOfSound.function(fieldL, temperatureL, &aL, advectionData->computeSpeedOfSound.context.get())); + PetscCall(advectionData->computePressure.function(fieldL, temperatureL, &pL, advectionData->computePressure.context.get())); + } + + PetscReal densityR; + PetscReal normalVelocityR; + PetscReal velocityR[3]; + PetscReal internalEnergyR; + PetscReal aR; + PetscReal pR; + + { // decode right state + densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + PetscReal temperatureR; + + PetscCall(advectionData->computeTemperature.function(fieldR, auxR[aOff[0]] * .67 + .33 * auxL[aOff[0]], &temperatureR, advectionData->computeTemperature.context.get())); + + // Get the velocity in this direction + normalVelocityR = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velocityR[d] = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityR; + normalVelocityR += velocityR[d] * norm[d]; + } + + PetscCall(advectionData->computeInternalEnergy.function(fieldR, temperatureR, &internalEnergyR, advectionData->computeInternalEnergy.context.get())); + PetscCall(advectionData->computeSpeedOfSound.function(fieldR, temperatureR, &aR, advectionData->computeSpeedOfSound.context.get())); + PetscCall(advectionData->computePressure.function(fieldR, temperatureR, &pR, advectionData->computePressure.context.get())); + } + + // get the face values + PetscReal massFlux; + PetscReal p12; + + fluxCalculator::Direction direction = + advectionData->fluxCalculatorFunction(advectionData->fluxCalculatorCtx, normalVelocityL, aL, densityL, pL, normalVelocityR, aR, densityR, pR, &massFlux, &p12); + + if (direction == fluxCalculator::LEFT) { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); + PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; + flux[CompressibleFlowFields::RHOE] = HL * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = velocityL[n] * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldL[uOff[RHOYI_FIELD] + ns] / densityL * areaMag; + // EV's + for (PetscInt ev = 0; ev < advectionData->numberEV; ev++) flux[uOff[RHOEV_FIELD] + ev] = (massFlux * fieldL[uOff[RHOEV_FIELD] + ev] / densityL) * areaMag; + } else if (direction == fluxCalculator::RIGHT) { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); + PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; + flux[CompressibleFlowFields::RHOE] = HR * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = velocityR[n] * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldR[uOff[RHOYI_FIELD] + ns] / densityR * areaMag; + // EV's + for (PetscInt ev = 0; ev < advectionData->numberEV; ev++) flux[uOff[RHOEV_FIELD] + ev] = (massFlux * fieldR[uOff[RHOEV_FIELD] + ev] / densityR) * areaMag; + } else { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + + PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); + PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; + + PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); + PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; + + flux[CompressibleFlowFields::RHOE] = 0.5 * (HL + HR) * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = 0.5 * (velocityL[n] + velocityR[n]) * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) + flux[uOff[RHOYI_FIELD] + ns] = massFlux * 0.5 * (fieldR[uOff[RHOYI_FIELD] + ns] + fieldL[uOff[RHOYI_FIELD] + ns]) / (0.5 * (densityL + densityR)) * areaMag; + // EV's + for (PetscInt ev = 0; ev < advectionData->numberEV; ev++) flux[uOff[RHOEV_FIELD] + ev] = (massFlux * fieldR[uOff[RHOEV_FIELD] + ev] / densityR) * areaMag; + // flux[uOff[RHOEV_FIELD]+ev] = massFlux * 0.5 * (fieldR[uOff[RHOEV_FIELD] + ev] + fieldL[uOff[RHOEV_FIELD] + ev] ) / ( 0.5 * (densityR + densityL)) * areaMag; + } + + PetscFunctionReturn(0); +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver &flow, void *ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto timeStepData = (CflTimeStepData *)ctx; + auto advectionData = timeStepData->advectionData; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar *locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + const PetscScalar *x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for euler and densityYi + auto eulerId = flow.GetSubDomain().GetField("euler").id; + + // Get alpha if provided + PetscReal pgsAlpha = 1.0; + if (timeStepData->pgs) { + pgsAlpha = timeStepData->pgs->GetAlpha(); + } + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal *euler; + const PetscReal *conserved = NULL; + const PetscReal *cellCharacteristics = NULL; + DMPlexPointGlobalFieldRead(dm, cell, eulerId, x, &euler) >> utilities::PetscUtilities::checkError; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (euler) { // must be real cell and not ghost + PetscReal rho = euler[CompressibleFlowFields::RHO]; + + // Get the speed of sound from the eos + // TODO:: Replace this with a better temperature guess (see compute conduction Time Step below) + PetscReal temperature; + advectionData->computeTemperature.function(conserved, 300, &temperature, advectionData->computeTemperature.context.get()) >> utilities::PetscUtilities::checkError; + PetscReal a; + advectionData->computeSpeedOfSound.function(conserved, temperature, &a, advectionData->computeSpeedOfSound.context.get()) >> utilities::PetscUtilities::checkError; + + PetscReal dx = 2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]; + + PetscReal velSum = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velSum += PetscAbsReal(euler[CompressibleFlowFields::RHOU + d]) / rho; + } + PetscReal dt = advectionData->cfl * dx / (a / pgsAlpha + velSum); + + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + return dtMin; +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver &flow, void *ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData *)ctx; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar *locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar *x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar *aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto kFunction = diffusionData->kFunction.function; + auto kFunctionContext = diffusionData->kFunction.context.get(); + auto cvFunction = diffusionData->specificHeat.function; + auto cvFunctionContext = diffusionData->specificHeat.context.get(); + auto density = diffusionData->density.function; + auto densityContext = diffusionData->density.context.get(); + auto stabFactor = diffusionData->conductionStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal *conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal *temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + const PetscReal *cellCharacteristics = NULL; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal k; + kFunction(conserved, *temperature, &k, kFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal cv; + cvFunction(conserved, *temperature, &cv, cvFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal rho; + density(conserved, *temperature, &rho, densityContext) >> utilities::PetscUtilities::checkError; + + // Compute alpha + PetscReal alpha = k / (rho * cv); + + PetscReal dx2 = PetscSqr(2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]); + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / alpha); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + + return dtMin; +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver &flow, void *ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData *)ctx; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar *locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar *x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar *aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto muFunction = diffusionData->muFunction.function; + auto muFunctionContext = diffusionData->muFunction.context.get(); + auto density = diffusionData->density.function; + auto densityContext = diffusionData->density.context.get(); + auto stabFactor = diffusionData->viscousStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal *conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal *temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + const PetscReal *cellCharacteristics = NULL; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal mu; + muFunction(conserved, *temperature, &mu, muFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal rho; + density(conserved, *temperature, &rho, densityContext) >> utilities::PetscUtilities::checkError; + + // Compute nu + PetscReal nu = mu / rho; + + PetscReal dx2 = PetscSqr(2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]); + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / nu); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + return dtMin; +} +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::ComputeViscousSpeciesDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver &flow, void *ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData *)ctx; + + // Get the fv geom + PetscReal minCellRadius; + DMPlexGetGeometryFVM(dm, NULL, NULL, &minCellRadius) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar *x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar *aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // assume the smallest cell is the limiting factor for now + const PetscReal dx2 = PetscSqr(2.0 * minCellRadius); + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto diffFunction = diffusionData->diffFunction.function; + auto diffFunctionContext = diffusionData->diffFunction.context.get(); + auto stabFactor = diffusionData->diffusiveStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal *conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal *temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal diff; + diffFunction(conserved, *temperature, &diff, diffFunctionContext) >> utilities::PetscUtilities::checkError; + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / diff); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + return dtMin; +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionFlux(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int T = 0; + const int VEL = 1; + + auto flowParameters = (DiffusionData *)ctx; + + // Compute mu and k + PetscReal mu = 0.0; + flowParameters->muFunction.function(field, aux[aOff[T]], &mu, flowParameters->muFunction.context.get()); + PetscReal k = 0.0; + flowParameters->kFunction.function(field, aux[aOff[T]], &k, flowParameters->kFunction.context.get()); + + // Compute the stress tensor tau + PetscReal tau[9]; // Maximum size without symmetry + PetscCall(ablate::finiteVolume::processes::NavierStokesTransport::CompressibleFlowComputeStressTensor(dim, mu, gradAux + aOff_x[VEL], tau)); + + // for each velocity component + for (PetscInt c = 0; c < dim; ++c) { + PetscReal viscousFlux = 0.0; + + // March over each direction + for (PetscInt d = 0; d < dim; ++d) { + viscousFlux += -fg->normal[d] * tau[c * dim + d]; // This is tau[c][d] + } + + // add in the contribution + flux[CompressibleFlowFields::RHOU + c] = viscousFlux; + } + + // energy equation + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; ++d) { + PetscReal heatFlux = 0.0; + // add in the contributions for this viscous terms + for (PetscInt c = 0; c < dim; ++c) { + heatFlux += aux[aOff[VEL] + c] * tau[d * dim + c]; + } + + // heat conduction (-k dT/dx - k dT/dy - k dT/dz) . n A + heatFlux += k * gradAux[aOff_x[T] + d]; + + // Multiply by the area normal + heatFlux *= -fg->normal[d]; + + flux[CompressibleFlowFields::RHOE] += heatFlux; + } + + // zero out the density flux + flux[CompressibleFlowFields::RHO] = 0.0; + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionEnergyFlux(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + // compute the temperature in this volume + const PetscReal temperature = aux[aOff[temp]]; + PetscCall(flowParameters->computeSpeciesSensibleEnthalpyFunction.function( + field, temperature, flowParameters->speciesSpeciesSensibleEnthalpy.data(), flowParameters->computeSpeciesSensibleEnthalpyFunction.context.get())); + + // set the non rho E fluxes to zero + flux[CompressibleFlowFields::RHO] = 0.0; + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; d++) { + flux[CompressibleFlowFields::RHOU + d] = 0.0; + } + + // compute diff, this can be constant or variable + PetscReal diff = 0.0; + flowParameters->diffFunction.function(field, temperature, &diff, flowParameters->diffFunction.context.get()); + + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * diff * flowParameters->speciesSpeciesSensibleEnthalpy[sp] * gradAux[offset]; + flux[CompressibleFlowFields::RHOE] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionEnergyFluxVariableDiffusionCoefficient( + PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], + const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + // compute the temperature in this volume + const PetscReal temperature = aux[aOff[temp]]; + PetscCall(flowParameters->computeSpeciesSensibleEnthalpyFunction.function( + field, temperature, flowParameters->speciesSpeciesSensibleEnthalpy.data(), flowParameters->computeSpeciesSensibleEnthalpyFunction.context.get())); + + // set the non rho E fluxes to zero + flux[CompressibleFlowFields::RHO] = 0.0; + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; d++) { + flux[CompressibleFlowFields::RHOU + d] = 0.0; + } + + // compute diff, this can be constant or variable + flowParameters->diffFunction.function(field, temperature, flowParameters->speciesDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * flowParameters->speciesDiffusionCoefficient[sp] * flowParameters->speciesSpeciesSensibleEnthalpy[sp] * gradAux[offset]; + flux[CompressibleFlowFields::RHOE] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionSpeciesFlux(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + const PetscReal temperature = aux[aOff[temp]]; + + // compute diff + PetscReal diff = 0.0; + flowParameters->diffFunction.function(field, temperature, &diff, flowParameters->diffFunction.context.get()); + + // species equations + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + flux[sp] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * diff * gradAux[offset]; + flux[sp] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionSpeciesFluxVariableDiffusionCoefficient( + PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], + const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + const PetscReal temperature = aux[aOff[temp]]; + + // compute diff + flowParameters->diffFunction.function(field, temperature, flowParameters->speciesDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + + // species equations + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + flux[sp] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * flowParameters->speciesDiffusionCoefficient[sp] * gradAux[offset]; + flux[sp] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionEVFlux(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int EULER_FIELD = 0; + const int EV_FIELD = 0; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + const PetscReal temperature = aux[aOff[1]]; // Temperature is second aux in + // compute diff + PetscReal diff = 0.0; + flowParameters->evDiffFunction.function(field, temperature, &diff, flowParameters->evDiffFunction.context.get()); + + // species equations + for (PetscInt ev = 0; ev < flowParameters->numberEV; ++ev) { + flux[ev] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[EV_FIELD] + (ev * dim) + d; + PetscReal evFlux = -fg->normal[d] * density * diff * gradAux[offset]; + flux[ev] += evFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport::DiffusionEVFluxVariableDiffusionCoefficient( + PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], + const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void *ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int EULER_FIELD = 0; + const int EV_FIELD = 0; + + auto flowParameters = (DiffusionData *)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + const PetscReal temperature = aux[aOff[1]]; // Temperature is second aux in + // compute diff + flowParameters->diffFunction.function(field, temperature, flowParameters->evDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + + // species equations + for (PetscInt ev = 0; ev < flowParameters->numberEV; ++ev) { + flux[ev] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[EV_FIELD] + (ev * dim) + d; + PetscReal evFlux = -fg->normal[d] * density * flowParameters->evDiffusionCoefficient[ev] * gradAux[offset]; + flux[ev] += evFlux; + } + } + PetscFunctionReturn(0); +} diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.hpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.hpp new file mode 100644 index 000000000..56d674705 --- /dev/null +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.hpp @@ -0,0 +1,213 @@ +#ifndef ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESSingleProgressTRANSPORT_H +#define ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESSingleProgressTRANSPORT_H + +#include +#include "eos/transport/transportModel.hpp" +#include "finiteVolume/fluxCalculator/fluxCalculator.hpp" +#include "flowProcess.hpp" +#include "pressureGradientScaling.hpp" + +namespace ablate::finiteVolume::processes { + +class CompactCompressibleNSSpeciesSingleProgressTransport : public FlowProcess { + private: + // store ctx needed for function advection function that is passed into Petsc + struct AdvectionData { + // flow CFL + PetscReal cfl; + + /* number of gas species and extra species */ + PetscInt numberSpecies; + PetscInt numberEV; + + // EOS function calls (For Temperature it's better to just use the TemperatureTemperature function so we can guess Temperature from t + eos::ThermodynamicTemperatureFunction computeTemperature; + eos::ThermodynamicTemperatureFunction computeInternalEnergy; + eos::ThermodynamicTemperatureFunction computeSpeedOfSound; + eos::ThermodynamicTemperatureFunction computePressure; + + /* store method used for flux calculator */ + ablate::finiteVolume::fluxCalculator::FluxCalculatorFunction fluxCalculatorFunction; + void* fluxCalculatorCtx; + }; + AdvectionData advectionData; + + struct DiffusionData { + /* thermal conductivity Diffusivity, and Dynamic Viscosity*/ + eos::ThermodynamicTemperatureFunction kFunction; + eos::ThermodynamicTemperatureFunction muFunction; + eos::ThermodynamicTemperatureFunction diffFunction; + eos::ThermodynamicTemperatureFunction evDiffFunction; + + /* number of gas species and components in the SingleProgress */ + PetscInt numberSpecies; + PetscInt numberEV; + + /* functions to compute species enthalpy */ + eos::ThermodynamicTemperatureFunction computeSpeciesSensibleEnthalpyFunction; + /* store a scratch space for speciesSpeciesSensibleEnthalpy */ + std::vector speciesSpeciesSensibleEnthalpy; + /* store an optional scratch space for individual species diffusion */ + std::vector speciesDiffusionCoefficient; + /* store a scratch space for evDiffusionCoefficient */ + std::vector evDiffusionCoefficient; + }; + DiffusionData diffusionData; + + // Store the required ctx for time stepping + struct CflTimeStepData { + /* thermal conductivity*/ + AdvectionData* advectionData; + /* pressure gradient scaling */ + std::shared_ptr pgs; + }; + CflTimeStepData timeStepData; + + //! methods and functions to compute diffusion based time stepping + struct DiffusionTimeStepData { + /* number of gas species */ + PetscInt numberSpecies; + /* store an optional scratch space for individual species diffusion */ + std::vector speciesDiffusionCoefficient; + PetscReal SingleProgressDiffusionCoefficient; + + //! stability factor for condition time step. 0 (default) does not compute factor + PetscReal diffusiveStabilityFactor; + //! stability factor for condition time step. 0 (default) does not compute factor + PetscReal conductionStabilityFactor; + //! stability factor for viscous diffusion time step. 0 (default) does not compute factor + PetscReal viscousStabilityFactor; + + /* diffusivity */ + eos::ThermodynamicTemperatureFunction diffFunction; + /* thermal conductivity*/ + eos::ThermodynamicTemperatureFunction kFunction; + /* dynamic viscosity*/ + eos::ThermodynamicTemperatureFunction muFunction; + /* specific heat*/ + eos::ThermodynamicTemperatureFunction specificHeat; + /* density */ + eos::ThermodynamicTemperatureFunction density; + }; + DiffusionTimeStepData diffusionTimeStepData; + + const std::shared_ptr fluxCalculator; + const std::shared_ptr eos; + const std::shared_ptr transportModel; + const std::shared_ptr SingleProgressTransportModel; + + eos::ThermodynamicTemperatureFunction computeTemperatureFunction; + eos::ThermodynamicFunction computePressureFunction; + + public: + explicit CompactCompressibleNSSpeciesSingleProgressTransport(const std::shared_ptr& parameters, std::shared_ptr eos, + std::shared_ptr fluxCalcIn = {}, std::shared_ptr baseTransport = {}, + std::shared_ptr evTransport = {}, + std::shared_ptr = {}); + /** + * public function to link this process with the flow + * @param flow + */ + void Setup(ablate::finiteVolume::FiniteVolumeSolver& flow) override; + + /** + * This Computes the Advective Flow for rho, rhoE, and rhoVel, rhoYi, and rhoEV. + * u = {"euler"} or {"euler", "densityYi"} if species are tracked + * a = {} + * ctx = FlowData_CompressibleFlow + * @return + */ + static PetscErrorCode AdvectionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscScalar fieldL[], const PetscScalar fieldR[], const PetscInt aOff[], + const PetscScalar auxL[], const PetscScalar auxR[], PetscScalar* flux, void* ctx); + + /** + * This Computes the diffusion flux for euler rhoE, rhoVel + * u = {"euler", "densityYi"} + * a = {"temperature", "velocity"} + * ctx = FlowData_CompressibleFlow + * @return + */ + static PetscErrorCode DiffusionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the energy transfer for species diffusion flux for rhoE + * f = "euler" + * u = {"euler", "densityYi"} + * a = {"yi"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEnergyFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the energy transfer for species diffusion flux for rhoE for variable diffusion coefficient + * f = "euler" + * u = {"euler", "densityYi"} + * a = {"yi"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEnergyFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the species transfer for species diffusion flux + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionSpeciesFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the species transfer for species diffusion flux for variable diffusion coefficient + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionSpeciesFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + /** + * This computes the species transfer for species diffusion flux + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEVFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the species transfer for species diffusion flux for variable diffusion coefficient + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEVFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], + PetscScalar flux[], void* ctx); + + // static function to compute time step for euler advection + static double ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeViscousSpeciesDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); +}; + +} // namespace ablate::finiteVolume::processes + +#endif // ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESSingleProgressTRANSPORT_H diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp new file mode 100644 index 000000000..bff9a76a2 --- /dev/null +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp @@ -0,0 +1,795 @@ +#include "compactCompressibleNSSpeciesTransport.hpp" +#include +#include "finiteVolume/compressibleFlowFields.hpp" +#include "finiteVolume/fluxCalculator/ausm.hpp" +#include "finiteVolume/processes/navierStokesTransport.hpp" +#include "finiteVolume/processes/speciesTransport.hpp" +#include "parameters/emptyParameters.hpp" +#include "utilities/constants.hpp" +#include "utilities/mathUtilities.hpp" +#include "utilities/petscUtilities.hpp" + +ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::CompactCompressibleNSSpeciesTransport(const std::shared_ptr& parametersIn, + std::shared_ptr eosIn, + std::shared_ptr fluxCalcIn, + std::shared_ptr baseTransport, + std::shared_ptr pgs) + : advectionData(), fluxCalculator(fluxCalcIn), eos(std::move(eosIn)), transportModel(std::move(baseTransport)) { + auto parameters = ablate::parameters::EmptyParameters::Check(parametersIn); + + if (fluxCalculator) { + // cfl + advectionData.cfl = parameters->Get("cfl", 0.5); + + // extract the difference function from fluxDifferencer object + advectionData.fluxCalculatorFunction = fluxCalculator->GetFluxCalculatorFunction(); + advectionData.fluxCalculatorCtx = fluxCalculator->GetFluxCalculatorContext(); + + advectionData.numberSpecies = (PetscInt)eos->GetSpeciesVariables().size(); + } + + timeStepData.advectionData = &advectionData; + timeStepData.pgs = std::move(pgs); + + if (transportModel) { + // Add in the time stepping + diffusionTimeStepData.conductionStabilityFactor = parameters->Get("conductionStabilityFactor", 0.0); + diffusionTimeStepData.viscousStabilityFactor = parameters->Get("viscousStabilityFactor", 0.0); + diffusionTimeStepData.diffusiveStabilityFactor = parameters->Get("speciesStabilityFactor", 0.0); + diffusionTimeStepData.speciesDiffusionCoefficient.resize(eos->GetSpeciesVariables().size()); + + diffusionData.numberSpecies = (PetscInt)eos->GetSpeciesVariables().size(); + diffusionData.speciesSpeciesSensibleEnthalpy.resize(eos->GetSpeciesVariables().size()); + diffusionData.speciesDiffusionCoefficient.resize(eos->GetSpeciesVariables().size()); + } +} + +void ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::Setup(ablate::finiteVolume::FiniteVolumeSolver& flow) { + // Register the euler,eulerYi, and species source terms + if (fluxCalculator) { + // I don't know why we wouldn't push through the old temperature fields, maybe slower for perfect gas/idealized gas's but when there is a temperature iterative method this should be better + // If it is worse for perfect gas's, going to need to add in an option switch -klb + flow.RegisterRHSFunction(AdvectionFlux, + &advectionData, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::TEMPERATURE_FIELD}); + + // Set the ComputeCFLTimestepFrom flow Process through + flow.RegisterComputeTimeStepFunction(ComputeCflTimeStep, &timeStepData, "cfl"); + + advectionData.computeTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + advectionData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, flow.GetSubDomain().GetFields()); + advectionData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); + advectionData.computePressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); + } + + // if there are any coefficients for diffusion, compute diffusion, for here we will follow suit in allowing multiple diffusion functions + + if (transportModel) { + // Store the required data for the low level c functions + diffusionData.muFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Viscosity, flow.GetSubDomain().GetFields()); + diffusionData.kFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Conductivity, flow.GetSubDomain().GetFields()); + diffusionData.diffFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Diffusivity, flow.GetSubDomain().GetFields()); + diffusionData.computeSpeciesSensibleEnthalpyFunction = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeciesSensibleEnthalpy, flow.GetSubDomain().GetFields()); + + /* For now we will do 3 different diffusive flux calls just since it doesn't seem like their splitting costs too much time in redundant calculations */ + if (diffusionData.muFunction.function || diffusionData.kFunction.function) { + // Register the Diffusion Source term + flow.RegisterRHSFunction(DiffusionFlux, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::TEMPERATURE_FIELD, CompressibleFlowFields::VELOCITY_FIELD}); + } + if (diffusionData.diffFunction.function) { + // Specify a different rhs function depending on if the diffusion flux is constant + if (diffusionData.diffFunction.propertySize == 1) { + flow.RegisterRHSFunction(DiffusionEnergyFlux, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + flow.RegisterRHSFunction(DiffusionSpeciesFlux, + &diffusionData, + {CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else if (diffusionData.diffFunction.propertySize == diffusionData.numberSpecies) { + flow.RegisterRHSFunction(DiffusionEnergyFluxVariableDiffusionCoefficient, + &diffusionData, + {CompressibleFlowFields::EULER_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + flow.RegisterRHSFunction(DiffusionSpeciesFluxVariableDiffusionCoefficient, + &diffusionData, + {CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); + } else { + throw std::invalid_argument("The diffusion property size must be 1 or number of species in ablate::finiteVolume::processes::SpeciesTransport."); + } + } + + // Check to see if time step calculations should be added for viscosity or conduction + if (diffusionTimeStepData.conductionStabilityFactor > 0 || diffusionTimeStepData.viscousStabilityFactor > 0) { + diffusionTimeStepData.kFunction = diffusionData.kFunction; + diffusionTimeStepData.muFunction = diffusionData.muFunction; + diffusionTimeStepData.specificHeat = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpecificHeatConstantVolume, flow.GetSubDomain().GetFields()); + diffusionTimeStepData.density = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Density, flow.GetSubDomain().GetFields()); + + diffusionTimeStepData.numberSpecies = eos->GetSpeciesVariables().size(); + + if (diffusionTimeStepData.conductionStabilityFactor > 0) { + flow.RegisterComputeTimeStepFunction(ComputeConductionTimeStep, &diffusionTimeStepData, "cond"); + } + if (diffusionTimeStepData.viscousStabilityFactor > 0) { + flow.RegisterComputeTimeStepFunction(ComputeViscousDiffusionTimeStep, &diffusionTimeStepData, "visc"); + } + } + if (diffusionTimeStepData.diffusiveStabilityFactor > 0) { + diffusionTimeStepData.numberSpecies = diffusionData.numberSpecies; + diffusionTimeStepData.diffFunction = diffusionData.diffFunction; + flow.RegisterComputeTimeStepFunction(ComputeViscousSpeciesDiffusionTimeStep, &diffusionTimeStepData, "spec"); + } + } + + // Setup up aux updates and normalizations + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::VELOCITY_FIELD)) { + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxVelocityField, + nullptr, + std::vector{CompressibleFlowFields::VELOCITY_FIELD}, + {CompressibleFlowFields::EULER_FIELD}); + } + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::TEMPERATURE_FIELD)) { + computeTemperatureFunction = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + flow.RegisterAuxFieldUpdate( + ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxTemperatureField, &computeTemperatureFunction, std::vector{CompressibleFlowFields::TEMPERATURE_FIELD}, {}); + } + if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::PRESSURE_FIELD)) { + computePressureFunction = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); + flow.RegisterAuxFieldUpdate( + ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxPressureField, &computePressureFunction, std::vector{CompressibleFlowFields::PRESSURE_FIELD}, {}); + } + flow.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::SpeciesTransport::UpdateAuxMassFractionField, + &advectionData.numberSpecies, + std::vector{CompressibleFlowFields::YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}); + // clean up the species + flow.RegisterPostEvaluate(ablate::finiteVolume::processes::SpeciesTransport::NormalizeSpecies); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::AdvectionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt* uOff, const PetscScalar* fieldL, + const PetscScalar* fieldR, const PetscInt* aOff, const PetscScalar* auxL, const PetscScalar* auxR, + PetscScalar* flux, void* ctx) { + PetscFunctionBeginUser; + + auto advectionData = (AdvectionData*)ctx; + + const int EULER_FIELD = 0; + const int RHOYI_FIELD = 1; + + // Compute the norm + PetscReal norm[3]; + utilities::MathUtilities::NormVector(dim, fg->normal, norm); + const PetscReal areaMag = utilities::MathUtilities::MagVector(dim, fg->normal); + + // Decode the left and right states + PetscReal densityL; + PetscReal normalVelocityL; + PetscReal velocityL[3]; + PetscReal internalEnergyL; + PetscReal aL; + PetscReal pL; + + // decode the left side + { + densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + PetscReal temperatureL; + + PetscCall(advectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + .33 * auxR[aOff[0]], &temperatureL, advectionData->computeTemperature.context.get())); + + // Get the velocity in this direction + normalVelocityL = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velocityL[d] = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityL; + normalVelocityL += velocityL[d] * norm[d]; + } + PetscCall(advectionData->computeInternalEnergy.function(fieldL, temperatureL, &internalEnergyL, advectionData->computeInternalEnergy.context.get())); + PetscCall(advectionData->computeSpeedOfSound.function(fieldL, temperatureL, &aL, advectionData->computeSpeedOfSound.context.get())); + PetscCall(advectionData->computePressure.function(fieldL, temperatureL, &pL, advectionData->computePressure.context.get())); + } + + PetscReal densityR; + PetscReal normalVelocityR; + PetscReal velocityR[3]; + PetscReal internalEnergyR; + PetscReal aR; + PetscReal pR; + + { // decode right state + densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + PetscReal temperatureR; + + PetscCall(advectionData->computeTemperature.function(fieldR, auxR[aOff[0]] * .67 + .33 * auxL[aOff[0]], &temperatureR, advectionData->computeTemperature.context.get())); + + // Get the velocity in this direction + normalVelocityR = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velocityR[d] = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityR; + normalVelocityR += velocityR[d] * norm[d]; + } + + PetscCall(advectionData->computeInternalEnergy.function(fieldR, temperatureR, &internalEnergyR, advectionData->computeInternalEnergy.context.get())); + PetscCall(advectionData->computeSpeedOfSound.function(fieldR, temperatureR, &aR, advectionData->computeSpeedOfSound.context.get())); + PetscCall(advectionData->computePressure.function(fieldR, temperatureR, &pR, advectionData->computePressure.context.get())); + } + + // get the face values + PetscReal massFlux; + PetscReal p12; + + fluxCalculator::Direction direction = + advectionData->fluxCalculatorFunction(advectionData->fluxCalculatorCtx, normalVelocityL, aL, densityL, pL, normalVelocityR, aR, densityR, pR, &massFlux, &p12); + + if (direction == fluxCalculator::LEFT) { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); + PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; + flux[CompressibleFlowFields::RHOE] = HL * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = velocityL[n] * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldL[uOff[RHOYI_FIELD] + ns] / densityL * areaMag; + } else if (direction == fluxCalculator::RIGHT) { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); + PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; + flux[CompressibleFlowFields::RHOE] = HR * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = velocityR[n] * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldR[uOff[RHOYI_FIELD] + ns] / densityR * areaMag; + } else { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; + + PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); + PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; + + PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); + PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; + + flux[CompressibleFlowFields::RHOE] = 0.5 * (HL + HR) * massFlux * areaMag; + for (PetscInt n = 0; n < dim; n++) { + flux[CompressibleFlowFields::RHOU + n] = 0.5 * (velocityL[n] + velocityR[n]) * massFlux * areaMag + p12 * fg->normal[n]; + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) + flux[uOff[RHOYI_FIELD] + ns] = massFlux * 0.5 * (fieldR[uOff[RHOYI_FIELD] + ns] + fieldL[uOff[RHOYI_FIELD] + ns]) / (0.5 * (densityL + densityR)) * areaMag; + } + + PetscFunctionReturn(0); +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto timeStepData = (CflTimeStepData*)ctx; + auto advectionData = timeStepData->advectionData; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar* locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + const PetscScalar* x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for euler and densityYi + auto eulerId = flow.GetSubDomain().GetField("euler").id; + + // Get alpha if provided + PetscReal pgsAlpha = 1.0; + if (timeStepData->pgs) { + pgsAlpha = timeStepData->pgs->GetAlpha(); + } + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal* euler; + const PetscReal* conserved = NULL; + const PetscReal* cellCharacteristics = NULL; + DMPlexPointGlobalFieldRead(dm, cell, eulerId, x, &euler) >> utilities::PetscUtilities::checkError; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (euler) { // must be real cell and not ghost + PetscReal rho = euler[CompressibleFlowFields::RHO]; + + // Get the speed of sound from the eos + // TODO:: Replace this with a better temperature guess (see compute conduction Time Step below) + PetscReal temperature; + advectionData->computeTemperature.function(conserved, 300, &temperature, advectionData->computeTemperature.context.get()) >> utilities::PetscUtilities::checkError; + PetscReal a; + advectionData->computeSpeedOfSound.function(conserved, temperature, &a, advectionData->computeSpeedOfSound.context.get()) >> utilities::PetscUtilities::checkError; + + PetscReal dx = 2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]; + + PetscReal velSum = 0.0; + for (PetscInt d = 0; d < dim; d++) { + velSum += PetscAbsReal(euler[CompressibleFlowFields::RHOU + d]) / rho; + } + PetscReal dt = advectionData->cfl * dx / (a / pgsAlpha + velSum); + + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + return dtMin; +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData*)ctx; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar* locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar* x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar* aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto kFunction = diffusionData->kFunction.function; + auto kFunctionContext = diffusionData->kFunction.context.get(); + auto cvFunction = diffusionData->specificHeat.function; + auto cvFunctionContext = diffusionData->specificHeat.context.get(); + auto density = diffusionData->density.function; + auto densityContext = diffusionData->density.context.get(); + auto stabFactor = diffusionData->conductionStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal* conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal* temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + const PetscReal* cellCharacteristics = NULL; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal k; + kFunction(conserved, *temperature, &k, kFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal cv; + cvFunction(conserved, *temperature, &cv, cvFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal rho; + density(conserved, *temperature, &rho, densityContext) >> utilities::PetscUtilities::checkError; + + // Compute alpha + PetscReal alpha = k / (rho * cv); + + PetscReal dx2 = PetscSqr(2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]); + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / alpha); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + + return dtMin; +} + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData*)ctx; + + // Get the fv geom + Vec locCharacteristicsVec; + DM characteristicsDm; + const PetscScalar* locCharacteristicsArray; + flow.GetMeshCharacteristics(characteristicsDm, locCharacteristicsVec); + VecGetArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar* x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar* aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto muFunction = diffusionData->muFunction.function; + auto muFunctionContext = diffusionData->muFunction.context.get(); + auto density = diffusionData->density.function; + auto densityContext = diffusionData->density.context.get(); + auto stabFactor = diffusionData->viscousStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal* conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal* temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + const PetscReal* cellCharacteristics = NULL; + DMPlexPointLocalRead(characteristicsDm, cell, locCharacteristicsArray, &cellCharacteristics) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal mu; + muFunction(conserved, *temperature, &mu, muFunctionContext) >> utilities::PetscUtilities::checkError; + PetscReal rho; + density(conserved, *temperature, &rho, densityContext) >> utilities::PetscUtilities::checkError; + + // Compute nu + PetscReal nu = mu / rho; + + PetscReal dx2 = PetscSqr(2.0 * cellCharacteristics[FiniteVolumeSolver::MIN_CELL_RADIUS]); + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / nu); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(locCharacteristicsVec, &locCharacteristicsArray) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + return dtMin; +} +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::ComputeViscousSpeciesDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx) { + // Get the dm and current solution vector + DM dm; + TSGetDM(ts, &dm) >> utilities::PetscUtilities::checkError; + Vec v; + TSGetSolution(ts, &v) >> utilities::PetscUtilities::checkError; + + // Get the flow param + auto diffusionData = (DiffusionTimeStepData*)ctx; + + // Get the fv geom + PetscReal minCellRadius; + DMPlexGetGeometryFVM(dm, NULL, NULL, &minCellRadius) >> utilities::PetscUtilities::checkError; + + // Get the valid cell range over this region + ablate::domain::Range cellRange; + flow.GetCellRangeWithoutGhost(cellRange); + + // Get the solution data + const PetscScalar* x; + VecGetArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + + // Get the auxData + const PetscScalar* aux; + const DM auxDM = flow.GetSubDomain().GetAuxDM(); + VecGetArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + + // Get the dim from the dm + PetscInt dim; + DMGetDimension(dm, &dim) >> utilities::PetscUtilities::checkError; + + // assume the smallest cell is the limiting factor for now + const PetscReal dx2 = PetscSqr(2.0 * minCellRadius); + + // Get field location for temperature + auto temperatureField = flow.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD).id; + + // get functions + auto diffFunction = diffusionData->diffFunction.function; + auto diffFunctionContext = diffusionData->diffFunction.context.get(); + auto stabFactor = diffusionData->diffusiveStabilityFactor; + + // March over each cell + PetscReal dtMin = ablate::utilities::Constants::large; + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + auto cell = cellRange.GetPoint(c); + + const PetscReal* conserved = NULL; + DMPlexPointGlobalRead(dm, cell, x, &conserved) >> utilities::PetscUtilities::checkError; + + const PetscReal* temperature = NULL; + DMPlexPointLocalFieldRead(auxDM, cell, temperatureField, aux, &temperature) >> utilities::PetscUtilities::checkError; + + if (conserved) { // must be real cell and not ghost + PetscReal diff; + diffFunction(conserved, *temperature, &diff, diffFunctionContext) >> utilities::PetscUtilities::checkError; + + // compute dt + double dt = PetscAbs(stabFactor * dx2 / diff); + dtMin = PetscMin(dtMin, dt); + } + } + VecRestoreArrayRead(v, &x) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(flow.GetSubDomain().GetAuxGlobalVector(), &aux) >> utilities::PetscUtilities::checkError; + flow.RestoreRange(cellRange); + return dtMin; +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::DiffusionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], + const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], + const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], + void* ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int T = 0; + const int VEL = 1; + + auto flowParameters = (DiffusionData*)ctx; + + // Compute mu and k + PetscReal mu = 0.0; + flowParameters->muFunction.function(field, aux[aOff[T]], &mu, flowParameters->muFunction.context.get()); + PetscReal k = 0.0; + flowParameters->kFunction.function(field, aux[aOff[T]], &k, flowParameters->kFunction.context.get()); + + // Compute the stress tensor tau + PetscReal tau[9]; // Maximum size without symmetry + PetscCall(ablate::finiteVolume::processes::NavierStokesTransport::CompressibleFlowComputeStressTensor(dim, mu, gradAux + aOff_x[VEL], tau)); + + // for each velocity component + for (PetscInt c = 0; c < dim; ++c) { + PetscReal viscousFlux = 0.0; + + // March over each direction + for (PetscInt d = 0; d < dim; ++d) { + viscousFlux += -fg->normal[d] * tau[c * dim + d]; // This is tau[c][d] + } + + // add in the contribution + flux[CompressibleFlowFields::RHOU + c] = viscousFlux; + } + + // energy equation + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; ++d) { + PetscReal heatFlux = 0.0; + // add in the contributions for this viscous terms + for (PetscInt c = 0; c < dim; ++c) { + heatFlux += aux[aOff[VEL] + c] * tau[d * dim + c]; + } + + // heat conduction (-k dT/dx - k dT/dy - k dT/dz) . n A + heatFlux += k * gradAux[aOff_x[T] + d]; + + // Multiply by the area normal + heatFlux *= -fg->normal[d]; + + flux[CompressibleFlowFields::RHOE] += heatFlux; + } + + // zero out the density flux + flux[CompressibleFlowFields::RHO] = 0.0; + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::DiffusionEnergyFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], + const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], + const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], + PetscScalar flux[], void* ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData*)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + // compute the temperature in this volume + const PetscReal temperature = aux[aOff[temp]]; + PetscCall(flowParameters->computeSpeciesSensibleEnthalpyFunction.function( + field, temperature, flowParameters->speciesSpeciesSensibleEnthalpy.data(), flowParameters->computeSpeciesSensibleEnthalpyFunction.context.get())); + + // set the non rho E fluxes to zero + flux[CompressibleFlowFields::RHO] = 0.0; + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; d++) { + flux[CompressibleFlowFields::RHOU + d] = 0.0; + } + + // compute diff, this can be constant or variable + PetscReal diff = 0.0; + flowParameters->diffFunction.function(field, temperature, &diff, flowParameters->diffFunction.context.get()); + + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * diff * flowParameters->speciesSpeciesSensibleEnthalpy[sp] * gradAux[offset]; + flux[CompressibleFlowFields::RHOE] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::DiffusionEnergyFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], + const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData*)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + // compute the temperature in this volume + const PetscReal temperature = aux[aOff[temp]]; + PetscCall(flowParameters->computeSpeciesSensibleEnthalpyFunction.function( + field, temperature, flowParameters->speciesSpeciesSensibleEnthalpy.data(), flowParameters->computeSpeciesSensibleEnthalpyFunction.context.get())); + + // set the non rho E fluxes to zero + flux[CompressibleFlowFields::RHO] = 0.0; + flux[CompressibleFlowFields::RHOE] = 0.0; + for (PetscInt d = 0; d < dim; d++) { + flux[CompressibleFlowFields::RHOU + d] = 0.0; + } + + // compute diff, this can be constant or variable + flowParameters->diffFunction.function(field, temperature, flowParameters->speciesDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * flowParameters->speciesDiffusionCoefficient[sp] * flowParameters->speciesSpeciesSensibleEnthalpy[sp] * gradAux[offset]; + flux[CompressibleFlowFields::RHOE] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::DiffusionSpeciesFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], + const PetscScalar field[], const PetscScalar grad[], const PetscInt aOff[], + const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], + PetscScalar flux[], void* ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData*)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + + const PetscReal temperature = aux[aOff[temp]]; + + // compute diff + PetscReal diff = 0.0; + flowParameters->diffFunction.function(field, temperature, &diff, flowParameters->diffFunction.context.get()); + + // species equations + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + flux[sp] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * diff * gradAux[offset]; + flux[sp] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::DiffusionSpeciesFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], + const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], + const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx) { + PetscFunctionBeginUser; + // this order is based upon the order that they are passed into RegisterRHSFunction + const int yi = 0; + const int euler = 0; + const int temp = 1; + + auto flowParameters = (DiffusionData*)ctx; + + // get the current density from euler + const PetscReal density = field[uOff[euler] + CompressibleFlowFields::RHO]; + const PetscReal temperature = aux[aOff[temp]]; + + // compute diff + flowParameters->diffFunction.function(field, temperature, flowParameters->speciesDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + + // species equations + for (PetscInt sp = 0; sp < flowParameters->numberSpecies; ++sp) { + flux[sp] = 0; + for (PetscInt d = 0; d < dim; ++d) { + // speciesFlux(-rho Di dYi/dx - rho Di dYi/dy - rho Di dYi//dz) . n A + const int offset = aOff_x[yi] + (sp * dim) + d; + PetscReal speciesFlux = -fg->normal[d] * density * flowParameters->speciesDiffusionCoefficient[sp] * gradAux[offset]; + flux[sp] += speciesFlux; + } + } + + PetscFunctionReturn(0); +} diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp new file mode 100644 index 000000000..91cdb47ac --- /dev/null +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp @@ -0,0 +1,181 @@ +#ifndef ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESTRANSPORT_H +#define ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESTRANSPORT_H + +#include +#include "eos/transport/transportModel.hpp" +#include "finiteVolume/fluxCalculator/fluxCalculator.hpp" +#include "flowProcess.hpp" +#include "pressureGradientScaling.hpp" + +namespace ablate::finiteVolume::processes { + +class CompactCompressibleNSSpeciesTransport : public FlowProcess { + private: + // store ctx needed for function advection function that is passed into Petsc + struct AdvectionData { + // flow CFL + PetscReal cfl; + + /* number of gas species and extra species */ + PetscInt numberSpecies; + + // EOS function calls (For Temperature it's better to just use the TemperatureTemperature function so we can guess Temperature from t + eos::ThermodynamicTemperatureFunction computeTemperature; + eos::ThermodynamicTemperatureFunction computeInternalEnergy; + eos::ThermodynamicTemperatureFunction computeSpeedOfSound; + eos::ThermodynamicTemperatureFunction computePressure; + + /* store method used for flux calculator */ + ablate::finiteVolume::fluxCalculator::FluxCalculatorFunction fluxCalculatorFunction; + void* fluxCalculatorCtx; + }; + AdvectionData advectionData; + + struct DiffusionData { + /* thermal conductivity Diffusivity, and Dynamic Viscosity*/ + eos::ThermodynamicTemperatureFunction kFunction; + eos::ThermodynamicTemperatureFunction muFunction; + eos::ThermodynamicTemperatureFunction diffFunction; + /* number of gas species */ + PetscInt numberSpecies; + /* functions to compute species enthalpy */ + eos::ThermodynamicTemperatureFunction computeSpeciesSensibleEnthalpyFunction; + /* store a scratch space for speciesSpeciesSensibleEnthalpy */ + std::vector speciesSpeciesSensibleEnthalpy; + /* store an optional scratch space for individual species diffusion */ + std::vector speciesDiffusionCoefficient; + }; + DiffusionData diffusionData; + + // Store the required ctx for time stepping + struct CflTimeStepData { + /* thermal conductivity*/ + AdvectionData* advectionData; + /* pressure gradient scaling */ + std::shared_ptr pgs; + }; + CflTimeStepData timeStepData; + + //! methods and functions to compute diffusion based time stepping + struct DiffusionTimeStepData { + /* number of gas species */ + PetscInt numberSpecies; + /* store an optional scratch space for individual species diffusion */ + std::vector speciesDiffusionCoefficient; + + //! stability factor for condition time step. 0 (default) does not compute factor + PetscReal diffusiveStabilityFactor; + //! stability factor for condition time step. 0 (default) does not compute factor + PetscReal conductionStabilityFactor; + //! stability factor for viscous diffusion time step. 0 (default) does not compute factor + PetscReal viscousStabilityFactor; + + /* diffusivity */ + eos::ThermodynamicTemperatureFunction diffFunction; + /* thermal conductivity*/ + eos::ThermodynamicTemperatureFunction kFunction; + /* dynamic viscosity*/ + eos::ThermodynamicTemperatureFunction muFunction; + /* specific heat*/ + eos::ThermodynamicTemperatureFunction specificHeat; + /* density */ + eos::ThermodynamicTemperatureFunction density; + }; + DiffusionTimeStepData diffusionTimeStepData; + + const std::shared_ptr fluxCalculator; + const std::shared_ptr eos; + const std::shared_ptr transportModel; + + eos::ThermodynamicTemperatureFunction computeTemperatureFunction; + eos::ThermodynamicFunction computePressureFunction; + + public: + explicit CompactCompressibleNSSpeciesTransport(const std::shared_ptr& parameters, std::shared_ptr eos, + std::shared_ptr fluxCalcIn = {}, std::shared_ptr baseTransport = {}, + std::shared_ptr = {}); + /** + * public function to link this process with the flow + * @param flow + */ + void Setup(ablate::finiteVolume::FiniteVolumeSolver& flow) override; + + /** + * This Computes the Advective Flow for rho, rhoE, and rhoVel, rhoYi, and rhoEV. + * u = {"euler"} or {"euler", "densityYi"} if species are tracked + * a = {} + * ctx = FlowData_CompressibleFlow + * @return + */ + static PetscErrorCode AdvectionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscScalar fieldL[], const PetscScalar fieldR[], const PetscInt aOff[], + const PetscScalar auxL[], const PetscScalar auxR[], PetscScalar* flux, void* ctx); + + /** + * This Computes the diffusion flux for euler rhoE, rhoVel + * u = {"euler", "densityYi"} + * a = {"temperature", "velocity"} + * ctx = FlowData_CompressibleFlow + * @return + */ + static PetscErrorCode DiffusionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the energy transfer for species diffusion flux for rhoE + * f = "euler" + * u = {"euler", "densityYi"} + * a = {"yi"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEnergyFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the energy transfer for species diffusion flux for rhoE for variable diffusion coefficient + * f = "euler" + * u = {"euler", "densityYi"} + * a = {"yi"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionEnergyFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the species transfer for species diffusion flux + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionSpeciesFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], const PetscScalar grad[], + const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + /** + * This computes the species transfer for species diffusion flux for variable diffusion coefficient + * f = "densityYi" + * u = {"euler"} + * a = {"yi", "T"} + * ctx = SpeciesDiffusionData + * @return + */ + static PetscErrorCode DiffusionSpeciesFluxVariableDiffusionCoefficient(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar field[], + const PetscScalar grad[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar aux[], + const PetscScalar gradAux[], PetscScalar flux[], void* ctx); + + // static function to compute time step for euler advection + static double ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeViscousSpeciesDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); +}; + +} // namespace ablate::finiteVolume::processes + +#endif // ABLATELIBRARY_COMPACTCOMPRESSIBLENSSPECIESTRANSPORT_H diff --git a/src/finiteVolume/processes/evTransport.cpp b/src/finiteVolume/processes/evTransport.cpp index ee2de66cf..03dee5d3a 100644 --- a/src/finiteVolume/processes/evTransport.cpp +++ b/src/finiteVolume/processes/evTransport.cpp @@ -34,24 +34,29 @@ void ablate::finiteVolume::processes::EVTransport::Setup(ablate::finiteVolume::F advectionData.fluxCalculatorCtx = fluxCalculator->GetFluxCalculatorContext(); // set decode state functions - advectionData.computeTemperature = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + advectionData.computeTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); advectionData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, flow.GetSubDomain().GetFields()); advectionData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); advectionData.computePressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); - flow.RegisterRHSFunction(AdvectionFlux, &advectionData, evConservedField.name, {CompressibleFlowFields::EULER_FIELD, evConservedField.name}, {}); + flow.RegisterRHSFunction(AdvectionFlux, &advectionData, {evConservedField.name}, {CompressibleFlowFields::EULER_FIELD, evConservedField.name}, {CompressibleFlowFields::TEMPERATURE_FIELD}); } if (transportModel) { diffusionData.evDiffusionCoefficient.resize(diffusionData.numberEV); - diffusionData.diffFunction = transportModel->GetTransportFunction(eos::transport::TransportProperty::Diffusivity, flow.GetSubDomain().GetFields()); + diffusionData.diffFunction = transportModel->GetTransportTemperatureFunction(eos::transport::TransportProperty::Diffusivity, flow.GetSubDomain().GetFields()); if (diffusionData.diffFunction.function) { if (diffusionData.diffFunction.propertySize == 1) { - flow.RegisterRHSFunction(DiffusionEVFlux, &diffusionData, evConservedField.name, {CompressibleFlowFields::EULER_FIELD}, {nonConserved}); + flow.RegisterRHSFunction( + DiffusionEVFlux, &diffusionData, {evConservedField.name}, {CompressibleFlowFields::EULER_FIELD}, {nonConserved, CompressibleFlowFields::TEMPERATURE_FIELD}); } else if (diffusionData.diffFunction.propertySize == diffusionData.numberEV) { - flow.RegisterRHSFunction(DiffusionEVFluxVariableDiffusionCoefficient, &diffusionData, evConservedField.name, {CompressibleFlowFields::EULER_FIELD}, {nonConserved}); + flow.RegisterRHSFunction(DiffusionEVFluxVariableDiffusionCoefficient, + &diffusionData, + {evConservedField.name}, + {CompressibleFlowFields::EULER_FIELD}, + {nonConserved, CompressibleFlowFields::TEMPERATURE_FIELD}); } else { throw std::invalid_argument("The diffusion property size must be 1 or number of ev in ablate::finiteVolume::processes::EVTransport."); } @@ -160,7 +165,7 @@ PetscErrorCode ablate::finiteVolume::processes::EVTransport::AdvectionFlux(Petsc densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureL; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldL, &temperatureL, eulerAdvectionData->computeTemperature.context.get())); + PetscCall(eulerAdvectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + auxR[aOff[0]] * .33, &temperatureL, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityL = 0.0; @@ -184,7 +189,7 @@ PetscErrorCode ablate::finiteVolume::processes::EVTransport::AdvectionFlux(Petsc densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureR; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldR, &temperatureR, eulerAdvectionData->computeTemperature.context.get())); + PetscCall(eulerAdvectionData->computeTemperature.function(fieldR, auxL[aOff[0]] * .33 + auxR[aOff[0]] * .67, &temperatureR, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityR = 0.0; @@ -231,10 +236,10 @@ PetscErrorCode ablate::finiteVolume::processes::EVTransport::DiffusionEVFlux(Pet // get the current density from euler const PetscReal density = field[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; - + const PetscReal temperature = aux[aOff[1]]; // Temperature is second aux in // compute diff PetscReal diff = 0.0; - flowParameters->diffFunction.function(field, &diff, flowParameters->diffFunction.context.get()); + flowParameters->diffFunction.function(field, temperature, &diff, flowParameters->diffFunction.context.get()); // species equations for (PetscInt ev = 0; ev < flowParameters->numberEV; ++ev) { @@ -263,9 +268,9 @@ PetscErrorCode ablate::finiteVolume::processes::EVTransport::DiffusionEVFluxVari // get the current density from euler const PetscReal density = field[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; - + const PetscReal temperature = aux[aOff[1]]; // Temperature is second aux in // compute diff - flowParameters->diffFunction.function(field, flowParameters->evDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); + flowParameters->diffFunction.function(field, temperature, flowParameters->evDiffusionCoefficient.data(), flowParameters->diffFunction.context.get()); // species equations for (PetscInt ev = 0; ev < flowParameters->numberEV; ++ev) { diff --git a/src/finiteVolume/processes/evTransport.hpp b/src/finiteVolume/processes/evTransport.hpp index 1d7283924..159e2ed42 100644 --- a/src/finiteVolume/processes/evTransport.hpp +++ b/src/finiteVolume/processes/evTransport.hpp @@ -25,7 +25,7 @@ class EVTransport : public FlowProcess { PetscInt numberEV; // EOS function calls - eos::ThermodynamicFunction computeTemperature; + eos::ThermodynamicTemperatureFunction computeTemperature; eos::ThermodynamicTemperatureFunction computeInternalEnergy; eos::ThermodynamicTemperatureFunction computeSpeedOfSound; eos::ThermodynamicTemperatureFunction computePressure; @@ -40,7 +40,7 @@ class EVTransport : public FlowProcess { PetscInt numberEV; /* functions to compute diffusion */ - eos::ThermodynamicFunction diffFunction; + eos::ThermodynamicTemperatureFunction diffFunction; /* store a scratch space for evDiffusionCoefficient */ std::vector evDiffusionCoefficient; diff --git a/src/finiteVolume/processes/les.cpp b/src/finiteVolume/processes/les.cpp index 7db5010f0..64dfa33b6 100644 --- a/src/finiteVolume/processes/les.cpp +++ b/src/finiteVolume/processes/les.cpp @@ -24,21 +24,21 @@ void ablate::finiteVolume::processes::LES::Setup(ablate::finiteVolume::FiniteVol // Register the euler/Energy LESdiffusion source term tke flow.RegisterRHSFunction(LesEnergyFlux, nullptr, - CompressibleFlowFields::EULER_FIELD, + {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::VELOCITY_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); // Register the N_S LESdiffusion source terms - flow.RegisterRHSFunction(LesMomentumFlux, nullptr, CompressibleFlowFields::EULER_FIELD, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::VELOCITY_FIELD}); + flow.RegisterRHSFunction(LesMomentumFlux, nullptr, {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::VELOCITY_FIELD}); // Register the tke LESdiffusion source term - flow.RegisterRHSFunction(LesTkeFlux, nullptr, conservedFieldName, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::VELOCITY_FIELD}); + flow.RegisterRHSFunction(LesTkeFlux, nullptr, {conservedFieldName}, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::VELOCITY_FIELD}); // Register the Species LESdiffusion source term if (flow.GetSubDomain().ContainsField(CompressibleFlowFields::DENSITY_YI_FIELD)) { // the species are treated like any other transported ev auto& numberComponent = numberComponents.emplace_back(flow.GetSubDomain().GetField(CompressibleFlowFields::YI_FIELD).numberComponents); - flow.RegisterRHSFunction(LesEvFlux, &numberComponent, CompressibleFlowFields::DENSITY_YI_FIELD, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::YI_FIELD}); + flow.RegisterRHSFunction(LesEvFlux, &numberComponent, {CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::EULER_FIELD}, {tkeField, CompressibleFlowFields::YI_FIELD}); } // March over any ev @@ -48,7 +48,7 @@ void ablate::finiteVolume::processes::LES::Setup(ablate::finiteVolume::FiniteVol // the species are treated like any other transported ev auto& numberComponent = numberComponents.emplace_back(evConservedField.numberComponents); - flow.RegisterRHSFunction(LesEvFlux, &numberComponent, evConservedField.name, {CompressibleFlowFields::EULER_FIELD}, {tkeField, nonConservedEvName}); + flow.RegisterRHSFunction(LesEvFlux, &numberComponent, {evConservedField.name}, {CompressibleFlowFields::EULER_FIELD}, {tkeField, nonConservedEvName}); } } diff --git a/src/finiteVolume/processes/navierStokesTransport.cpp b/src/finiteVolume/processes/navierStokesTransport.cpp index 27eed9f1a..185600be2 100644 --- a/src/finiteVolume/processes/navierStokesTransport.cpp +++ b/src/finiteVolume/processes/navierStokesTransport.cpp @@ -35,12 +35,14 @@ ablate::finiteVolume::processes::NavierStokesTransport::NavierStokesTransport(co void ablate::finiteVolume::processes::NavierStokesTransport::Setup(ablate::finiteVolume::FiniteVolumeSolver& flow) { // Register the euler source terms if (fluxCalculator) { - flow.RegisterRHSFunction(AdvectionFlux, &advectionData, CompressibleFlowFields::EULER_FIELD, {CompressibleFlowFields::EULER_FIELD}, {}); + // I don't know why we wouldn't push through the old temperature fields, maybe slower for perfect gas/idealized gas's but when there is a temperature iterative method this should be better + // If it is worse for perfect gas's, going to need to add in an option switch -klb + flow.RegisterRHSFunction(AdvectionFlux, &advectionData, {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::TEMPERATURE_FIELD}); // PetscErrorCode PetscOptionsGetBool(PetscOptions options,const char pre[],const char name[],PetscBool *ivalue,PetscBool *set) flow.RegisterComputeTimeStepFunction(ComputeCflTimeStep, &timeStepData, "cfl"); - advectionData.computeTemperature = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + advectionData.computeTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); advectionData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, flow.GetSubDomain().GetFields()); advectionData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); advectionData.computePressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); @@ -56,7 +58,7 @@ void ablate::finiteVolume::processes::NavierStokesTransport::Setup(ablate::finit // Register the euler diffusion source terms flow.RegisterRHSFunction(DiffusionFlux, &diffusionData, - CompressibleFlowFields::EULER_FIELD, + {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::TEMPERATURE_FIELD, CompressibleFlowFields::VELOCITY_FIELD}); } @@ -102,6 +104,7 @@ PetscErrorCode ablate::finiteVolume::processes::NavierStokesTransport::Advection auto eulerAdvectionData = (AdvectionData*)ctx; const int EULER_FIELD = 0; + const int TEMP_FIELD = 0; // Compute the norm PetscReal norm[3]; @@ -121,7 +124,8 @@ PetscErrorCode ablate::finiteVolume::processes::NavierStokesTransport::Advection densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureL; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldL, &temperatureL, eulerAdvectionData->computeTemperature.context.get())); + PetscCall( + eulerAdvectionData->computeTemperature.function(fieldL, auxL[aOff[TEMP_FIELD]] * .67 + .33 * auxR[aOff[TEMP_FIELD]], &temperatureL, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityL = 0.0; @@ -146,7 +150,8 @@ PetscErrorCode ablate::finiteVolume::processes::NavierStokesTransport::Advection densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureR; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldR, &temperatureR, eulerAdvectionData->computeTemperature.context.get())); + PetscCall( + eulerAdvectionData->computeTemperature.function(fieldR, auxL[aOff[TEMP_FIELD]] * .33 + .67 * auxR[aOff[TEMP_FIELD]], &temperatureR, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityR = 0.0; @@ -255,8 +260,9 @@ double ablate::finiteVolume::processes::NavierStokesTransport::ComputeCflTimeSte PetscReal rho = euler[CompressibleFlowFields::RHO]; // Get the speed of sound from the eos + // TODO:: Replace this with a better temperature guess (see compute conduction Time Step below) PetscReal temperature; - advectionData->computeTemperature.function(conserved, &temperature, advectionData->computeTemperature.context.get()) >> utilities::PetscUtilities::checkError; + advectionData->computeTemperature.function(conserved, 300, &temperature, advectionData->computeTemperature.context.get()) >> utilities::PetscUtilities::checkError; PetscReal a; advectionData->computeSpeedOfSound.function(conserved, temperature, &a, advectionData->computeSpeedOfSound.context.get()) >> utilities::PetscUtilities::checkError; diff --git a/src/finiteVolume/processes/navierStokesTransport.hpp b/src/finiteVolume/processes/navierStokesTransport.hpp index b957f2285..abdfa9355 100644 --- a/src/finiteVolume/processes/navierStokesTransport.hpp +++ b/src/finiteVolume/processes/navierStokesTransport.hpp @@ -20,7 +20,7 @@ class NavierStokesTransport : public FlowProcess { PetscInt numberSpecies; // EOS function calls - eos::ThermodynamicFunction computeTemperature; + eos::ThermodynamicTemperatureFunction computeTemperature; eos::ThermodynamicTemperatureFunction computeInternalEnergy; eos::ThermodynamicTemperatureFunction computeSpeedOfSound; eos::ThermodynamicTemperatureFunction computePressure; @@ -38,6 +38,13 @@ class NavierStokesTransport : public FlowProcess { eos::ThermodynamicTemperatureFunction muFunction; }; + // static function to compute time step for euler advection + static double ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + // static function to compute the conduction based time step + static double ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); + private: const std::shared_ptr fluxCalculator; const std::shared_ptr eos; @@ -81,15 +88,6 @@ class NavierStokesTransport : public FlowProcess { }; DiffusionTimeStepData diffusionTimeStepData; - // static function to compute time step for euler advection - static double ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); - - // static function to compute the conduction based time step - static double ComputeConductionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); - - // static function to compute the conduction based time step - static double ComputeViscousDiffusionTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx); - public: /** * Function to compute the temperature field. This function assumes that the input values will be {"euler", "densityYi"} diff --git a/src/finiteVolume/processes/speciesTransport.cpp b/src/finiteVolume/processes/speciesTransport.cpp index c4f741967..911ca6b95 100644 --- a/src/finiteVolume/processes/speciesTransport.cpp +++ b/src/finiteVolume/processes/speciesTransport.cpp @@ -35,8 +35,12 @@ ablate::finiteVolume::processes::SpeciesTransport::SpeciesTransport(std::shared_ void ablate::finiteVolume::processes::SpeciesTransport::Setup(ablate::finiteVolume::FiniteVolumeSolver &flow) { if (!eos->GetSpeciesVariables().empty()) { if (fluxCalculator) { - flow.RegisterRHSFunction(AdvectionFlux, &advectionData, CompressibleFlowFields::DENSITY_YI_FIELD, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {}); - advectionData.computeTemperature = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); + flow.RegisterRHSFunction(AdvectionFlux, + &advectionData, + {CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, + {CompressibleFlowFields::TEMPERATURE_FIELD}); + advectionData.computeTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Temperature, flow.GetSubDomain().GetFields()); advectionData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, flow.GetSubDomain().GetFields()); advectionData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); advectionData.computePressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::Pressure, flow.GetSubDomain().GetFields()); @@ -49,23 +53,23 @@ void ablate::finiteVolume::processes::SpeciesTransport::Setup(ablate::finiteVolu if (diffusionData.diffFunction.propertySize == 1) { flow.RegisterRHSFunction(DiffusionEnergyFlux, &diffusionData, - CompressibleFlowFields::EULER_FIELD, + {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); flow.RegisterRHSFunction(DiffusionSpeciesFlux, &diffusionData, - CompressibleFlowFields::DENSITY_YI_FIELD, + {CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); } else if (diffusionData.diffFunction.propertySize == numberSpecies) { flow.RegisterRHSFunction(DiffusionEnergyFluxVariableDiffusionCoefficient, &diffusionData, - CompressibleFlowFields::EULER_FIELD, + {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); flow.RegisterRHSFunction(DiffusionSpeciesFluxVariableDiffusionCoefficient, &diffusionData, - CompressibleFlowFields::DENSITY_YI_FIELD, + {CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::YI_FIELD, CompressibleFlowFields::TEMPERATURE_FIELD}); } else { @@ -280,7 +284,7 @@ PetscErrorCode ablate::finiteVolume::processes::SpeciesTransport::AdvectionFlux( densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureL; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldL, &temperatureL, eulerAdvectionData->computeTemperature.context.get())); + PetscCall(eulerAdvectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + auxR[aOff[0]] * .33, &temperatureL, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityL = 0.0; @@ -302,7 +306,7 @@ PetscErrorCode ablate::finiteVolume::processes::SpeciesTransport::AdvectionFlux( densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureR; - PetscCall(eulerAdvectionData->computeTemperature.function(fieldR, &temperatureR, eulerAdvectionData->computeTemperature.context.get())); + PetscCall(eulerAdvectionData->computeTemperature.function(fieldR, auxL[aOff[0]] * .33 + auxR[aOff[0]] * .67, &temperatureR, eulerAdvectionData->computeTemperature.context.get())); // Get the velocity in this direction normalVelocityR = 0.0; diff --git a/src/finiteVolume/processes/speciesTransport.hpp b/src/finiteVolume/processes/speciesTransport.hpp index 8f2863124..3fabd2575 100644 --- a/src/finiteVolume/processes/speciesTransport.hpp +++ b/src/finiteVolume/processes/speciesTransport.hpp @@ -19,7 +19,7 @@ class SpeciesTransport : public FlowProcess { PetscInt numberSpecies; // EOS function calls - eos::ThermodynamicFunction computeTemperature; + eos::ThermodynamicTemperatureFunction computeTemperature; eos::ThermodynamicTemperatureFunction computeInternalEnergy; eos::ThermodynamicTemperatureFunction computeSpeedOfSound; eos::ThermodynamicTemperatureFunction computePressure; @@ -87,7 +87,6 @@ class SpeciesTransport : public FlowProcess { */ static void NormalizeSpecies(TS ts, ablate::solver::Solver&); - private: /** * This computes the energy transfer for species diffusion flux for rhoE * f = "euler" diff --git a/src/finiteVolume/processes/thermophoreticDiffusion.cpp b/src/finiteVolume/processes/thermophoreticDiffusion.cpp index f1aebceff..4a5ea334d 100644 --- a/src/finiteVolume/processes/thermophoreticDiffusion.cpp +++ b/src/finiteVolume/processes/thermophoreticDiffusion.cpp @@ -9,17 +9,17 @@ ablate::finiteVolume::processes::ThermophoreticDiffusion::ThermophoreticDiffusio void ablate::finiteVolume::processes::ThermophoreticDiffusion::Setup(ablate::finiteVolume::FiniteVolumeSolver &flow) { flow.RegisterRHSFunction(ThermophoreticDiffusionEnergyFlux, &viscosityTemperatureFunction, - CompressibleFlowFields::EULER_FIELD, + {CompressibleFlowFields::EULER_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD, CompressibleFlowFields::DENSITY_PROGRESS_FIELD}, {CompressibleFlowFields::TEMPERATURE_FIELD}); flow.RegisterRHSFunction(ThermophoreticDiffusionVariableFlux, &viscosityTemperatureFunction, - CompressibleFlowFields::DENSITY_YI_FIELD, + {CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_YI_FIELD}, {CompressibleFlowFields::TEMPERATURE_FIELD}); flow.RegisterRHSFunction(ThermophoreticDiffusionVariableFlux, &viscosityTemperatureFunction, - CompressibleFlowFields::DENSITY_PROGRESS_FIELD, + {CompressibleFlowFields::DENSITY_PROGRESS_FIELD}, {CompressibleFlowFields::EULER_FIELD, CompressibleFlowFields::DENSITY_PROGRESS_FIELD}, {CompressibleFlowFields::TEMPERATURE_FIELD}); diff --git a/src/finiteVolume/processes/twoPhaseEulerAdvection.cpp b/src/finiteVolume/processes/twoPhaseEulerAdvection.cpp index b7247b0ad..8f045c18e 100644 --- a/src/finiteVolume/processes/twoPhaseEulerAdvection.cpp +++ b/src/finiteVolume/processes/twoPhaseEulerAdvection.cpp @@ -232,8 +232,8 @@ void ablate::finiteVolume::processes::TwoPhaseEulerAdvection::Setup(ablate::fini decoder = CreateTwoPhaseDecoder(flow.GetSubDomain().GetDimensions(), eosGas, eosLiquid); // Currently, no option for species advection - flow.RegisterRHSFunction(CompressibleFlowComputeEulerFlux, this, CompressibleFlowFields::EULER_FIELD, {VOLUME_FRACTION_FIELD, DENSITY_VF_FIELD, CompressibleFlowFields::EULER_FIELD}, {}); - flow.RegisterRHSFunction(CompressibleFlowComputeVFFlux, this, DENSITY_VF_FIELD, {VOLUME_FRACTION_FIELD, DENSITY_VF_FIELD, CompressibleFlowFields::EULER_FIELD}, {}); + flow.RegisterRHSFunction(CompressibleFlowComputeEulerFlux, this, {CompressibleFlowFields::EULER_FIELD}, {VOLUME_FRACTION_FIELD, DENSITY_VF_FIELD, CompressibleFlowFields::EULER_FIELD}, {}); + flow.RegisterRHSFunction(CompressibleFlowComputeVFFlux, this, {DENSITY_VF_FIELD}, {VOLUME_FRACTION_FIELD, DENSITY_VF_FIELD, CompressibleFlowFields::EULER_FIELD}, {}); flow.RegisterComputeTimeStepFunction(ComputeCflTimeStep, &timeStepData, "cfl"); timeStepData.computeSpeedOfSound = eosTwoPhase->GetThermodynamicFunction(eos::ThermodynamicProperty::SpeedOfSound, flow.GetSubDomain().GetFields()); diff --git a/tests/integrationTests/inputs/reactingFlow/sampleSootDiffusionFlameCompact2.yaml b/tests/integrationTests/inputs/reactingFlow/sampleSootDiffusionFlameCompact2.yaml new file mode 100644 index 000000000..47becb47c --- /dev/null +++ b/tests/integrationTests/inputs/reactingFlow/sampleSootDiffusionFlameCompact2.yaml @@ -0,0 +1,132 @@ +# Example 1D diffusion flame with soot and thermophoretic diffusion. The diffusion flame only uses +# the transport processes of the Navier–Stokes equations. +--- +test: + # a unique test name for this integration tests + name: sampleSootDiffusionFlameCompact2 + # create a default assert that compares the log file + assert: "inputs/reactingFlow/sampleSootDiffusionFlame.txt" + +environment: + title: _1DSampleSootDiffusionFlameTchemCompact + tagDirectory: false +arguments: { } +timestepper: + arguments: + ts_type: rk + ts_max_time: 1 + ts_max_steps: 50 + ts_dt: 1.0e-10 + ts_adapt_safety: 0.95 + ts_adapt_type: physicsConstrained + domain: !ablate::domain::BoxMeshBoundaryCells + name: simpleBoxField + faces: [ 50 ] + lower: [ 0.0 ] + upper: [ 0.01 ] + options: + dm_plex_hash_location: true + preModifiers: + # distribute the mesh across the mpi rank with ghost cells + - !ablate::domain::modifiers::DistributeWithGhostCells + ghostCellDepth: 2 + postModifiers: + - !ablate::domain::modifiers::GhostBoundaryCells + fields: + - !ablate::finiteVolume::CompressibleFlowFields + eos: !ablate::eos::TChemSoot &eos + mechFile: ../mechanisms/gri30.yml + options: + # set a minimum temperature for the chemical kinetics ode integration + thresholdTemperature: 560 + conservedFieldOptions: + petscfv_type: leastsquares + region: + name: domain + - !ablate::domain::FieldDescription + name: pressure + type: FV + location: aux + region: + name: domain + initialization: + - !ablate::finiteVolume::fieldFunctions::Euler + state: &flowFieldState + eos: *eos + temperature: !ablate::mathFunctions::Peak + startValues: [ 300 ] + peakValues: [ 2710 ] + endValues: [ 300 ] + start: 0.0 + peak: 0.007 + end: 0.01 + pressure: 101325.0 + velocity: "0.0" + other: !ablate::finiteVolume::fieldFunctions::MassFractions + &massFracs + eos: *eos + values: + - fieldName: O2 + field: !ablate::mathFunctions::Linear + startValues: [ 0.0 ] + endValues: [ 1.0 ] + end: .01 + - fieldName: C2H4 + field: !ablate::mathFunctions::Linear + startValues: [ 1.0 ] + endValues: [ 0.0 ] + end: .01 + - !ablate::finiteVolume::fieldFunctions::DensityMassFractions + state: *flowFieldState + + # Set the number density to zero + - fieldName: densityProgress + field: "0.0" + + +solvers: + - !ablate::finiteVolume::CompressibleFlowSolver + id: flowField + compact: 2 + region: + name: interiorCells + eos: *eos + transport: !ablate::eos::tChemSoot::SootSpeciesTransportModel + transport: !ablate::eos::transport::Sutherland + eos: *eos + evTransport: !ablate::eos::tChemSoot::SootProgressTransportModel + transport: !ablate::eos::transport::Sutherland + eos: *eos + monitors: + - !ablate::monitors::TimeStepMonitor + interval: 10 + - !ablate::monitors::CurveMonitor + interval: 25 + additionalProcesses: + - !ablate::finiteVolume::processes::Chemistry + eos: *eos + - !ablate::finiteVolume::processes::ConstantPressureFix + eos: *eos + pressure: 101325.0 + - !ablate::finiteVolume::processes::ThermophoreticDiffusion + transport: !ablate::eos::transport::Sutherland + eos: *eos + - !ablate::boundarySolver::BoundarySolver + id: walls + region: + name: boundaryCellsRight + fieldBoundary: + name: boundaryFaces + mergeFaces: true + processes: + - !ablate::boundarySolver::lodi::Inlet + eos: *eos + - !ablate::boundarySolver::BoundarySolver + id: slab boundary + region: + name: boundaryCellsLeft + fieldBoundary: + name: boundaryFaces + processes: + - !ablate::boundarySolver::lodi::Inlet + eos: *eos diff --git a/tests/integrationTests/inputs/reactingFlow/simpleReactingFlow_zerork_Compact.yaml b/tests/integrationTests/inputs/reactingFlow/simpleReactingFlow_zerork_Compact.yaml new file mode 100644 index 000000000..917122d70 --- /dev/null +++ b/tests/integrationTests/inputs/reactingFlow/simpleReactingFlow_zerork_Compact.yaml @@ -0,0 +1,73 @@ +--- +test: + # a unique test name for this integration tests + name: simpleReactingFlowZerorkCompact1 + # create a default assert that compares the log file + assert: "inputs/reactingFlow/simpleReactingFlow.zerork.txt" + +environment: + title: _simpleReactingFlowZerorkCompact + tagDirectory: false +arguments: + dm_plex_separate_marker: "" +timestepper: + name: theMainTimeStepper + io: + interval: 0 + arguments: + ts_type: rk + ts_max_time: 0.2 + ts_max_steps: 25 + ts_dt: 1E-6 + ts_adapt_type: none + domain: !ablate::domain::BoxMesh + name: simpleBoxField + faces: [ 10, 10 ] + lower: [ -0.1, -0.1 ] + upper: [ .1, .1 ] + boundary: [ "NONE", "NONE" ] + simplex: false + modifiers: + - !ablate::domain::modifiers::GhostBoundaryCells + - !ablate::domain::modifiers::DistributeWithGhostCells + fields: + - !ablate::finiteVolume::CompressibleFlowFields + eos: !ablate::eos::zerorkEOS &eos + reactionFile: ../mechanisms/gri30.inp + thermoFile: ../mechanisms/gri30.dat + options: + reactorType: ConstantPressure #zerork deafult is constant volume + initialization: + - &eulerField + fieldName: "euler" #for euler all components are in a single field + field: >- + 1.0, + sqrt(x*x+y*y) <.05 ? 1498029.067485712: -58970.06564527616, + 0.0, + 0.0 + - &densityYiField + fieldName: "densityYi" #H2,H,O,O2,OH,H2O,HO2,H2O2,C,CH,CH2,CH2(S),CH3,CH4,CO,CO2,HCO,CH2O,CH2OH,CH3O,CH3OH,C2H,C2H2,C2H3,C2H4,C2H5,C2H6,HCCO,CH2CO,HCCOH,N,NH,NH2,NH3,NNH,NO,NO2,N2O,HNO,CN,HCN,H2CN,HCNN,HCNO,HOCN,HNCO,NCO,N2,AR,C3H7,C3H8,CH2CHO,CH3CHO + field: 0,0,0,0.2,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,0 +solver: !ablate::finiteVolume::CompressibleFlowSolver + id: vortexFlowField + compact: 1 + parameters: {} + transport: + k: 0.0 + mu: 0.0 + additionalProcesses: + # add in the reaction processes + - !ablate::finiteVolume::processes::Chemistry + eos: *eos + boundaryConditions: + - !ablate::finiteVolume::boundaryConditions::EssentialGhost + boundaryName: "walls" + labelIds: [1, 2, 3, 4] + boundaryValue: *eulerField + - !ablate::finiteVolume::boundaryConditions::EssentialGhost + boundaryName: "walls" + labelIds: [1, 2, 3, 4] + boundaryValue: *densityYiField + monitors: + - !ablate::monitors::TimeStepMonitor + eos: *eos diff --git a/tests/unitTests/finiteVolume/processes/navierStokesTransportTests.cpp b/tests/unitTests/finiteVolume/processes/navierStokesTransportTests.cpp index b7bd3c303..10c160b9b 100644 --- a/tests/unitTests/finiteVolume/processes/navierStokesTransportTests.cpp +++ b/tests/unitTests/finiteVolume/processes/navierStokesTransportTests.cpp @@ -30,7 +30,7 @@ TEST_P(NavierStokesTransportFluxTestFixture, ShouldComputeCorrectFlux) { // set a perfect gas for testing auto eos = std::make_shared(std::make_shared()); auto eulerFieldMock = ablateTesting::domain::MockField::Create("euler", 3); - eulerFlowData.computeTemperature = eos->GetThermodynamicFunction(ablate::eos::ThermodynamicProperty::Temperature, {eulerFieldMock}); + eulerFlowData.computeTemperature = eos->GetThermodynamicTemperatureFunction(ablate::eos::ThermodynamicProperty::Temperature, {eulerFieldMock}); eulerFlowData.computeInternalEnergy = eos->GetThermodynamicTemperatureFunction(ablate::eos::ThermodynamicProperty::InternalSensibleEnergy, {eulerFieldMock}); eulerFlowData.computeSpeedOfSound = eos->GetThermodynamicTemperatureFunction(ablate::eos::ThermodynamicProperty::SpeedOfSound, {eulerFieldMock}); eulerFlowData.computePressure = eos->GetThermodynamicTemperatureFunction(ablate::eos::ThermodynamicProperty::Pressure, {eulerFieldMock}); @@ -42,7 +42,10 @@ TEST_P(NavierStokesTransportFluxTestFixture, ShouldComputeCorrectFlux) { // act std::vector computedFlux(params.expectedFlux.size()); PetscInt uOff[1] = {0}; - ablate::finiteVolume::processes::NavierStokesTransport::AdvectionFlux(params.area.size(), &faceGeom, uOff, ¶ms.xLeft[0], ¶ms.xRight[0], NULL, NULL, NULL, &computedFlux[0], &eulerFlowData); + PetscInt aOff[1] = {0}; + PetscReal TempGuess[1] = {300}; + ablate::finiteVolume::processes::NavierStokesTransport::AdvectionFlux( + params.area.size(), &faceGeom, uOff, ¶ms.xLeft[0], ¶ms.xRight[0], aOff, TempGuess, TempGuess, &computedFlux[0], &eulerFlowData); // assert for (std::size_t i = 0; i < params.expectedFlux.size(); i++) {