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) 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..06598b8fe 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,11 @@ 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..dd811f465 --- /dev/null +++ b/src/boundarySolver/lodi/massFluxInlet.cpp @@ -0,0 +1,208 @@ +#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..d31d7ec22 --- /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 diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 9f9702f5c..4f69fd6c2 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -2,8 +2,9 @@ #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; @@ -285,7 +286,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 +296,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 +424,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 +454,29 @@ 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 */ + PetscBool cancel = PETSC_FALSE; + for (PetscInt pd = 0; pd < dof; ++pd) { + if (cellPhi[pd] == 0) { + for (PetscInt d = 0; d < dim; d++) { + // 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; + } + } + } + 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 +709,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/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/compressibleFlowSolver.cpp b/src/finiteVolume/compressibleFlowSolver.cpp index f5607d6a4..9bbeacd4b 100644 --- a/src/finiteVolume/compressibleFlowSolver.cpp +++ b/src/finiteVolume/compressibleFlowSolver.cpp @@ -6,6 +6,7 @@ #include "finiteVolume/processes/evTransport.hpp" #include "finiteVolume/processes/navierStokesTransport.hpp" #include "finiteVolume/processes/speciesTransport.hpp" +#include "utilities/constants.hpp" #include "utilities/vectorUtilities.hpp" /** @@ -20,7 +21,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 +42,10 @@ ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string std::make_shared(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport), }, additionalProcesses), - std::move(boundaryConditions)) {} + 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; +} 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, @@ -60,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(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/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..6c21ed1be 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..88352e50a 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -11,6 +11,7 @@ #include "solver/cellSolver.hpp" #include "solver/solver.hpp" #include "solver/timeStepper.hpp" +#include "utilities/constants.hpp" #include "utilities/vectorUtilities.hpp" namespace ablate::finiteVolume { @@ -79,6 +80,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 = ablate::utilities::Constants::large; + public: FiniteVolumeSolver(std::string solverId, std::shared_ptr, std::shared_ptr options, std::vector> flowProcesses, std::vector> boundaryConditions); diff --git a/src/solver/timeStepper.cpp b/src/solver/timeStepper.cpp index 803d89089..2ce8ebb7a 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) {