From f027dcd0733dcf099a262b3964becfd968264dc7 Mon Sep 17 00:00:00 2001 From: kolosret Date: Sun, 27 Apr 2025 12:30:06 -0700 Subject: [PATCH 1/6] added limiter --- src/domain/subDomain.hpp | 3 + src/finiteVolume/cellInterpolant.cpp | 72 +++++++++++++++------ src/finiteVolume/compressibleFlowFields.cpp | 5 +- src/finiteVolume/compressibleFlowFields.hpp | 4 +- src/finiteVolume/compressibleFlowSolver.cpp | 6 +- 5 files changed, 65 insertions(+), 25 deletions(-) diff --git a/src/domain/subDomain.hpp b/src/domain/subDomain.hpp index b41284c3e..1fdb1d523 100644 --- a/src/domain/subDomain.hpp +++ b/src/domain/subDomain.hpp @@ -36,6 +36,9 @@ class SubDomain : public io::Serializable { //! contains the DM field numbers for the fields in this DS, or NULL IS fieldMap; + //! contains the DM field numbers for the fields in this DS, or NULL + PetscReal maxGradient; + //! Keep track of fields that live in this subDomain; std::map fieldsByName; diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 9f9702f5c..4c13efdf0 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -285,7 +285,7 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter PetscFunctionBegin; PetscCall(DMPlexGetTreeChildren(dm, face, &numChildren, &children)); - if (numChildren) { + if (numChildren) { //if the tree contains children PetscInt c; for (c = 0; c < numChildren; c++) { @@ -295,27 +295,30 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter PetscCall(DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field, childFace, fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad)); } } - } else { - PetscScalar* ncx; - PetscFVCellGeom* ncg; - const PetscInt* fcells; - PetscInt ncell, d; - PetscReal v[3]; + } else { //if the tree doesn't contain children + PetscScalar* ncx; //neighbor cell centered values + PetscFVCellGeom* ncg; // neighbor cell geometry + const PetscInt* fcells; // cells attached to this face + PetscInt ncell, d; //neighbor cell and for loop index + PetscReal v[3]; //centr PetscCall(DMPlexGetSupport(dm, face, &fcells)); - ncell = cell == fcells[0] ? fcells[1] : fcells[0]; + ncell = cell == fcells[0] ? fcells[1] : fcells[0]; //figure out which cell is the neighbor and not this cell + // Read in the neighbor cell information if (field >= 0) { PetscCall(DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx)); } else { PetscCall(DMPlexPointLocalRead(dm, ncell, x, &ncx)); } PetscCall(DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg)); - DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v); + // Calculate the distance between the neighbor cell and this cell + DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v); // v_i = NeighborCentroid_i - ThisCentroid_i = dx_i for (d = 0; d < dof; ++d) { /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */ - PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v); - PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom; - + PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v); // denominator = \grad u \cdot dx + PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom; // f = 1/2 * (u_i+1-u_i)/(\Delta u from cell gradient dot with \delta x) + //What the above means is that if any cell face has no change, but there was ample enough change close to it such that + // the cell gradient dot with delta x is not 0, there is no limiting... i.e f = 0 and all limiters = 0 PetscCall(PetscLimiterLimit(lim, flim, &phi)); cellPhi[d] = PetscMin(cellPhi[d], phi); } @@ -420,17 +423,17 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: VecGetArrayRead(cellGeomVec, &cellGeometryArray); // create a temp work array - PetscReal* cellPhi; - DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi) >> utilities::PetscUtilities::checkError; + PetscReal* cellPhi; //Actual Limiter + DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi) >> utilities::PetscUtilities::checkError; //Size it up to be dof length of reals for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { PetscInt cell = cellRange.points ? cellRange.points[c] : c; - const PetscInt* cellFaces; - PetscScalar* cx; - PetscFVCellGeom* cg; - PetscScalar* cgrad; - PetscInt coneSize; + const PetscInt* cellFaces; //faces belonging to the cell + PetscScalar* cx; //cell solution/field values + PetscFVCellGeom* cg; //cell geometry + PetscScalar* cgrad; //cell gradient + PetscInt coneSize; //cell conectivity DMPlexGetConeSize(dm, cell, &coneSize) >> utilities::PetscUtilities::checkError; DMPlexGetCone(dm, cell, &cellFaces) >> utilities::PetscUtilities::checkError; @@ -450,11 +453,39 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field.id, cellFaces[f], faceRange.start, faceRange.end, cellPhi, xLocalArray, cellGeometryArray, cg, cx, cgrad) >> utilities::PetscUtilities::checkError; } + +// /* Apply limiter to gradient. By default, use the minimum limiter. */ +// PetscBool cancel = PETSC_FALSE; +// for (PetscInt pd = 0; pd < dof; ++pd) { +// if (cellPhi[2] == 0) cancel = PETSC_TRUE; +// } +// for (PetscInt pd = 0; pd < dof; ++pd) { +// /* Scalar limiter applied to each component separately */ +// for (PetscInt d = 0; d < dim; ++d) { +// if (cancel) cgrad[pd * dim + d] *= 0; +// else cgrad[pd * dim + d] *= cellPhi[pd]; +// } +// } + +// Gradlim; + /* Apply limiter to gradient */ + PetscBool cancel = PETSC_FALSE; + for (PetscInt pd = 0; pd < dof; ++pd) { + if (cellPhi[pd] ==0) { + for (PetscInt d =0; d < dim; d++) { + if (cgrad[pd * dim + d] > 10000) cancel = PETSC_TRUE; + } + } + } + for (PetscInt pd = 0; pd < dof; ++pd) { /* Scalar limiter applied to each component separately */ for (PetscInt d = 0; d < dim; ++d) { - cgrad[pd * dim + d] *= cellPhi[pd]; + if (cancel) + cgrad[pd * dim + d] *= 0; + else + cgrad[pd * dim + d] *= cellPhi[pd]; } } } @@ -687,6 +718,7 @@ static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, DMLabel region ++usedFaces; } } + //Free the memory allocated earlier with PetscMalloc3 PetscCall(PetscFree3(dx, grad, gref)); PetscFunctionReturn(0); } diff --git a/src/finiteVolume/compressibleFlowFields.cpp b/src/finiteVolume/compressibleFlowFields.cpp index 91c460af9..74abae3df 100644 --- a/src/finiteVolume/compressibleFlowFields.cpp +++ b/src/finiteVolume/compressibleFlowFields.cpp @@ -5,7 +5,7 @@ #include "utilities/vectorUtilities.hpp" ablate::finiteVolume::CompressibleFlowFields::CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region, - std::shared_ptr conservedFieldParameters) + std::shared_ptr conservedFieldParameters,double maxgradient) : eos(std::move(eos)), region(std::move(region)), conservedFieldOptions(std::move(conservedFieldParameters)) {} std::vector> ablate::finiteVolume::CompressibleFlowFields::GetFields() { @@ -78,4 +78,5 @@ std::istream& ablate::finiteVolume::operator>>(std::istream& is, ablate::finiteV #include "registrar.hpp" REGISTER(ablate::domain::FieldDescriptor, ablate::finiteVolume::CompressibleFlowFields, "FVM fields need for compressible flow", ARG(ablate::eos::EOS, "eos", "the equation of state to be used for the flow"), OPT(ablate::domain::Region, "region", "the region for the compressible flow (defaults to entire domain)"), - OPT(ablate::parameters::Parameters, "conservedFieldOptions", "petsc options used for the conserved fields. Common options would be petscfv_type and petsclimiter_type")); + OPT(ablate::parameters::Parameters, "conservedFieldOptions", "petsc options used for the conserved fields. Common options would be petscfv_type and petsclimiter_type"), + OPT(double, "maxgradient", "Maximum value of the gradient for the limiter in cellinterpolant")); diff --git a/src/finiteVolume/compressibleFlowFields.hpp b/src/finiteVolume/compressibleFlowFields.hpp index 4951ad7b4..5b415ba29 100644 --- a/src/finiteVolume/compressibleFlowFields.hpp +++ b/src/finiteVolume/compressibleFlowFields.hpp @@ -48,6 +48,8 @@ class CompressibleFlowFields : public domain::FieldDescriptor { inline const static std::string VELOCITY_FIELD = "velocity"; inline const static std::string PRESSURE_FIELD = "pressure"; + const std::shared_ptr limiterOptions; + protected: const std::shared_ptr eos; const std::shared_ptr region; @@ -61,7 +63,7 @@ class CompressibleFlowFields : public domain::FieldDescriptor { * @param region the region for all of the fields * @param conservedFieldParameters override the default field parameters for the conserved field */ - explicit CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region = {}, std::shared_ptr conservedFieldParameters = {}); + explicit CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region = {}, std::shared_ptr conservedFieldParameters = {},double maxgradient=0); /** * override and return the compressible flow fields diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index f5607d6a4..fc09a0dab 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -53,8 +53,10 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string #include "registrar.hpp" REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, "compressible finite volume flow", ARG(std::string, "id", "the name of the flow field"), - OPT(ablate::domain::Region, "region", "the region to apply this solver. Default is entire domain"), OPT(ablate::parameters::Parameters, "options", "the options passed to PETSc"), - ARG(ablate::eos::EOS, "eos", "the equation of state used to describe the flow"), OPT(ablate::parameters::Parameters, "parameters", "the parameters used for field values"), + OPT(ablate::domain::Region, "region", "the region to apply this solver. Default is entire domain"), + OPT(ablate::parameters::Parameters, "options", "the options passed to PETSc"), + ARG(ablate::eos::EOS, "eos", "the equation of state used to describe the flow"), + OPT(ablate::parameters::Parameters, "parameters", "the parameters used for field values"), OPT(ablate::eos::transport::TransportModel, "transport", "the diffusion transport model"), OPT(ablate::finiteVolume::fluxCalculator::FluxCalculator, "fluxCalculator", "the flux calculators (defaults to none)"), OPT(std::vector, "additionalProcesses", "any additional processes besides euler/yi/ev transport"), From 1d3830d31bc54f9a6cd00c5a81154dca202f365b Mon Sep 17 00:00:00 2001 From: kolosret Date: Mon, 28 Apr 2025 09:40:11 -0700 Subject: [PATCH 2/6] Maximum gradient for limiter is set from finitevolume solver, from the compressible flow solver register. --- src/finiteVolume/cellInterpolant.cpp | 22 ++++----------------- src/finiteVolume/cellInterpolant.hpp | 5 ++++- src/finiteVolume/compressibleFlowFields.cpp | 5 ++--- src/finiteVolume/compressibleFlowFields.hpp | 5 +++-- src/finiteVolume/compressibleFlowSolver.cpp | 9 ++++++--- src/finiteVolume/compressibleFlowSolver.hpp | 2 +- src/finiteVolume/finiteVolumeSolver.cpp | 4 ++-- src/finiteVolume/finiteVolumeSolver.hpp | 3 +++ 8 files changed, 25 insertions(+), 30 deletions(-) diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 4c13efdf0..6bdbfd52a 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -2,8 +2,8 @@ #include #include -ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptr subDomainIn, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec) - : subDomain(std::move(std::move(subDomainIn))) { +ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptr subDomainIn, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec,double maxGradIn) + : subDomain(std::move(std::move(subDomainIn))), maxLimGrad(maxGradIn) { auto getGradientDm = [this, solverRegion, faceGeomVec, cellGeomVec](const domain::Field& fieldInfo, std::vector& gradDMs) { auto petscField = subDomain->GetPetscFieldObject(fieldInfo); auto petscFieldFV = (PetscFV)petscField; @@ -454,27 +454,13 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: utilities::PetscUtilities::checkError; } -// /* Apply limiter to gradient. By default, use the minimum limiter. */ -// PetscBool cancel = PETSC_FALSE; -// for (PetscInt pd = 0; pd < dof; ++pd) { -// if (cellPhi[2] == 0) cancel = PETSC_TRUE; -// } -// for (PetscInt pd = 0; pd < dof; ++pd) { -// /* Scalar limiter applied to each component separately */ -// for (PetscInt d = 0; d < dim; ++d) { -// if (cancel) cgrad[pd * dim + d] *= 0; -// else cgrad[pd * dim + d] *= cellPhi[pd]; -// } -// } - -// Gradlim; - /* Apply limiter to gradient */ PetscBool cancel = PETSC_FALSE; for (PetscInt pd = 0; pd < dof; ++pd) { if (cellPhi[pd] ==0) { for (PetscInt d =0; d < dim; d++) { - if (cgrad[pd * dim + d] > 10000) cancel = PETSC_TRUE; + //limit the maxium gradient, this is required for strong shocks in rocket simulations + if (cgrad[pd * dim + d] > maxLimGrad) cancel = PETSC_TRUE; } } } diff --git a/src/finiteVolume/cellInterpolant.hpp b/src/finiteVolume/cellInterpolant.hpp index 6bfbdc7c9..15c8025ac 100644 --- a/src/finiteVolume/cellInterpolant.hpp +++ b/src/finiteVolume/cellInterpolant.hpp @@ -50,6 +50,9 @@ class CellInterpolant { //! store the dmGrad, these are specific to this finite volume solver std::vector gradientCellDms; + // Maximum value for gradients for the multi-direction flux limiter + const double maxLimGrad; + /** * Function to compute the flux source terms */ @@ -99,7 +102,7 @@ class CellInterpolant { * @param faceGeomVec * @param cellGeomVec */ - CellInterpolant(std::shared_ptr subDomain, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec); + CellInterpolant(std::shared_ptr subDomain, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec, double maxGradIn); ~CellInterpolant(); /** diff --git a/src/finiteVolume/compressibleFlowFields.cpp b/src/finiteVolume/compressibleFlowFields.cpp index 74abae3df..91c460af9 100644 --- a/src/finiteVolume/compressibleFlowFields.cpp +++ b/src/finiteVolume/compressibleFlowFields.cpp @@ -5,7 +5,7 @@ #include "utilities/vectorUtilities.hpp" ablate::finiteVolume::CompressibleFlowFields::CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region, - std::shared_ptr conservedFieldParameters,double maxgradient) + std::shared_ptr conservedFieldParameters) : eos(std::move(eos)), region(std::move(region)), conservedFieldOptions(std::move(conservedFieldParameters)) {} std::vector> ablate::finiteVolume::CompressibleFlowFields::GetFields() { @@ -78,5 +78,4 @@ std::istream& ablate::finiteVolume::operator>>(std::istream& is, ablate::finiteV #include "registrar.hpp" REGISTER(ablate::domain::FieldDescriptor, ablate::finiteVolume::CompressibleFlowFields, "FVM fields need for compressible flow", ARG(ablate::eos::EOS, "eos", "the equation of state to be used for the flow"), OPT(ablate::domain::Region, "region", "the region for the compressible flow (defaults to entire domain)"), - OPT(ablate::parameters::Parameters, "conservedFieldOptions", "petsc options used for the conserved fields. Common options would be petscfv_type and petsclimiter_type"), - OPT(double, "maxgradient", "Maximum value of the gradient for the limiter in cellinterpolant")); + OPT(ablate::parameters::Parameters, "conservedFieldOptions", "petsc options used for the conserved fields. Common options would be petscfv_type and petsclimiter_type")); diff --git a/src/finiteVolume/compressibleFlowFields.hpp b/src/finiteVolume/compressibleFlowFields.hpp index 5b415ba29..9cf6db5b9 100644 --- a/src/finiteVolume/compressibleFlowFields.hpp +++ b/src/finiteVolume/compressibleFlowFields.hpp @@ -48,7 +48,7 @@ class CompressibleFlowFields : public domain::FieldDescriptor { inline const static std::string VELOCITY_FIELD = "velocity"; inline const static std::string PRESSURE_FIELD = "pressure"; - const std::shared_ptr limiterOptions; + protected: const std::shared_ptr eos; @@ -56,6 +56,7 @@ class CompressibleFlowFields : public domain::FieldDescriptor { const std::shared_ptr conservedFieldOptions; const std::shared_ptr auxFieldOptions = ablate::parameters::MapParameters::Create({{"petscfv_type", "leastsquares"}, {"petsclimiter_type", "none"}}); + public: /** * Create a helper class that produces the required compressible flow fields based upon the eos and specifed region @@ -63,7 +64,7 @@ class CompressibleFlowFields : public domain::FieldDescriptor { * @param region the region for all of the fields * @param conservedFieldParameters override the default field parameters for the conserved field */ - explicit CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region = {}, std::shared_ptr conservedFieldParameters = {},double maxgradient=0); + explicit CompressibleFlowFields(std::shared_ptr eos, std::shared_ptr region = {}, std::shared_ptr conservedFieldParameters = {}); /** * override and return the compressible flow fields diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index fc09a0dab..f490a643e 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -20,7 +20,7 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string const std::shared_ptr& fluxCalculatorIn, std::vector> additionalProcesses, std::vector> boundaryConditions, - const std::shared_ptr& evTransport, int compact) + const std::shared_ptr& evTransport, int compact, double maxLimGrad) : FiniteVolumeSolver( std::move(solverId), std::move(region), std::move(options), compact ? (compact == 1 ? utilities::VectorUtilities::Merge({std::make_shared( @@ -41,7 +41,9 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string std::make_shared(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport), }, additionalProcesses), - std::move(boundaryConditions)) {} + std::move(boundaryConditions)) { + ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit=maxLimGrad; +} 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, @@ -62,4 +64,5 @@ REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, " 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"), - 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 + 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)"), + OPT(int, "maxLimGrad", "Maximum gradient in any direction. If Cellinterpolant sees larger gradient, the limiter is set to 0.")); \ No newline at end of file diff --git a/src/finiteVolume/compressibleFlowSolver.hpp b/src/finiteVolume/compressibleFlowSolver.hpp index 696e70765..4d1110ed7 100644 --- a/src/finiteVolume/compressibleFlowSolver.hpp +++ b/src/finiteVolume/compressibleFlowSolver.hpp @@ -32,7 +32,7 @@ class CompressibleFlowSolver : public FiniteVolumeSolver { 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 = {}, - int compact = 0); + int compact = 0, double maxLimGrad = 1E8); /** * Constructor without ev or additional processes diff --git a/src/finiteVolume/finiteVolumeSolver.cpp b/src/finiteVolume/finiteVolumeSolver.cpp index 9883e0e42..eacf0ef1c 100644 --- a/src/finiteVolume/finiteVolumeSolver.cpp +++ b/src/finiteVolume/finiteVolumeSolver.cpp @@ -239,7 +239,7 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets StartEvent("FiniteVolumeSolver::ComputeRHSFunction::discontinuousFluxFunction"); if (!discontinuousFluxFunctionDescriptions.empty()) { if (cellInterpolant == nullptr) { - cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec); + cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit); } cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), discontinuousFluxFunctionDescriptions, faceRange, cellRange, cellGeomVec, faceGeomVec); @@ -253,7 +253,7 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets StartEvent("FiniteVolumeSolver::ComputeRHSFunction::pointFunction"); if (!pointFunctionDescriptions.empty()) { if (cellInterpolant == nullptr) { - cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec); + cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec,maxlimit); } cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), pointFunctionDescriptions, cellRange, cellGeomVec); diff --git a/src/finiteVolume/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index fd40dc2cf..9c5b57813 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -79,6 +79,9 @@ class FiniteVolumeSolver : public solver::CellSolver, //! Store a dm, vec and array for mesh characteristics specific to the fvm Vec meshCharacteristicsLocalVec = nullptr; + protected: + double maxlimit; + public: FiniteVolumeSolver(std::string solverId, std::shared_ptr, std::shared_ptr options, std::vector> flowProcesses, std::vector> boundaryConditions); From cba4b8b4389fd99c5b31d77d62c603df05db8484 Mon Sep 17 00:00:00 2001 From: kolosret Date: Mon, 28 Apr 2025 10:07:52 -0700 Subject: [PATCH 3/6] Updated timestepper for accurate profiling results. --- src/solver/timeStepper.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/solver/timeStepper.cpp b/src/solver/timeStepper.cpp index 803d89089..92b573766 100644 --- a/src/solver/timeStepper.cpp +++ b/src/solver/timeStepper.cpp @@ -363,6 +363,9 @@ PetscErrorCode ablate::solver::TimeStepper::SolverComputeRHSFunction(TS ts, Pets VecZeroEntries(locF); CHKMEMQ; + //Sync all the ranks here, to avoid imbalance form PreRHS calls (eg. radiative gain) to distort computeRHS profiling results + MPI_Barrier(PETSC_COMM_WORLD); + // Call each of the provided RHS functions timeStepper->StartEvent("SolverComputeRHSFunction::ComputeRHSFunction"); for (auto& solver : timeStepper->rhsFunctionSolvers) { From 7c108e32bff2fcd11b556f4d40a84960af59bcf1 Mon Sep 17 00:00:00 2001 From: kolosret Date: Mon, 28 Apr 2025 10:19:13 -0700 Subject: [PATCH 4/6] Patch from Kenny containing the massflux boundary condition. --- src/boundarySolver/lodi/CMakeLists.txt | 3 + src/boundarySolver/lodi/lodiBoundary.cpp | 12 +- src/boundarySolver/lodi/massFluxInlet.cpp | 209 ++++++++++++++++++++++ src/boundarySolver/lodi/massFluxInlet.hpp | 27 +++ 4 files changed, 249 insertions(+), 2 deletions(-) create mode 100644 src/boundarySolver/lodi/massFluxInlet.cpp create mode 100644 src/boundarySolver/lodi/massFluxInlet.hpp diff --git a/src/boundarySolver/lodi/CMakeLists.txt b/src/boundarySolver/lodi/CMakeLists.txt index 23afd5e8d..695ffe748 100644 --- a/src/boundarySolver/lodi/CMakeLists.txt +++ b/src/boundarySolver/lodi/CMakeLists.txt @@ -4,10 +4,13 @@ target_sources(ablateLibrary isothermalWall.cpp openBoundary.cpp inlet.cpp + massFluxInlet.cpp + PUBLIC lodiBoundary.hpp isothermalWall.hpp openBoundary.hpp inlet.hpp + massFluxInlet.hpp ) diff --git a/src/boundarySolver/lodi/lodiBoundary.cpp b/src/boundarySolver/lodi/lodiBoundary.cpp index ee8da5081..d7ce65186 100644 --- a/src/boundarySolver/lodi/lodiBoundary.cpp +++ b/src/boundarySolver/lodi/lodiBoundary.cpp @@ -72,6 +72,7 @@ void ablate::boundarySolver::lodi::LODIBoundary::GetmdFdn(const PetscInt sOff[], } mdFdn[sOff[eulerId] + RHO] = -d[0]; mdFdn[sOff[eulerId] + RHOVELN] = -(vel[0] * d[0] + rho * d[2]); // Wall normal component momentum, not really rho u + double KE = vel[0] * vel[0]; double dvelterm = vel[0] * d[2]; for (int ndim = 1; ndim < dims; ndim++) { // Tangential components for momentum @@ -81,11 +82,11 @@ void ablate::boundarySolver::lodi::LODIBoundary::GetmdFdn(const PetscInt sOff[], } KE = 0.5e+0 * KE; mdFdn[sOff[eulerId] + RHOE] = -(d[0] * (KE + Enth - Cp * T) + d[1] / (Cp / Cv - 1.e+0 + 1.0E-30) + rho * dvelterm); + for (int ns = 0; ns < nSpecEqs; ns++) { const PetscReal *rhoYi = conserved + uOff[speciesId]; mdFdn[sOff[speciesId] + ns] = -(rhoYi[ns] / rho * d[0] + rho * d[2 + dims + ns]); // species } - int ne = 0; for (std::size_t ev = 0; ev < evIds.size(); ++ev) { const PetscReal *rhoEV = conserved + uOff[evIds[ev]]; @@ -117,7 +118,7 @@ void ablate::boundarySolver::lodi::LODIBoundary::Setup(ablate::boundarySolver::B // Compute the number of equations that need to be solved dims = bSolver.GetSubDomain().GetDimensions(); - // check if the eos need to do any updates + // check if the eos need to do any updates (I.e Chemtab and converting P.V's to Yi's auto chemModel = std::dynamic_pointer_cast(eos); if (chemModel) { for (auto &updateFunction : chemModel->GetSolutionFieldUpdates()) { @@ -144,6 +145,13 @@ void ablate::boundarySolver::lodi::LODIBoundary::Setup(ablate::boundarySolver::B std::vector{finiteVolume::CompressibleFlowFields::TEMPERATURE_FIELD}, {}); } + if (bSolver.GetSubDomain().ContainsField(finiteVolume::CompressibleFlowFields::PRESSURE_FIELD)) { + // add in aux update variables + bSolver.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxPressureField, + &computePressure, + std::vector{finiteVolume::CompressibleFlowFields::PRESSURE_FIELD}, + {}); + } } if (bSolver.GetSubDomain().ContainsField(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD)) { nSpecEqs = bSolver.GetSubDomain().GetField(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD).numberComponents; diff --git a/src/boundarySolver/lodi/massFluxInlet.cpp b/src/boundarySolver/lodi/massFluxInlet.cpp new file mode 100644 index 000000000..522bf30ac --- /dev/null +++ b/src/boundarySolver/lodi/massFluxInlet.cpp @@ -0,0 +1,209 @@ +#include "massFluxInlet.hpp" +#include +#include "finiteVolume/compressibleFlowFields.hpp" +#include "mathFunctions/functionFactory.hpp" + +using fp = ablate::finiteVolume::CompressibleFlowFields; + +ablate::boundarySolver::lodi::MassFluxInlet::MassFluxInlet(std::shared_ptr eos, std::shared_ptr pressureGradientScaling, + std::shared_ptr prescribedMassFlux) + : LODIBoundary(std::move(eos), std::move(pressureGradientScaling)), prescribedMassFlux(std::move(prescribedMassFlux)) {} + +void ablate::boundarySolver::lodi::MassFluxInlet::Setup(ablate::boundarySolver::BoundarySolver &bSolver) { + ablate::boundarySolver::lodi::LODIBoundary::Setup(bSolver); + + bSolver.RegisterFunction(InletFunction, this, fieldNames, fieldNames, {}); + + // Register a pre function step to update mass flux over this solver if specified + if (prescribedMassFlux) { + // define an update field function + auto updateFieldFunction = + std::make_shared(finiteVolume::CompressibleFlowFields::EULER_FIELD, ablate::mathFunctions::Create(UpdateMassFluxFunction, prescribedMassFlux.get())); + + bSolver.RegisterPreStep([&bSolver, updateFieldFunction](auto ts, auto &solver) { + // Get the current time + PetscReal time; + TSGetTime(ts, &time) >> utilities::PetscUtilities::checkError; + + bSolver.InsertFieldFunctions({updateFieldFunction}, time); + }); + } +} +PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::InletFunction(PetscInt dim, const ablate::boundarySolver::BoundarySolver::BoundaryFVFaceGeom *fg, const PetscFVCellGeom *boundaryCell, + const PetscInt *uOff, const PetscScalar *boundaryValues, const PetscScalar **stencilValues, const PetscInt *aOff, + const PetscScalar *auxValues, const PetscScalar **stencilAuxValues, PetscInt stencilSize, const PetscInt *stencil, + const PetscScalar *stencilWeights, const PetscInt *sOff, PetscScalar *source, void *ctx) { + PetscFunctionBeginUser; + auto inletBoundary = (MassFluxInlet *)ctx; + + // Compute the transformation matrix + PetscReal transformationMatrix[3][3]; + utilities::MathUtilities::ComputeTransformationMatrix(dim, fg->normal, transformationMatrix); + + // Compute the pressure/values on the boundary + PetscReal boundaryDensity; + PetscReal boundaryTemperature; + PetscReal boundaryVel[3]; + PetscReal boundaryNormalVelocity = 0.0; + PetscReal boundarySpeedOfSound; + PetscReal boundaryPressure; + + // Get the densityYi pointer if available + const PetscScalar *boundaryDensityYi = inletBoundary->nSpecEqs > 0 ? boundaryValues + uOff[inletBoundary->speciesId] : nullptr; + + // Get the velocity and pressure on the surface + { + boundaryDensity = boundaryValues[uOff[inletBoundary->eulerId] + finiteVolume::CompressibleFlowFields::RHO]; + for (PetscInt d = 0; d < dim; d++) { + boundaryVel[d] = boundaryValues[uOff[inletBoundary->eulerId] + finiteVolume::CompressibleFlowFields::RHOU + d] / boundaryDensity; + boundaryNormalVelocity += boundaryVel[d] * fg->normal[d]; + } + PetscCall(inletBoundary->computeTemperature.function(boundaryValues, &boundaryTemperature, inletBoundary->computeTemperature.context.get())); + PetscCall(inletBoundary->computeSpeedOfSound.function(boundaryValues, boundaryTemperature, &boundarySpeedOfSound, inletBoundary->computeSpeedOfSound.context.get())); + PetscCall(inletBoundary->computePressureFromTemperature.function(boundaryValues, boundaryTemperature, &boundaryPressure, inletBoundary->computePressureFromTemperature.context.get())); + } + + // Map the boundary velocity into the normal coord system + PetscReal boundaryVelNormCord[3]; + utilities::MathUtilities::Multiply(dim, transformationMatrix, boundaryVel, boundaryVelNormCord); + + // Compute each stencil point + std::vector stencilDensity(stencilSize); + std::vector> stencilVel(stencilSize, std::vector(dim)); + std::vector stencilNormalVelocity(stencilSize); + std::vector stencilPressure(stencilSize); + + for (PetscInt s = 0; s < stencilSize; s++) { + stencilDensity[s] = stencilValues[s][uOff[inletBoundary->eulerId] + finiteVolume::CompressibleFlowFields::RHO]; + for (PetscInt d = 0; d < dim; d++) { + stencilVel[s][d] = stencilValues[s][uOff[inletBoundary->eulerId] + finiteVolume::CompressibleFlowFields::RHOU + d] / stencilDensity[s]; + stencilNormalVelocity[s] += stencilVel[s][d] * fg->normal[d]; + } + PetscCall(inletBoundary->computePressure.function(stencilValues[s], &stencilPressure[s], inletBoundary->computePressure.context.get())); + } + + // Interpolate the normal velocity gradient to the surface + PetscScalar dVeldNorm; + BoundarySolver::ComputeGradientAlongNormal(dim, fg, boundaryNormalVelocity, stencilSize, &stencilNormalVelocity[0], stencilWeights, dVeldNorm); + PetscScalar dPdNorm; + BoundarySolver::ComputeGradientAlongNormal(dim, fg, boundaryPressure, stencilSize, &stencilPressure[0], stencilWeights, dPdNorm); + + // Compute the cp, cv from the eos + std::vector boundaryYi(inletBoundary->nSpecEqs); + for (PetscInt i = 0; i < inletBoundary->nSpecEqs; i++) { + boundaryYi[i] = boundaryDensityYi[i] / boundaryDensity; + } + PetscReal boundaryCp, boundaryCv; + inletBoundary->computeSpecificHeatConstantPressure.function(boundaryValues, boundaryTemperature, &boundaryCp, inletBoundary->computeSpecificHeatConstantPressure.context.get()); + inletBoundary->computeSpecificHeatConstantVolume.function(boundaryValues, boundaryTemperature, &boundaryCv, inletBoundary->computeSpecificHeatConstantVolume.context.get()); + + // Compute the enthalpy + PetscReal boundarySensibleEnthalpy; + inletBoundary->computeSensibleEnthalpyFunction.function(boundaryValues, boundaryTemperature, &boundarySensibleEnthalpy, inletBoundary->computeSensibleEnthalpyFunction.context.get()); + + // get_vel_and_c_prims(PGS, velwall[0], C, Cp, Cv, velnprm, Cprm); + PetscReal velNormPrim, speedOfSoundPrim; + inletBoundary->GetVelAndCPrims(boundaryNormalVelocity, boundarySpeedOfSound, boundaryCp, boundaryCv, velNormPrim, speedOfSoundPrim); + + // get_eigenvalues + std::vector lambda(inletBoundary->nEqs); + inletBoundary->GetEigenValues(boundaryNormalVelocity, boundarySpeedOfSound, velNormPrim, speedOfSoundPrim, &lambda[0]); + + // Get alpha + PetscReal pgsAlpha = inletBoundary->pressureGradientScaling ? inletBoundary->pressureGradientScaling->GetAlpha() : 1.0; + + // Get scriptL + std::vector scriptL(inletBoundary->nEqs); + // Outgoing acoustic wave + scriptL[1 + dim] = lambda[1 + dim] * (dPdNorm - boundaryDensity * PetscSqr(pgsAlpha) * dVeldNorm * (velNormPrim - boundaryNormalVelocity - speedOfSoundPrim)); + + //The Following comes from Simit's Massflux_NSCBC sideset, but i assume it's just what falls out from the constant + //mass flux isothermal lodi BC w/ pgs scaling + PetscReal gam,F; + gam = boundaryCp/boundaryCv; + F = gam*boundaryNormalVelocity*speedOfSoundPrim; + F = F/(boundarySpeedOfSound*boundarySpeedOfSound/PetscSqr(pgsAlpha) + + gam*boundaryNormalVelocity * (velNormPrim-boundaryNormalVelocity) ); + F = (1.+F)/(1.-F); + // Incoming acoustic wave + scriptL[0] = F*scriptL[1 + dim]; + // Entropy wave + scriptL[1] = 0.5e+0 * (boundaryCp / boundaryCv - 1.e+0) * (scriptL[1 + dim] + scriptL[0]) - + 0.5 * (boundaryCp / boundaryCv + 1.e+0) * (scriptL[0] - scriptL[1 + dim]) * (velNormPrim - boundaryNormalVelocity) / speedOfSoundPrim; + + // Tangential velocities + for (int d = 1; d < dim; d++) { + scriptL[1 + d] = 0.e+0; + } + for (int ns = 0; ns < inletBoundary->nSpecEqs; ns++) { + // Species + scriptL[2 + dim + ns] = 0.e+0; + } + for (int ne = 0; ne < inletBoundary->nEvEqs; ne++) { + // Extra variables + scriptL[2 + dim + inletBoundary->nSpecEqs + ne] = 0.e+0; + } + + // Directly compute the source terms, note that this may be problem in the future with multiple source terms on the same boundary cell + inletBoundary->GetmdFdn(sOff, + boundaryVelNormCord, + boundaryDensity, + boundaryTemperature, + boundaryCp, + boundaryCv, + boundarySpeedOfSound, + boundarySensibleEnthalpy, + velNormPrim, + speedOfSoundPrim, + boundaryValues, + uOff, + scriptL.data(), + transformationMatrix, + source); + + PetscFunctionReturn(0); +} + +PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::UpdateMassFluxFunction(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, void *ctx) { + PetscFunctionBeginUser; + + auto massFluxFunction = (ablate::mathFunctions::MathFunction *)ctx; + + // Get the current velocity + PetscScalar massFluxArray[3]; + PetscScalar kineticEnergy = 0.0; + PetscScalar density = u[fp::RHO]; + for (PetscInt d = 0; d < dim; d++) { + massFluxArray[d] = u[fp::RHOU + d]; + kineticEnergy += PetscSqr(massFluxArray[d]/ density); + } + kineticEnergy *= 0.5; + + // Get the internal energy + PetscScalar internalEnergy = u[fp::RHOE] / density; + + // Compute the sensible energy + PetscScalar sensibleEnergy = internalEnergy - kineticEnergy; + + // Update velocity + massFluxFunction->GetPetscFunction()(dim, time, x, dim, massFluxArray, massFluxFunction->GetContext()) >> utilities::PetscUtilities::checkError; + + // Update the momentum terms + kineticEnergy = 0.0; + for (PetscInt d = 0; d < dim; d++) { + u[fp::RHOU + d] = massFluxArray[d]; + kineticEnergy += PetscSqr(massFluxArray[d]/density); + } + kineticEnergy *= 0.5; + + // Update the new internal energy + u[fp::RHOE] = (sensibleEnergy + kineticEnergy) * density; + + PetscFunctionReturn(0); +} + +#include "registrar.hpp" +REGISTER(ablate::boundarySolver::BoundaryProcess, ablate::boundarySolver::lodi::MassFluxInlet, "Enforces an inlet with specified velocity", + ARG(ablate::eos::EOS, "eos", "The EOS describing the flow field at the wall"), + OPT(ablate::finiteVolume::processes::PressureGradientScaling, "pgs", "Pressure gradient scaling is used to scale the acoustic propagation speed and increase time step for low speed flows"), + OPT(ablate::mathFunctions::MathFunction, "massFlux", "optional velocity function that can change over time")); diff --git a/src/boundarySolver/lodi/massFluxInlet.hpp b/src/boundarySolver/lodi/massFluxInlet.hpp new file mode 100644 index 000000000..101d7eb69 --- /dev/null +++ b/src/boundarySolver/lodi/massFluxInlet.hpp @@ -0,0 +1,27 @@ +#ifndef ABLATELIBRARY_LODI_MASSFLUXINLET_HPP +#define ABLATELIBRARY_LODI_MASSFLUXINLET_HPP + +#include "lodiBoundary.hpp" + +namespace ablate::boundarySolver::lodi { +//This class assumes Isothermal inlet with constant mass flux +class MassFluxInlet : public LODIBoundary { + private: + //! The prescribed velocity in normal cartesian components (vx, vy, vz) + std::shared_ptr prescribedMassFlux; + static PetscErrorCode UpdateMassFluxFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar* u, void* ctx); + + public: + explicit MassFluxInlet(std::shared_ptr eos, std::shared_ptr pressureGradientScaling = {}, + std::shared_ptr prescribedVelocity = {}); + + void Setup(ablate::boundarySolver::BoundarySolver& bSolver) override; + + static PetscErrorCode InletFunction(PetscInt dim, const boundarySolver::BoundarySolver::BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[], + const PetscScalar* boundaryValues, const PetscScalar* stencilValues[], const PetscInt aOff[], const PetscScalar* auxValues, + const PetscScalar* stencilAuxValues[], PetscInt stencilSize, const PetscInt stencil[], const PetscScalar stencilWeights[], const PetscInt sOff[], + PetscScalar source[], void* ctx); +}; + +} // namespace ablate::boundarySolver::lodi +#endif // ABLATELIBRARY_MASSFLUXINLET_HPP From c2f4742d2610990d9eca10dbf09f6f94bcb8cbe2 Mon Sep 17 00:00:00 2001 From: kolosret Date: Mon, 28 Apr 2025 11:10:28 -0700 Subject: [PATCH 5/6] Version bump --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d3e7a48c4..e27377f5c 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.02) +project(ablateLibrary VERSION 0.13.03) # Load the Required 3rd Party Libaries pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc) From 326edca3c481a3406ed0bdce250c6892f07061ac Mon Sep 17 00:00:00 2001 From: kolosret Date: Wed, 30 Apr 2025 13:34:02 -0700 Subject: [PATCH 6/6] Correct formatting --- src/boundarySolver/lodi/lodiBoundary.cpp | 6 +-- src/boundarySolver/lodi/massFluxInlet.cpp | 29 ++++++------ src/boundarySolver/lodi/massFluxInlet.hpp | 4 +- src/finiteVolume/cellInterpolant.cpp | 51 +++++++++++---------- src/finiteVolume/compressibleFlowFields.hpp | 3 -- src/finiteVolume/compressibleFlowSolver.cpp | 8 ++-- src/finiteVolume/finiteVolumeSolver.cpp | 2 +- src/solver/timeStepper.cpp | 2 +- 8 files changed, 49 insertions(+), 56 deletions(-) diff --git a/src/boundarySolver/lodi/lodiBoundary.cpp b/src/boundarySolver/lodi/lodiBoundary.cpp index d7ce65186..06598b8fe 100644 --- a/src/boundarySolver/lodi/lodiBoundary.cpp +++ b/src/boundarySolver/lodi/lodiBoundary.cpp @@ -147,10 +147,8 @@ void ablate::boundarySolver::lodi::LODIBoundary::Setup(ablate::boundarySolver::B } if (bSolver.GetSubDomain().ContainsField(finiteVolume::CompressibleFlowFields::PRESSURE_FIELD)) { // add in aux update variables - bSolver.RegisterAuxFieldUpdate(ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxPressureField, - &computePressure, - std::vector{finiteVolume::CompressibleFlowFields::PRESSURE_FIELD}, - {}); + bSolver.RegisterAuxFieldUpdate( + ablate::finiteVolume::processes::NavierStokesTransport::UpdateAuxPressureField, &computePressure, std::vector{finiteVolume::CompressibleFlowFields::PRESSURE_FIELD}, {}); } } if (bSolver.GetSubDomain().ContainsField(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD)) { diff --git a/src/boundarySolver/lodi/massFluxInlet.cpp b/src/boundarySolver/lodi/massFluxInlet.cpp index 522bf30ac..dd811f465 100644 --- a/src/boundarySolver/lodi/massFluxInlet.cpp +++ b/src/boundarySolver/lodi/massFluxInlet.cpp @@ -6,7 +6,7 @@ using fp = ablate::finiteVolume::CompressibleFlowFields; ablate::boundarySolver::lodi::MassFluxInlet::MassFluxInlet(std::shared_ptr eos, std::shared_ptr pressureGradientScaling, - std::shared_ptr prescribedMassFlux) + std::shared_ptr prescribedMassFlux) : LODIBoundary(std::move(eos), std::move(pressureGradientScaling)), prescribedMassFlux(std::move(prescribedMassFlux)) {} void ablate::boundarySolver::lodi::MassFluxInlet::Setup(ablate::boundarySolver::BoundarySolver &bSolver) { @@ -30,9 +30,9 @@ void ablate::boundarySolver::lodi::MassFluxInlet::Setup(ablate::boundarySolver:: } } PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::InletFunction(PetscInt dim, const ablate::boundarySolver::BoundarySolver::BoundaryFVFaceGeom *fg, const PetscFVCellGeom *boundaryCell, - const PetscInt *uOff, const PetscScalar *boundaryValues, const PetscScalar **stencilValues, const PetscInt *aOff, - const PetscScalar *auxValues, const PetscScalar **stencilAuxValues, PetscInt stencilSize, const PetscInt *stencil, - const PetscScalar *stencilWeights, const PetscInt *sOff, PetscScalar *source, void *ctx) { + const PetscInt *uOff, const PetscScalar *boundaryValues, const PetscScalar **stencilValues, const PetscInt *aOff, + const PetscScalar *auxValues, const PetscScalar **stencilAuxValues, PetscInt stencilSize, const PetscInt *stencil, + const PetscScalar *stencilWeights, const PetscInt *sOff, PetscScalar *source, void *ctx) { PetscFunctionBeginUser; auto inletBoundary = (MassFluxInlet *)ctx; @@ -117,16 +117,15 @@ PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::InletFunction(PetscI // Outgoing acoustic wave scriptL[1 + dim] = lambda[1 + dim] * (dPdNorm - boundaryDensity * PetscSqr(pgsAlpha) * dVeldNorm * (velNormPrim - boundaryNormalVelocity - speedOfSoundPrim)); - //The Following comes from Simit's Massflux_NSCBC sideset, but i assume it's just what falls out from the constant - //mass flux isothermal lodi BC w/ pgs scaling - PetscReal gam,F; - gam = boundaryCp/boundaryCv; - F = gam*boundaryNormalVelocity*speedOfSoundPrim; - F = F/(boundarySpeedOfSound*boundarySpeedOfSound/PetscSqr(pgsAlpha) + - gam*boundaryNormalVelocity * (velNormPrim-boundaryNormalVelocity) ); - F = (1.+F)/(1.-F); + // The Following comes from Simit's Massflux_NSCBC sideset, but i assume it's just what falls out from the constant + // mass flux isothermal lodi BC w/ pgs scaling + PetscReal gam, F; + gam = boundaryCp / boundaryCv; + F = gam * boundaryNormalVelocity * speedOfSoundPrim; + F = F / (boundarySpeedOfSound * boundarySpeedOfSound / PetscSqr(pgsAlpha) + gam * boundaryNormalVelocity * (velNormPrim - boundaryNormalVelocity)); + F = (1. + F) / (1. - F); // Incoming acoustic wave - scriptL[0] = F*scriptL[1 + dim]; + scriptL[0] = F * scriptL[1 + dim]; // Entropy wave scriptL[1] = 0.5e+0 * (boundaryCp / boundaryCv - 1.e+0) * (scriptL[1 + dim] + scriptL[0]) - 0.5 * (boundaryCp / boundaryCv + 1.e+0) * (scriptL[0] - scriptL[1 + dim]) * (velNormPrim - boundaryNormalVelocity) / speedOfSoundPrim; @@ -175,7 +174,7 @@ PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::UpdateMassFluxFuncti PetscScalar density = u[fp::RHO]; for (PetscInt d = 0; d < dim; d++) { massFluxArray[d] = u[fp::RHOU + d]; - kineticEnergy += PetscSqr(massFluxArray[d]/ density); + kineticEnergy += PetscSqr(massFluxArray[d] / density); } kineticEnergy *= 0.5; @@ -192,7 +191,7 @@ PetscErrorCode ablate::boundarySolver::lodi::MassFluxInlet::UpdateMassFluxFuncti kineticEnergy = 0.0; for (PetscInt d = 0; d < dim; d++) { u[fp::RHOU + d] = massFluxArray[d]; - kineticEnergy += PetscSqr(massFluxArray[d]/density); + kineticEnergy += PetscSqr(massFluxArray[d] / density); } kineticEnergy *= 0.5; diff --git a/src/boundarySolver/lodi/massFluxInlet.hpp b/src/boundarySolver/lodi/massFluxInlet.hpp index 101d7eb69..d31d7ec22 100644 --- a/src/boundarySolver/lodi/massFluxInlet.hpp +++ b/src/boundarySolver/lodi/massFluxInlet.hpp @@ -4,7 +4,7 @@ #include "lodiBoundary.hpp" namespace ablate::boundarySolver::lodi { -//This class assumes Isothermal inlet with constant mass flux +// This class assumes Isothermal inlet with constant mass flux class MassFluxInlet : public LODIBoundary { private: //! The prescribed velocity in normal cartesian components (vx, vy, vz) @@ -13,7 +13,7 @@ class MassFluxInlet : public LODIBoundary { public: explicit MassFluxInlet(std::shared_ptr eos, std::shared_ptr pressureGradientScaling = {}, - std::shared_ptr prescribedVelocity = {}); + std::shared_ptr prescribedVelocity = {}); void Setup(ablate::boundarySolver::BoundarySolver& bSolver) override; diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 6bdbfd52a..120e346f0 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -2,7 +2,8 @@ #include #include -ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptr subDomainIn, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec,double maxGradIn) +ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptr subDomainIn, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec, + double maxGradIn) : subDomain(std::move(std::move(subDomainIn))), maxLimGrad(maxGradIn) { auto getGradientDm = [this, solverRegion, faceGeomVec, cellGeomVec](const domain::Field& fieldInfo, std::vector& gradDMs) { auto petscField = subDomain->GetPetscFieldObject(fieldInfo); @@ -285,7 +286,7 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter PetscFunctionBegin; PetscCall(DMPlexGetTreeChildren(dm, face, &numChildren, &children)); - if (numChildren) { //if the tree contains children + if (numChildren) { // if the tree contains children PetscInt c; for (c = 0; c < numChildren; c++) { @@ -295,15 +296,15 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter PetscCall(DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field, childFace, fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad)); } } - } else { //if the tree doesn't contain children - PetscScalar* ncx; //neighbor cell centered values - PetscFVCellGeom* ncg; // neighbor cell geometry - const PetscInt* fcells; // cells attached to this face - PetscInt ncell, d; //neighbor cell and for loop index - PetscReal v[3]; //centr + } else { // if the tree doesn't contain children + PetscScalar* ncx; // neighbor cell centered values + PetscFVCellGeom* ncg; // neighbor cell geometry + const PetscInt* fcells; // cells attached to this face + PetscInt ncell, d; // neighbor cell and for loop index + PetscReal v[3]; // centr PetscCall(DMPlexGetSupport(dm, face, &fcells)); - ncell = cell == fcells[0] ? fcells[1] : fcells[0]; //figure out which cell is the neighbor and not this cell + ncell = cell == fcells[0] ? fcells[1] : fcells[0]; // figure out which cell is the neighbor and not this cell // Read in the neighbor cell information if (field >= 0) { PetscCall(DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx)); @@ -312,13 +313,13 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter } PetscCall(DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg)); // Calculate the distance between the neighbor cell and this cell - DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v); // v_i = NeighborCentroid_i - ThisCentroid_i = dx_i + DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v); // v_i = NeighborCentroid_i - ThisCentroid_i = dx_i for (d = 0; d < dof; ++d) { /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */ - PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v); // denominator = \grad u \cdot dx - PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom; // f = 1/2 * (u_i+1-u_i)/(\Delta u from cell gradient dot with \delta x) - //What the above means is that if any cell face has no change, but there was ample enough change close to it such that - // the cell gradient dot with delta x is not 0, there is no limiting... i.e f = 0 and all limiters = 0 + PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v); // denominator = \grad u \cdot dx + PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom; // f = 1/2 * (u_i+1-u_i)/(\Delta u from cell gradient dot with \delta x) + // What the above means is that if any cell face has no change, but there was ample enough change close to it such that + // the cell gradient dot with delta x is not 0, there is no limiting... i.e f = 0 and all limiters = 0 PetscCall(PetscLimiterLimit(lim, flim, &phi)); cellPhi[d] = PetscMin(cellPhi[d], phi); } @@ -423,17 +424,17 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: VecGetArrayRead(cellGeomVec, &cellGeometryArray); // create a temp work array - PetscReal* cellPhi; //Actual Limiter - DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi) >> utilities::PetscUtilities::checkError; //Size it up to be dof length of reals + PetscReal* cellPhi; // Actual Limiter + DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi) >> utilities::PetscUtilities::checkError; // Size it up to be dof length of reals for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { PetscInt cell = cellRange.points ? cellRange.points[c] : c; - const PetscInt* cellFaces; //faces belonging to the cell - PetscScalar* cx; //cell solution/field values - PetscFVCellGeom* cg; //cell geometry - PetscScalar* cgrad; //cell gradient - PetscInt coneSize; //cell conectivity + const PetscInt* cellFaces; // faces belonging to the cell + PetscScalar* cx; // cell solution/field values + PetscFVCellGeom* cg; // cell geometry + PetscScalar* cgrad; // cell gradient + PetscInt coneSize; // cell conectivity DMPlexGetConeSize(dm, cell, &coneSize) >> utilities::PetscUtilities::checkError; DMPlexGetCone(dm, cell, &cellFaces) >> utilities::PetscUtilities::checkError; @@ -457,9 +458,9 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: /* Apply limiter to gradient */ PetscBool cancel = PETSC_FALSE; for (PetscInt pd = 0; pd < dof; ++pd) { - if (cellPhi[pd] ==0) { - for (PetscInt d =0; d < dim; d++) { - //limit the maxium gradient, this is required for strong shocks in rocket simulations + if (cellPhi[pd] == 0) { + for (PetscInt d = 0; d < dim; d++) { + // limit the maxium gradient, this is required for strong shocks in rocket simulations if (cgrad[pd * dim + d] > maxLimGrad) cancel = PETSC_TRUE; } } @@ -704,7 +705,7 @@ static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, DMLabel region ++usedFaces; } } - //Free the memory allocated earlier with PetscMalloc3 + // Free the memory allocated earlier with PetscMalloc3 PetscCall(PetscFree3(dx, grad, gref)); PetscFunctionReturn(0); } diff --git a/src/finiteVolume/compressibleFlowFields.hpp b/src/finiteVolume/compressibleFlowFields.hpp index 9cf6db5b9..4951ad7b4 100644 --- a/src/finiteVolume/compressibleFlowFields.hpp +++ b/src/finiteVolume/compressibleFlowFields.hpp @@ -48,15 +48,12 @@ class CompressibleFlowFields : public domain::FieldDescriptor { inline const static std::string VELOCITY_FIELD = "velocity"; inline const static std::string PRESSURE_FIELD = "pressure"; - - protected: const std::shared_ptr eos; const std::shared_ptr region; const std::shared_ptr conservedFieldOptions; const std::shared_ptr auxFieldOptions = ablate::parameters::MapParameters::Create({{"petscfv_type", "leastsquares"}, {"petsclimiter_type", "none"}}); - public: /** * Create a helper class that produces the required compressible flow fields based upon the eos and specifed region diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index f490a643e..e947c5d44 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -42,7 +42,7 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string }, additionalProcesses), std::move(boundaryConditions)) { - ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit=maxLimGrad; + ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit = maxLimGrad; } ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, @@ -55,10 +55,8 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string #include "registrar.hpp" REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, "compressible finite volume flow", ARG(std::string, "id", "the name of the flow field"), - OPT(ablate::domain::Region, "region", "the region to apply this solver. Default is entire domain"), - OPT(ablate::parameters::Parameters, "options", "the options passed to PETSc"), - ARG(ablate::eos::EOS, "eos", "the equation of state used to describe the flow"), - OPT(ablate::parameters::Parameters, "parameters", "the parameters used for field values"), + OPT(ablate::domain::Region, "region", "the region to apply this solver. Default is entire domain"), OPT(ablate::parameters::Parameters, "options", "the options passed to PETSc"), + ARG(ablate::eos::EOS, "eos", "the equation of state used to describe the flow"), OPT(ablate::parameters::Parameters, "parameters", "the parameters used for field values"), OPT(ablate::eos::transport::TransportModel, "transport", "the diffusion transport model"), OPT(ablate::finiteVolume::fluxCalculator::FluxCalculator, "fluxCalculator", "the flux calculators (defaults to none)"), OPT(std::vector, "additionalProcesses", "any additional processes besides euler/yi/ev transport"), diff --git a/src/finiteVolume/finiteVolumeSolver.cpp b/src/finiteVolume/finiteVolumeSolver.cpp index eacf0ef1c..6c21ed1be 100644 --- a/src/finiteVolume/finiteVolumeSolver.cpp +++ b/src/finiteVolume/finiteVolumeSolver.cpp @@ -253,7 +253,7 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets StartEvent("FiniteVolumeSolver::ComputeRHSFunction::pointFunction"); if (!pointFunctionDescriptions.empty()) { if (cellInterpolant == nullptr) { - cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec,maxlimit); + cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit); } cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), pointFunctionDescriptions, cellRange, cellGeomVec); diff --git a/src/solver/timeStepper.cpp b/src/solver/timeStepper.cpp index 92b573766..2ce8ebb7a 100644 --- a/src/solver/timeStepper.cpp +++ b/src/solver/timeStepper.cpp @@ -363,7 +363,7 @@ PetscErrorCode ablate::solver::TimeStepper::SolverComputeRHSFunction(TS ts, Pets VecZeroEntries(locF); CHKMEMQ; - //Sync all the ranks here, to avoid imbalance form PreRHS calls (eg. radiative gain) to distort computeRHS profiling results + // Sync all the ranks here, to avoid imbalance form PreRHS calls (eg. radiative gain) to distort computeRHS profiling results MPI_Barrier(PETSC_COMM_WORLD); // Call each of the provided RHS functions