diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 0e0c2e624bd0..3c1ded682430 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -27,6 +27,7 @@ #pragma once +#include #include #include @@ -94,6 +95,27 @@ class CFluidModel { public: virtual ~CFluidModel() {} + // Fully incoming flow needs inlet data; subsonic reverse flow uses extrapolated entropy and tangential velocity. + enum class PressureOutletStatus { Valid, InvalidState, Backflow }; + + /*! + * \brief Reconstruct a characteristic pressure outlet along this fluid's isentrope. + * Overwrites the thermodynamic scratch state. See the implementation for wave relations and references. + * \param[in] nDim - Spatial dimension (2 or 3). + * \param[in] interior - Conservative state [rho, rho*u_1, ..., rho*u_nDim, rho*E]. + * \param[in] normal - Outward unit normal. + * \param[in] gridNormal - Grid velocity dotted with normal. + * \param[in] targetPressure - Positive static pressure, in the interior state's units. + * \param[out] boundary - Conservative outlet state; use only when the result is Valid. + * \param[in] turbulentK - Frozen SST specific turbulent kinetic energy. + * \param[in] idealGamma - Constant ideal-gas gamma; zero selects general EOS integration. + * \return Valid, InvalidState (invalid input or EOS/integration failure), or Backflow. + */ + PressureOutletStatus ComputePressureOutletState(unsigned short nDim, const su2double* interior, + const su2double* normal, su2double gridNormal, + su2double targetPressure, su2double* boundary, + su2double turbulentK = 0, su2double idealGamma = 0); + /*! * \brief Get fluid pressure. */ @@ -403,4 +425,29 @@ class CFluidModel { * \return Newton solver iteration count at termination. */ virtual unsigned long GetnIter_Newton() { return 0; } + + private: + using PressureOutletState = std::array; // Density, thermal energy, relative normal velocity. + + /*! \brief Reconstruct the ideal-gas isentropic state, limited by sonic outflow. */ + PressureOutletStatus ComputeIdealGasPressureOutletState(su2double targetPressure, su2double idealGamma, + PressureOutletState& state); + + /*! \brief Integrate the general EOS isentrope up to the requested pressure or the sonic point. */ + PressureOutletStatus IntegratePressureOutletIsentrope(su2double targetPressure, PressureOutletState& state); + + /*! \brief Locate the sonic crossing and retain its sensitivity to the interior state. */ + bool FindPressureOutletSonicState(su2double& logPressure, PressureOutletState& state, su2double step, + su2double referenceSound2); + + /*! \brief Take one RK4 step, leaving the EOS at the resulting state. */ + bool StepPressureOutletIsentrope(su2double logPressure, const PressureOutletState& state, su2double step, + PressureOutletState& result, su2double referenceSound2); + + /*! \brief Evaluate isentropic derivatives with respect to log pressure at an admissible EOS state. */ + bool EvaluatePressureOutletIsentrope(su2double logPressure, const PressureOutletState& state, + PressureOutletState& derivative, su2double referenceSound2); + + /*! \brief Check the current EOS pressure, temperature and sound speed before using them. */ + bool IsPressureOutletStateValid() const; }; diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp index a5139f3e557b..ee70493abc0d 100644 --- a/SU2_CFD/src/fluid/CFluidModel.cpp +++ b/SU2_CFD/src/fluid/CFluidModel.cpp @@ -29,6 +29,8 @@ #include +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + #include "../../include/fluid/CConstantConductivity.hpp" #include "../../include/fluid/CConstantConductivityRANS.hpp" #include "../../include/fluid/CConstantDiffusivity.hpp" @@ -153,3 +155,243 @@ unique_ptr CFluidModel::MakeMassDiffusivityModel(const CConfi void CFluidModel::SetMassDiffusivityModel(const CConfig* config) { MassDiffusivity = MakeMassDiffusivityModel(config, 0); } + +/*! + * \brief Build a pressure-outlet state while preserving interior entropy and tangential velocity. + * + * The acoustic wave speeds are q-a and q+a, where q is outward normal velocity relative to the + * grid and a is sound speed. For 0 < q < a, only q-a enters the domain, so pressure supplies the + * one missing condition. The outgoing acoustic relation is dq = -dp/(rho*a) at constant entropy. + * Expansion stops at q=a: lower external pressure can no longer send information upstream. + * + * Characteristic boundary theory: K. W. Thompson, "Time dependent boundary conditions for + * hyperbolic systems", J. Comput. Phys. 68 (1987), 1-24, https://doi.org/10.1016/0021-9991(87)90041-6; + * T. J. Poinsot and S. K. Lele, "Boundary conditions for direct simulations of compressible viscous + * flows", J. Comput. Phys. 101 (1992), 104-129, https://doi.org/10.1016/0021-9991(92)90046-2. + * These references cover characteristic wave counting. The finite-change integration along the + * equation of state (EOS) isentrope and the sonic-root derivative below are implementation choices. + * + * Evaluations overwrite this fluid model's thermodynamic state. + * \param[in] nDim - Spatial dimension (2 or 3). + * \param[in] interior - Conservative state [rho, rho*u_1, ..., rho*u_nDim, rho*E]. + * \param[in] normal - Outward unit normal; tangential velocity is unchanged. + * \param[in] gridNormal - Grid velocity dotted with normal, in the same units as interior velocity. + * \param[in] targetPressure - Positive static thermodynamic pressure in the interior state's units. + * \param[out] boundary - Conservative outlet state; use only when the returned status is Valid. + * \param[in] turbulentK - Specific SST turbulent kinetic energy, held fixed during reconstruction. + * \param[in] idealGamma - Constant ideal-gas gamma; pass zero to integrate the general EOS. + * \return Valid, InvalidState (invalid input, EOS state or failed integration), or Backflow. + */ +CFluidModel::PressureOutletStatus CFluidModel::ComputePressureOutletState(unsigned short nDim, + const su2double* interior, + const su2double* normal, su2double gridNormal, + su2double targetPressure, su2double* boundary, + su2double turbulentK, su2double idealGamma) { + // Reject non-finite inputs and nonpositive density or pressure before using the EOS. + for (auto v = 0u; v < nDim + 2; ++v) + if (!std::isfinite(SU2_TYPE::GetValue(interior[v]))) return PressureOutletStatus::InvalidState; + const auto rho = interior[0]; + if (!(rho > 0) || !(targetPressure > 0) || !std::isfinite(SU2_TYPE::GetValue(targetPressure))) + return PressureOutletStatus::InvalidState; + // Recover velocity and thermal energy from the interior conservative state. + su2double velocity[3] = {}; + for (auto d = 0u; d < nDim; ++d) velocity[d] = interior[d + 1] / rho; + // SST total energy contains frozen turbulent kinetic energy; the EOS needs thermal energy only. + const su2double energy = interior[nDim + 1] / rho - .5 * GeometryToolbox::SquaredNorm(nDim, velocity) - turbulentK; + SetTDState_rhoe(rho, energy); + const auto sound2 = GetSoundSpeed2(); + if (!IsPressureOutletStateValid()) return PressureOutletStatus::InvalidState; + const su2double sound = sqrt(sound2); + // Wave direction depends on normal speed relative to the boundary, not total flow speed. + const su2double normalVelocity = GeometryToolbox::DotProduct(nDim, velocity, normal) - gridNormal; + // Subsonic reverse flow uses interior entropy and tangential velocity as inflow assumptions. + // Fully incoming flow needs inlet data that a pressure-only boundary cannot provide. + if (normalVelocity <= -sound) return PressureOutletStatus::Backflow; + if (normalVelocity >= sound) { + // No wave enters the domain, so the prescribed pressure cannot change this state. + for (auto v = 0u; v < nDim + 2; ++v) boundary[v] = interior[v]; + return PressureOutletStatus::Valid; + } + + PressureOutletState state = {rho, energy, normalVelocity}; + const auto status = idealGamma > 1 ? ComputeIdealGasPressureOutletState(targetPressure, idealGamma, state) + : IntegratePressureOutletIsentrope(targetPressure, state); + if (status != PressureOutletStatus::Valid) return status; + // Convert the outlet state back to lab-frame momentum, keeping tangential velocity unchanged. + boundary[0] = state[0]; + for (auto d = 0u; d < nDim; ++d) { + velocity[d] += (state[2] - normalVelocity) * normal[d]; + boundary[d + 1] = state[0] * velocity[d]; + } + // Restore kinetic and SST energy to form the total conservative energy. + boundary[nDim + 1] = state[0] * (GetStaticEnergy() + .5 * GeometryToolbox::SquaredNorm(nDim, velocity) + turbulentK); + return PressureOutletStatus::Valid; +} + +// Ideal gases have a closed-form solution, so no numerical integration is needed. +CFluidModel::PressureOutletStatus CFluidModel::ComputeIdealGasPressureOutletState(su2double targetPressure, + su2double idealGamma, + PressureOutletState& state) { + const auto rho = state[0], normalVelocity = state[2]; + const auto pressure = GetPressure(); + const su2double sound = sqrt(GetSoundSpeed2()); + const su2double gm1 = idealGamma - 1; + // Compute the sound speed at the sonic limit and at the requested pressure. + const su2double sonicSound = (gm1 * normalVelocity + 2 * sound) / (idealGamma + 1); + const su2double targetSound = sound * pow(targetPressure / pressure, gm1 / (2 * idealGamma)); + // Limit expansion to sonic speed, then recover density and normal velocity. + const su2double boundarySound = fmax(targetSound, sonicSound); + const su2double density = rho * pow(boundarySound / sound, 2 / gm1); + const su2double q = normalVelocity + 2 * (sound - boundarySound) / gm1; + if (q <= -boundarySound) return PressureOutletStatus::Backflow; + // Evaluate the EOS at the actual outlet pressure, which may be limited by the sonic state. + const su2double p = pressure * pow(density / rho, idealGamma); + SetTDState_Prho(p, density); + if (!IsPressureOutletStateValid()) return PressureOutletStatus::InvalidState; + state = {density, GetStaticEnergy(), q}; + return PressureOutletStatus::Valid; +} + +// Use log pressure so each integration step represents a relative pressure change. +CFluidModel::PressureOutletStatus CFluidModel::IntegratePressureOutletIsentrope(su2double targetPressure, + PressureOutletState& state) { + const auto rho = state[0], energy = state[1]; + const auto pressure = GetPressure(), sound2 = GetSoundSpeed2(); + const su2double sound = sqrt(sound2); + // Start at the interior pressure and choose the direction toward the requested pressure. + su2double x = log(pressure), end = log(targetPressure); + const auto direction = end >= x ? 1.0 : -1.0; + su2double h = direction * .05; + bool finished = false; + for (auto iteration = 0u; iteration < 10000; ++iteration) { + if (fabs(end - x) < 1e-13) { + // Finish the tiny remaining pressure change while retaining its sensitivity to the inputs. + PressureOutletState derivative; + if (!EvaluatePressureOutletIsentrope(x, state, derivative, sound2)) return PressureOutletStatus::InvalidState; + for (auto k = 0u; k < 3; ++k) state[k] += (end - x) * derivative[k]; + x = end; + finished = true; + break; + } + // Shorten the last step so it cannot pass the requested pressure. + if (fabs(h) > fabs(end - x)) h = end - x; + // Compare one full step with two half steps and retain the finer solution when accepted. + PressureOutletState coarse, half, fine; + bool valid = StepPressureOutletIsentrope(x, state, h, coarse, sound2) && + StepPressureOutletIsentrope(x, state, .5 * h, half, sound2) && + StepPressureOutletIsentrope(x + .5 * h, half, .5 * h, fine, sound2); + su2double error = 1; + // Scale each error by a typical density, energy or speed so their sizes can be compared. + if (valid) + error = fmax(fabs(fine[0] - coarse[0]) / (rho + fabs(fine[0])), + fmax(fabs(fine[1] - coarse[1]) / (fabs(energy) + sound2 + fabs(fine[1])), + fabs(fine[2] - coarse[2]) / (sound + fabs(fine[2])))); + // Retry with a smaller step if the EOS rejects a stage or the integration error is too large. + if (!valid || error > 1e-10) { + h *= .5; + if (fabs(h) < 1e-12) return PressureOutletStatus::InvalidState; + continue; + } + const su2double nextSound = sqrt(GetSoundSpeed2()); + // If this step reaches sonic speed, stop expansion at the crossing instead of at the target pressure. + if (fine[2] >= nextSound) { + if (!FindPressureOutletSonicState(x, state, h, sound2)) return PressureOutletStatus::InvalidState; + finished = true; + break; + } + // Accept the two-half-step result; grow the next step if the error was comfortably small. + x += h; + state = fine; + if (error < 1e-12) h *= 2; + } + // Reject an unfinished integration and check the final EOS state before building the output. + if (!finished) return PressureOutletStatus::InvalidState; + SetTDState_rhoe(state[0], state[1]); + if (!IsPressureOutletStateValid()) return PressureOutletStatus::InvalidState; + // Check the accumulated pressure error at the target or sonic stop, not at RK predictors. + if (fabs(GetPressure() / exp(x) - 1) > 1e-6) return PressureOutletStatus::InvalidState; + if (state[2] <= -GetSoundSpeed()) return PressureOutletStatus::Backflow; + return PressureOutletStatus::Valid; +} + +bool CFluidModel::FindPressureOutletSonicState(su2double& x, PressureOutletState& state, su2double h, + su2double referenceSound2) { + PressureOutletState half, fine; + su2double lower = 0, upper = 1; + // Bisect the fraction of this step, reintegrating from the last accepted state each time. + for (auto k = 0u; k < 40; ++k) { + const su2double fraction = .5 * (lower + upper); + if (!StepPressureOutletIsentrope(x, state, fraction * h * .5, half, referenceSound2) || + !StepPressureOutletIsentrope(x + fraction * h * .5, half, fraction * h * .5, fine, referenceSound2)) + return false; + if (fine[2] >= sqrt(GetSoundSpeed2())) + upper = fraction; + else + lower = fraction; + } + // Keep the state and pressure at the sonic crossing. + x += .5 * (lower + upper) * h; + state = fine; + // Bisection finds the sonic point; a Newton correction recovers how it moves with the inputs. + PressureOutletState derivative; + if (!EvaluatePressureOutletIsentrope(x, state, derivative, referenceSound2)) return false; + const su2double sonicResidual = state[2] - sqrt(GetSoundSpeed2()); + // Sample sound speed on either side to estimate how q-a changes along the pressure path. + constexpr double delta = 1e-5; + PressureOutletState shifted = state, unused; + for (auto k = 0u; k < 3; ++k) shifted[k] += delta * derivative[k]; + if (!EvaluatePressureOutletIsentrope(x + delta, shifted, unused, referenceSound2)) return false; + const double soundPlus = SU2_TYPE::GetValue(GetSoundSpeed()); + shifted = state; + for (auto k = 0u; k < 3; ++k) shifted[k] -= delta * derivative[k]; + if (!EvaluatePressureOutletIsentrope(x - delta, shifted, unused, referenceSound2)) return false; + const double soundMinus = SU2_TYPE::GetValue(GetSoundSpeed()); + // Only the residual needs automatic differentiation here; the slope's derivative has negligible weight. + const double slope = SU2_TYPE::GetValue(derivative[2]) - (soundPlus - soundMinus) / (2 * delta); + if (!std::isfinite(slope) || fabs(slope) < 1e-12) return false; + // Apply the correction to both the state and its stopping pressure. + const su2double correction = -sonicResidual / slope; + for (auto k = 0u; k < 3; ++k) state[k] += correction * derivative[k]; + x += correction; + return true; +} + +// Classical fourth-order Runge-Kutta; each stage must remain inside the EOS validity range. +bool CFluidModel::StepPressureOutletIsentrope(su2double x, const PressureOutletState& y, su2double h, + PressureOutletState& out, su2double referenceSound2) { + PressureOutletState k1, k2, k3, k4, trial; + // Sample the slope at the start and at two estimates of the midpoint. + if (!EvaluatePressureOutletIsentrope(x, y, k1, referenceSound2)) return false; + for (auto k = 0u; k < 3; ++k) trial[k] = y[k] + .5 * h * k1[k]; + if (!EvaluatePressureOutletIsentrope(x + .5 * h, trial, k2, referenceSound2)) return false; + for (auto k = 0u; k < 3; ++k) trial[k] = y[k] + .5 * h * k2[k]; + if (!EvaluatePressureOutletIsentrope(x + .5 * h, trial, k3, referenceSound2)) return false; + // Sample the endpoint, then combine all four slopes to advance the state. + for (auto k = 0u; k < 3; ++k) trial[k] = y[k] + h * k3[k]; + if (!EvaluatePressureOutletIsentrope(x + h, trial, k4, referenceSound2)) return false; + for (auto k = 0u; k < 3; ++k) out[k] = y[k] + h * (k1[k] + 2 * k2[k] + 2 * k3[k] + k4[k]) / 6; + // Validate the combined state and leave the EOS at that endpoint. + return EvaluatePressureOutletIsentrope(x + h, out, trial, referenceSound2); +} + +// Return how density, thermal energy and relative normal velocity change with log pressure. +bool CFluidModel::EvaluatePressureOutletIsentrope(su2double logPressure, const PressureOutletState& y, + PressureOutletState& derivative, su2double referenceSound2) { + if (!(y[0] > 0) || !std::isfinite(SU2_TYPE::GetValue(y[0] + y[1] + y[2]))) return false; + const su2double p = exp(logPressure); + SetTDState_rhoe(y[0], y[1]); + // Intermediate stages can have pressure error, but the EOS must not silently change their density or energy. + if (!IsPressureOutletStateValid() || fabs(GetDensity() / y[0] - 1) > 1e-6 || + fabs(GetStaticEnergy() - y[1]) > 1e-6 * (referenceSound2 + fabs(y[1]))) + return false; + const auto a2 = GetSoundSpeed2(); + // Constant-entropy derivatives: d(rho)/dp=1/a^2, de/dp=p/(rho^2*a^2), dq/dp=-1/(rho*a). + derivative = {p / a2, p * p / (y[0] * y[0] * a2), -p / (y[0] * sqrt(a2))}; + return true; +} + +bool CFluidModel::IsPressureOutletStateValid() const { + const auto p = SU2_TYPE::GetValue(GetPressure()), t = SU2_TYPE::GetValue(GetTemperature()); + const auto a2 = SU2_TYPE::GetValue(GetSoundSpeed2()); + return p > 0 && t > 0 && a2 > 0 && std::isfinite(p + t + a2); +} diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..0c6427207f0f 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -5308,7 +5308,9 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Density_i = nodes->GetDensity(iPoint); const auto Energy_i = nodes->GetEnergy(iPoint); - const su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + const bool pressureOutflow = config->GetKind_Data_Riemann(Marker_Tag) == STATIC_PRESSURE; + const auto tke_i = pressureOutflow && tkeNeeded ? solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0) : su2double(0.0); + const su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i - tke_i; GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5436,18 +5438,15 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, case STATIC_PRESSURE: - /*--- Retrieve the static pressure for this boundary. ---*/ + /*--- Only pressure is prescribed; CFluidModel::ComputePressureOutletState derives density and energy on the isentrope. ---*/ Pressure_e = config->GetRiemann_Var1(Marker_Tag); Pressure_e /= config->GetPressure_Ref(); Density_e = Density_i; /* --- Compute the boundary state u_e --- */ - GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e); for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Velocity_i[iDim]; - - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); - Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + Energy_e = Energy_i; // The nonlinear outflow construction below supplies the pressure state. break; default: @@ -5455,11 +5454,16 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, break; } - /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + // A linear characteristic jump can leave the EOS validity range for a large pressure change; + // static-pressure outlets instead follow the nonlinear isentrope in CFluidModel::ComputePressureOutletState. + if (!pressureOutflow || implicit) { + /*--- Compute P (matrix of right eigenvectors) ---*/ + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + + /*--- Compute inverse P (matrix of left eigenvectors)---*/ + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); - /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + } /*--- eigenvalues contribution due to grid motion ---*/ su2double ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal); @@ -5484,32 +5488,69 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, u_e[nVar-1] = Energy_e*Density_e; u_i[nVar-1] = Energy_i*Density_i; - /*--- Compute the characteristic jumps ---*/ - for (auto iVar = 0u; iVar < nVar; iVar++){ - dw[iVar] = 0; - for (auto jVar = 0u; jVar < nVar; jVar++) - dw[iVar] += invP_Tensor[iVar][jVar] * (u_e[jVar] - u_i[jVar]); + if (!pressureOutflow) { + /*--- Compute the characteristic jumps ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++){ + dw[iVar] = 0; + for (auto jVar = 0u; jVar < nVar; jVar++) + dw[iVar] += invP_Tensor[iVar][jVar] * (u_e[jVar] - u_i[jVar]); + } + + /*--- Compute the boundary state u_b using characteristics ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) { + u_b[iVar] = u_i[iVar]; + for (auto jVar = 0u; jVar < nVar; jVar++) { + if (Lambda_i[jVar] < 0) + u_b[iVar] += P_Tensor[iVar][jVar]*dw[jVar]; + } + } + } - /*--- Compute the boundary state u_b using characteristics ---*/ - for (auto iVar = 0u; iVar < nVar; iVar++) { - u_b[iVar] = u_i[iVar]; - for (auto jVar = 0u; jVar < nVar; jVar++) { - if (Lambda_i[jVar] < 0) - u_b[iVar] += P_Tensor[iVar][jVar]*dw[jVar]; + const su2double gridNormal = dynamic_grid ? GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal) : su2double(0.0); + if (pressureOutflow) { + const auto status = GetFluidModel()->ComputePressureOutletState(nDim, u_i, UnitNormal, gridNormal, Pressure_e, u_b, tke_i, ideal_gas ? Gamma : su2double(0.0)); + if (status != CFluidModel::PressureOutletStatus::Valid) { + // MPI error reporting uses shared state; only one OpenMP thread may enter it. + SU2_OMP_CRITICAL + { + if (status == CFluidModel::PressureOutletStatus::Backflow) + SU2_MPI::Error("Static-pressure outflow requires an inflow state for backflow with all characteristics incoming on marker " + Marker_Tag, CURRENT_FUNCTION); + else + SU2_MPI::Error("Unable to construct an admissible static-pressure outflow state on marker " + Marker_Tag + + ". Check the prescribed pressure and the EOS validity range.", CURRENT_FUNCTION); + } END_SU2_OMP_CRITICAL } } /*--- Compute the thermodynamic state in u_b ---*/ su2double Velocity_b[MAXNDIM] = {0}; const auto Density_b = u_b[0]; + // Reject invalid density before dividing momentum or energy by it. + if (!(Density_b > 0.0) || !std::isfinite(SU2_TYPE::GetValue(Density_b))) { + SU2_MPI::Error("Nonphysical characteristic boundary density on marker " + Marker_Tag + + ". Check the prescribed state and initialization; pressure continuation may be needed for a steady solve.", + CURRENT_FUNCTION); + } + for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_b[iDim] = u_b[iDim+1]/Density_b; const auto Velocity2_b = GeometryToolbox::SquaredNorm(nDim, Velocity_b); const su2double Energy_b = u_b[nVar-1]/Density_b; - const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b - tke_i; GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); + const auto boundaryPressure = SU2_TYPE::GetValue(GetFluidModel()->GetPressure()); + const auto boundaryTemperature = SU2_TYPE::GetValue(GetFluidModel()->GetTemperature()); + const auto boundarySoundSpeed2 = SU2_TYPE::GetValue(GetFluidModel()->GetSoundSpeed2()); + // Diagnose invalid EOS states here, before a flux or Jacobian turns them into a NaN residual. + if (!(boundaryPressure > 0.0 && boundaryTemperature > 0.0 && boundarySoundSpeed2 > 0.0) || + !std::isfinite(boundaryPressure) || !std::isfinite(boundaryTemperature) || !std::isfinite(boundarySoundSpeed2)) { + SU2_MPI::Error("Nonphysical characteristic boundary thermodynamic state on marker " + Marker_Tag + + ". Check the prescribed state and initialization; pressure continuation may be needed for a steady solve.", + CURRENT_FUNCTION); + } + /*--- Store number of Newton iterations at BC ---*/ if(config->GetKind_FluidModel() == DATADRIVEN_FLUID) @@ -5519,7 +5560,17 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Temperature_b = GetFluidModel()->GetTemperature(); const su2double Enthalpy_b = Energy_b + Pressure_b/Density_b; const su2double Kappa_b = GetFluidModel()->GetdPde_rho() / Density_b; - const su2double Chi_b = GetFluidModel()->GetdPdrho_e() - Kappa_b * StaticEnergy_b; + // With SST, E=e+|u|^2/2+k; differentiating pressure at fixed k requires e+k in chi. + const su2double Chi_b = GetFluidModel()->GetdPdrho_e() - Kappa_b * (StaticEnergy_b + tke_i); + + /*--- Scalar boundary solvers also need the reconstructed primitive state in Euler runs. ---*/ + auto* V_boundary = GetCharacPrimVar(val_marker, iVertex); + V_boundary[prim_idx.Temperature()] = Temperature_b; + for (auto d = 0u; d < nDim; ++d) V_boundary[prim_idx.Velocity() + d] = Velocity_b[d]; + V_boundary[prim_idx.Pressure()] = Pressure_b; + V_boundary[prim_idx.Density()] = Density_b; + V_boundary[prim_idx.Enthalpy()] = Enthalpy_b; + V_boundary[prim_idx.SoundSpeed()] = GetFluidModel()->GetSoundSpeed(); /*--- Compute the residuals ---*/ su2double Residual[MAXNVAR] = {0}; @@ -5551,7 +5602,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, DubDu[iVar][iVar]=1; } - /*--- Compute DubDu -= RNL---*/ + // Keep outgoing characteristic increments; freeze the wave basis in this approximate matrix. for (auto iVar=0u; iVarGetLaminarViscosity(); @@ -5706,6 +5748,9 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain bool viscous = config->GetViscous(); bool gravity = (config->GetGravityForce()); bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + // TurboRiemann reads MARKER_RIEMANN; a matching MARKER_GILES entry is not required or consulted. + const auto dataKind = config->GetKind_Data_Riemann(Marker_Tag); + const bool pressureOutflow = dataKind == STATIC_PRESSURE || dataKind == RADIAL_EQUILIBRIUM; const bool ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS; @@ -5780,7 +5825,10 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Density_i = nodes->GetDensity(iPoint); Energy_i = nodes->GetEnergy(iPoint); - StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + // Keep SST kinetic energy in the conserved state, but exclude it from the thermal EOS input. + const su2double tke_i = pressureOutflow && tkeNeeded + ? solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0) : su2double(0.0); + StaticEnergy_i = Energy_i - 0.5*Velocity2_i - tke_i; GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5790,29 +5838,32 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain SoundSpeed_i = GetFluidModel()->GetSoundSpeed(); Kappa_i = GetFluidModel()->GetdPde_rho() / Density_i; - Chi_i = GetFluidModel()->GetdPdrho_e() - Kappa_i * StaticEnergy_i; + Chi_i = GetFluidModel()->GetdPdrho_e() - Kappa_i * (StaticEnergy_i + tke_i); ProjVelocity_i = 0.0; for (iDim = 0; iDim < nDim; iDim++) ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; - su2double donorAverages[5] = {0.0}; - switch (config->GetKind_Data_Giles(Marker_Tag)){ - case MIXING_IN: case MIXING_IN_1D: case MIXING_OUT: case MIXING_OUT_1D: - for (auto mixVar = 0u; mixVar < 5; mixVar++) donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar); - break; - default: - break; + // Static and radial-equilibrium pressure outlets need no donor mixing state. + if (!pressureOutflow) { + su2double donorAverages[5] = {0.0}; + switch (dataKind){ + case MIXING_IN: case MIXING_IN_1D: case MIXING_OUT: case MIXING_OUT_1D: + for (auto mixVar = 0u; mixVar < 5; mixVar++) donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar); + break; + default: + break; + } + ExtAverageDensity = donorAverages[0]; + ExtAveragePressure = donorAverages[1]; + ExtAverageTurboVelocity[0] = donorAverages[2]; + ExtAverageTurboVelocity[1] = donorAverages[3]; + if (nDim == 3) ExtAverageTurboVelocity[2] = donorAverages[4]; } - ExtAverageDensity = donorAverages[0]; - ExtAveragePressure = donorAverages[1]; - ExtAverageTurboVelocity[0] = donorAverages[2]; - ExtAverageTurboVelocity[1] = donorAverages[3]; - if (nDim == 3) ExtAverageTurboVelocity[2] = donorAverages[4]; /*--- Build the external state u_e from boundary data and internal node ---*/ - switch(config->GetKind_Data_Riemann(Marker_Tag)) + switch(dataKind) { //TODO(turbo), generalize for 3D case //TODO(turbo), generalize for Inlet and Outlet in for backflow treatment @@ -5899,38 +5950,15 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; break; - case STATIC_PRESSURE: + case STATIC_PRESSURE: case RADIAL_EQUILIBRIUM: - /*--- Retrieve the static pressure for this boundary. ---*/ - Pressure_e = config->GetRiemann_Var1(Marker_Tag); - Pressure_e /= config->GetPressure_Ref(); + // The radial profile is already nondimensional; a configured uniform pressure is dimensional. + Pressure_e = dataKind == RADIAL_EQUILIBRIUM + ? RadialEquilibriumPressure[val_marker][iSpan] + : config->GetRiemann_Var1(Marker_Tag) / config->GetPressure_Ref(); Density_e = Density_i; - - /* --- Compute the boundary state u_e --- */ - GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e); - Velocity2_e = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_e[iDim] = Velocity_i[iDim]; - Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; - } - Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; - break; - - - case RADIAL_EQUILIBRIUM: - - /*--- Retrieve the static pressure for this boundary. ---*/ - Pressure_e = RadialEquilibriumPressure[val_marker][iSpan]; - Density_e = Density_i; - - /* --- Compute the boundary state u_e --- */ - GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e); - Velocity2_e = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_e[iDim] = Velocity_i[iDim]; - Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; - } - Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + Energy_e = Energy_i; // The shared pressure-outflow reconstruction supplies the boundary state. + for (iDim = 0; iDim < nDim; iDim++) Velocity_e[iDim] = Velocity_i[iDim]; break; default: @@ -5938,11 +5966,13 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain break; } - /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + if (!pressureOutflow || implicit) { + /*--- Compute P (matrix of right eigenvectors) ---*/ + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); - /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + /*--- Compute inverse P (matrix of left eigenvectors)---*/ + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + } /*--- eigenvalues contribution due to grid motion ---*/ if (dynamic_grid){ @@ -5972,31 +6002,39 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain u_i[iDim+1] = Velocity_i[iDim]*Density_i; u_i[nVar-1] = Energy_i*Density_i; - /*--- Compute the characteristic jumps ---*/ - for (iVar = 0; iVar < nVar; iVar++) - { - dw[iVar] = 0; - for (jVar = 0; jVar < nVar; jVar++) - dw[iVar] += invP_Tensor[iVar][jVar] * (u_e[jVar] - u_i[jVar]); - - } - - /*--- Compute the boundary state u_b using characteristics ---*/ - for (iVar = 0; iVar < nVar; iVar++) - { - u_b[iVar] = u_i[iVar]; - - for (jVar = 0; jVar < nVar; jVar++) - { - if (Lambda_i[jVar] < 0) + if (pressureOutflow) { + // Apply each span's pressure to the local state; using a span average here would discard local waves. + const su2double gridNormal = dynamic_grid ? GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal) : su2double(0.0); + const auto status = GetFluidModel()->ComputePressureOutletState(nDim, u_i, UnitNormal, gridNormal, Pressure_e, + u_b, tke_i, ideal_gas ? Gamma : su2double(0.0)); + if (status != CFluidModel::PressureOutletStatus::Valid) { + // MPI error reporting uses shared state; only one OpenMP thread may enter it. + SU2_OMP_CRITICAL { - u_b[iVar] += P_Tensor[iVar][jVar]*dw[jVar]; + if (status == CFluidModel::PressureOutletStatus::Backflow) + SU2_MPI::Error("Static-pressure outflow requires an inflow state for backflow with all characteristics incoming on marker " + Marker_Tag, CURRENT_FUNCTION); + else + SU2_MPI::Error("Unable to construct an admissible static-pressure outflow state on marker " + Marker_Tag + + ". Check the prescribed pressure and the EOS validity range.", CURRENT_FUNCTION); + } END_SU2_OMP_CRITICAL + } + } else { + /*--- Compute the characteristic jumps ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + dw[iVar] = 0; + for (jVar = 0; jVar < nVar; jVar++) + dw[iVar] += invP_Tensor[iVar][jVar] * (u_e[jVar] - u_i[jVar]); + } + /*--- Compute the boundary state u_b using characteristics ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + u_b[iVar] = u_i[iVar]; + for (jVar = 0; jVar < nVar; jVar++) { + if (Lambda_i[jVar] < 0) u_b[iVar] += P_Tensor[iVar][jVar]*dw[jVar]; } } } - /*--- Compute the thermodynamic state in u_b ---*/ Density_b = u_b[0]; Velocity2_b = 0; @@ -6006,13 +6044,21 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_b += Velocity_b[iDim]*Velocity_b[iDim]; } Energy_b = u_b[nVar-1]/Density_b; - StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + StaticEnergy_b = Energy_b - 0.5*Velocity2_b - tke_i; GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); Pressure_b = GetFluidModel()->GetPressure(); Temperature_b = GetFluidModel()->GetTemperature(); Enthalpy_b = Energy_b + Pressure_b/Density_b; Kappa_b = GetFluidModel()->GetdPde_rho() / Density_b; - Chi_b = GetFluidModel()->GetdPdrho_e() - Kappa_b * StaticEnergy_b; + Chi_b = GetFluidModel()->GetdPdrho_e() - Kappa_b * (StaticEnergy_b + tke_i); + + // Turbulence and other scalar boundary solvers read this state even when viscous fluxes are disabled. + V_boundary[prim_idx.Temperature()] = Temperature_b; + for (auto d = 0u; d < nDim; ++d) V_boundary[prim_idx.Velocity() + d] = Velocity_b[d]; + V_boundary[prim_idx.Pressure()] = Pressure_b; + V_boundary[prim_idx.Density()] = Density_b; + V_boundary[prim_idx.Enthalpy()] = Enthalpy_b; + V_boundary[prim_idx.SoundSpeed()] = GetFluidModel()->GetSoundSpeed(); /*--- Compute the residuals ---*/ conv_numerics->GetInviscidProjFlux(&Density_b, Velocity_b, &Pressure_b, &Enthalpy_b, Normal, Residual); @@ -6048,7 +6094,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain DubDu[iVar][iVar]= 1; } - /*--- Compute DubDu -= RNL---*/ + // Keep outgoing characteristic increments; freeze the wave basis in this approximate matrix. for (iVar=0; iVarGetLaminarViscosity(); @@ -7448,175 +7486,103 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_ZONE_SCOPED - unsigned short iVar, iDim; - unsigned long iVertex, iPoint; - su2double Pressure, P_Exit, Velocity[3], - Velocity2, Entropy, Density, Energy, Riemann, Vn, SoundSpeed, Mach_Exit, Vn_Exit, - Area, UnitNormal[3]; - su2double *V_outlet, *V_domain; - - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - su2double Gas_Constant = config->GetGas_ConstantND(); - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - bool gravity = (config->GetGravityForce()); - bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); - - auto *Normal = new su2double[nDim]; - - /*--- Loop over all the vertices on this boundary marker ---*/ + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + const bool tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST; + const string Marker_Tag = config->GetMarker_All_TagBound(val_marker); SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - /*--- Allocate the value at the outlet ---*/ - V_outlet = GetCharacPrimVar(val_marker, iVertex); - - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); - - Area = GeometryToolbox::Norm(nDim, Normal); - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; - - /*--- Current solution at this boundary node ---*/ - V_domain = nodes->GetPrimitive(iPoint); - - /*--- Build the fictitious inlet state based on characteristics ---*/ - - /*--- Retrieve the specified back pressure for this outlet. ---*/ - if (gravity) P_Exit = config->GetOutlet_Pressure(Marker_Tag) - geometry->nodes->GetCoord(iPoint, nDim-1)*STANDARD_GRAVITY; - else P_Exit = config->GetOutlet_Pressure(Marker_Tag); - - /*--- Non-dim. the inputs if necessary. ---*/ - P_Exit = P_Exit/config->GetPressure_Ref(); - - /*--- Check whether the flow is supersonic at the exit. The type - of boundary update depends on this. ---*/ - Density = V_domain[nDim+2]; - Velocity2 = 0.0; Vn = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity[iDim] = V_domain[iDim+1]; - Velocity2 += Velocity[iDim]*Velocity[iDim]; - Vn += Velocity[iDim]*UnitNormal[iDim]; - } - Pressure = V_domain[nDim+1]; - SoundSpeed = sqrt(Gamma*Pressure/Density); - Mach_Exit = sqrt(Velocity2)/SoundSpeed; - - if (Mach_Exit >= 1.0) { - - /*--- Supersonic exit flow: there are no incoming characteristics, - so no boundary condition is necessary. Set outlet state to current - state so that upwinding handles the direction of propagation. ---*/ - for (iVar = 0; iVar < nPrimVar; iVar++) V_outlet[iVar] = V_domain[iVar]; + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; ++iVertex) { + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + if (!geometry->nodes->GetDomain(iPoint)) continue; - } else { + su2double Normal[MAXNDIM] = {}, UnitNormal[MAXNDIM] = {}; + // Reconstruction uses an outward unit normal; the numerical flux uses the outward area normal. + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (auto d = 0u; d < nDim; ++d) Normal[d] = -Normal[d]; + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + for (auto d = 0u; d < nDim; ++d) UnitNormal[d] = Normal[d] / Area; + conv_numerics->SetNormal(Normal); - /*--- Subsonic exit flow: there is one incoming characteristic, - therefore one variable can be specified (back pressure) and is used - to update the conservative variables. Compute the entropy and the - acoustic Riemann variable. These invariants, as well as the - tangential velocity components, are extrapolated. Adapted from an - original implementation in the Stanford University multi-block - (SUmb) solver in the routine bcSubsonicOutflow.f90 by Edwin van - der Weide, last modified 09-10-2007. ---*/ + su2double pressure = config->GetOutlet_Pressure(Marker_Tag); + if (config->GetGravityForce()) pressure -= geometry->nodes->GetCoord(iPoint, nDim - 1) * STANDARD_GRAVITY; + pressure /= config->GetPressure_Ref(); + // The helper classifies flow in the grid's frame while keeping conservative energy in the lab frame. + const su2double gridNormal = dynamic_grid ? GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal) : su2double(0.0); + const su2double tke = tkeNeeded ? solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0) : su2double(0.0); + const auto* interior = nodes->GetSolution(iPoint); + su2double boundary[MAXNVAR] = {}; + // Ordinary compressible outlets use the ideal-gas path; real-gas pressure outlets use BC_Riemann. + const auto status = GetFluidModel()->ComputePressureOutletState(nDim, interior, UnitNormal, gridNormal, pressure, + boundary, tke, Gamma); + if (status != CFluidModel::PressureOutletStatus::Valid) { + // MPI error reporting uses shared state; only one OpenMP thread may enter it. + SU2_OMP_CRITICAL + { + if (status == CFluidModel::PressureOutletStatus::Backflow) + SU2_MPI::Error("Static-pressure outflow requires an inflow state for backflow with all characteristics incoming on marker " + Marker_Tag, CURRENT_FUNCTION); + else + SU2_MPI::Error("Unable to construct an admissible static-pressure outflow state on marker " + Marker_Tag + + ". Check the prescribed pressure and the EOS validity range.", CURRENT_FUNCTION); + } END_SU2_OMP_CRITICAL + } - Entropy = Pressure*pow(1.0/Density, Gamma); - Riemann = Vn + 2.0*SoundSpeed/Gamma_Minus_One; + const su2double density = boundary[0], energy = boundary[nVar - 1] / density; + const auto* V_domain = nodes->GetPrimitive(iPoint); + auto* V_outlet = GetCharacPrimVar(val_marker, iVertex); + // Preserve primitive fields outside the reconstructed thermodynamic and velocity entries. + for (auto variable = 0u; variable < nPrimVar; ++variable) V_outlet[variable] = V_domain[variable]; + for (auto d = 0u; d < nDim; ++d) V_outlet[prim_idx.Velocity() + d] = boundary[d + 1] / density; + // Recover the primitive thermodynamic fields from the reconstructed conservative state. + GetFluidModel()->SetTDState_rhoe(density, energy - .5 * GeometryToolbox::SquaredNorm(nDim, V_outlet + prim_idx.Velocity()) - tke); + V_outlet[prim_idx.Temperature()] = GetFluidModel()->GetTemperature(); + V_outlet[prim_idx.Pressure()] = GetFluidModel()->GetPressure(); + V_outlet[prim_idx.Density()] = density; + V_outlet[prim_idx.Enthalpy()] = energy + V_outlet[prim_idx.Pressure()] / density; + V_outlet[prim_idx.SoundSpeed()] = GetFluidModel()->GetSoundSpeed(); + conv_numerics->SetPrimitive(V_domain, V_outlet); + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Compute the new fictious state at the outlet ---*/ - Density = pow(P_Exit/Entropy,1.0/Gamma); - Pressure = P_Exit; - SoundSpeed = sqrt(Gamma*P_Exit/Density); - Vn_Exit = Riemann - 2.0*SoundSpeed/Gamma_Minus_One; - Velocity2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity[iDim] = Velocity[iDim] + (Vn_Exit-Vn)*UnitNormal[iDim]; - Velocity2 += Velocity[iDim]*Velocity[iDim]; + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); + if (implicit) { + su2double Pdata[MAXNVAR][MAXNVAR] = {}, Idata[MAXNVAR][MAXNVAR] = {}; + su2double projection[MAXNVAR][MAXNVAR] = {}, matrix[MAXNVAR][MAXNVAR] = {}; + su2double *P[MAXNVAR], *Pinv[MAXNVAR], *rows[MAXNVAR]; + for (auto row = 0u; row < nVar; ++row) { + P[row] = Pdata[row]; + Pinv[row] = Idata[row]; + rows[row] = matrix[row]; + } + const auto* velocity = V_domain + prim_idx.Velocity(); + const su2double sound = V_domain[prim_idx.SoundSpeed()]; + // Freeze the wave basis in the implicit update; the residual keeps the nonlinear outlet state. + // The general eigenvectors include frozen SST energy through enthalpy and chi. + const su2double rho = interior[0], enthalpy = V_domain[prim_idx.Enthalpy()]; + const su2double kappa = Gamma_Minus_One, chi = -kappa * tke; + conv_numerics->GetPMatrix(&rho, velocity, &sound, &enthalpy, &chi, &kappa, UnitNormal, P); + conv_numerics->GetPMatrix_inv(Pinv, &rho, velocity, &sound, &chi, &kappa, UnitNormal); + const su2double q = GeometryToolbox::DotProduct(nDim, velocity, UnitNormal) - gridNormal; + // Exclude incoming waves from this approximate state derivative. + for (auto row = 0u; row < nVar; ++row) + for (auto col = 0u; col < nVar; ++col) { + projection[row][col] = row == col ? 1.0 : 0.0; + for (auto mode = 0u; mode < nVar; ++mode) { + su2double speed = q; + if (mode == nVar - 2) speed += sound; + if (mode == nVar - 1) speed -= sound; + if (speed < 0) projection[row][col] -= P[row][mode] * Pinv[mode][col]; + } } - Energy = P_Exit/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); - - /*--- Conservative variables, using the derived quantities ---*/ - V_outlet[0] = Pressure / ( Gas_Constant * Density); - for (iDim = 0; iDim < nDim; iDim++) - V_outlet[iDim+1] = Velocity[iDim]; - V_outlet[nDim+1] = Pressure; - V_outlet[nDim+2] = Density; - V_outlet[nDim+3] = Energy + Pressure/Density; - - } - - /*--- Set various quantities in the solver class ---*/ - conv_numerics->SetPrimitive(V_domain, V_outlet); - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - - /*--- Compute the residual using an upwind scheme ---*/ - - auto residual = conv_numerics->ComputeResidual(config); - - /*--- Add Residuals and Jacobians ---*/ - - LinSysRes.AddBlock(iPoint, residual); - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - -// /*--- Viscous contribution, commented out because serious convergence problems ---*/ -// -// if (viscous) { -// -// /*--- Set laminar and eddy viscosity at the infinity ---*/ -// -// V_outlet[nDim+5] = nodes->GetLaminarViscosity(iPoint); -// V_outlet[nDim+6] = nodes->GetEddyViscosity(iPoint); -// -// /*--- Set the normal vector and the coordinates ---*/ -// -// visc_numerics->SetNormal(Normal); -// su2double Coord_Reflected[MAXNDIM]; -// GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), -// geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// -// /*--- Primitive variables, and gradient ---*/ -// -// visc_numerics->SetPrimitive(V_domain, V_outlet); -// visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); -// -// /*--- Turbulent kinetic energy ---*/ -// -// if (config->GetKind_Turb_Model() == TURB_MODEL::SST) -// visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), -// solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); -// -// /*--- Compute and update residual ---*/ -// -// auto residual = visc_numerics->ComputeResidual(config); -// LinSysRes.SubtractBlock(iPoint, residual); -// -// /*--- Jacobian contribution for implicit integration ---*/ -// -// if (implicit) -// Jacobian.SubtractBlock2Diag(iPoint, residual.acobian_i); -// -// } - + for (auto row = 0u; row < nVar; ++row) + for (auto col = 0u; col < nVar; ++col) { + matrix[row][col] = residual.jacobian_i[row][col]; + for (auto k = 0u; k < nVar; ++k) matrix[row][col] += residual.jacobian_j[row][k] * projection[k][col]; + } + Jacobian.AddBlock2Diag(iPoint, rows); } } END_SU2_OMP_FOR - - /*--- Free locally allocated memory ---*/ - delete [] Normal; - } void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_container, diff --git a/TestCases/euler/wedge/inv_wedge_HLLC.cfg b/TestCases/euler/wedge/inv_wedge_HLLC.cfg index c89cdbd4b276..dd43193bb4d8 100644 --- a/TestCases/euler/wedge/inv_wedge_HLLC.cfg +++ b/TestCases/euler/wedge/inv_wedge_HLLC.cfg @@ -47,7 +47,8 @@ CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) ITER= 11000 -LINEAR_SOLVER= BCGSTAB +% FGMRES remains effective with short mixed-precision solves at this high CFL. +LINEAR_SOLVER= FGMRES LINEAR_SOLVER_ERROR= 1E-1 LINEAR_SOLVER_ITER= 3 diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 2475709cf540..f16c5e22295b 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -80,7 +80,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [ 0.036140, 6.976887, -0.219795, 0.055341] + wedge.test_vals = [-3.366758, 2.363084, -0.249525, 0.043953] test_list.append(wedge) # ONERA M6 Wing @@ -116,7 +116,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 100 - flatplate.test_vals = [-8.265203, -2.790380, 0.001084, 0.036280, 2.361500, -2.325300, 0.000000, 0.000000] + flatplate.test_vals = [-8.270397, -2.795629, 0.001084, 0.036280, 2.361500, -2.325300, 0.000000, 0.000000] test_list.append(flatplate) # Laminar cylinder (steady) @@ -149,7 +149,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.004080, -7.638116, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.004071, -7.637239, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -158,7 +158,7 @@ def main(): periodic2d.cfg_dir = "navierstokes/periodic2D" periodic2d.cfg_file = "config.cfg" periodic2d.test_iter = 1400 - periodic2d.test_vals = [-10.314545, -8.056213, -7.708669, -4.833054, -1.059806, -2953.400000] + periodic2d.test_vals = [-10.314271, -8.055939, -7.708395, -4.832780, -1.059805, -2953.400000] test_list.append(periodic2d) ########################## @@ -194,7 +194,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.443902] + turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.444207] test_list.append(turb_flatplate) # ONERA M6 Wing @@ -265,7 +265,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-2.663059, 2.911787, -2.522191, 1.990732, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.663033, 2.911808, -2.517853, 1.992380, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000] test_list.append(axi_rans_air_nozzle_restart) @@ -711,7 +711,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.455981, -1.541069] + bars_SST_2D.test_vals = [13.000000, -0.455981, -1.525412] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -773,7 +773,7 @@ def main(): fsi_cht_restart.cfg_dir = "fea_fsi/stat_fsi" fsi_cht_restart.cfg_file = "config_restart.cfg" fsi_cht_restart.test_iter = 0 - fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327033, -9.644570, -9.319618, 608.380000, -0.012974, 0.000000, 20.000000] + fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327047, -9.640535, -9.319270, 608.380000, -0.012974, 0.000000, 20.000000] fsi_cht_restart.multizone = True test_list.append(fsi_cht_restart) diff --git a/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg b/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg index f7798c39f785..fb0bbfd18360 100644 --- a/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg +++ b/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg @@ -20,7 +20,8 @@ MESH_FORMAT= CGNS TABULAR_FORMAT= TECPLOT CONV_FILENAME= history -RESTART_FILENAME= restart_flow +% Write a separate output so repeated runs use the same original restart input. +RESTART_FILENAME= restart_supersonic VOLUME_FILENAME= volume_flow SURFACE_FILENAME= surface_flow OUTPUT_WRT_FREQ= 500 diff --git a/TestCases/navierstokes/flatplate/lam_flatplate.cfg b/TestCases/navierstokes/flatplate/lam_flatplate.cfg index 77c342091f84..557ff6c0556b 100644 --- a/TestCases/navierstokes/flatplate/lam_flatplate.cfg +++ b/TestCases/navierstokes/flatplate/lam_flatplate.cfg @@ -88,7 +88,8 @@ MESH_FORMAT= SU2 SOLUTION_FILENAME= solution_flow SOLUTION_ADJ_FILENAME= solution_adj TABULAR_FORMAT= CSV -RESTART_FILENAME= restart_flow +% Keep the supersonic restart input intact when both tests run in this directory. +RESTART_FILENAME= restart_flatplate RESTART_ADJ_FILENAME= restart_adj GRAD_OBJFUNC_FILENAME= of_grad OUTPUT_WRT_FREQ= 250 diff --git a/TestCases/navierstokes/periodic2D/no_limiter_riemann.cfg b/TestCases/navierstokes/periodic2D/no_limiter_riemann.cfg new file mode 100644 index 000000000000..12f6b1e65333 --- /dev/null +++ b/TestCases/navierstokes/periodic2D/no_limiter_riemann.cfg @@ -0,0 +1,81 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Rotational periodic sector without slope limiter, % +% testing the Riemann outlet implicit matrix % +% Author: SU2 Contributors % +% Date: Oct 2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.1 +AOA= 22.5 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_PRESSURE= 99000 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1 +REF_AREA= 1 +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.87 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,45, 0,0,0 ) +MARKER_RIEMANN= ( outlet, STATIC_PRESSURE, 99000, 0, 0, 0, 0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300, 100000, 1, 0, 0 ) +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +% +MARKER_ANALYZE= ( outlet, inlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 100 +CFL_ADAPT= NO +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 0.1 +LINEAR_SOLVER_ITER= 4 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -11 +CONV_STARTITER= 0 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= sector.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_RES, LINSOL_RESIDUAL, SURFACE_PRESSURE_DROP ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ddb1c8049b58..88e9e7a00c78 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -274,7 +274,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-4.142624, 1.585638, -0.249527, 0.043953] + wedge.test_vals = [-4.188807, 1.541867, -0.249526, 0.043953] test_list.append(wedge) # ONERA M6 Wing @@ -328,7 +328,7 @@ def main(): ramp.cfg_dir = "euler/ramp" ramp.cfg_file = "inv_ramp.cfg" ramp.test_iter = 10 - ramp.test_vals = [-13.648343, -8.010490, -0.076277, 0.054839] + ramp.test_vals = [-13.648309, -8.010607, -0.076277, 0.054839] ramp.test_vals_aarch64 = [-13.648406, -8.014579, -0.076277, 0.054839] test_list.append(ramp) @@ -357,7 +357,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 100 - flatplate.test_vals = [-8.165135, -2.690978, 0.001084, 0.036277, 2.361500, -2.325300, 0.000000, 0.000000] + flatplate.test_vals = [-8.160319, -2.686084, 0.001084, 0.036277, 2.361500, -2.325300, 0.000000, 0.000000] test_list.append(flatplate) # Supersonic laminar flat plate @@ -365,7 +365,7 @@ def main(): flatplate_supersonic.cfg_dir = "navierstokes/flatplate" flatplate_supersonic.cfg_file = "flatplate_supersonic.cfg" flatplate_supersonic.test_iter = 100 - flatplate_supersonic.test_vals = [100.000000, -2.940836, -0.677210, -0.570631, 2.495392, 0.001469] + flatplate_supersonic.test_vals = [100.000000, -9.034577, -8.611322, -9.029955, -3.700944, 0.002511] test_list.append(flatplate_supersonic) # Custom objective function @@ -373,7 +373,7 @@ def main(): flatplate_udobj.cfg_dir = "user_defined_functions" flatplate_udobj.cfg_file = "lam_flatplate.cfg" flatplate_udobj.test_iter = 20 - flatplate_udobj.test_vals = [-0.001015, 596.450000, 299.550000, 296.900000, 21.318000, 0.586640, 36.553000, 2.188800] + flatplate_udobj.test_vals = [-0.001015, 596.450000, 299.560000, 296.900000, 21.317000, 0.586750, 36.553000, 2.188800] test_list.append(flatplate_udobj) # Probe performance test (11 probes, ADT path) @@ -381,7 +381,7 @@ def main(): probe_performance_11.cfg_dir = "user_defined_functions" probe_performance_11.cfg_file = "test_11_probes.cfg" probe_performance_11.test_iter = 4 - probe_performance_11.test_vals = [-6.290748, 101020, 101050, 99123] # RMS_DENSITY, probe1, probe6, probe11 + probe_performance_11.test_vals = [-6.290996, 101020.000000, 101050.000000, 99125.000000] # RMS_DENSITY, probe1, probe6, probe11 # Tolerances are typically 0.001 in TestCase.py test_list.append(probe_performance_11) @@ -435,15 +435,23 @@ def main(): periodic2d_no_limiter.cfg_dir = "navierstokes/periodic2D" periodic2d_no_limiter.cfg_file = "no_limiter.cfg" periodic2d_no_limiter.test_iter = 100 - periodic2d_no_limiter.test_vals = [-3.216881, -0.582813, -0.629049, 2.266141, -1.056877, -811.470000] + periodic2d_no_limiter.test_vals = [-3.217331, -0.583709, -0.629299, 2.265684, -1.057572, -810.150000] test_list.append(periodic2d_no_limiter) + # The fully coupled nonlinear boundary derivative used to diverge before iteration 100 at CFL 100. + periodic2d_no_limiter_riemann = TestCase('periodic2d_no_limiter_riemann') + periodic2d_no_limiter_riemann.cfg_dir = "navierstokes/periodic2D" + periodic2d_no_limiter_riemann.cfg_file = "no_limiter_riemann.cfg" + periodic2d_no_limiter_riemann.test_iter = 100 + periodic2d_no_limiter_riemann.test_vals = [-3.216878, -0.582314, -0.628636, 2.266530, -1.057965, -811.500000] + test_list.append(periodic2d_no_limiter_riemann) + # 2D rotational periodic sector with multigrid periodic2d_multigrid = TestCase('periodic2d_multigrid') periodic2d_multigrid.cfg_dir = "navierstokes/periodic2D" periodic2d_multigrid.cfg_file = "multigrid.cfg" periodic2d_multigrid.test_iter = 300 - periodic2d_multigrid.test_vals = [-4.422778, -1.615251, -1.453736, 1.047350, -1.410225, -2060.600000] + periodic2d_multigrid.test_vals = [-4.422875, -1.615366, -1.453862, 1.047230, -1.410284, -2060.900000] test_list.append(periodic2d_multigrid) # 2D rotational periodic sector with limiter @@ -451,7 +459,7 @@ def main(): periodic2d.cfg_dir = "navierstokes/periodic2D" periodic2d.cfg_file = "config.cfg" periodic2d.test_iter = 100 - periodic2d.test_vals = [-3.266248, -0.617784, -0.620252, 2.216476, -1.007761, -1035.000000] + periodic2d.test_vals = [-3.266432, -0.617956, -0.620044, 2.216317, -1.007813, -1035.000000] test_list.append(periodic2d) ########################## @@ -495,7 +503,7 @@ def main(): turb_flatplate_species.cfg_dir = "rans/flatplate" turb_flatplate_species.cfg_file = "turb_SA_flatplate_species.cfg" turb_flatplate_species.test_iter = 20 - turb_flatplate_species.test_vals = [-5.359043, -1.881674, -3.075839, 0.124416, -3.602253, 2.000000, -1.239286, 1.000000, -1.158672, 3.000000, -0.849003, 0.999876, 0.999876] + turb_flatplate_species.test_vals = [-5.359058, -1.881677, -3.075574, 0.124434, -3.602252, 2.000000, -1.239891, 1.000000, -1.158675, 3.000000, -0.849018, 0.999876, 0.999876] test_list.append(turb_flatplate_species) # Flat plate SST compressibility correction Wilcox @@ -519,7 +527,7 @@ def main(): turb_flatplate_sst_roughBCKnopp.cfg_dir = "rans/flatplate/roughness/bc_knopp" turb_flatplate_sst_roughBCKnopp.cfg_file = "turb_SST_flatplate_roughBCKnopp.cfg" turb_flatplate_sst_roughBCKnopp.test_iter = 10 - turb_flatplate_sst_roughBCKnopp.test_vals = [-5.058634, -2.460850, -2.847064, 0.447200, -2.595042, 1.497149, -0.188079, 0.004571] + turb_flatplate_sst_roughBCKnopp.test_vals = [-5.066141, -2.470493, -2.851854, 0.439238, -2.595926, 1.496827, -0.188074, 0.004571] test_list.append(turb_flatplate_sst_roughBCKnopp) # FLAT PLATE, ROUGHNESS BC AUPOIX SST @@ -527,7 +535,7 @@ def main(): turb_flatplate_sst_roughBCAupoix.cfg_dir = "rans/flatplate/roughness/bc_aupoix" turb_flatplate_sst_roughBCAupoix.cfg_file = "turb_SST_flatplate_roughBCAupoix.cfg" turb_flatplate_sst_roughBCAupoix.test_iter = 10 - turb_flatplate_sst_roughBCAupoix.test_vals = [-5.278728, -2.302982, -2.890419, 0.227585, -1.393786, 3.192601, -0.188729, 0.0071821] + turb_flatplate_sst_roughBCAupoix.test_vals = [-5.290263, -2.305835, -2.892120, 0.214802, -1.393793, 3.192602, -0.188752, 0.007182] test_list.append(turb_flatplate_sst_roughBCAupoix) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SA, restart from a solution without wall functions @@ -535,7 +543,7 @@ def main(): turb_wallfunction_flatplate_sa_restart.cfg_dir = "wallfunctions/flatplate/compressible_SA" turb_wallfunction_flatplate_sa_restart.cfg_file = "turb_SA_flatplate_restart.cfg" turb_wallfunction_flatplate_sa_restart.test_iter = 10 - turb_wallfunction_flatplate_sa_restart.test_vals = [-3.603550, -1.991111, -1.848187, 1.875995, -5.284151, 10.000000, -2.255938, -0.064941, 0.002054] + turb_wallfunction_flatplate_sa_restart.test_vals = [-3.603551, -1.991100, -1.848182, 1.875994, -5.284153, 10.000000, -2.255956, -0.064942, 0.002054] test_list.append(turb_wallfunction_flatplate_sa_restart) # ONERA M6 Wing @@ -660,7 +668,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-2.662880, 2.912015, -2.712454, 1.890642, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.662833, 2.912056, -2.707804, 1.889981, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -1292,7 +1300,7 @@ def main(): coolprop_fluidModel.cfg_dir = "nicf/coolprop" coolprop_fluidModel.cfg_file = "fluidModel.cfg" coolprop_fluidModel.test_iter = 5 - coolprop_fluidModel.test_vals = [-4.666146, -0.172864, 4.492817, 0.000000, 0.000000] + coolprop_fluidModel.test_vals = [-4.666868, -0.172864, 4.492817, 0.000000, 0.000000] coolprop_fluidModel.enabled_on_cpu_arch = ["x86_64"] test_list.append(coolprop_fluidModel) @@ -1301,7 +1309,7 @@ def main(): coolprop_transportModel.cfg_dir = "nicf/coolprop" coolprop_transportModel.cfg_file = "transportModel.cfg" coolprop_transportModel.test_iter = 5 - coolprop_transportModel.test_vals = [-4.667122, -0.172865, 4.687652, 0.000000, 0.000000] + coolprop_transportModel.test_vals = [-4.667851, -0.172865, 4.687645, 0.000000, 0.000000] coolprop_transportModel.enabled_on_cpu_arch = ["x86_64"] test_list.append(coolprop_transportModel) @@ -1310,7 +1318,7 @@ def main(): datadriven_fluidModel.cfg_dir = "nicf/datadriven" datadriven_fluidModel.cfg_file = "datadriven_nozzle.cfg" datadriven_fluidModel.test_iter = 50 - datadriven_fluidModel.test_vals = [-4.879421, -2.353319, -2.344004, 0.390613, -1.559800, 1.322761] + datadriven_fluidModel.test_vals = [-4.690176, -1.714239, -2.218210, 0.764338, -1.558140, 1.350242] test_list.append(datadriven_fluidModel) ###################################### @@ -1439,7 +1447,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047] + bars_SST_2D.test_vals = [13.000000, -0.457912, -1.525399] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1690,7 +1698,7 @@ def main(): pywrapper_unsteadyFSI.cfg_dir = "py_wrapper/dyn_fsi" pywrapper_unsteadyFSI.cfg_file = "config.cfg" pywrapper_unsteadyFSI.test_iter = 4 - pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 47, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156] + pywrapper_unsteadyFSI.test_vals = [0.000000, 31.000000, 5.000000, 48.000000, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156] pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") pywrapper_unsteadyFSI.unsteady = True pywrapper_unsteadyFSI.multizone = True diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index b80601e77efe..eb9cf6556295 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -150,7 +150,7 @@ def main(): discadj_axisymmetric_rans_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" discadj_axisymmetric_rans_nozzle.cfg_file = "air_nozzle_restart.cfg" discadj_axisymmetric_rans_nozzle.test_iter = 10 - discadj_axisymmetric_rans_nozzle.test_vals = [9.911310, 5.072379, 7.118130, 2.484452] + discadj_axisymmetric_rans_nozzle.test_vals = [9.911225, 5.072252, 7.118033, 2.484389] discadj_axisymmetric_rans_nozzle.no_restart = True test_list.append(discadj_axisymmetric_rans_nozzle) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 2f256854d5c0..9273b90e38f7 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -119,7 +119,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-4.399745, 1.331048, -0.249527, 0.043953] + wedge.test_vals = [-4.354414, 1.376401, -0.249526, 0.043953] test_list.append(wedge) # ONERA M6 Wing @@ -176,7 +176,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 20 - flatplate.test_vals = [-5.623336, -0.144687, 0.001202, 0.028141, 2.361500, -2.333400, 0.000000, 0.000000] + flatplate.test_vals = [-5.561710, -0.083360, 0.001203, 0.028157, 2.361500, -2.333400, 0.000000, 0.000000] test_list.append(flatplate) # Laminar cylinder (steady) @@ -216,7 +216,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003115, -7.626023, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003110, -7.625378, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -260,7 +260,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-5.286997, -6.483481, -0.187358, 0.003712, 10.000000, -2.181754] + turb_flatplate.test_vals = [-5.287013, -6.483448, -0.187358, 0.003711, 10.000000, -2.181875] test_list.append(turb_flatplate) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SST @@ -268,7 +268,7 @@ def main(): turb_wallfunction_flatplate_sst.cfg_dir = "wallfunctions/flatplate/compressible_SST" turb_wallfunction_flatplate_sst.cfg_file = "turb_SST_flatplate.cfg" turb_wallfunction_flatplate_sst.test_iter = 10 - turb_wallfunction_flatplate_sst.test_vals = [-4.036300, -1.916264, -1.821790, 1.443878, -1.579554, 1.521535, 10.000000, -2.351460, 0.030011, 0.002409] + turb_wallfunction_flatplate_sst.test_vals = [-4.036186, -1.915130, -1.821604, 1.444009, -1.579554, 1.521877, 10.000000, -2.348138, 0.029998, 0.002409] test_list.append(turb_wallfunction_flatplate_sst) # FLAT PLATE, ROUGHNESS BC WILCOX2006 SST @@ -276,7 +276,7 @@ def main(): turb_flatplate_sst_roughBCWilcox2006.cfg_dir = "rans/flatplate/roughness/bc_wilcox2006" turb_flatplate_sst_roughBCWilcox2006.cfg_file = "turb_SST_flatplate_roughBCWilcox2006.cfg" turb_flatplate_sst_roughBCWilcox2006.test_iter = 10 - turb_flatplate_sst_roughBCWilcox2006.test_vals = [-5.117900, -2.534224, -2.904279, 0.381710, -3.100346, 1.180161, -0.188797, 0.004029] + turb_flatplate_sst_roughBCWilcox2006.test_vals = [-5.120194, -2.539726, -2.906124, 0.378835, -3.104047, 1.179642, -0.188795, 0.004029] test_list.append(turb_flatplate_sst_roughBCWilcox2006) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SA @@ -284,7 +284,7 @@ def main(): turb_wallfunction_flatplate_sa.cfg_dir = "wallfunctions/flatplate/compressible_SA" turb_wallfunction_flatplate_sa.cfg_file = "turb_SA_flatplate.cfg" turb_wallfunction_flatplate_sa.test_iter = 10 - turb_wallfunction_flatplate_sa.test_vals = [-3.580017, -1.860982, -1.506109, 1.898317, -5.258602, 10.000000, -2.557123, -0.295793, 0.002533] + turb_wallfunction_flatplate_sa.test_vals = [-3.580020, -1.860973, -1.506099, 1.898314, -5.258586, 10.000000, -2.556814, -0.295784, 0.002533] test_list.append(turb_wallfunction_flatplate_sa) # ONERA M6 Wing @@ -375,7 +375,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-2.662936, 2.911899, -2.338516, 2.161131, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.662911, 2.911919, -2.337545, 2.161293, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -385,7 +385,7 @@ def main(): axi_rans_air_nozzle_species.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_species.cfg_file = "air_nozzle_species.cfg" axi_rans_air_nozzle_species.test_iter = 10 - axi_rans_air_nozzle_species.test_vals = [-1.676190, 3.896581, -2.912396, 5.762790, -3.565909, 0.000000] + axi_rans_air_nozzle_species.test_vals = [-1.666920, 3.899160, -2.874820, 5.767107, -3.524094, 0.000000] axi_rans_air_nozzle_species.tol = 0.0001 test_list.append(axi_rans_air_nozzle_species) @@ -928,7 +928,7 @@ def main(): ls89_sa.cfg_dir = "nicf/LS89" ls89_sa.cfg_file = "turb_SA_PR.cfg" ls89_sa.test_iter = 20 - ls89_sa.test_vals = [-5.041500, -13.379523, 0.179019, 0.420831] + ls89_sa.test_vals = [-5.041528, -13.379614, 0.179019, 0.420831] test_list.append(ls89_sa) # Rarefaction shock wave edge_VW @@ -1086,7 +1086,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051] + bars_SST_2D.test_vals = [13.000000, -0.456143, -1.525401] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1164,7 +1164,7 @@ def main(): fsi_cht.cfg_dir = "fea_fsi/stat_fsi" fsi_cht.cfg_file = "config.cfg" fsi_cht.test_iter = 20 - fsi_cht.test_vals = [5.000000, -5.077002, -5.379450, -9.247804, -9.320014, -9.185034, 608.350000, -0.012973, 0.000000, 30.000000] + fsi_cht.test_vals = [5.000000, -5.077000, -5.379479, -9.247837, -9.319619, -9.184898, 608.350000, -0.012973, 0.000000, 30.000000] fsi_cht.multizone = True test_list.append(fsi_cht) @@ -1759,7 +1759,7 @@ def main(): pywrapper_custom_inlet.cfg_dir = "py_wrapper/custom_inlet" pywrapper_custom_inlet.cfg_file = "lam_flatplate.cfg" pywrapper_custom_inlet.test_iter = 20 - pywrapper_custom_inlet.test_vals = [-4.123921, -1.544230, -3.886906, 1.338713, -0.688132, 0.160844, -0.009845, 0.520880, -0.530730] + pywrapper_custom_inlet.test_vals = [-4.123923, -1.544231, -3.887532, 1.338711, -0.688139, 0.160847, -0.009845, 0.520890, -0.530730] pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index d59aad17bf96..c4fcd01712a8 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -249,7 +249,7 @@ def main(): tutorial_trans_flatplate_T3A.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A" tutorial_trans_flatplate_T3A.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3A.test_iter = 20 - tutorial_trans_flatplate_T3A.test_vals = [-5.790137, -2.054834, -3.894659, -0.255074, -1.747220, 5.119341, -3.493237, 0.393262] + tutorial_trans_flatplate_T3A.test_vals = [-5.790325, -2.054834, -3.896287, -0.255197, -1.746976, 5.119339, -3.493229, 0.393180] tutorial_trans_flatplate_T3A.test_vals_aarch64 = [-5.808996, -2.070606, -3.969765, -0.277943, -1.953289, 1.708472, -3.514943, 0.357411] tutorial_trans_flatplate_T3A.no_restart = True test_list.append(tutorial_trans_flatplate_T3A) @@ -259,7 +259,7 @@ def main(): tutorial_trans_flatplate_T3Am.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A-" tutorial_trans_flatplate_T3Am.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3Am.test_iter = 20 - tutorial_trans_flatplate_T3Am.test_vals = [-5.587389, -1.700868, -3.093935, -0.102834, -3.750523, 3.287643, -2.394575, 1.119623] + tutorial_trans_flatplate_T3Am.test_vals = [-5.589599, -1.700867, -3.097326, -0.104799, -3.750523, 3.287643, -2.394576, 1.119623] tutorial_trans_flatplate_T3Am.test_vals_aarch64 = [-5.540938, -1.681627, -2.878831, -0.058224, -3.695533, 3.413628, -2.385345, 1.103633] tutorial_trans_flatplate_T3Am.no_restart = True test_list.append(tutorial_trans_flatplate_T3Am) @@ -295,7 +295,7 @@ def main(): tutorial_nicfd_nozzle.cfg_dir = "../Tutorials/compressible_flow/NICFD_nozzle" tutorial_nicfd_nozzle.cfg_file = "NICFD_nozzle.cfg" tutorial_nicfd_nozzle.test_iter = 20 - tutorial_nicfd_nozzle.test_vals = [-1.800072, -6.002691, 4.635204, 0.000000, 0.000000] + tutorial_nicfd_nozzle.test_vals = [-1.794899, -6.051798, 4.639477, 0.000000, 0.000000] tutorial_nicfd_nozzle.no_restart = True test_list.append(tutorial_nicfd_nozzle) diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 516f51cc5ebd..899baba67968 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -63,7 +63,7 @@ def main(): flatplate_sst1994m.cfg_dir = "vandv/rans/flatplate" flatplate_sst1994m.cfg_file = "turb_flatplate_sst.cfg" flatplate_sst1994m.test_iter = 5 - flatplate_sst1994m.test_vals = [-13.040943, -10.136971, -10.942764, -7.985181, -10.323857, -4.732487, 0.002801] + flatplate_sst1994m.test_vals = [-6.757337, -3.881283, -4.244894, -1.293607, -4.474870, -0.678315, 0.002801] flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807] test_list.append(flatplate_sst1994m) @@ -72,7 +72,7 @@ def main(): bump_sst1994m.cfg_dir = "vandv/rans/bump_in_channel" bump_sst1994m.cfg_file = "turb_bump_sst.cfg" bump_sst1994m.test_iter = 5 - bump_sst1994m.test_vals = [-11.928292, -10.095796, -9.512953, -6.445652, -11.774088, -6.988752, 0.004931] + bump_sst1994m.test_vals = [-8.032540, -5.141490, -6.166538, -2.601015, -6.571019, -2.377920, 0.004931] bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911] test_list.append(bump_sst1994m) @@ -81,7 +81,7 @@ def main(): swbli_sa.cfg_dir = "vandv/rans/swbli" swbli_sa.cfg_file = "config_sa.cfg" swbli_sa.test_iter = 5 - swbli_sa.test_vals = [-11.502718, -10.939184, -12.034284, -10.581169, -16.088844, 0.002242, -1.664946, 1.258100] + swbli_sa.test_vals = [-8.692514, -8.335554, -8.419342, -7.670662, -12.947166, 0.002242, -1.240847, 1.257600] swbli_sa.test_vals_aarch64 = [-11.504424, -10.941741, -12.049925, -10.586263, -16.090385, 0.002242, -1.614365, 1.340100] test_list.append(swbli_sa) @@ -91,7 +91,7 @@ def main(): swbli_sst.cfg_dir = "vandv/rans/swbli" swbli_sst.cfg_file = "config_sst.cfg" swbli_sst.test_iter = 5 - swbli_sst.test_vals = [-11.319578, -10.641523, -11.224600, -10.150215, -11.407538, -2.637660, 0.001816, -1.839484, -3.514593, 11.136000] + swbli_sst.test_vals = [-8.817539, -8.137335, -8.606270, -7.611652, -9.007191, -1.508190, 0.001816, -2.802542, -4.644043, 8.615500] test_list.append(swbli_sst) # DSMA661 - SA diff --git a/UnitTests/SU2_CFD/fluid/MLP_nitrogen_outlet.mlp b/UnitTests/SU2_CFD/fluid/MLP_nitrogen_outlet.mlp new file mode 100644 index 000000000000..33acbec8dbc6 --- /dev/null +++ b/UnitTests/SU2_CFD/fluid/MLP_nitrogen_outlet.mlp @@ -0,0 +1,78 @@ +
+ +[number of layers] +4 + +[neurons per layer] +2 +12 +12 +1 + +[activation function] +linear +exponential +exponential +linear + +[input names] +Density +Energy + +[input regularization method] +minmax + +[input normalization] ++2.1090497188197315e-02 +1.3149764159282715e+03 +-8.7440483449298932e+04 +1.9651855451196916e+06 + +[output names] +s + +[output regularization method] +minmax + +[output normalization] ++2.7967895102447687e+03 +9.7307708439237176e+03 + +
+ +[weights per layer] + +-3.6931488653013372e-01 -3.1474377747118677e+02 -7.8849044681531866e-01 -3.6267711872705242e+01 -9.4693133112296881e-04 -2.4759101201954605e+00 +3.6859588904967044e-01 -5.0599594919629425e+00 +9.8363564142113091e-01 -9.5973334878440486e-01 +1.9619213003259452e+00 -5.3480266381775685e+00 ++5.4575757021893469e-02 +2.2748894127143549e-03 +3.8016598848103861e-01 +1.8052223383547446e-02 -1.2347916434431229e-01 -7.7692862405548332e-02 +3.6960457731459249e-02 -1.0148725116734734e-01 -2.2166557551333326e+00 -6.6492982296293741e+00 -5.3598673815060827e-01 +3.9361084436121579e-01 + + +-6.0289903346074269e-01 -3.9951172103553416e+00 -2.2040250557871737e+00 -4.8092917866087673e-02 -3.5019717606530099e-01 -2.1178003722755276e-01 +3.3610626255824066e-01 -1.7572402599106512e-01 +1.3795282848988730e-01 +5.3542074183822723e-02 -1.1060525758624411e+00 -8.9126619859982104e-01 ++1.0999659702233211e+00 +3.9482740845858855e+00 +2.1715624069187149e+00 +1.4315443142891550e-02 -4.3513271093960628e-03 -4.4748521946946160e-02 -2.8331497763081557e-02 -1.9574676138448106e-02 -1.5964493668808994e-02 +1.8515158367985804e-01 +3.1523943980199180e-03 -1.0347253629184561e-01 +-1.8450855176577385e-01 -4.7305893876638888e-01 -2.2939956576919096e+00 +2.7972185892274029e-02 -2.5174517196989155e-01 +4.1674646590455910e-01 +2.4389381929091177e-02 -1.1506031345005816e+00 +1.4593750594266303e-01 -1.3341116852452784e-02 -1.6773663446073919e+00 -6.4540800321549163e-01 +-1.0786752489344005e+00 -1.5968939022876310e+00 -1.0681376729491869e+00 +1.5749813923225323e-01 +2.4718027905077058e-01 +2.6907670615987791e+00 -4.6742758504205534e-02 -3.5565514260358899e-01 +2.9143199538146930e-01 -1.8549853942551744e+00 +2.8711314462584575e-01 -1.6208994609294269e+00 +-1.2264765482666773e+00 -2.4674919552514600e+00 -1.3265176065977726e+00 -7.9141424729064969e-01 -4.0141907114642217e-01 -2.1563100057410614e-01 +1.3420268405649921e-01 -7.6076476638588331e-01 +6.3262617904080032e-01 +3.2728408194119374e-01 -3.4376106978785187e-01 +9.5964078326406876e-01 +-5.6554083347952928e-02 -2.8561417230466311e+00 -1.7262656034757942e+00 -7.9731568302412625e-02 -1.4701690023700852e-02 +6.1537787079645445e-01 +5.8596314725702459e-01 +5.2372641773803164e-01 -9.3584202734540056e-01 +1.9430032120400176e-01 +1.6558824068472527e-01 -3.3813471747014906e-01 +-1.6404418380135983e+00 -3.6256806347533748e+00 -1.9510896963160578e+00 +1.4047705287010601e-01 -7.3224816960376093e-01 -1.1816450873146218e+00 -7.4837642269864391e-01 +8.1551255144823184e-02 -5.5181791290108151e-01 -1.1824504105148357e+00 -1.3971674262310756e+00 -8.2542852042334747e-01 +-6.7607155586922496e-01 -3.3400299198433361e+00 -2.6375974912330991e+00 +1.6416639579357875e-01 -1.4220151573121304e-01 +3.7772129606488758e-01 -4.7228940157279992e-01 -7.4399289762720047e-02 -3.2290813450396694e-01 +4.3858239632455737e-01 +1.0164701735185362e-01 -3.3222358111657461e-02 +-3.4774715255690947e-01 +4.8950766001144164e-01 -3.4378531837692144e-01 -1.1679338423937842e+00 +1.3849313394711733e+00 +3.2393331547257442e-01 -1.2946005564391727e+00 +5.1115287970759704e-01 -2.5785233904399696e-01 +6.6021164623134965e-01 -4.7439733068883655e-01 -1.0192935577129223e+00 ++6.7849424196606736e-02 +5.5134370817853919e-02 +2.5319641231188278e-01 -1.6444969219945896e-01 +1.4176720779837590e+00 -3.0121917899206711e-02 +7.6592053680163308e-01 +3.4498207782656359e-01 +4.8360343008648565e-01 -7.3603980293366553e-02 +4.9794189079347193e+00 -6.1550558220266955e+00 +-4.4536951029960153e-01 -5.1374824869859192e-01 -1.5830799746340365e+00 -3.1825549970147937e-01 -1.4138701094696484e+00 -1.7480260116648119e+00 -1.7371258820972322e+00 +2.1697392732932461e-01 -1.8760676798491716e+00 -1.5271356949826242e+00 +1.3567944311103606e+00 +3.3936225827318928e-01 +-2.0340884443227991e-01 -1.9875192710887564e+00 -1.3690603946085547e+00 +2.3618086548763309e-01 +9.5711434428773612e-02 +2.3086929427470457e-01 -4.3063327262015250e-01 -8.5660991011086729e-01 +1.2072540292660496e-04 +9.4544676622201873e-01 +2.7549871149128580e-01 -2.8035823287386563e-01 + + ++6.0375568451703152e-02 ++2.7623189372863434e-02 +-6.4576408824258202e-01 ++3.4502521675902814e-01 +-1.1677097032696395e+00 ++4.4896276611858071e-01 ++7.0784230816897908e-01 +-3.7543360726436248e-01 ++8.2313797915585174e-01 ++5.5378236116127455e-01 +-4.4541146539877080e-01 ++1.2768696085188918e+00 + + +[biases per layer] ++0.0000000000000000e+00 +0.0000000000000000e+00 ++1.1207495814336845e+00 +2.5533329348987270e+00 -7.0659576209784636e-01 -3.8615060918494654e-01 -5.9507558268617333e-01 +7.8985446383846680e-01 +1.2916322810880783e+00 +1.1508396571030062e+00 -6.2407788653699292e-01 +6.9955173074603794e-01 -8.5587392374930227e-01 -1.0952945666386391e+00 +-1.9636295603087286e+00 -4.0044919187392631e+00 -8.8524618313204073e-01 -5.9839400291935296e-01 +1.8819710713781790e-01 -7.3343736858168218e-01 -2.7315727055497185e-01 -6.3204585891050125e-01 -2.4026293759253928e-01 -6.5513899674899456e-01 -5.7887459135893582e-01 +1.5254915040858916e-01 ++4.3424817896574419e-01 diff --git a/UnitTests/SU2_CFD/solvers/characteristic_pressure_outlet_tests.cpp b/UnitTests/SU2_CFD/solvers/characteristic_pressure_outlet_tests.cpp new file mode 100644 index 000000000000..b22088603f89 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/characteristic_pressure_outlet_tests.cpp @@ -0,0 +1,264 @@ +/*! + * \file characteristic_pressure_outlet_tests.cpp + * \brief Regression for characteristic pressure boundary admissibility. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "pressure_outlet_test_utils.hpp" +#include "../../UnitQuadTestCase.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +#include "../../../SU2_CFD/include/fluid/CFluidModel.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../../SU2_CFD/include/fluid/CVanDerWaalsGas.hpp" +#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp" +#include "../../../SU2_CFD/include/fluid/CDataDrivenFluid.hpp" + +namespace { +void CharacteristicPressureFixture(su2double pressure, bool backflow = false) { + UnitQuadTestCase test; + std::ostringstream options; + options << "SOLVER= EULER\nMESH_FORMAT= RECTANGLE\nMESH_BOX_SIZE= 5,129,5\n" + << "REF_DIMENSIONALIZATION= DIMENSIONAL\nMACH_NUMBER= 0.1\n" + << "FREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 230\nFREESTREAM_PRESSURE= 1500112\n" + << "FLUID_MODEL= PR_GAS\nGAMMA_VALUE= 1.0125\nGAS_CONSTANT= 35.152\n" + << "CRITICAL_PRESSURE= 1415200\nCRITICAL_TEMPERATURE= 564.1\nACENTRIC_FACTOR= 0.529\n" + << "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n" + << "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 20,0,0\n" + << "MARKER_EULER= (x_minus,y_minus,y_plus)\nMARKER_MONITORING= (x_minus)\n" + << "MARKER_RIEMANN= (x_plus,STATIC_PRESSURE," << pressure << ",0,0,0,0)\n"; + test.config_options = options.str(); + test.InitConfig(); + test.InitGeometry(); + test.InitSolver(); + auto* solver = test.solver[FLOW_SOL]; + auto* nodes = static_cast(solver->GetNodes()); + CPengRobinson gas(1.0125, 35.152, 1415200, 564.1, 0.529); + // Captured conservative startup state: low back pressure previously produced an invalid EOS state. + su2double state[4] = {139.72158881490299, 11271.637336535776, -2338.5053086132339, 221781581.74179047}; + if (backflow) { + const auto oldVelocity = state[1] / state[0]; + state[1] = state[0] * -200; + // Change only normal velocity; adjust total energy to keep the original thermal state. + state[3] += .5 * state[0] * (40000 - oldVelocity * oldVelocity); + } + const su2double grid[2] = {20, 0}; + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) { + nodes->SetSolution(p, state); + nodes->SetSolution_Old(p, state); + REQUIRE(nodes->SetPrimVar(p, &gas)); + nodes->SetSecondaryVar(p, &gas); + test.geometry->nodes->SetGridVel(p, grid); + for (auto k = 0u; k < 4; ++k) solver->LinSysRes(p, k) = 0; + su2double zero[16] = {}; + solver->Jacobian.SetBlock(p, p, zero); + } + CUpwGeneralHLLC_Flow numerics(2, 4, test.config.get()); + for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker) + if (test.config->GetMarker_All_TagBound(marker) == "x_plus") SU2_OMP_PARALLEL { + CUpwGeneralHLLC_Flow threadNumerics(2, 4, test.config.get()); + solver->BC_Riemann(test.geometry.get(), test.solver, &threadNumerics, nullptr, test.config.get(), marker); + } + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) + for (auto k = 0u; k < 4; ++k) CHECK(std::isfinite(SU2_TYPE::GetValue(solver->LinSysRes(p, k)))); + for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker) { + if (test.config->GetMarker_All_TagBound(marker) != "x_plus") continue; + for (auto v = 0ul; v < test.geometry->nVertex[marker]; ++v) { + const auto point = test.geometry->vertex[marker][v]->GetNode(); + if (!test.geometry->nodes->GetDomain(point)) continue; + su2double normal[2], unit[2]; + test.geometry->vertex[marker][v]->GetNormal(normal); + const auto area = GeometryToolbox::Norm(2, normal); + for (auto d = 0u; d < 2; ++d) unit[d] = (normal[d] = -normal[d]) / area; + // BC_Riemann uses the physical boundary flux, including transport by the moving grid. + su2double boundary[4]; + REQUIRE(gas.ComputePressureOutletState(2, state, unit, GeometryToolbox::DotProduct(2, grid, unit), pressure, + boundary) == CFluidModel::PressureOutletStatus::Valid); + auto flux = [&](const su2double* boundary, su2double* output) { + su2double velocity[2]; + for (auto d = 0u; d < 2; ++d) velocity[d] = boundary[d + 1] / boundary[0]; + const auto energy = boundary[3] / boundary[0]; + gas.SetTDState_rhoe(boundary[0], energy - .5 * GeometryToolbox::SquaredNorm(2, velocity)); + const auto p = gas.GetPressure(), h = energy + p / boundary[0]; + numerics.GetInviscidProjFlux(&boundary[0], velocity, &p, &h, normal, output); + const auto speed = GeometryToolbox::DotProduct(2, grid, normal); + for (auto k = 0u; k < 4; ++k) output[k] -= speed * boundary[k]; + }; + su2double reference[4]; + flux(boundary, reference); + const auto* primitive = solver->GetCharacPrimVar(marker, v); + CHECK(SU2_TYPE::GetValue(primitive[0]) == Approx(SU2_TYPE::GetValue(gas.GetTemperature())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[3]) == Approx(SU2_TYPE::GetValue(gas.GetPressure())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[4]) == Approx(SU2_TYPE::GetValue(gas.GetDensity())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[6]) == Approx(SU2_TYPE::GetValue(gas.GetSoundSpeed())).epsilon(1e-8)); + for (auto k = 0u; k < 4; ++k) + CHECK(SU2_TYPE::GetValue(solver->LinSysRes(point, k)) == + Approx(SU2_TYPE::GetValue(reference[k])).epsilon(1e-8)); + for (auto column = 0u; column < 4; ++column) { + // Freeze the wave basis and differentiate the physical boundary flux along outgoing increments. + su2double plus[4], minus[4], fPlus[4], fMinus[4], increment[4]; + OutgoingOutletIncrement(gas, 2, state, unit, GeometryToolbox::DotProduct(2, grid, unit), 0, column, increment); + const auto delta = 2e-5 * fabs(state[column]); + for (auto k = 0u; k < 4; ++k) { + plus[k] = boundary[k] + delta * increment[k]; + minus[k] = boundary[k] - delta * increment[k]; + } + flux(plus, fPlus); + flux(minus, fMinus); + for (auto row = 0u; row < 4; ++row) { + const auto expected = SU2_TYPE::GetValue((fPlus[row] - fMinus[row]) / (2 * delta)); + CHECK(solver->Jacobian.GetBlock(point, point)[row * 4 + column] == + Approx(expected).epsilon(3e-4).margin(1e-5)); + } + } + } + } +} +} // namespace + +TEST_CASE("Characteristic pressure boundary reconstructs the captured startup states", "[Pressure-outlet]") { + for (const auto pressure : {1400000., 100000., 1000.}) CharacteristicPressureFixture(pressure); +} + +TEST_CASE("Nonlinear pressure outflow preserves the acoustic invariant and frame", "[Pressure-outlet]") { + for (auto dim : {2u, 3u}) { + for (auto kind : {0u, 1u, 2u}) { + std::unique_ptr gas; + if (kind == 0) gas.reset(new CIdealGas(1.4, 35.152)); + if (kind == 1) gas.reset(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, .529)); + if (kind == 2) gas.reset(new CVanDerWaalsGas(1.0125, 35.152, 1415200, 564.1)); + const su2double p = 1500112, rho = 230, normal[3] = {1, 0, 0}, boost[3] = {20, -7, 3}; + gas->SetTDState_Prho(p, rho); + const auto e = gas->GetStaticEnergy(), sound = gas->GetSoundSpeed(), entropy = gas->GetEntropy(); + // The smallest ratio is 1000/1500112, reproducing the startup expansion past the sonic limit. + for (const auto ratio : {1., .9, 1.05, .0006666168959384}) { + su2double out[2][5] = {}; + // Boost both flow and grid: the relative wave speeds and reconstructed thermodynamics must agree. + for (auto frame = 0u; frame < 2; ++frame) { + su2double velocity[3] = {.3 * sound, 5, -2}, state[5] = {rho}; + for (auto d = 0u; d < dim; ++d) { + if (frame) velocity[d] += boost[d]; + state[d + 1] = rho * velocity[d]; + } + state[dim + 1] = rho * (e + .5 * GeometryToolbox::SquaredNorm(dim, velocity)); + REQUIRE(gas->ComputePressureOutletState(dim, state, normal, frame ? boost[0] : 0, ratio * p, out[frame]) == + CFluidModel::PressureOutletStatus::Valid); + const auto density = out[frame][0]; + su2double v[3] = {}; + for (auto d = 0u; d < dim; ++d) v[d] = out[frame][d + 1] / density; + gas->SetTDState_rhoe(density, out[frame][dim + 1] / density - .5 * GeometryToolbox::SquaredNorm(dim, v)); + const auto q = v[0] - (frame ? boost[0] : 0); + CHECK(density > 0); + CHECK(gas->GetTemperature() > 0); + CHECK(gas->GetSoundSpeed2() > 0); + CHECK(SU2_TYPE::GetValue(gas->GetEntropy()) == Approx(SU2_TYPE::GetValue(entropy)).margin(1e-5)); + if (ratio > .1) { + CHECK(SU2_TYPE::GetValue(gas->GetPressure()) == Approx(SU2_TYPE::GetValue(ratio * p)).epsilon(1e-7)); + CHECK(q < gas->GetSoundSpeed()); + } else { + CHECK(gas->GetPressure() > ratio * p); + CHECK(SU2_TYPE::GetValue(q) == Approx(SU2_TYPE::GetValue(gas->GetSoundSpeed())).epsilon(1e-7)); + } + if (kind == 0) { + // Gamma=1.4 gives 2/(gamma-1)=5, so q+5*a provides an analytic integration check. + const su2double aRatio = ratio > .1 ? pow(ratio, .4 / 2.8) : (.4 * .3 + 2) / 2.4; + CHECK(SU2_TYPE::GetValue(density) == Approx(SU2_TYPE::GetValue(rho * pow(aRatio, 5))).epsilon(1e-7)); + CHECK(SU2_TYPE::GetValue(q + 5 * gas->GetSoundSpeed()) == + Approx(SU2_TYPE::GetValue(5.3 * sound)).epsilon(1e-7)); + } + } + CHECK(SU2_TYPE::GetValue(out[1][0]) == Approx(SU2_TYPE::GetValue(out[0][0])).epsilon(1e-9)); + // In the boosted frame, rho*E gains boost dot momentum + rho*|boost|^2/2. + su2double energyBoost = 0, boost2 = 0; + for (auto d = 0u; d < dim; ++d) { + CHECK(SU2_TYPE::GetValue(out[1][d + 1]) == + Approx(SU2_TYPE::GetValue(out[0][d + 1] + boost[d] * out[0][0])).epsilon(1e-8)); + energyBoost += boost[d] * out[0][d + 1]; + boost2 += boost[d] * boost[d]; + } + CHECK(SU2_TYPE::GetValue(out[1][dim + 1]) == + Approx(SU2_TYPE::GetValue(out[0][dim + 1] + energyBoost + .5 * boost2 * out[0][0])).epsilon(1e-8)); + } + su2double state[5] = {rho}, unchanged[5] = {}; + state[1] = 2 * sound * rho; + state[dim + 1] = rho * (e + 2 * sound * sound); + REQUIRE(gas->ComputePressureOutletState(dim, state, normal, 0, 1000, unchanged) == + CFluidModel::PressureOutletStatus::Valid); + for (auto v = 0u; v < dim + 2; ++v) CHECK(unchanged[v] == state[v]); + // Frozen SST k must remain in total energy without entering the EOS. + state[1] = .3 * sound * rho; + state[dim + 1] = rho * (e + .5 * .3 * .3 * sound * sound); + su2double plain[5], withK[5]; + REQUIRE(gas->ComputePressureOutletState(dim, state, normal, 0, .9 * p, plain) == + CFluidModel::PressureOutletStatus::Valid); + state[dim + 1] += 17.5 * rho; + REQUIRE(gas->ComputePressureOutletState(dim, state, normal, 0, .9 * p, withK, 17.5) == + CFluidModel::PressureOutletStatus::Valid); + for (auto v = 0u; v < dim + 1; ++v) CHECK(withK[v] == Approx(SU2_TYPE::GetValue(plain[v])).epsilon(1e-8)); + CHECK(withK[dim + 1] == Approx(SU2_TYPE::GetValue(plain[dim + 1] + 17.5 * plain[0])).epsilon(1e-8)); + // Reversing subsonic flow retains the pressure-only closure; fully incoming flow needs inlet data. + state[1] = -state[1]; + CHECK(gas->ComputePressureOutletState(dim, state, normal, 0, p, unchanged, 17.5) == + CFluidModel::PressureOutletStatus::Valid); + state[1] = -2 * sound * rho; + state[dim + 1] = rho * (e + 2 * sound * sound + 17.5); + CHECK(gas->ComputePressureOutletState(dim, state, normal, 0, p, unchanged, 17.5) == + CFluidModel::PressureOutletStatus::Backflow); + } + } +} + +// Hidden Catch cases terminate the process; the Python checker runs them separately to inspect the error. +TEST_CASE("Characteristic pressure boundary rejects nonpositive pressure", "[.][Boundary-invalid-pressure]") { + CharacteristicPressureFixture(0); +} +TEST_CASE("Characteristic pressure boundary requires data for backflow", "[.][Boundary-backflow]") { + CharacteristicPressureFixture(1400000, true); +} + +#ifdef HAVE_MLPCPP +TEST_CASE("Pressure outlet follows the learned EOS without pressure-density inversion", "[Pressure-outlet]") { + // MLP_nitrogen_outlet.mlp maps (density, thermal energy) to entropy; its derivatives define the EOS. + // This fixture checks that outlet integration can use that native input pair without a (p,rho) solve. + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= ROE\nFLUID_MODEL= DATADRIVEN_FLUID\nINTERPOLATION_METHOD= " + "MLP\nUSE_PINN= YES\n" + "FILENAMES_INTERPOLATOR= (src/SU2/UnitTests/SU2_CFD/fluid/MLP_nitrogen_outlet.mlp)\n"; + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CDataDrivenFluid gas(&config, false); + gas.SetTDState_rhoe(16, 130000); + const su2double pressure = gas.GetPressure(), sound = gas.GetSoundSpeed(), entropy = gas.GetEntropy(); + const su2double normal[2] = {1, 0}; + for (const auto mach : {-.3, 0., .3}) { + su2double state[4] = {16, 16 * mach * sound, 0, 16 * (130000 + .5 * mach * mach * sound * sound)}, boundary[4]; + REQUIRE(gas.ComputePressureOutletState(2, state, normal, 0, .98 * pressure, boundary) == + CFluidModel::PressureOutletStatus::Valid); + const su2double velocity = boundary[1] / boundary[0]; + gas.SetTDState_rhoe(boundary[0], boundary[3] / boundary[0] - .5 * velocity * velocity); + CHECK(SU2_TYPE::GetValue(gas.GetPressure()) == Approx(SU2_TYPE::GetValue(.98 * pressure)).epsilon(1e-7)); + CHECK(SU2_TYPE::GetValue(gas.GetEntropy()) == Approx(SU2_TYPE::GetValue(entropy)).margin(1e-5)); + CHECK(gas.GetTemperature() > 0); + CHECK(gas.GetSoundSpeed2() > 0); + } +} +#endif diff --git a/UnitTests/SU2_CFD/solvers/check_characteristic_boundary_admissibility.py b/UnitTests/SU2_CFD/solvers/check_characteristic_boundary_admissibility.py new file mode 100644 index 000000000000..91e68fac088c --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/check_characteristic_boundary_admissibility.py @@ -0,0 +1,36 @@ +import subprocess +import sys + +# SU2_MPI::Error exits the test process, so run each hidden Catch case in isolation. +# Check the boundary diagnostic as well as the exit code: a later flux NaN is a different failure. +for tag, expected in [ + ( + "Turbo-invalid-pressure", + "Unable to construct an admissible static-pressure outflow state", + ), + ("Turbo-backflow", "requires an inflow state for backflow"), + ( + "Ordinary-invalid-pressure", + "Unable to construct an admissible static-pressure outflow state", + ), + ("Ordinary-backflow", "requires an inflow state for backflow"), + ( + "Boundary-invalid-pressure", + "Unable to construct an admissible static-pressure outflow state", + ), + ( + "Boundary-backflow", + "requires an inflow state for backflow", + ), +]: + result = subprocess.run( + [sys.argv[1], f"[{tag}]", "-r", "compact"], + capture_output=True, + text=True, + timeout=60, + ) + output = result.stdout + result.stderr + assert result.returncode != 0, output + assert expected in output, output + assert "NaN detected" not in output, output +print("Invalid characteristic boundary states are rejected before flux assembly.") diff --git a/UnitTests/SU2_CFD/solvers/pressure_outlet_ad_tests.cpp b/UnitTests/SU2_CFD/solvers/pressure_outlet_ad_tests.cpp new file mode 100644 index 000000000000..548c2700f96f --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/pressure_outlet_ad_tests.cpp @@ -0,0 +1,96 @@ +/*! + * \file pressure_outlet_ad_tests.cpp + * \brief Differentiation of subsonic, sonic-limited and supersonic pressure outlets. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "../../../SU2_CFD/include/fluid/CFluidModel.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../../SU2_CFD/include/fluid/CIdealGas.hpp" +#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp" +#include "../../../SU2_CFD/include/fluid/CVanDerWaalsGas.hpp" + +TEST_CASE("Pressure outlet state derivatives match finite differences", "[Pressure-outlet-AD]") { + // Check the analytic ideal-gas path and integrated Peng-Robinson/van der Waals paths with the same inputs. + for (auto kind : {0u, 1u, 2u}) { + std::unique_ptr gas; + if (kind == 0) gas.reset(new CIdealGas(1.4, 35.152)); + if (kind == 1) gas.reset(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, .529)); + if (kind == 2) gas.reset(new CVanDerWaalsGas(1.0125, 35.152, 1415200, 564.1)); + const double rho = 230, p = 1500112; + gas->SetTDState_Prho(p, rho); + const double e = SU2_TYPE::GetValue(gas->GetStaticEnergy()), a = SU2_TYPE::GetValue(gas->GetSoundSpeed()); + const su2double normal[2] = {1, 0}; + // Cover reverse subsonic flow, ordinary outflow, near-sonic flow and supersonic extrapolation; + // the lowest pressure ratio makes subsonic outward expansion stop at the moving sonic limit. + for (const auto mach : {-.3, .3, .99, 1.4}) + for (const auto ratio : {1., .95, .001}) { + const double values[6] = { + rho, rho * (mach * a + 20), rho * 5, rho * (e + .5 * (pow(mach * a + 20, 2) + 25)), ratio * p, 20}; + // The first four inputs are conservative variables; the last two are pressure and grid speed. + auto evaluate = [&](su2double* input, su2double* output) { + return gas->ComputePressureOutletState(2, input, normal, input[5], input[4], output, 0, kind == 0 ? 1.4 : 0.); + }; + for (auto column = 0u; column < 6; ++column) { + const double scale[6] = {rho, rho * a, rho * a, rho * a * a, p, a}; + const double delta = 2e-5 * fmax(fabs(values[column]), scale[column]); + su2double plus[6], minus[6], fp[4], fm[4]; + for (auto k = 0u; k < 6; ++k) plus[k] = minus[k] = values[k]; + plus[column] += delta; + minus[column] -= delta; + REQUIRE(evaluate(plus, fp) == CFluidModel::PressureOutletStatus::Valid); + REQUIRE(evaluate(minus, fm) == CFluidModel::PressureOutletStatus::Valid); + // Finite differences rerun the complete map, so the reference includes movement of the sonic root. + for (auto row = 0u; row < 4; ++row) { + const double expected = SU2_TYPE::GetValue((fp[row] - fm[row]) / (2 * delta)); + su2double input[6], output[4]; + for (auto k = 0u; k < 6; ++k) input[k] = values[k]; +#ifdef CODI_REVERSE_TYPE + // Record one scalar input/output pair per tape to keep sensitivities independent. + AD::Reset(); + AD::StartRecording(); + AD::RegisterInput(input[column]); +#else + // Forward mode propagates a unit seed for the selected input column. + SU2_TYPE::SetDerivative(input[column], 1); +#endif + REQUIRE(evaluate(input, output) == CFluidModel::PressureOutletStatus::Valid); +#ifdef CODI_REVERSE_TYPE + AD::RegisterOutput(output[row]); + AD::StopRecording(); + SU2_TYPE::SetDerivative(output[row], 1); + AD::ComputeAdjoint(); + const double derivative = SU2_TYPE::GetDerivative(input[column]); + AD::Reset(); +#else + const double derivative = SU2_TYPE::GetDerivative(output[row]); +#endif + INFO("gas=" << kind << " Mach=" << mach << " pressure ratio=" << ratio << " row=" << row + << " column=" << column); + CHECK(derivative == Approx(expected).epsilon(3e-4).margin(2e-4)); + } + } + } + } +} diff --git a/UnitTests/SU2_CFD/solvers/pressure_outlet_boundary_tests.cpp b/UnitTests/SU2_CFD/solvers/pressure_outlet_boundary_tests.cpp new file mode 100644 index 000000000000..b1c5a0fd1f23 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/pressure_outlet_boundary_tests.cpp @@ -0,0 +1,172 @@ +/*! + * \file pressure_outlet_boundary_tests.cpp + * \brief Pressure outlet classification and boundary Jacobian regression. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include "pressure_outlet_test_utils.hpp" +#include "../../UnitQuadTestCase.hpp" +#include "../../../SU2_CFD/include/fluid/CIdealGas.hpp" +#include "../../../SU2_CFD/include/variables/CNSVariable.hpp" + +namespace { +// Unit flux matrices expose I + the outgoing-wave projection in the assembled implicit block. +class OutletStateProbe final : public CNumerics { + su2double state[5] = {}, identity[5][5] = {}; + const su2double* rows[5]; + + public: + OutletStateProbe(unsigned short dim, const CConfig* config) : CNumerics(dim, dim + 2, config) { + for (auto k = 0u; k < nVar; ++k) { + rows[k] = identity[k]; + identity[k][k] = 1; + } + } + ResidualType<> ComputeResidual(const CConfig*) override { + const su2double rho = V_j[nDim + 2], p = V_j[nDim + 1]; + state[0] = rho; + for (auto d = 0u; d < nDim; ++d) state[d + 1] = rho * V_j[d + 1]; + state[nDim + 1] = rho * V_j[nDim + 3] - p; + return {state, rows, rows}; + } +}; + +void IdealOutlet(unsigned short dim, double normalMach, double pressureRatio, bool moving, bool threaded = false, + double turbulentK = 0) { + UnitQuadTestCase test; + std::ostringstream options; + options << "SOLVER= " << (turbulentK ? "RANS" : "EULER") << "\nMESH_FORMAT= " << (dim == 2 ? "RECTANGLE" : "BOX") + << "\nMESH_BOX_SIZE= " << (threaded ? "5,17,17" : "5,5,5") << "\nREF_DIMENSIONALIZATION= DIMENSIONAL\n" + << "MACH_NUMBER= 0.3\nFLUID_MODEL= IDEAL_GAS\nGAMMA_VALUE= 1.4\nGAS_CONSTANT= 287\n" + << "FREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 1.2\nFREESTREAM_PRESSURE= 100000\n" + << "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n" + << "MARKER_EULER= (x_minus,y_minus,y_plus" << (dim == 3 ? ",z_minus,z_plus" : "") + << ")\nMARKER_MONITORING= (x_minus)\nMARKER_OUTLET= (x_plus," << pressureRatio * 100000 << ")\n"; + if (moving) options << "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 500,240,0\n"; + if (turbulentK) + options << "KIND_TURB_MODEL= SST\nINIT_OPTION= TD_CONDITIONS\nREYNOLDS_NUMBER= 1E6\n" + "VISCOSITY_MODEL= CONSTANT_VISCOSITY\nCONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY\n"; + test.config_options = options.str(); + test.InitConfig(); + test.InitGeometry(); + test.InitSolver(); + auto* solver = test.solver[FLOW_SOL]; + auto* nodes = static_cast(solver->GetNodes()); + CIdealGas gas(1.4, 287); + if (turbulentK) { + gas.SetLaminarViscosityModel(test.config.get()); + gas.SetThermalConductivityModel(test.config.get()); + } + const double rho = 1.2, p = 100000, a = sqrt(1.4 * p / rho); + const su2double grid[3] = {moving ? 500. : 0., moving ? 240. : 0., 0.}; + // Tangential Mach 1.5 makes total speed supersonic even when normalMach is subsonic; + // the grid boost also distinguishes lab-frame speed from the relative normal wave speed. + su2double state[5] = {rho}, velocity[3] = {normalMach * a + grid[0], 1.5 * a + grid[1], 0}; + for (auto d = 0u; d < dim; ++d) state[d + 1] = rho * velocity[d]; + state[dim + 1] = p / .4 + rho * (.5 * GeometryToolbox::SquaredNorm(dim, velocity) + turbulentK); + for (auto point = 0ul; point < test.geometry->GetnPoint(); ++point) { + nodes->SetSolution(point, state); + if (turbulentK) { + test.solver[TURB_SOL]->GetNodes()->SetSolution(point, 0, turbulentK); + REQUIRE(static_cast(nodes)->SetPrimVar(point, 0, turbulentK, &gas)); + } else { + REQUIRE(nodes->SetPrimVar(point, &gas)); + } + if (moving) test.geometry->nodes->SetGridVel(point, grid); + for (auto k = 0u; k < dim + 2; ++k) solver->LinSysRes(point, k) = 0; + su2double zero[25] = {}; + solver->Jacobian.SetBlock(point, point, zero); + } + for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker) { + if (test.config->GetMarker_All_TagBound(marker) != "x_plus") continue; + SU2_OMP_PARALLEL { + OutletStateProbe probe(dim, test.config.get()); + solver->BC_Outlet(test.geometry.get(), test.solver, &probe, nullptr, test.config.get(), marker); + } + // Independent ideal-gas reference: gamma=1.4 gives q+5*a=constant and rho proportional to a^5. + auto reference = [&](const su2double* u, su2double* out) { + const su2double density = u[0], q = u[1] / density - grid[0]; + su2double v[3] = {}; + for (auto d = 0u; d < dim; ++d) v[d] = u[d + 1] / density; + const su2double pressure = .4 * (u[dim + 1] - density * (.5 * GeometryToolbox::SquaredNorm(dim, v) + turbulentK)); + const su2double sound = sqrt(1.4 * pressure / density); + if (q >= sound) { + for (auto k = 0u; k < dim + 2; ++k) out[k] = u[k]; + return; + } + const su2double sonic = (.4 * q + 2 * sound) / 2.4; + const su2double target = sound * pow(pressureRatio * p / pressure, .4 / 2.8); + const su2double boundarySound = fmax(target, sonic); + out[0] = density * pow(boundarySound / sound, 5); + v[0] += 5 * (sound - boundarySound); + for (auto d = 0u; d < dim; ++d) out[d + 1] = out[0] * v[d]; + out[dim + 1] = + out[0] * (boundarySound * boundarySound / .56 + .5 * GeometryToolbox::SquaredNorm(dim, v) + turbulentK); + }; + su2double expected[5]; + reference(state, expected); + for (auto vertex = 0ul; vertex < test.geometry->nVertex[marker]; ++vertex) { + const auto point = test.geometry->vertex[marker][vertex]->GetNode(); + if (!test.geometry->nodes->GetDomain(point)) continue; + for (auto row = 0u; row < dim + 2; ++row) { + CHECK(SU2_TYPE::GetValue(solver->LinSysRes(point, row)) == + Approx(SU2_TYPE::GetValue(expected[row])).epsilon(1e-8).margin(1e-8)); + } + const su2double normal[3] = {1, 0, 0}; + for (auto column = 0u; column < dim + 2; ++column) { + su2double increment[5]; + OutgoingOutletIncrement(gas, dim, state, normal, grid[0], turbulentK, column, increment); + for (auto row = 0u; row < dim + 2; ++row) { + CHECK(solver->Jacobian.GetBlock(point, point)[row * (dim + 2) + column] == + Approx(SU2_TYPE::GetValue(increment[row]) + (row == column ? 1. : 0.)).epsilon(1e-8).margin(1e-8)); + } + } + } + } +} +} // namespace + +TEST_CASE("Ordinary pressure outlet uses relative normal characteristics and an outgoing-wave matrix", + "[Pressure-outlet]") { + for (auto dim : {2u, 3u}) + for (const auto moving : {false, true}) + for (const auto normalMach : {-.3, 0., .3, .99, 1.4}) + for (const auto ratio : {1., .95, 1.02, .001}) IdealOutlet(dim, normalMach, ratio, moving); +} + +TEST_CASE("Ordinary pressure outlet rejects nonpositive pressure", "[.][Ordinary-invalid-pressure]") { + IdealOutlet(2, .3, 0, false); +} +TEST_CASE("Ordinary pressure outlet requires data for backflow", "[.][Ordinary-backflow]") { + IdealOutlet(2, -1.4, .95, false); +} + +TEST_CASE("Pressure outlet reconstruction is thread local", "[Pressure-outlet]") { + IdealOutlet(3, .3, .95, true, true); +} + +TEST_CASE("Ordinary outlet projection keeps SST energy in the characteristic basis", "[Pressure-outlet]") { + for (auto dim : {2u, 3u}) + for (auto mach : {-.3, .3, 1.4}) IdealOutlet(dim, mach, .95, true, false, 1500); +} diff --git a/UnitTests/SU2_CFD/solvers/pressure_outlet_test_utils.hpp b/UnitTests/SU2_CFD/solvers/pressure_outlet_test_utils.hpp new file mode 100644 index 000000000000..251455a74886 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/pressure_outlet_test_utils.hpp @@ -0,0 +1,67 @@ +/*! + * \file pressure_outlet_test_utils.hpp + * \brief Independent wave increments for pressure-outlet matrix tests. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../../../SU2_CFD/include/fluid/CFluidModel.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +// Independent primitive-variable wave decomposition for checking the implicit outlet projection. +inline void OutgoingOutletIncrement(CFluidModel& gas, unsigned short dim, const su2double* state, + const su2double* normal, su2double gridNormal, su2double turbulentK, + unsigned short column, su2double* increment) { + const auto rho = state[0]; + su2double velocity[3] = {}, dv[3] = {}; + const su2double drho = column == 0 ? 1.0 : 0.0; + for (auto d = 0u; d < dim; ++d) { + velocity[d] = state[d + 1] / rho; + dv[d] = ((column == d + 1 ? 1.0 : 0.0) - velocity[d] * drho) / rho; + } + const auto energy = state[dim + 1] / rho; + gas.SetTDState_rhoe(rho, energy - .5 * GeometryToolbox::SquaredNorm(dim, velocity) - turbulentK); + const auto sound = gas.GetSoundSpeed(), dpde = gas.GetdPde_rho(), dpdrho = gas.GetdPdrho_e(); + const su2double de = + ((column == dim + 1 ? 1.0 : 0.0) - energy * drho) / rho - GeometryToolbox::DotProduct(dim, velocity, dv); + const su2double dp = dpdrho * drho + dpde * de; + const su2double q = GeometryToolbox::DotProduct(dim, velocity, normal) - gridNormal; + const su2double dq = GeometryToolbox::DotProduct(dim, dv, normal); + su2double projectedRho = 0, projectedP = 0, projectedVelocity[3] = {}; + if (q >= 0) { + projectedRho = drho - dp / (sound * sound); + for (auto d = 0u; d < dim; ++d) projectedVelocity[d] = dv[d] - dq * normal[d]; + } + for (const auto sign : {-1, 1}) { + if (q + sign * sound < 0) continue; + const su2double amplitude = .5 * (dp + sign * rho * sound * dq); + projectedRho += amplitude / (sound * sound); + projectedP += amplitude; + for (auto d = 0u; d < dim; ++d) projectedVelocity[d] += sign * amplitude / (rho * sound) * normal[d]; + } + increment[0] = projectedRho; + for (auto d = 0u; d < dim; ++d) increment[d + 1] = velocity[d] * projectedRho + rho * projectedVelocity[d]; + increment[dim + 1] = energy * projectedRho + rho * (projectedP - dpdrho * projectedRho) / dpde + + rho * GeometryToolbox::DotProduct(dim, velocity, projectedVelocity); +} diff --git a/UnitTests/SU2_CFD/solvers/turbo_pressure_outlet_tests.cpp b/UnitTests/SU2_CFD/solvers/turbo_pressure_outlet_tests.cpp new file mode 100644 index 000000000000..5ae36c5d4552 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/turbo_pressure_outlet_tests.cpp @@ -0,0 +1,287 @@ +/*! + * \file turbo_pressure_outlet_tests.cpp + * \brief Turbomachinery pressure-outlet reconstruction and ALE regression. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include +#include "catch.hpp" +#include "pressure_outlet_test_utils.hpp" +#include "../../UnitQuadTestCase.hpp" +#include "../../../SU2_CFD/include/variables/CNSVariable.hpp" +#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp" +#include "../../../SU2_CFD/include/fluid/CVanDerWaalsGas.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +#include "../../../SU2_CFD/include/fluid/CFluidModel.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +namespace { +// Supply a known spanwise pressure profile to test how the boundary consumes it, independently +// of the radial-equilibrium profile update elsewhere in the solver. +class TurboOutletSolver final : public CEulerSolver { + public: + TurboOutletSolver(CGeometry* geometry, CConfig* config, double pressure, bool radial) + : CEulerSolver(geometry, config, MESH_0, config->GetViscous()) { + const auto spans = config->GetnSpanWiseSections(); + RadialEquilibriumPressure.resize(nMarker, spans); + ExtAverageTurboVelocity.resize(nDim); + for (auto marker = 0u; marker < nMarker; ++marker) + for (auto span = 0u; span < spans; ++span) + RadialEquilibriumPressure[marker][span] = pressure * (radial ? .98 + .02 * span : 1.); + } +}; + +// Exercise the viscous boundary path with zero diffusion so the checks measure only convective fluxes. +class ZeroOutletDiffusion final : public CNumerics { + su2double residual[5] = {}, matrix[5][5] = {}; + const su2double* rows[5]; + + public: + ZeroOutletDiffusion(unsigned short dim, const CConfig* config) : CNumerics(dim, dim + 2, config) { + for (auto row = 0u; row < nVar; ++row) rows[row] = matrix[row]; + } + ResidualType<> ComputeResidual(const CConfig*) override { return {residual, rows, rows}; } +}; + +// Supply fixed SST k to verify thermal-energy bookkeeping without advancing a turbulence model. +class FrozenOutletTurbulence final : public CSolver { + CVariable variables; + CVariable* GetBaseClassPointerToNodes() override { return &variables; } + + public: + FrozenOutletTurbulence(unsigned long points, const CConfig* config, double k) : variables(points, 2, config) { + SetBaseClassPointerToNodes(); + for (auto point = 0ul; point < points; ++point) variables.SetSolution(point, 0, k); + } +}; + +void TurboPressureOutlet(unsigned short dim, unsigned short kind, double mach, double ratio, bool moving, bool radial, + bool implicit = true, bool totalConditions = false, bool threaded = false, + double turbulentK = 0) { + CAPTURE(dim, kind, mach, ratio, moving, radial, implicit, totalConditions, threaded, turbulentK); + // Axial turbo normals point along x in the 2D fixture and z in the 3D fixture. + const unsigned short normalDim = dim == 2 ? 0 : 2; + const std::string inlet = dim == 2 ? "x_minus" : "z_minus", outlet = dim == 2 ? "x_plus" : "z_plus"; + const double p = 1500112, rho = 230; + UnitQuadTestCase test; + std::ostringstream options; + options + << std::setprecision(17) << "SOLVER= " << (turbulentK ? "RANS" : "EULER") + << "\nMESH_FORMAT= " << (dim == 2 ? "RECTANGLE" : "BOX") + << "\nMESH_BOX_SIZE= " << (threaded ? "5,129,5" : "5,5,5") + << "\nMESH_BOX_OFFSET= 1,0,0\nREF_DIMENSIONALIZATION= DIMENSIONAL\n" + << "MACH_NUMBER= 0.3\nFREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 230\n" + << "FREESTREAM_PRESSURE= 1500112\nFLUID_MODEL= " + << (kind == 0 ? "IDEAL_GAS" + : kind == 1 ? "PR_GAS" + : "VW_GAS") + << "\nGAMMA_VALUE= " << (kind == 0 ? "1.4" : "1.0125") + << "\nGAS_CONSTANT= 35.152\nCRITICAL_PRESSURE= 1415200\nCRITICAL_TEMPERATURE= 564.1\nACENTRIC_FACTOR= 0.529\n" + << "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= " << (implicit ? "EULER_IMPLICIT" : "RUNGE-KUTTA_EXPLICIT") + << "\nMARKER_EULER= " << (dim == 2 ? "(y_minus,y_plus)" : "(x_minus,x_plus)") << "\nMARKER_MONITORING= (" << inlet + << ")\nMARKER_TURBOMACHINERY= (" << inlet << ',' << outlet << ")\n" + << "TURBOMACHINERY_KIND= AXIAL\nTURBO_PERF_KIND= TURBINE\nSPANWISE_KIND= EQUISPACED\nNUM_SPANWISE_SECTIONS= 3\n" + << "MARKER_RIEMANN= (" << inlet << ",TOTAL_CONDITIONS_PT,1600000,600,1,0,0," << outlet << ',' + << (totalConditions ? "TOTAL_CONDITIONS_PT" + : radial ? "RADIAL_EQUILIBRIUM" + : "STATIC_PRESSURE") + << ',' << ratio * p << (totalConditions ? ",600,1,0,0)\n" : ",0,0,0,0)\n"); + if (moving) options << "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 20,-7,30\n"; + if (dim == 3) options << "MARKER_PERIODIC= (y_minus,y_plus,0,0,0,0,0,0,0,1,0)\n"; + if (turbulentK) + options << "KIND_TURB_MODEL= SST\nINIT_OPTION= TD_CONDITIONS\nREYNOLDS_NUMBER= 1E6\n" + "VISCOSITY_MODEL= CONSTANT_VISCOSITY\nCONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY\n"; + test.config_options = options.str(); + test.InitConfig(); + test.InitGeometry(); + auto* config = test.config.get(); + std::ostringstream quiet; + std::cout.rdbuf(quiet.rdbuf()); + if (dim == 3) { + test.geometry->MatchPeriodic(config, 1); + test.geometry->PreprocessPeriodicComms(test.geometry.get(), config); + } + for (auto flag : {INFLOW, OUTFLOW}) test.geometry->ComputeNSpan(config, 0, flag, true); + for (auto flag : {INFLOW, OUTFLOW}) test.geometry->SetTurboVertex(config, 0, flag, true); + config->SetnSpanMaxAllZones(config->GetnSpanWiseSections()); + test.solver = new CSolver*[MAX_SOLS](); + auto* solver = new TurboOutletSolver(test.geometry.get(), config, ratio * p, radial); + test.solver[FLOW_SOL] = solver; + std::cout.rdbuf(test.orig_buf); + std::unique_ptr gas; + if (kind == 0) gas.reset(new CIdealGas(1.4, 35.152)); + if (kind == 1) gas.reset(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, .529)); + if (kind == 2) gas.reset(new CVanDerWaalsGas(1.0125, 35.152, 1415200, 564.1)); + if (turbulentK) { + gas->SetLaminarViscosityModel(config); + gas->SetThermalConductivityModel(config); + } + gas->SetTDState_Prho(p, rho); + const su2double sound = gas->GetSoundSpeed(), e = gas->GetStaticEnergy(); + const su2double grid[3] = {moving ? 20. : 0., moving ? -7. : 0., moving ? 30. : 0.}; + su2double state[5] = {rho}, velocity[3] = {3, -2, 1}; + for (auto d = 0u; d < dim; ++d) { + velocity[d] = (d == normalDim ? mach * sound : velocity[d]) + grid[d]; + state[d + 1] = rho * velocity[d]; + } + state[dim + 1] = rho * (e + .5 * GeometryToolbox::SquaredNorm(dim, velocity) + turbulentK); + std::unique_ptr turbulence; + if (turbulentK) { + turbulence.reset(new FrozenOutletTurbulence(test.geometry->GetnPoint(), config, turbulentK)); + test.solver[TURB_SOL] = turbulence.get(); + } + auto* nodes = static_cast(solver->GetNodes()); + for (auto point = 0ul; point < test.geometry->GetnPoint(); ++point) { + nodes->SetSolution(point, state); + if (turbulentK) { + auto* viscousNodes = static_cast(nodes); + REQUIRE(viscousNodes->SetPrimVar(point, 0, turbulentK, gas.get())); + viscousNodes->SetSecondaryVar(point, gas.get()); + } else { + REQUIRE(nodes->SetPrimVar(point, gas.get())); + nodes->SetSecondaryVar(point, gas.get()); + } + if (moving) test.geometry->nodes->SetGridVel(point, grid); + for (auto k = 0u; k < dim + 2; ++k) solver->LinSysRes(point, k) = 0; + su2double zero[25] = {}; + solver->Jacobian.SetBlock(point, point, zero); + } + for (auto marker = 0u; marker < config->GetnMarker_All(); ++marker) { + if (config->GetMarker_All_TagBound(marker) != outlet) continue; + const auto spans = test.geometry->GetnSpanWiseSections(OUTFLOW); + auto target = [&](unsigned short span) { return ratio * p * (radial ? .98 + .02 * span : 1.); }; + SU2_OMP_PARALLEL { + CUpwGeneralHLLC_Flow numerics(dim, dim + 2, config); + ZeroOutletDiffusion diffusion(dim, config); + solver->BC_TurboRiemann(test.geometry.get(), test.solver, &numerics, &diffusion, config, marker); + } + CHECK(solver->GetFluidModel()->GetPressure() > 0); + CUpwGeneralHLLC_Flow numerics(dim, dim + 2, config); + for (auto span = 0u; span < spans; ++span) + for (auto vertex = 0ul; vertex < test.geometry->GetnVertexSpan(marker, span); ++vertex) { + const auto* turbo = test.geometry->turbovertex[marker][span][vertex]; + const auto point = turbo->GetNode(), oldVertex = turbo->GetOldVertex(); + su2double normal[3], unit[3]; + test.geometry->vertex[marker][oldVertex]->GetNormal(normal); + turbo->GetNormal(unit); + for (auto d = 0u; d < dim; ++d) normal[d] = -normal[d]; + su2double boundary[5]; + REQUIRE(gas->ComputePressureOutletState(dim, state, unit, GeometryToolbox::DotProduct(dim, grid, unit), + target(span), boundary, turbulentK, + kind == 0 ? 1.4 : 0) == CFluidModel::PressureOutletStatus::Valid); + auto flux = [&](const su2double* boundary, su2double* result) { + su2double v[3]; + for (auto d = 0u; d < dim; ++d) v[d] = boundary[d + 1] / boundary[0]; + const su2double energy = boundary[dim + 1] / boundary[0]; + gas->SetTDState_rhoe(boundary[0], energy - .5 * GeometryToolbox::SquaredNorm(dim, v) - turbulentK); + const su2double pressure = gas->GetPressure(), h = energy + pressure / boundary[0]; + numerics.GetInviscidProjFlux(&boundary[0], v, &pressure, &h, normal, result); + // The moving-grid transport term occurs once, in both explicit and implicit boundary paths. + const su2double speed = GeometryToolbox::DotProduct(dim, grid, normal); + for (auto k = 0u; k < dim + 2; ++k) result[k] -= speed * boundary[k]; + }; + su2double expected[5]; + flux(boundary, expected); + for (auto row = 0u; row < dim + 2; ++row) { + CHECK(std::isfinite(SU2_TYPE::GetValue(solver->LinSysRes(point, row)))); + CHECK(SU2_TYPE::GetValue(solver->LinSysRes(point, row)) == + Approx(SU2_TYPE::GetValue(expected[row])).epsilon(1e-8).margin(1e-7)); + } + const auto* primitive = solver->GetCharacPrimVar(marker, oldVertex); + CHECK(SU2_TYPE::GetValue(primitive[0]) == Approx(SU2_TYPE::GetValue(gas->GetTemperature())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[dim + 1]) == Approx(SU2_TYPE::GetValue(gas->GetPressure())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[dim + 2]) == Approx(SU2_TYPE::GetValue(gas->GetDensity())).epsilon(1e-8)); + CHECK(SU2_TYPE::GetValue(primitive[dim + 4]) == Approx(SU2_TYPE::GetValue(gas->GetSoundSpeed())).epsilon(1e-8)); + if (mach == .3 && ratio > .1) + CHECK(SU2_TYPE::GetValue(primitive[dim + 1]) == Approx(target(span)).epsilon(1e-8)); + if (!implicit) continue; + // The implicit matrix uses a frozen characteristic projection, with local span pressure and SST k. + for (auto column = 0u; column < dim + 2; ++column) { + su2double plus[5], minus[5], fp[5], fm[5], increment[5]; + OutgoingOutletIncrement(*gas, dim, state, unit, GeometryToolbox::DotProduct(dim, grid, unit), turbulentK, + column, increment); + const double delta = + 2e-5 * fmax(fabs(SU2_TYPE::GetValue(state[column])), column == 0 ? rho : rho * SU2_TYPE::GetValue(sound)); + for (auto k = 0u; k < dim + 2; ++k) { + plus[k] = boundary[k] + delta * increment[k]; + minus[k] = boundary[k] - delta * increment[k]; + } + flux(plus, fp); + flux(minus, fm); + for (auto row = 0u; row < dim + 2; ++row) + CHECK(solver->Jacobian.GetBlock(point, point)[row * (dim + 2) + column] == + Approx(SU2_TYPE::GetValue((fp[row] - fm[row]) / (2 * delta))).epsilon(5e-4).margin(2e-5)); + } + } + } +} +} // namespace + +TEST_CASE("TurboRiemann pressure outlets preserve admissibility and the ALE flux", "[Pressure-outlet][Turbo-outlet]") { + for (auto dim : {2u, 3u}) + for (auto kind : {0u, 1u, 2u}) + for (const auto radial : {false, true}) + for (const auto moving : {false, true}) + for (const auto mach : {-.3, 0., .3, .99, 1.4}) + for (const auto ratio : {.95, .0006666168959384}) + TurboPressureOutlet(dim, kind, mach, ratio, moving, radial); +} + +TEST_CASE("Explicit TurboRiemann pressure outlet uses one grid flux correction", "[Pressure-outlet][Turbo-outlet]") { + TurboPressureOutlet(2, 0, 1.4, .95, true, false, false); +} + +TEST_CASE("TurboRiemann real-gas startup pressure remains admissible", "[Pressure-outlet][Turbo-startup]") { + // This ratio prescribes 1000 Pa from an interior pressure of 1500112 Pa. + TurboPressureOutlet(2, 1, .3, .0006666168959384, false, false); +} + +TEST_CASE("Implicit TurboRiemann pressure outlet uses one grid flux correction", "[Pressure-outlet][Turbo-moving]") { + TurboPressureOutlet(2, 0, 1.4, .95, true, false); +} + +TEST_CASE("TurboRiemann pressure reconstruction is thread local", "[Pressure-outlet][Turbo-outlet]") { + // Use enough outlet vertices to enter the OpenMP work-sharing loop with separate EOS instances. + TurboPressureOutlet(2, 1, .3, .95, true, false, true, false, true); +} + +TEST_CASE("Other TurboRiemann modes retain one ALE correction", "[Pressure-outlet][Turbo-outlet]") { + for (const auto implicit : {false, true}) TurboPressureOutlet(2, 0, 1.4, .95, true, false, implicit, true); +} + +TEST_CASE("TurboRiemann pressure outlet rejects nonpositive pressure", "[.][Turbo-invalid-pressure]") { + TurboPressureOutlet(2, 1, .3, 0, false, false); +} + +TEST_CASE("TurboRiemann radial pressure outlet requires data for backflow", "[.][Turbo-backflow]") { + TurboPressureOutlet(3, 1, -1.4, .95, true, true); +} + +TEST_CASE("TurboRiemann real-gas radial startup remains admissible", "[Pressure-outlet][Turbo-radial-startup]") { + TurboPressureOutlet(3, 1, .3, .0006666168959384, false, true); +} + +TEST_CASE("TurboRiemann pressure outlet separates frozen turbulence energy from the EOS", + "[Pressure-outlet][Turbo-outlet]") { + for (auto kind : {0u, 1u}) TurboPressureOutlet(3, kind, .3, .95, true, true, true, false, false, 13); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 1d1a97f3b827..b703dce88ccf 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -21,6 +21,9 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/periodic_limiters.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', + 'SU2_CFD/solvers/characteristic_pressure_outlet_tests.cpp', + 'SU2_CFD/solvers/pressure_outlet_boundary_tests.cpp', + 'SU2_CFD/solvers/turbo_pressure_outlet_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', @@ -29,11 +32,13 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', # Reverse-mode (algorithmic differentiation) tests: su2_cfd_tests_ad = files(['Common/simple_ad_test.cpp', + 'SU2_CFD/solvers/pressure_outlet_ad_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests_AD.cpp', 'Common/toolboxes/multilayer_perceptron/MLP_Jacobian_tests.cpp']) # Forward-mode (direct differentiation) tests: -su2_cfd_tests_dd = files(['Common/simple_directdiff_test.cpp']) +su2_cfd_tests_dd = files(['Common/simple_directdiff_test.cpp', + 'SU2_CFD/solvers/pressure_outlet_ad_tests.cpp']) # ------------------------------------------------------------------------- # End of unit test listings @@ -51,6 +56,10 @@ if get_option('enable-tests') cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) + # Fatal boundary cases need separate processes; they cannot run inside the normal Catch suite. + test('Characteristic boundary admissibility', python, + args: files('SU2_CFD/solvers/check_characteristic_boundary_admissibility.py') + [test_driver.full_path()], + timeout: 180) test('NEMO convergence field migration', python, args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()], timeout: 300) diff --git a/config_template.cfg b/config_template.cfg index 5e611f22ed48..2ba054ae22b8 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1214,6 +1214,13 @@ MARKER_INLET= ( NONE ) % % Outlet boundary marker(s) (NONE = no marker) % Compressible: ( outlet marker, back pressure (static thermodynamic), ... ) +% Ideal-gas finite-volume outlets use normal velocity relative to the grid. +% Tangential speed does not determine whether pressure can affect the interior. +% Subsonic outflow preserves interior entropy and tangential velocity while imposing pressure. +% If the requested pressure would accelerate flow past sonic, use the sonic state instead; +% its pressure can be higher than the configured back pressure. +% Already-supersonic outward flow is extrapolated. Subsonic reverse flow uses interior +% entropy and tangential velocity as inflow assumptions; fully incoming flow requires inlet data. % Inc. Pressure: ( outlet marker, back pressure (static gauge in Pa), ... ) % Inc. Mass Flow: ( outlet marker, mass flow target (kg/s), ... ) MARKER_OUTLET= ( NONE ) @@ -1301,6 +1308,11 @@ MARKER_PRESSURE= ( NONE ) % % Riemann boundary marker(s) (NONE = no marker) % Format: (marker, data kind flag, list of data) +% For finite-volume STATIC_PRESSURE outlets (including turbomachinery) and +% TurboRiemann RADIAL_EQUILIBRIUM outlets, the boundary state follows the EOS +% isentrope (constant entropy) for ideal/real gases, using normal speed relative to the grid. +% Apply the same sonic limit and reverse-flow assumptions described for MARKER_OUTLET above. +% RADIAL_EQUILIBRIUM supplies a pressure for each span; reconstruction uses the local flow state. MARKER_RIEMANN= ( NONE ) % % Shroud boundary marker(s) (NONE = no marker)