From f027dcd0733dcf099a262b3964becfd968264dc7 Mon Sep 17 00:00:00 2001 From: kolosret Date: Sun, 27 Apr 2025 12:30:06 -0700 Subject: [PATCH 1/9] 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/9] 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/9] 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/9] 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/9] 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/9] 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 From e37d94607697265bd9b2577fa5ed5ee57f3b85ec Mon Sep 17 00:00:00 2001 From: kolosret Date: Wed, 30 Apr 2025 13:56:53 -0700 Subject: [PATCH 7/9] removing extra variable.. --- src/domain/subDomain.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/domain/subDomain.hpp b/src/domain/subDomain.hpp index 1fdb1d523..b41284c3e 100644 --- a/src/domain/subDomain.hpp +++ b/src/domain/subDomain.hpp @@ -36,9 +36,6 @@ 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; From e49c2a683bb0986184b0b427ce633bb16302054d Mon Sep 17 00:00:00 2001 From: kolosret Date: Tue, 6 May 2025 12:08:07 -0700 Subject: [PATCH 8/9] Bugfix: serializer returns 0 by default for optional arguments. Using the absolute value of the gradient... Setting maxgradient to ablate::utilities::constants::large in finitevolume solver --- src/finiteVolume/cellInterpolant.cpp | 8 ++++++-- src/finiteVolume/compressibleFlowSolver.cpp | 4 +++- src/finiteVolume/finiteVolumeSolver.hpp | 3 ++- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 120e346f0..4f69fd6c2 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -460,8 +460,12 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: 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 (cgrad[pd * dim + d] > maxLimGrad) cancel = PETSC_TRUE; + // Due to directional limiting being difficult in unstructured grids, + // a strict gradient limiter is introduced here to revert back to cell centered + // reconstructions if a component limiter is 0 even though there is a strong + // component gradient in a direction (Usually happens at strong shocks seen in + // high pressured rocket simulations) + if (PetscAbsReal(cgrad[pd * dim + d]) > maxLimGrad) cancel = PETSC_TRUE; } } } diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index e947c5d44..3d1505287 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -7,6 +7,7 @@ #include "finiteVolume/processes/navierStokesTransport.hpp" #include "finiteVolume/processes/speciesTransport.hpp" #include "utilities/vectorUtilities.hpp" +#include "utilities/constants.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. @@ -42,7 +43,8 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string }, additionalProcesses), std::move(boundaryConditions)) { - ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit = maxLimGrad; + // The user can input something small like 1E-10 but > 0 to turn of projections... + if(maxLimGrad!=0) ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit = maxLimGrad; } ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, diff --git a/src/finiteVolume/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index 9c5b57813..7716a41a2 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -12,6 +12,7 @@ #include "solver/solver.hpp" #include "solver/timeStepper.hpp" #include "utilities/vectorUtilities.hpp" +#include "utilities/constants.hpp" namespace ablate::finiteVolume { @@ -80,7 +81,7 @@ class FiniteVolumeSolver : public solver::CellSolver, Vec meshCharacteristicsLocalVec = nullptr; protected: - double maxlimit; + double maxlimit=ablate::utilities::Constants::large; public: FiniteVolumeSolver(std::string solverId, std::shared_ptr, std::shared_ptr options, std::vector> flowProcesses, From 5db80c427b746a03c9375dac92e232fd8c5ef8f7 Mon Sep 17 00:00:00 2001 From: kolosret Date: Wed, 7 May 2025 10:46:30 -0700 Subject: [PATCH 9/9] Final formatting, changing input type --- src/finiteVolume/compressibleFlowSolver.cpp | 6 +++--- src/finiteVolume/finiteVolumeSolver.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/finiteVolume/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index 3d1505287..9bbeacd4b 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -6,8 +6,8 @@ #include "finiteVolume/processes/evTransport.hpp" #include "finiteVolume/processes/navierStokesTransport.hpp" #include "finiteVolume/processes/speciesTransport.hpp" -#include "utilities/vectorUtilities.hpp" #include "utilities/constants.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. @@ -44,7 +44,7 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string additionalProcesses), std::move(boundaryConditions)) { // The user can input something small like 1E-10 but > 0 to turn of projections... - if(maxLimGrad!=0) ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit = maxLimGrad; + if (maxLimGrad != 0) ablate::finiteVolume::FiniteVolumeSolver::FiniteVolumeSolver::maxlimit = maxLimGrad; } ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr region, std::shared_ptr options, @@ -65,4 +65,4 @@ REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, " 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)"), - 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 + OPT(double, "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/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index 7716a41a2..88352e50a 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -11,8 +11,8 @@ #include "solver/cellSolver.hpp" #include "solver/solver.hpp" #include "solver/timeStepper.hpp" -#include "utilities/vectorUtilities.hpp" #include "utilities/constants.hpp" +#include "utilities/vectorUtilities.hpp" namespace ablate::finiteVolume { @@ -81,7 +81,7 @@ class FiniteVolumeSolver : public solver::CellSolver, Vec meshCharacteristicsLocalVec = nullptr; protected: - double maxlimit=ablate::utilities::Constants::large; + double maxlimit = ablate::utilities::Constants::large; public: FiniteVolumeSolver(std::string solverId, std::shared_ptr, std::shared_ptr options, std::vector> flowProcesses,