diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 550729877fba..81fb736ce20d 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4607,6 +4607,36 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Rotating_Frame = (Kind_GridMovement == ROTATING_FRAME); + /*--- The finite volume flow solvers apply these source terms with a single numerics object, which is built + for only one of them: combinations would apply one source several times and ignore the others. Only these + solvers (and their discrete adjoints) are checked: e.g. elasticity uses BODY_FORCE and GRAVITY_FORCE together. ---*/ + + if (GetFluidProblem() && !GetNEMOProblem()) { + const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE); + const int nSources = int(Body_Force) + int(Rotating_Frame) + int(Axisymmetric) + int(GravityForce && !inc) + + int(inc && Kind_DensityModel == INC_DENSITYMODEL::BOUSSINESQ) + + int(inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE); + if (nSources > 1) { + SU2_MPI::Error("Only one of BODY_FORCE, GRID_MOVEMENT= ROTATING_FRAME, AXISYMMETRIC, GRAVITY_FORCE (compressible), " + "INC_DENSITY_MODEL= BOUSSINESQ and KIND_STREAMWISE_PERIODIC (incompressible) can be used at a time.", + CURRENT_FUNCTION); + } + } + + /*--- The second source object, for radiation, the energy correction of incompressible streamwise periodic flow + (without periodic temperature) or vorticity confinement, is also built for only one of them. ---*/ + + if (GetFluidProblem() && !GetNEMOProblem()) { + const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE); + const bool periodicOutlet = inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE && Energy_Equation && + !Streamwise_Periodic_Temperature; + if (int(Radiation) + int(periodicOutlet) + int(VorticityConfinement) > 1) { + SU2_MPI::Error("Only one of RADIATION, the energy correction of KIND_STREAMWISE_PERIODIC (incompressible, " + "STREAMWISE_PERIODIC_TEMPERATURE= NO) and VORTICITY_CONFINEMENT can be used at a time.", + CURRENT_FUNCTION); + } + } + /*--- In case the grid movement parameters have not been declared in the config file, set them equal to zero for safety. Also check to make sure that for each option, a value has been declared for each moving marker. ---*/ @@ -4966,6 +4996,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Restart || ((Kind_MGCycle == MG_CYCLE::FULL) && ContinuousAdjoint)) Kind_MGCycle = MG_CYCLE::V; + /*--- The multigrid arrays hold at most 10 coarse levels (CMultiGridIntegration::MAX_MG_LEVELS). ---*/ + if (nMGLevels > 10) { + SU2_MPI::Error("MGLEVEL must be between 0 and 10; larger multigrid hierarchies are not supported.", CURRENT_FUNCTION); + } + FinestMesh = MESH_0; if (Kind_MGCycle == MG_CYCLE::FULL) FinestMesh = nMGLevels; diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp index f60c0422f0cc..89eb028c6aeb 100644 --- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp +++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp @@ -69,7 +69,14 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER const auto chunkSize = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_CHUNK); #endif - static constexpr size_t MAXNVAR = 20; + static constexpr size_t MAXNVAR = 32; + + /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/ + if (varEnd > MAXNVAR) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS + SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION); + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } /*--- For each (non-halo) volume integrate over its faces (edges). ---*/ diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp index 7d8beafb03e0..0cdc60207a27 100644 --- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -80,8 +80,12 @@ void computeLimiters_impl(CSolver* solver, { constexpr size_t MAXNVAR = 32; - if (varEnd > MAXNVAR) + /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/ + if (varEnd > MAXNVAR) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION); + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } const size_t nPointDomain = geometry.GetnPointDomain(); const size_t nPoint = geometry.GetnPoint(); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..650b26b91e7f 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1155,6 +1155,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS); const auto iVel = prim_idx.Velocity(); + const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2; /*--- Blazek chapter 8.: * The components of the momentum residual normal to the symmetry plane are zeroed out. @@ -1184,8 +1185,8 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area; } - /*--- Energy terms due to grid movement (aka work of pressure forces). ---*/ - if (dynamic_grid) { + /*--- The classical compressible projection below already supplies wall pressure work. ---*/ + if (dynamic_grid && !compressibleProjection) { su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex); conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -1217,8 +1218,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve /*--- Compute the residual using an upwind scheme. ---*/ auto residual = conv_numerics->ComputeResidual(config); - /*--- Use just the energy fluxes to update the residual, adding the others would - * increase numerical diffusion which we wish to avoid if possible. ---*/ + /*--- Energy fluxes for the other conservative layouts. ---*/ for (auto iVar = iVel + nDim; iVar < nVar; iVar++) { LinSysRes(iPoint, iVar) += residual.residual[iVar]; } @@ -1243,14 +1243,16 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve su2double* solutionOld = nodes->GetSolution_Old(iPoint); su2double gridVel[MAXNDIM] = {}; + su2double gridNormalVel = 0.0; if (dynamic_grid) { for (auto iDim = 0u; iDim < nDim; iDim++) { gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim]; } + gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal); if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) { for(auto iDim = 0u; iDim < nDim; iDim++) { /*--- Multiply by density since we are correcting conservative variables. ---*/ - gridVel[iDim] *= nodes->GetDensity(iPoint); + gridVel[iDim] *= compressibleProjection ? solutionOld[0] : nodes->GetDensity(iPoint); } } } @@ -1261,20 +1263,27 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) { solutionOld[iVel + iDim] -= vp * UnitNormal[iDim]; } + if (compressibleProjection) { + solutionOld[iVel + nDim] -= gridNormalVel * vp; + } /*--- Keep only the tangential part of the momentum residuals. ---*/ su2double normalRes = 0.0; for (auto iDim = 0u; iDim < nDim; iDim++) { normalRes += LinSysRes(iPoint, iVel + iDim) * UnitNormal[iDim]; } + if (compressibleProjection) { + /*--- Project in the wall frame: normal momentum follows density, with consistent pressure work. ---*/ + normalRes -= gridNormalVel * LinSysRes(iPoint, 0); + LinSysRes(iPoint, iVel + nDim) -= gridNormalVel * normalRes; + } for (auto iDim = 0u; iDim < nDim; iDim++) { LinSysRes(iPoint, iVel + iDim) -= normalRes * UnitNormal[iDim]; } /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) { - /*--- Modify the Jacobians according to the modification of the residual - * J_new = (I - n * n^T) * J where n = {0, nx, ny, nz, 0, ...} ---*/ + /*--- Apply the same wall-frame projection to the Jacobians. ---*/ su2double mat[MAXNVAR * MAXNVAR] = {}; for (auto iVar = 0u; iVar < nVar; iVar++) @@ -1282,6 +1291,13 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) for (auto jDim = 0u; jDim < nDim; jDim++) mat[(iDim + iVel) * nVar + jDim + iVel] -= UnitNormal[iDim] * UnitNormal[jDim]; + if (compressibleProjection) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + mat[(iDim + iVel) * nVar] += gridNormalVel * UnitNormal[iDim]; + mat[(iVel + nDim) * nVar + iVel + iDim] -= gridNormalVel * UnitNormal[iDim]; + } + mat[(iVel + nDim) * nVar] += gridNormalVel * gridNormalVel; + } auto ModifyJacobian = [&](const unsigned long jPoint) { su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR]; @@ -1306,6 +1322,10 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolve for (auto iDim = 0u; iDim < nDim; iDim++) { normalRes += Res_TruncError[iVel + iDim] * UnitNormal[iDim]; } + if (compressibleProjection) { + normalRes -= gridNormalVel * Res_TruncError[0]; + Res_TruncError[iVel + nDim] -= gridNormalVel * normalRes; + } for (auto iDim = 0u; iDim < nDim; iDim++) { Res_TruncError[iVel + iDim] -= normalRes * UnitNormal[iDim]; } @@ -2974,6 +2994,7 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo bool use_solution_old) { const auto iVel = prim_idx.Velocity(); const auto nDim = geometry->GetnDim(); + const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) != EULER_WALL) continue; @@ -3004,12 +3025,15 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo * For dynamic grids subtract the grid velocity to enforce (v - v_grid).n = 0, * multiplying by density for compressible flow (conservative variables). ---*/ su2double gridVel[MAXNDIM] = {}; - if (dynamic_grid && !use_solution_old) { + su2double gridNormalVel = 0.0; + if (dynamic_grid && (compressibleProjection || !use_solution_old)) { for (auto iDim = 0u; iDim < nDim; iDim++) gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim]; + gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal); if constexpr (FlowRegime == ENUM_REGIME::COMPRESSIBLE) { - for (auto iDim = 0u; iDim < nDim; iDim++) - gridVel[iDim] *= nodes->GetDensity(iPoint); + /*--- Corrections carry their own density increment, rather than the current solution density. ---*/ + const auto density = compressibleProjection ? sol[0] : nodes->GetDensity(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) gridVel[iDim] *= density; } } @@ -3019,6 +3043,7 @@ void CFVMFlowSolverBase::MultigridProjectEulerWall(CGeometry* geo /*--- Project to tangent plane. ---*/ for (auto iDim = 0u; iDim < nDim; iDim++) sol[iVel + iDim] -= momentum_n * UnitNormal[iDim]; + if (compressibleProjection) sol[iVel + nDim] -= gridNormalVel * momentum_n; } END_SU2_OMP_FOR } diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index b0130279f76f..713ca2e86662 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -31,7 +31,7 @@ CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - variable_density = ((config->GetKind_DensityModel() == INC_DENSITYMODEL::VARIABLE) || (config->GetKind_DensityModel() == INC_DENSITYMODEL::FLAMELET)); + variable_density = config->GetVariable_Density_Model(); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); fix_factor = config->GetCent_Inc_Jac_Fix_Factor(); @@ -161,17 +161,25 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim+1; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, + su2double cp) { + const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j); } } @@ -250,7 +258,7 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi CCentJSTInc_Flow::CCentJSTInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - variable_density = (config->GetVariable_Density_Model()); + variable_density = config->GetVariable_Density_Model(); energy = config->GetEnergy_Equation(); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); @@ -401,17 +409,25 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim+1; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T, + su2double cp) { + const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0); + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j); } } diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp index e806127da203..57c7ed9f6f5a 100644 --- a/SU2_CFD/src/numerics/flow/convection/fds.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp @@ -249,19 +249,23 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; /*--- Residual contributions ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*--- Jacobian contributions ---*/ - /*--- Implicit terms ---*/ - if (implicit) { - for (iDim = 0; iDim < nDim; iDim++){ - Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added + once. The density depends on h (variable density models) and not on the pressure. ---*/ + if (implicit) { + auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double dRhodh) { + const su2double factor = 0.5*ProjVelocity; + Jac[0][nDim+1] -= factor*dRhodh; + for (iDim = 0; iDim < nDim; iDim++) { + Jac[iDim+1][iDim+1] -= factor*rho; + Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh; } - Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j; - } + Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh); + }; + addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, dRhodh_i); + addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, dRhodh_j); } } diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index e3797c7c5a5f..87784e73f39f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -138,7 +138,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[3][0] = -Gamma*U_i[2]*U_i[3]/(U_i[0]*U_i[0]) + (Gamma-1)*U_i[2]*(U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]*U_i[0]); jacobian[3][1] = -(Gamma-1)*U_i[2]*U_i[1]/(U_i[0]*U_i[0]); - jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 1/2*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); + jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 0.5*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); jacobian[3][3] = Gamma*U_i[2]/U_i[0]; for (iVar=0; iVar < nVar; iVar++) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index b1c61ad377de..c05a097c4a2f 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1947,27 +1947,39 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } - /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ + /*--- Check for non-physical solutions after reconstruction. If found, use the + cell-average value of the solution. This is a locally 1st order approximation, + which is typically only active during the start-up of a calculation. + Pressure and density are checked before the fluid model is evaluated with them. ---*/ - if (!ideal_gas || low_mach_corr) { - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); - } + auto nonPositive = [](su2double value) { return !(value > 0.0) || !std::isfinite(SU2_TYPE::GetValue(value)); }; + bool neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]); + bool neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]); - /*--- Low-Mach number correction. ---*/ + if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) { - if (low_mach_corr) { - LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); - } + /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ - /*--- Check for non-physical solutions after reconstruction. If found, use the - cell-average value of the solution. This is a locally 1st order approximation, - which is typically only active during the start-up of a calculation. ---*/ + if (!ideal_gas || low_mach_corr) { + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); + } + + /*--- Low-Mach number correction. ---*/ + + if (low_mach_corr) { + LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); + } + + /*--- The correction changes the face states, check them again. ---*/ - bool neg_pres_or_rho_i = (Primitive_i[prim_idx.Pressure()] < 0.0) || (Primitive_i[prim_idx.Density()] < 0.0); - bool neg_pres_or_rho_j = (Primitive_j[prim_idx.Pressure()] < 0.0) || (Primitive_j[prim_idx.Density()] < 0.0); + neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]); + neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]); + } - su2double R = sqrt(fabs(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()])); + su2double R = 1.0; + if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) + R = sqrt(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()]); su2double sq_vel = 0.0; for (auto iDim = 0u; iDim < nDim; iDim++) { su2double RoeVelocity = (R * Primitive_j[iDim + prim_idx.Velocity()] + @@ -5801,7 +5813,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; su2double donorAverages[5] = {0.0}; - switch (config->GetKind_Data_Giles(Marker_Tag)){ + switch (config->GetKind_Data_Riemann(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; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 512a99a97dc6..8c6ee1719a01 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2235,7 +2235,7 @@ void CIncEulerSolver::SetRangePressure(CGeometry *geometry, CSolver **solver_con BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { minP = MinP; - SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); maxP = MaxP; SU2_MPI::Allreduce(&maxP, &MaxP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index a499a075fa1d..2dfe8b65cf70 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -88,6 +88,20 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") { CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT); } +TEST_CASE("Supported multigrid levels are retained", "[Config][MultigridLimits]") { + for (const auto levels : {0, 10}) { + std::stringstream options(base_options + "MGLEVEL= " + std::to_string(levels) + "\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + CHECK(config.GetnMGLevels() == levels); + } +} + +TEST_CASE("Excessive multigrid levels are rejected", "[.MultigridLimitError]") { + std::stringstream options(base_options + "MGLEVEL= 11\n"); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + FAIL("An unsupported multigrid hierarchy was accepted."); +} + TEST_CASE("HLLC and AUSM SLAU Jacobian options are independent", "[Config][HLLCOptions]") { for (const auto accurate : {false, true}) { for (const auto exactHLLC : {false, true}) { diff --git a/UnitTests/Common/check_multigrid_limits.py b/UnitTests/Common/check_multigrid_limits.py new file mode 100644 index 000000000000..5070fdaed3cc --- /dev/null +++ b/UnitTests/Common/check_multigrid_limits.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python3 +"""Check the fatal configuration diagnostic in its own process.""" +import subprocess +import sys + +result = subprocess.run( + [sys.argv[1], "[.MultigridLimitError]"], capture_output=True, text=True, timeout=30 +) +expected = "MGLEVEL must be between 0 and 10" +if result.returncode == 0 or expected not in result.stdout + result.stderr: + raise RuntimeError( + "Unsupported MGLEVEL did not fail with its diagnostic:\n" + + result.stdout + + result.stderr + ) +print("PASS: MGLEVEL=11 is rejected with its supported range") diff --git a/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp new file mode 100644 index 000000000000..1baf0c271428 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp @@ -0,0 +1,182 @@ +/*! + * \file moving_wall_tests.cpp + * \brief Regression for compressible wall projection in fixed and uniformly translating frames. + * \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 +#include +#include "../../UnitQuadTestCase.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp" +#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp" + +TEST_CASE("Compressible wall projection commutes with a Galilean boost", "[HLLC-wall]") { + for (const auto nDim : {2u, 3u}) { + for (const bool realGas : {false, true}) { + UnitQuadTestCase fixed, moving; + const std::string markers = + nDim == 2 ? "(x_minus,x_plus,y_minus,y_plus)" : "(x_minus,x_plus,y_minus,y_plus,z_minus,z_plus)"; + const std::string options = "SOLVER= EULER\nMESH_FORMAT= " + std::string(nDim == 2 ? "RECTANGLE" : "BOX") + + "\nMESH_BOX_SIZE= 5,5,5\nREF_DIMENSIONALIZATION= DIMENSIONAL\n" + "MACH_NUMBER= 0.1\nFREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 230\n" + "FREESTREAM_PRESSURE= 1500112\nGAS_CONSTANT= 35.152\n" + "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n" + "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 0,0,0\nMARKER_EULER= " + + markers + "\nMARKER_MONITORING= (x_minus)\n" + + (realGas ? "FLUID_MODEL= PR_GAS\nGAMMA_VALUE= 1.0125\nCRITICAL_PRESSURE= " + "1415200\nCRITICAL_TEMPERATURE= 564.1\nACENTRIC_FACTOR= 0.529\n" + : "FLUID_MODEL= IDEAL_GAS\nGAMMA_VALUE= 1.4\n"); + fixed.config_options = moving.config_options = options; + fixed.InitConfig(); + fixed.InitGeometry(); + fixed.InitSolver(); + moving.InitConfig(); + moving.InitGeometry(); + moving.InitSolver(); + const auto nVar = nDim + 2; + const su2double boost[3] = {20, -7, 3}; + const su2double velocity[3] = {11, 4, -2}; + su2double B[5][5] = {}, T[5][5] = {}, baseJac[5][5] = {}, boostedJac[5][5] = {}; + su2double c2 = 0; + for (auto i = 0u; i < nVar; ++i) { + B[i][i] = T[i][i] = 1; + for (auto j = 0u; j < nVar; ++j) baseJac[i][j] = 1 + i + 2 * j; + } + for (auto i = 0u; i < nDim; ++i) { + c2 += boost[i] * boost[i]; + B[i + 1][0] = boost[i]; + T[i + 1][0] = -boost[i]; + B[nVar - 1][i + 1] = boost[i]; + T[nVar - 1][i + 1] = -boost[i]; + } + B[nVar - 1][0] = T[nVar - 1][0] = 0.5 * c2; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) + for (auto k = 0u; k < nVar; ++k) + for (auto l = 0u; l < nVar; ++l) boostedJac[i][j] += B[i][k] * baseJac[k][l] * T[l][j]; + std::unique_ptr fluid = + realGas ? std::unique_ptr(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, 0.529)) + : std::unique_ptr(new CIdealGas(1.4, 35.152)); + fluid->SetTDState_Prho(1500112, 230); + std::array baseState{}; + baseState[0] = 230; + su2double speed2 = 0; + for (auto i = 0u; i < nDim; ++i) { + baseState[i + 1] = 230 * velocity[i]; + speed2 += velocity[i] * velocity[i]; + } + baseState[nVar - 1] = 230 * (fluid->GetStaticEnergy() + 0.5 * speed2); + const su2double baseResidual[5] = {0.3, -0.5, 0.2, 0.7, -0.1}; + for (auto frame = 0u; frame < 2; ++frame) { + auto& test = frame ? moving : fixed; + auto* solver = test.solver[FLOW_SOL]; + auto* nodes = static_cast(solver->GetNodes()); + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) { + su2double state[5] = {}, residual[5] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) { + state[i] += (frame ? B[i][j] : (i == j)) * baseState[j]; + residual[i] += (frame ? B[i][j] : (i == j)) * baseResidual[j]; + } + nodes->SetSolution(p, state); + nodes->SetSolution_Old(p, state); + REQUIRE(nodes->SetPrimVar(p, fluid.get())); + nodes->SetSecondaryVar(p, fluid.get()); + for (auto i = 0u; i < nDim; ++i) test.geometry->nodes->SetGridVel(p, i, frame ? boost[i] : 0); + for (auto i = 0u; i < nVar; ++i) { + solver->LinSysRes(p, i) = residual[i]; + nodes->GetResTruncError(p)[i] = residual[i]; + } + su2double block[25] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) block[i * nVar + j] = frame ? boostedJac[i][j] : baseJac[i][j]; + solver->Jacobian.SetBlock(p, p, block); + for (auto k = 0u; k < test.geometry->nodes->GetnPoint(p); ++k) + solver->Jacobian.SetBlock(p, test.geometry->nodes->GetPoint(p, k), block); + } + std::unique_ptr numerics = + realGas ? std::unique_ptr(new CUpwGeneralHLLC_Flow(nDim, nVar, test.config.get())) + : std::unique_ptr(new CUpwHLLC_Flow(nDim, nVar, test.config.get())); + for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker) + if (test.config->GetMarker_All_TagBound(marker) == "x_minus") + solver->BC_Sym_Plane(test.geometry.get(), test.solver, numerics.get(), nullptr, test.config.get(), marker); + } + for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) { + auto* f = fixed.solver[FLOW_SOL]; + auto* m = moving.solver[FLOW_SOL]; + for (auto i = 0u; i < nVar; ++i) { + su2double state = 0, residual = 0, truncation = 0; + for (auto j = 0u; j < nVar; ++j) { + state += B[i][j] * f->GetNodes()->GetSolution_Old(p)[j]; + residual += B[i][j] * f->LinSysRes(p, j); + truncation += B[i][j] * f->GetNodes()->GetResTruncError(p)[j]; + } + CHECK(m->GetNodes()->GetSolution_Old(p)[i] == Approx(state).margin(1e-6)); + CHECK(m->LinSysRes(p, i) == Approx(residual).margin(1e-8)); + CHECK(m->GetNodes()->GetResTruncError(p)[i] == Approx(truncation).margin(1e-8)); + for (auto block = 0u; block <= fixed.geometry->nodes->GetnPoint(p); ++block) { + const auto neighbor = block ? fixed.geometry->nodes->GetPoint(p, block - 1) : p; + for (auto j = 0u; j < nVar; ++j) { + su2double expected = 0; + for (auto k = 0u; k < nVar; ++k) + for (auto l = 0u; l < nVar; ++l) + expected += B[i][k] * f->Jacobian.GetBlock(p, neighbor)[k * nVar + l] * T[l][j]; + CHECK(m->Jacobian.GetBlock(p, neighbor)[i * nVar + j] == Approx(expected).margin(1e-5)); + } + } + } + } + + for (const bool correction : {false, true}) { + const su2double delta[5] = {0.3, -0.5, 0.2, 0.7, -0.1}; + for (auto frame = 0u; frame < 2; ++frame) { + auto& test = frame ? moving : fixed; + auto* solver = test.solver[FLOW_SOL]; + for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) { + su2double state[5] = {}; + for (auto i = 0u; i < nVar; ++i) + for (auto j = 0u; j < nVar; ++j) + state[i] += (frame ? B[i][j] : (i == j)) * (correction ? delta[j] : baseState[j]); + if (correction) + solver->GetNodes()->SetSolution_Old(p, state); + else + solver->GetNodes()->SetSolution(p, state); + } + solver->MultigridProjectEulerWall(test.geometry.get(), test.config.get(), correction); + } + for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) { + const auto* f = correction ? fixed.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p) + : fixed.solver[FLOW_SOL]->GetNodes()->GetSolution(p); + const auto* m = correction ? moving.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p) + : moving.solver[FLOW_SOL]->GetNodes()->GetSolution(p); + for (auto i = 0u; i < nVar; ++i) { + su2double expected = 0; + for (auto j = 0u; j < nVar; ++j) expected += B[i][j] * f[j]; + CHECK(m[i] == Approx(expected).margin(1e-6)); + } + } + } + } + } +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 85b96280544e..c0b2f31f2f16 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -22,6 +22,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/output/COutput_convergence_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', + 'SU2_CFD/solvers/moving_wall_tests.cpp', 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', @@ -54,6 +55,8 @@ if get_option('enable-tests') cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) + test('Multigrid configuration limits', python, + args: files('Common/check_multigrid_limits.py') + [test_driver.full_path()]) test('NEMO convergence field migration', python, args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()], timeout: 300)