From 83e080c60d2d3f628be884a320623b6319d2cc6c Mon Sep 17 00:00:00 2001 From: bigfooted Date: Tue, 16 Apr 2024 22:33:09 +0200 Subject: [PATCH 01/20] fix corner node --- SU2_CFD/src/solvers/CEulerSolver.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 16cdfe59a4b..9385c9c874c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7111,6 +7111,23 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, break; } + /*--- Check if the inlet node is shared with a viscous wall. ---*/ + + if (geometry->nodes->GetViscousBoundary(iPoint)) { + + V_inlet[0] = nodes->GetTemperature(iPoint); + + /*--- Impose the wall velocity from the interior. ---*/ + for (iDim = 0; iDim < nDim; iDim++){ + V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + } + + /*--- Match the pressure and energy at the wall. ---*/ + V_inlet[nDim+1] = nodes->GetPressure(iPoint); + V_inlet[nDim+3] = nodes->GetPressure(iPoint)/(Density*Gamma_Minus_One) + nodes->GetPressure(iPoint)/Density; // + Pressure/Density; + } + + /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_inlet); From 8624e0df58e87694951b036375444a1e53eb095d Mon Sep 17 00:00:00 2001 From: bigfooted Date: Tue, 16 Apr 2024 23:38:05 +0200 Subject: [PATCH 02/20] fix corner node temp and dens --- SU2_CFD/src/solvers/CEulerSolver.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 9385c9c874c..feaf719b308 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7115,16 +7115,24 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (geometry->nodes->GetViscousBoundary(iPoint)) { - V_inlet[0] = nodes->GetTemperature(iPoint); - /*--- Impose the wall velocity from the interior. ---*/ + + Velocity2 = 0.0; for (iDim = 0; iDim < nDim; iDim++){ V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; } + Pressure = nodes->GetPressure(iPoint); + /*--- Match the pressure and energy at the wall. ---*/ - V_inlet[nDim+1] = nodes->GetPressure(iPoint); - V_inlet[nDim+3] = nodes->GetPressure(iPoint)/(Density*Gamma_Minus_One) + nodes->GetPressure(iPoint)/Density; // + Pressure/Density; + Density = Pressure/(Gas_Constant*Temperature); + Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; + if (tkeNeeded) Energy += GetTke_Inf(); + + V_inlet[nDim+1] = Pressure; + V_inlet[nDim+2] = Density; + V_inlet[nDim+3] = Energy + Pressure/Density; } From f4f8a2bf74b2b69627f3ee2692e446a1acb0f440 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 17 Apr 2024 00:09:54 +0200 Subject: [PATCH 03/20] fix warning for temperature maybe used uninitialized --- SU2_CFD/src/solvers/CEulerSolver.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index feaf719b308..7bac11d2e4d 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -6880,10 +6880,11 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short val_marker) { unsigned short iDim; unsigned long iVertex, iPoint; - su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Temperature, Riemann, + su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Riemann, Pressure, Density, Energy, Flow_Dir[MAXNDIM], Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM]; su2double *V_inlet, *V_domain; + su2double Temperature = 288.15; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One; @@ -7096,8 +7097,8 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ - - V_inlet[0] = Pressure / ( Gas_Constant * Density); + Temperature = Pressure / ( Gas_Constant * Density); + V_inlet[0] = Temperature; for (iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; V_inlet[nDim+1] = Pressure; From 4844ba99783a68a6f6789017f0023c3d22e540a6 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 17 Apr 2024 00:11:51 +0200 Subject: [PATCH 04/20] move comment --- SU2_CFD/src/solvers/CEulerSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 7bac11d2e4d..bf465b5cd14 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7124,9 +7124,9 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; } - Pressure = nodes->GetPressure(iPoint); + /*--- Match the pressure, density and energy at the wall. ---*/ - /*--- Match the pressure and energy at the wall. ---*/ + Pressure = nodes->GetPressure(iPoint); Density = Pressure/(Gas_Constant*Temperature); Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; if (tkeNeeded) Energy += GetTke_Inf(); From 0cb2f38c5dc923836e193b9c7dcf1a525a51434c Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Apr 2024 11:28:14 +0200 Subject: [PATCH 05/20] fix bounded scalar for compressible --- Common/src/CConfig.cpp | 6 +- SU2_CFD/include/solvers/CSpeciesSolver.hpp | 20 ++++ SU2_CFD/src/solvers/CEulerSolver.cpp | 114 ++++++++++++++------- SU2_CFD/src/solvers/CSpeciesSolver.cpp | 16 ++- TestCases/parallel_regression.py | 2 +- 5 files changed, 115 insertions(+), 43 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18da76556c4..e4a9f562c7a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5538,16 +5538,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i /*--- Define some variables for flamelet model. ---*/ if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) { /*--- The controlling variables are progress variable, total enthalpy, and optionally mixture fraction ---*/ - //n_control_vars = nSpecies - n_user_scalars; if (n_control_vars != (nSpecies - n_user_scalars)) - SU2_MPI::Error("Number of initial species incompatbile with number of controlling variables and user scalars.", CURRENT_FUNCTION); + SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION); /*--- We can have additional user defined transported scalars ---*/ n_scalars = n_control_vars + n_user_scalars; } - if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) { - SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION); - } } void CConfig::SetMarkers(SU2_COMPONENT val_software) { diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 1ee78d90286..717f855c59b 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -153,6 +153,26 @@ class CSpeciesSolver : public CScalarSolver { * (like BC_Isothermal_Wall) are implemented the respective CHeatSolver implementations can act as a starting * point ---*/ + /*! + * \brief Generic implementation of the isothermal wall also covering CHT cases, + * for which the wall temperature is given by GetConjugateHeatVariable. + */ + void BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker, + bool cht_mode = false); + + /*! + * \brief Impose the scalar wall boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override; + /*! * \brief Source term computation for axisymmetric flow. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index bf465b5cd14..f49b19ef5a4 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -341,6 +341,10 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, CommunicateInitialState(geometry, config); + /*--- Sizing edge mass flux array ---*/ + if (config->GetBounded_Scalar()) + EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0); + /*--- Add the solver name.. ---*/ SolverName = "C.FLOW"; @@ -1749,6 +1753,7 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta CConfig *config, unsigned short iMesh, unsigned short iRKStep) { EdgeFluxResidual(geometry, solver_container, config); + const bool bounded_scalar = config->GetBounded_Scalar(); } void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container, @@ -1772,6 +1777,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); + const bool bounded_scalar = config->GetBounded_Scalar(); /*--- Non-physical counter. ---*/ unsigned long counter_local = 0; @@ -1937,6 +1943,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain auto residual = numerics->ComputeResidual(config); + if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + /*--- Set the final value of the Roe dissipation coefficient ---*/ if ((kind_dissipation != NO_ROELOWDISS) && (MGLevel != MESH_0)) { @@ -5067,8 +5075,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, su2double Velocity_i[MAXNDIM]={0}; for (auto iDim=0u; iDim < nDim; iDim++) Velocity_i[iDim] = nodes->GetVelocity(iPoint,iDim); - - const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i); + + const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i); const auto Density_i = nodes->GetDensity(iPoint); @@ -5090,7 +5098,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, T_Total{0}, P_Total{0}, Density_e{0}, StaticEnthalpy_e{0}, StaticEnergy_e{0}; su2double Velocity2_e{0}, NormalVelocity{0}, TangVelocity{0}, VelMag_e{0}; - su2double Velocity_e[MAXNDIM] = {0}; + su2double Velocity_e[MAXNDIM] = {0}; const su2double * Flow_Dir, * Mach; switch(config->GetKind_Data_Riemann(Marker_Tag)) { @@ -5146,10 +5154,10 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, GetFluidModel()->SetTDState_PT(P_static, T_static); /* --- Compute the boundary state u_e --- */ - for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed(); - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; @@ -5176,7 +5184,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed(); - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; @@ -5208,10 +5216,10 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, /* --- Compute the boundary state u_e --- */ GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e); - for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Velocity_i[iDim]; - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; break; @@ -5235,7 +5243,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, } /*--- Flow eigenvalues, boundary state u_e and u_i ---*/ - su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0}, + su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0}, dw[MAXNVAR]={0}; u_e[0] = Density_e; u_i[0] = Density_i; @@ -5302,7 +5310,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - su2double **Jacobian_b = new su2double*[nVar], + su2double **Jacobian_b = new su2double*[nVar], **Jacobian_i = new su2double*[nVar]; su2double **DubDu = new su2double*[nVar]; for (auto iVar = 0u; iVar < nVar; iVar++){ @@ -5333,7 +5341,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, if (dynamic_grid){ const auto gridVel = geometry->nodes->GetGridVel(iPoint); const auto projVelocity = GeometryToolbox::DotProduct(nDim, gridVel, Normal); - for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) Jacobian_b[iVar][iVar] -= projVelocity; } @@ -5394,7 +5402,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); /*--- Secondary variables ---*/ - + auto S_domain = nodes->GetSecondary(iPoint); /*--- Compute secondary thermodynamic properties (partial derivatives...) ---*/ @@ -6922,7 +6930,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_domain = nodes->GetPrimitive(iPoint); - /*--- Build the fictitious intlet state based on characteristics ---*/ + /*--- Build the fictitious inlet state based on characteristics ---*/ /*--- Subsonic inflow: there is one outgoing characteristic (u-c), @@ -7116,27 +7124,61 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (geometry->nodes->GetViscousBoundary(iPoint)) { - /*--- Impose the wall velocity from the interior. ---*/ + switch (Kind_Inlet) { - Velocity2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++){ - V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); - Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; - } + /*--- Total properties have been specified at the inlet. ---*/ - /*--- Match the pressure, density and energy at the wall. ---*/ + case INLET_TYPE::TOTAL_CONDITIONS: { - Pressure = nodes->GetPressure(iPoint); - Density = Pressure/(Gas_Constant*Temperature); - Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- Impose the wall velocity from the interior wall node. ---*/ - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure/Density; - } + Velocity2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; + } + + /*--- Match the pressure, density and energy at the wall. ---*/ + + Pressure = nodes->GetPressure(iPoint); + Density = Pressure / (Gas_Constant * Temperature); + Energy = Pressure / (Density*Gamma_Minus_One) + 0.5 * Velocity2; + if (tkeNeeded) Energy += GetTke_Inf(); + + V_inlet[nDim+1] = Pressure; + V_inlet[nDim+2] = Density; + V_inlet[nDim+3] = Energy + Pressure/Density; + break; + } + case INLET_TYPE::MASS_FLOW: { + /*--- Impose the wall velocity from the interior wall node. ---*/ + + Velocity2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; + } + + Pressure = nodes->GetPressure(iPoint); + Density = nodes->GetDensity(iPoint); + + Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; + + if (tkeNeeded) Energy += GetTke_Inf(); + + V_inlet[nDim+3] = Energy + Pressure/Density; + + break; + } + + default: + SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); + break; + + } + } /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_inlet); @@ -8821,9 +8863,9 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon const auto nSpanWiseSections = config->GetnSpanWiseSections(); const auto iZone = config->GetiZone(); - + for (auto iSpan= 0u; iSpan < nSpanWiseSections; iSpan++){ - su2double TotalAreaVelocity[MAXNDIM]={0.0}, + su2double TotalAreaVelocity[MAXNDIM]={0.0}, TotalAreaPressure{0}, TotalAreaDensity{0}; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++){ @@ -9013,7 +9055,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC TotalAreaVelocity[iDim] += Area*Velocity[iDim]; TotalMassVelocity[iDim] += Area*(Density*TurboVelocity[0] )*Velocity[iDim]; } - + TotalFluxes[0] += Area*(Density*TurboVelocity[0]); TotalFluxes[1] += Area*(Density*TurboVelocity[0]*TurboVelocity[0] + Pressure); for (auto iDim = 2; iDim < nDim+1; iDim++) @@ -9062,7 +9104,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC UpdateTotalQuantities(iMarker, jSpan, iVertex); } } - } + } } // marker_flag match } // iMarkerTP match @@ -9277,7 +9319,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC OmegaOut[iMarkerTP - 1][iSpan] = AverageOmega[iMarker][iSpan]; NuOut[iMarkerTP - 1][iSpan] = AverageNu[iMarker][iSpan]; } - + auto TurboVel = (marker_flag == INFLOW) ? TurboVelocityIn[iMarkerTP - 1][iSpan] : TurboVelocityOut[iMarkerTP - 1][iSpan]; if (performance_average_process == MIXEDOUT) { @@ -9307,7 +9349,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC if(config->GetKind_Data_Giles(Marker_Tag) == RADIAL_EQUILIBRIUM){ RadialEquilibriumPressure[iMarker][nSpanWiseSections/2] = config->GetGiles_Var1(Marker_Tag)/config->GetPressure_Ref(); } - } + } for (auto iSpan= nSpanWiseSections/2; iSpan < nSpanWiseSections-1; iSpan++){ const auto Radius2 = geometry->GetTurboRadius(iMarker,iSpan+1); const auto Radius1 = geometry->GetTurboRadius(iMarker,iSpan); @@ -9351,9 +9393,9 @@ void CEulerSolver::MixedOut_Average(CConfig *config, su2double val_init_pressure su2double vel[MAXNDIM] = {0}; vel[0] = (val_Averaged_Flux[1] - pressure_mix) / val_Averaged_Flux[0]; - for (auto iDim = 1u; iDim < nDim; iDim++) + for (auto iDim = 1u; iDim < nDim; iDim++) vel[iDim] = val_Averaged_Flux[iDim+1] / val_Averaged_Flux[0]; - + const su2double velsq = GeometryToolbox::DotProduct(nDim, vel, vel); diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index f5f5ed16603..144c2afc4a4 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -345,7 +345,7 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C /*--- Identify the boundary by string name ---*/ string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - + if (config->GetMarker_StrongBC(Marker_Tag)==true) { nodes->SetSolution_Old(iPoint, Inlet_SpeciesVars[val_marker][iVertex]); @@ -406,6 +406,20 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C END_SU2_OMP_FOR } +void CSpeciesSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, + CNumerics* visc_numerics, CConfig* config, unsigned short val_marker, + bool cht_mode) { + +} + +void CSpeciesSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + +} + + void CSpeciesSolver::SetInletAtVertex(const su2double *val_inlet, unsigned short iMarker, unsigned long iVertex) { diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 51b18e48741..567e9adc016 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -302,7 +302,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 = [-6.653802, -1.181430, -0.794887, 0.000611, -0.000369, 0.000736, -0.001104, 596.690000, 299.800000, 296.890000, 21.492000, 0.563990, 37.148000, 2.278700] + flatplate_udobj.test_vals = [-6.653278, -1.180894, -0.794611, 0.000613, -0.000370, 0.000737, -0.001107, 596.680000, 299.790000, 296.890000, 21.488000, 0.564630, 37.153000, 2.284900] test_list.append(flatplate_udobj) # Laminar cylinder (steady) From 819b46ca14d2a3114f9957575bb6abcdb9294a36 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 20 Apr 2024 16:19:26 +0200 Subject: [PATCH 06/20] fix bounded scalar for compressible --- SU2_CFD/src/solvers/CEulerSolver.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index f49b19ef5a4..5e24a1ad039 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1753,7 +1753,6 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta CConfig *config, unsigned short iMesh, unsigned short iRKStep) { EdgeFluxResidual(geometry, solver_container, config); - const bool bounded_scalar = config->GetBounded_Scalar(); } void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container, @@ -1777,8 +1776,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); - const bool bounded_scalar = config->GetBounded_Scalar(); - + //const bool bounded_scalar = config->GetBounded_Scalar(); /*--- Non-physical counter. ---*/ unsigned long counter_local = 0; SU2_OMP_MASTER @@ -1943,7 +1941,9 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain auto residual = numerics->ComputeResidual(config); - if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; + // this does not seem to be necessary + cout << "bounded" << endl; + //if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; /*--- Set the final value of the Roe dissipation coefficient ---*/ From 7a5104d5627c001f2975b6e2e750cd8d796a58ea Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 10 Aug 2024 16:08:41 +0200 Subject: [PATCH 07/20] update regression tests --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 4 ++-- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index bd7635e1413..3b93ec68f68 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -139,7 +139,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.494728, -7.712546, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.494736, -7.712091, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.494717, -7.711274, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6ca9e503c0d..e2dd41e0c18 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -344,7 +344,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.492939, -7.672950, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.492865, -7.672526, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.492864, -7.671632, -0.000000, 2.085796] test_list.append(poiseuille_profile) @@ -1015,7 +1015,7 @@ def main(): flatplate_unsteady.cfg_dir = "navierstokes/flatplate" flatplate_unsteady.cfg_file = "lam_flatplate_unst.cfg" flatplate_unsteady.test_iter = 3 - flatplate_unsteady.test_vals = [7.9509e-06, -8.868859, -8.231652, -6.283262, -5.466675, -3.391163, 0.002078, -0.343642] + flatplate_unsteady.test_vals = [0.000008, -8.869318, -8.232471, -6.283031, -5.467151, -3.391524, 0.002078, -0.343471] flatplate_unsteady.unsteady = True test_list.append(flatplate_unsteady) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 59576fdc410..9c69ebf4792 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -524,7 +524,7 @@ def main(): pywrapper_wavy_wall_steady.cfg_dir = "py_wrapper/wavy_wall" pywrapper_wavy_wall_steady.cfg_file = "run_steady.py" pywrapper_wavy_wall_steady.test_iter = 100 - pywrapper_wavy_wall_steady.test_vals = [-1.360044, 2.580709, -2.892473] + pywrapper_wavy_wall_steady.test_vals = [-1.359409, 2.580816, -2.892697] pywrapper_wavy_wall_steady.command = TestCase.Command("mpirun -n 2", "python", "run_steady.py") pywrapper_wavy_wall_steady.timeout = 1600 pywrapper_wavy_wall_steady.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index fcdacb032a7..70e937959a9 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -198,7 +198,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.494681, -7.711642, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.494668, -7.710678, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.494684, -7.711379, -0.000000, 2.085796] #last 4 columns test_list.append(poiseuille_profile) @@ -1605,7 +1605,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.124164, -1.544359, -3.808866, 1.338411, -0.752679, 0.161436, -1.2391e-02, 5.1662e-01, -5.2901e-01] + pywrapper_custom_inlet.test_vals = [-4.124228, -1.544435, -3.806720, 1.338337, -0.752645, 0.162022, -0.012134, 0.516870, -0.529000] pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 From d179a4e116bd99ddcc6e22aaa9facc51a22933fd Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 10 Aug 2024 16:14:20 +0200 Subject: [PATCH 08/20] remove comments on compressible bounded scalar --- SU2_CFD/src/solvers/CEulerSolver.cpp | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 73c2ead44fa..789c068d894 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -341,10 +341,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, CommunicateInitialState(geometry, config); - /*--- Sizing edge mass flux array ---*/ - if (config->GetBounded_Scalar()) - EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0); - /*--- Add the solver name.. ---*/ SolverName = "C.FLOW"; @@ -1776,7 +1772,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); - //const bool bounded_scalar = config->GetBounded_Scalar(); + /*--- Non-physical counter. ---*/ unsigned long counter_local = 0; SU2_OMP_MASTER @@ -1941,10 +1937,6 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain auto residual = numerics->ComputeResidual(config); - // this does not seem to be necessary - cout << "bounded" << endl; - //if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0]; - /*--- Set the final value of the Roe dissipation coefficient ---*/ if ((kind_dissipation != NO_ROELOWDISS) && (MGLevel != MESH_0)) { From 1757060262add279f321123d48d8c3f298288f44 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 10 Aug 2024 18:13:34 +0200 Subject: [PATCH 09/20] update regression --- TestCases/parallel_regression.py | 4 ++-- TestCases/serial_regression.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 02815976ef2..9dc96e30afd 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -311,7 +311,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 = [-6.664134, -1.190073, -0.954366, 0.000641, -0.000633, 0.000548, -0.001181, 596.940000, 300.020000, 296.920000, 22.201000, 0.525750, 37.278000, 2.347900] + flatplate_udobj.test_vals = [-6.663514, -1.189436, -0.954938, 0.000644, -0.000635, 0.000549, -0.001183, 596.930000, 300.010000, 296.920000, 22.187000, 0.527090, 37.284000, 2.354200] test_list.append(flatplate_udobj) # Laminar cylinder (steady) @@ -1023,7 +1023,7 @@ def main(): flatplate_unsteady.cfg_dir = "navierstokes/flatplate" flatplate_unsteady.cfg_file = "lam_flatplate_unst.cfg" flatplate_unsteady.test_iter = 3 - flatplate_unsteady.test_vals = [0.000008, -8.869318, -8.232471, -6.283031, -5.467151, -3.391524, 0.002078, -0.343471] + flatplate_unsteady.test_vals = [-8.876870, -8.250064, -6.293918, -5.469416, -3.398952, 0.002075, -0.324126] flatplate_unsteady.unsteady = True test_list.append(flatplate_unsteady) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e1b0b85f0e8..6804176eac4 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1605,7 +1605,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.124228, -1.544435, -3.806720, 1.338337, -0.752645, 0.162022, -0.012134, 0.516870, -0.529000] + pywrapper_custom_inlet.test_vals = [-4.120437, -1.540129, -3.563086, 1.342567, -0.748783, 0.161976, -0.012959, 0.516320, -0.529280] pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 From 0faa67e9c9f5dd7a9210ce00abd0892c6a8faad5 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sat, 10 Aug 2024 23:15:41 +0200 Subject: [PATCH 10/20] update regression profile --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index a1ae9811bcb..ad6eb0b13ab 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -139,7 +139,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.494736, -7.712091, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.485991, -7.613038, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.494717, -7.711274, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e5b06814dd1..c21f7282fb0 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -344,7 +344,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.492865, -7.672526, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.483988, -7.577838, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.492864, -7.671632, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 6804176eac4..17350f42bae 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -198,7 +198,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.494668, -7.710678, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.485982, -7.613005, -0.000000, 2.085796] poiseuille_profile.test_vals_aarch64 = [-12.494684, -7.711379, -0.000000, 2.085796] #last 4 columns test_list.append(poiseuille_profile) From 48077f0c07e558542cf969f2e6416bfed8aa5749 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 27 Aug 2024 18:59:49 +0200 Subject: [PATCH 11/20] Update SU2_CFD/src/solvers/CEulerSolver.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/src/solvers/CEulerSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 29212a68775..b0253ed592a 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7080,7 +7080,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ - Temperature = Pressure / ( Gas_Constant * Density); + const su2double Temperature = Pressure / ( Gas_Constant * Density); V_inlet[0] = Temperature; for (iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; From 6ab0c5e98b91bb7175f2cf0e14f577860005e6c0 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 4 Oct 2024 10:08:12 +0200 Subject: [PATCH 12/20] add message, add P and rho --- Common/src/CConfig.cpp | 4 ++++ SU2_CFD/src/solvers/CEulerSolver.cpp | 4 +++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 27dcf568b79..f37f20cf9e5 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5594,6 +5594,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i n_scalars = n_control_vars + n_user_scalars; } + if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) { + SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION); + } + } void CConfig::SetMarkers(SU2_COMPONENT val_software) { diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 744240ebd15..419cbaa471c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -6868,7 +6868,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Pressure, Density, Energy, Flow_Dir[MAXNDIM], Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM]; su2double *V_inlet, *V_domain; - su2double Temperature = 288.15; + su2double Temperature; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One; @@ -7144,6 +7144,8 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (tkeNeeded) Energy += GetTke_Inf(); + V_inlet[nDim+1] = Pressure; + V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure/Density; break; From b8a970a643b1e053b66435d9a492bdec91887cf6 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 4 Oct 2024 22:17:53 +0200 Subject: [PATCH 13/20] update --- SU2_CFD/include/solvers/CSpeciesSolver.hpp | 20 -------------------- SU2_CFD/src/solvers/CEulerSolver.cpp | 11 ++++++----- SU2_CFD/src/solvers/CSpeciesSolver.cpp | 14 -------------- 3 files changed, 6 insertions(+), 39 deletions(-) diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 1ec74e85d3c..7236989afd4 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -152,26 +152,6 @@ class CSpeciesSolver : public CScalarSolver { * (like BC_Isothermal_Wall) are implemented the respective CHeatSolver implementations can act as a starting * point ---*/ - /*! - * \brief Generic implementation of the isothermal wall also covering CHT cases, - * for which the wall temperature is given by GetConjugateHeatVariable. - */ - void BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, - CNumerics* visc_numerics, CConfig* config, unsigned short val_marker, - bool cht_mode = false); - - /*! - * \brief Impose the scalar wall boundary condition. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] conv_numerics - Description of the numerical method. - * \param[in] visc_numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - * \param[in] val_marker - Surface marker where the boundary condition is applied. - */ - void BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, - CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override; - /*! * \brief Source term computation for axisymmetric flow. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index be6854cc955..7adac5cf4c3 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7096,7 +7096,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ - const su2double Temperature = Pressure / ( Gas_Constant * Density); + Temperature = Pressure / ( Gas_Constant * Density); V_inlet[0] = Temperature; for (iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; @@ -7134,11 +7134,12 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Pressure = nodes->GetPressure(iPoint); Density = Pressure / (Gas_Constant * Temperature); Energy = Pressure / (Density*Gamma_Minus_One) + 0.5 * Velocity2; + if (tkeNeeded) Energy += GetTke_Inf(); V_inlet[nDim+1] = Pressure; V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure/Density; + V_inlet[nDim+3] = Energy + Pressure / Density; break; } @@ -7154,14 +7155,14 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Pressure = nodes->GetPressure(iPoint); Density = nodes->GetDensity(iPoint); - + Temperature = Pressure / ( Gas_Constant * Density); Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + V_inlet[0] = Temperature; V_inlet[nDim+1] = Pressure; V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure/Density; + V_inlet[nDim+3] = Energy + Pressure / Density; break; } diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 3d2f7d04f5b..9d7c4494c51 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -406,20 +406,6 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C END_SU2_OMP_FOR } -void CSpeciesSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, - CNumerics* visc_numerics, CConfig* config, unsigned short val_marker, - bool cht_mode) { - -} - -void CSpeciesSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); - -} - - void CSpeciesSolver::SetInletAtVertex(const su2double *val_inlet, unsigned short iMarker, unsigned long iVertex) { From 843dc1f8c2379fbbbb18e9debf77028282fce573 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Mon, 7 Oct 2024 22:52:14 +0200 Subject: [PATCH 14/20] temperature variable --- SU2_CFD/src/solvers/CEulerSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 7adac5cf4c3..b9b7aa26b5d 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -6879,11 +6879,10 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short val_marker) { unsigned short iDim; unsigned long iVertex, iPoint; - su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Riemann, + su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Temperature, Riemann, Pressure, Density, Energy, Flow_Dir[MAXNDIM], Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM]; su2double *V_inlet, *V_domain; - su2double Temperature; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One; @@ -7137,6 +7136,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, if (tkeNeeded) Energy += GetTke_Inf(); + V_inlet[0] = Temperature; V_inlet[nDim+1] = Pressure; V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure / Density; From c9c2466584c55316aebd098ea80ff313349046d2 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 25 Oct 2024 13:54:36 +0200 Subject: [PATCH 15/20] merge w. develop, fix conflicts --- TestCases/parallel_regression.py | 8 -------- TestCases/serial_regression.py | 4 ---- 2 files changed, 12 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ee4db8097ac..d8e90de816e 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -311,11 +311,7 @@ def main(): flatplate_udobj.cfg_dir = "user_defined_functions" flatplate_udobj.cfg_file = "lam_flatplate.cfg" flatplate_udobj.test_iter = 20 -<<<<<<< HEAD flatplate_udobj.test_vals = [-6.663514, -1.189436, -0.954938, 0.000644, -0.000635, 0.000549, -0.001183, 596.930000, 300.010000, 296.920000, 22.187000, 0.527090, 37.284000, 2.354200] -======= - flatplate_udobj.test_vals = [-6.661084, -1.186974, -0.956148, 0.000640, -0.000642, 0.000539, -0.001181, 596.990000, 300.070000, 296.920000, 22.250000, 0.522670, 37.278000, 2.342900] ->>>>>>> develop test_list.append(flatplate_udobj) # Laminar cylinder (steady) @@ -1023,11 +1019,7 @@ def main(): flatplate_unsteady.cfg_dir = "navierstokes/flatplate" flatplate_unsteady.cfg_file = "lam_flatplate_unst.cfg" flatplate_unsteady.test_iter = 3 -<<<<<<< HEAD flatplate_unsteady.test_vals = [-8.876870, -8.250064, -6.293918, -5.469416, -3.398952, 0.002075, -0.324126] -======= - flatplate_unsteady.test_vals = [0.000008, -8.874953, -8.250030, -6.306009, -5.468992, -3.398179, 0.002075, -0.326075] ->>>>>>> develop flatplate_unsteady.unsteady = True test_list.append(flatplate_unsteady) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index d4fc8183064..d8797f27a27 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1605,11 +1605,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 -<<<<<<< HEAD pywrapper_custom_inlet.test_vals = [-4.120437, -1.540129, -3.563086, 1.342567, -0.748783, 0.161976, -0.012959, 0.516320, -0.529280] -======= - pywrapper_custom_inlet.test_vals = [-4.120585, -1.540325, -3.560392, 1.342419, -0.748768, 0.161248, -0.013208, 0.515780, -0.528990] ->>>>>>> develop pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 From 1a9776e0f2be22b4449f8ddbcb274069b68df700 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 30 Oct 2024 20:12:40 +0100 Subject: [PATCH 16/20] typo --- Common/src/geometry/CMultiGridGeometry.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 8c5bfddf39e..c3c600e8a12 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -36,8 +36,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- Create a queue system to do the agglomeration 1st) More than two markers ---> Vertices (never agglomerate) 2nd) Two markers ---> Edges (agglomerate if same BC, never agglomerate if different BC) - 3rd) One marker ---> Surface (always agglomarate) - 4th) No marker ---> Internal Volume (always agglomarate) ---*/ + 3rd) One marker ---> Surface (always agglomerate) + 4th) No marker ---> Internal Volume (always agglomerate) ---*/ /*--- Set a marker to indicate indirect agglomeration, for quads and hexs, i.e. consider up to neighbors of neighbors of neighbors. From 73e5e057e326675de4072397553fd312d9c8cf19 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 30 Oct 2024 21:04:54 +0100 Subject: [PATCH 17/20] using updated testcase repo --- .github/workflows/regression.yml | 6 +++--- TestCases/parallel_regression.py | 4 ++-- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index ca8bf820a84..1ee1b583e1d 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -209,7 +209,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_aachen_test_case -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: @@ -255,7 +255,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_aachen_test_case -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: @@ -300,7 +300,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_aachen_test_case -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index d8e90de816e..5452a94166f 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -311,7 +311,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 = [-6.663514, -1.189436, -0.954938, 0.000644, -0.000635, 0.000549, -0.001183, 596.930000, 300.010000, 296.920000, 22.187000, 0.527090, 37.284000, 2.354200] + flatplate_udobj.test_vals = [-6.660438, -1.186310, -0.956714, 0.000642, -0.000644, 0.000540, -0.001183, 596.990000, 300.060000, 296.920000, 22.236000, 0.524020, 37.284000, 2.349300] test_list.append(flatplate_udobj) # Laminar cylinder (steady) @@ -1019,7 +1019,7 @@ def main(): flatplate_unsteady.cfg_dir = "navierstokes/flatplate" flatplate_unsteady.cfg_file = "lam_flatplate_unst.cfg" flatplate_unsteady.test_iter = 3 - flatplate_unsteady.test_vals = [-8.876870, -8.250064, -6.293918, -5.469416, -3.398952, 0.002075, -0.324126] + flatplate_unsteady.test_vals = [-8.875126, -8.250188, -6.305789, -5.469452, -3.398228, 0.002075, -0.326018] flatplate_unsteady.unsteady = True test_list.append(flatplate_unsteady) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index f0d7bc6c416..7d0f0d10049 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -524,7 +524,7 @@ def main(): pywrapper_wavy_wall_steady.cfg_dir = "py_wrapper/wavy_wall" pywrapper_wavy_wall_steady.cfg_file = "run_steady.py" pywrapper_wavy_wall_steady.test_iter = 100 - pywrapper_wavy_wall_steady.test_vals = [-1.359409, 2.580816, -2.892697] + pywrapper_wavy_wall_steady.test_vals = [-1.352680, 2.579322, -2.898321] pywrapper_wavy_wall_steady.command = TestCase.Command("mpirun -n 2", "python", "run_steady.py") pywrapper_wavy_wall_steady.timeout = 1600 pywrapper_wavy_wall_steady.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index d8797f27a27..53b472c1434 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1605,7 +1605,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.120437, -1.540129, -3.563086, 1.342567, -0.748783, 0.161976, -0.012959, 0.516320, -0.529280] + pywrapper_custom_inlet.test_vals = [-4.120910, -1.540677, -3.567156, 1.342054, -0.748463, 0.161782, -0.012949, 0.516060, -0.529000] pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 From 787960ab5eedac7eb84bfa31b6f774644cfc1f07 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Wed, 30 Oct 2024 21:39:18 +0100 Subject: [PATCH 18/20] update poiseuille regressions --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index e9fc251fe1b..fe2634fb2a5 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -136,7 +136,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.485991, -7.613038, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.008990, -7.262477, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.494717, -7.711274, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 5452a94166f..44f11fe177b 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -344,7 +344,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.483988, -7.577838, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.007503, -7.227301, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.492864, -7.671632, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 53b472c1434..9fc971df334 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -198,7 +198,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.485982, -7.613005, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.009019, -7.262438, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.494684, -7.711379, -0.000000, 2.085796] #last 4 columns test_list.append(poiseuille_profile) From 9d1bfd091f2bca0348ca1dc9df0daf7208b726f2 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 8 Dec 2024 10:34:08 -0800 Subject: [PATCH 19/20] simplify --- SU2_CFD/src/solvers/CEulerSolver.cpp | 423 +++++++++++---------------- 1 file changed, 164 insertions(+), 259 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 25e9d8d892a..01626563d40 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -6882,7 +6882,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { unsigned short iDim; - unsigned long iVertex, iPoint; su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Temperature, Riemann, Pressure, Density, Energy, Flow_Dir[MAXNDIM], Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM]; @@ -6891,353 +6890,259 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One; const su2double Gas_Constant = config->GetGas_ConstantND(); - const auto Kind_Inlet = config->GetKind_Inlet(); + const auto Kind_Inlet_Cfg = config->GetKind_Inlet(); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the inlet ---*/ V_inlet = GetCharacPrimVar(val_marker, iVertex); - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto 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; - - /*--- Retrieve solution at this boundary node ---*/ - - V_domain = nodes->GetPrimitive(iPoint); - - /*--- Build the fictitious inlet state based on characteristics ---*/ - - - /*--- Subsonic inflow: there is one outgoing characteristic (u-c), - therefore we can specify all but one state variable at the inlet. - The outgoing Riemann invariant provides the final piece of info. - Adapted from an original implementation in the Stanford University - multi-block (SUmb) solver in the routine bcSubsonicInflow.f90 - written by Edwin van der Weide, last modified 04-20-2009. ---*/ - - switch (Kind_Inlet) { - - /*--- Total properties have been specified at the inlet. ---*/ - - case INLET_TYPE::TOTAL_CONDITIONS: { - - /*--- Retrieve the specified total conditions for this inlet. ---*/ - - P_Total = Inlet_Ptotal[val_marker][iVertex]; - T_Total = Inlet_Ttotal[val_marker][iVertex]; - const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; - const su2double mag = GeometryToolbox::Norm(nDim, dir); - for (iDim = 0; iDim < nDim; iDim++) { - Flow_Dir[iDim] = dir[iDim] / mag; - } - - /*--- Non-dim. the inputs if necessary. ---*/ - - P_Total /= config->GetPressure_Ref(); - T_Total /= config->GetTemperature_Ref(); - - /*--- Store primitives and set some variables for clarity. ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - Density = V_domain[nDim+2]; - Velocity2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity[iDim] = V_domain[iDim+1]; - Velocity2 += Velocity[iDim]*Velocity[iDim]; - } - Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; - Pressure = V_domain[nDim+1]; - H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; - SoundSpeed2 = Gamma*Pressure/Density; + /*--- Normal vector for this vertex (negate for outward convention) ---*/ - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); - Riemann = 2.0*sqrt(SoundSpeed2)/Gamma_Minus_One; - for (iDim = 0; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; + Area = GeometryToolbox::Norm(nDim, Normal); + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = Normal[iDim]/Area; - /*--- Total speed of sound ---*/ + /*--- Retrieve solution at this boundary node ---*/ - SoundSpeed_Total2 = Gamma_Minus_One*(H_Total - (Energy + Pressure/Density)+0.5*Velocity2) + SoundSpeed2; + V_domain = nodes->GetPrimitive(iPoint); - /*--- Dot product of normal and flow direction. This should - be negative due to outward facing boundary normal convention. ---*/ + /*--- On intersections with viscous walls we set the total conditions equal to + the static conditions of the wall point. This avoids a jump in pressure, + energy, etc. at this location. ---*/ - alpha = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - alpha += UnitNormal[iDim]*Flow_Dir[iDim]; + const auto Kind_Inlet = geometry->nodes->GetViscousBoundary(iPoint) + ? INLET_TYPE::TOTAL_CONDITIONS : Kind_Inlet_Cfg; - /*--- Coefficients in the quadratic equation for the velocity ---*/ + /*--- Build the fictitious inlet state based on characteristics ---*/ - aa = 1.0 + 0.5*Gamma_Minus_One*alpha*alpha; - bb = -1.0*Gamma_Minus_One*alpha*Riemann; - cc = 0.5*Gamma_Minus_One*Riemann*Riemann - -2.0*SoundSpeed_Total2/Gamma_Minus_One; + /*--- Subsonic inflow: there is one outgoing characteristic (u-c), + therefore we can specify all but one state variable at the inlet. + The outgoing Riemann invariant provides the final piece of info. + Adapted from an original implementation in the Stanford University + multi-block (SUmb) solver in the routine bcSubsonicInflow.f90 + written by Edwin van der Weide, last modified 04-20-2009. ---*/ - /*--- Solve quadratic equation for velocity magnitude. Value must - be positive, so the choice of root is clear. ---*/ + switch (Kind_Inlet) { - dd = bb*bb - 4.0*aa*cc; - dd = sqrt(max(0.0, dd)); - Vel_Mag = (-bb + dd)/(2.0*aa); - Vel_Mag = max(0.0, Vel_Mag); - Velocity2 = Vel_Mag*Vel_Mag; + /*--- Total properties have been specified at the inlet. ---*/ - /*--- Compute speed of sound from total speed of sound eqn. ---*/ + case INLET_TYPE::TOTAL_CONDITIONS: { - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; + /*--- Retrieve the specified total conditions for this inlet. ---*/ - /*--- Mach squared (cut between 0-1), use to adapt velocity ---*/ + P_Total = Inlet_Ptotal[val_marker][iVertex]; + T_Total = Inlet_Ttotal[val_marker][iVertex]; + const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; + const su2double mag = GeometryToolbox::Norm(nDim, dir); + for (iDim = 0; iDim < nDim; iDim++) { + Flow_Dir[iDim] = dir[iDim] / mag; + } - Mach2 = Velocity2/SoundSpeed2; - Mach2 = min(1.0, Mach2); - Velocity2 = Mach2*SoundSpeed2; - Vel_Mag = sqrt(Velocity2); - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; + /*--- Non-dim. the inputs if necessary. ---*/ - /*--- Compute new velocity vector at the inlet ---*/ + P_Total /= config->GetPressure_Ref(); + T_Total /= config->GetTemperature_Ref(); - for (iDim = 0; iDim < nDim; iDim++) - Velocity[iDim] = Vel_Mag*Flow_Dir[iDim]; + if (geometry->nodes->GetViscousBoundary(iPoint)) { + P_Total = nodes->GetPressure(iPoint); + T_Total = nodes->GetTemperature(iPoint); + } - /*--- Static temperature from the speed of sound relation ---*/ + /*--- Store primitives and set some variables for clarity. ---*/ - Temperature = SoundSpeed2/(Gamma*Gas_Constant); + Density = V_domain[nDim+2]; + Velocity2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity[iDim] = V_domain[iDim+1]; + Velocity2 += Velocity[iDim]*Velocity[iDim]; + } + Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + Pressure = V_domain[nDim+1]; + H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; + SoundSpeed2 = Gamma*Pressure/Density; - /*--- Static pressure using isentropic relation at a point ---*/ + /*--- Compute the acoustic Riemann invariant that is extrapolated + from the domain interior. ---*/ - Pressure = P_Total*pow((Temperature/T_Total), Gamma/Gamma_Minus_One); + Riemann = 2.0*sqrt(SoundSpeed2)/Gamma_Minus_One; + for (iDim = 0; iDim < nDim; iDim++) + Riemann += Velocity[iDim]*UnitNormal[iDim]; - /*--- Density at the inlet from the gas law ---*/ + /*--- Total speed of sound ---*/ - Density = Pressure/(Gas_Constant*Temperature); + SoundSpeed_Total2 = Gamma_Minus_One*(H_Total - (Energy + Pressure/Density)+0.5*Velocity2) + SoundSpeed2; - /*--- Using pressure, density, & velocity, compute the energy ---*/ + /*--- Dot product of normal and flow direction. This should + be negative due to outward facing boundary normal convention. ---*/ - Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + alpha = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + alpha += UnitNormal[iDim]*Flow_Dir[iDim]; - /*--- Primitive variables, using the derived quantities ---*/ + /*--- Coefficients in the quadratic equation for the velocity ---*/ - V_inlet[0] = Temperature; - for (iDim = 0; iDim < nDim; iDim++) - V_inlet[iDim+1] = Velocity[iDim]; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure/Density; + aa = 1.0 + 0.5*Gamma_Minus_One*alpha*alpha; + bb = -1.0*Gamma_Minus_One*alpha*Riemann; + cc = 0.5*Gamma_Minus_One*Riemann*Riemann + -2.0*SoundSpeed_Total2/Gamma_Minus_One; - break; - } - /*--- Mass flow has been specified at the inlet. ---*/ + /*--- Solve quadratic equation for velocity magnitude. Value must + be positive, so the choice of root is clear. ---*/ - case INLET_TYPE::MASS_FLOW: { + dd = bb*bb - 4.0*aa*cc; + dd = sqrt(max(0.0, dd)); + Vel_Mag = (-bb + dd)/(2.0*aa); + Vel_Mag = max(0.0, Vel_Mag); + Velocity2 = Vel_Mag*Vel_Mag; - /*--- Retrieve the specified mass flow for the inlet. ---*/ + /*--- Compute speed of sound from total speed of sound eqn. ---*/ - Density = Inlet_Ttotal[val_marker][iVertex]; - Vel_Mag = Inlet_Ptotal[val_marker][iVertex]; - const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; - const su2double mag = GeometryToolbox::Norm(nDim, dir); - for (iDim = 0; iDim < nDim; iDim++) { - Flow_Dir[iDim] = dir[iDim] / mag; - } + SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - /*--- Non-dim. the inputs if necessary. ---*/ + /*--- Mach squared (cut between 0-1), use to adapt velocity ---*/ - Density /= config->GetDensity_Ref(); - Vel_Mag /= config->GetVelocity_Ref(); + Mach2 = Velocity2/SoundSpeed2; + Mach2 = min(1.0, Mach2); + Velocity2 = Mach2*SoundSpeed2; + Vel_Mag = sqrt(Velocity2); + SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - /*--- Get primitives from current inlet state. ---*/ + /*--- Compute new velocity vector at the inlet ---*/ - for (iDim = 0; iDim < nDim; iDim++) - Velocity[iDim] = nodes->GetVelocity(iPoint,iDim); - Pressure = nodes->GetPressure(iPoint); - SoundSpeed2 = Gamma*Pressure/V_domain[nDim+2]; + for (iDim = 0; iDim < nDim; iDim++) + Velocity[iDim] = Vel_Mag*Flow_Dir[iDim]; - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ + /*--- Static temperature from the speed of sound relation ---*/ - Riemann = Two_Gamma_M1*sqrt(SoundSpeed2); - for (iDim = 0; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; + Temperature = SoundSpeed2/(Gamma*Gas_Constant); - /*--- Speed of sound squared for fictitious inlet state ---*/ + /*--- Static pressure using isentropic relation at a point ---*/ - SoundSpeed2 = Riemann; - for (iDim = 0; iDim < nDim; iDim++) - SoundSpeed2 -= Vel_Mag*Flow_Dir[iDim]*UnitNormal[iDim]; + Pressure = P_Total*pow((Temperature/T_Total), Gamma/Gamma_Minus_One); - SoundSpeed2 = max(0.0,0.5*Gamma_Minus_One*SoundSpeed2); - SoundSpeed2 = SoundSpeed2*SoundSpeed2; + /*--- Density at the inlet from the gas law ---*/ - /*--- Pressure for the fictitious inlet state ---*/ + Density = Pressure/(Gas_Constant*Temperature); - Pressure = SoundSpeed2*Density/Gamma; + /*--- Using pressure, density, & velocity, compute the energy ---*/ - /*--- Energy for the fictitious inlet state ---*/ + Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; + if (tkeNeeded) Energy += GetTke_Inf(); - Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Vel_Mag*Vel_Mag; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- Primitive variables, using the derived quantities ---*/ - /*--- Primitive variables, using the derived quantities ---*/ - Temperature = Pressure / ( Gas_Constant * Density); - V_inlet[0] = Temperature; - for (iDim = 0; iDim < nDim; iDim++) - V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure/Density; + V_inlet[0] = Temperature; + for (iDim = 0; iDim < nDim; iDim++) + V_inlet[iDim+1] = Velocity[iDim]; + V_inlet[nDim+1] = Pressure; + V_inlet[nDim+2] = Density; + V_inlet[nDim+3] = Energy + Pressure/Density; - break; - } - default: - SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); - break; + break; } + /*--- Mass flow has been specified at the inlet. ---*/ - /*--- Check if the inlet node is shared with a viscous wall. ---*/ + case INLET_TYPE::MASS_FLOW: { - if (geometry->nodes->GetViscousBoundary(iPoint)) { + /*--- Retrieve the specified mass flow for the inlet. ---*/ - switch (Kind_Inlet) { + Density = Inlet_Ttotal[val_marker][iVertex]; + Vel_Mag = Inlet_Ptotal[val_marker][iVertex]; + const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; + const su2double mag = GeometryToolbox::Norm(nDim, dir); + for (iDim = 0; iDim < nDim; iDim++) { + Flow_Dir[iDim] = dir[iDim] / mag; + } - /*--- Total properties have been specified at the inlet. ---*/ + /*--- Non-dim. the inputs if necessary. ---*/ - case INLET_TYPE::TOTAL_CONDITIONS: { + Density /= config->GetDensity_Ref(); + Vel_Mag /= config->GetVelocity_Ref(); - /*--- Impose the wall velocity from the interior wall node. ---*/ + /*--- Get primitives from current inlet state. ---*/ - Velocity2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); - Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; - } + for (iDim = 0; iDim < nDim; iDim++) + Velocity[iDim] = nodes->GetVelocity(iPoint,iDim); + Pressure = nodes->GetPressure(iPoint); + SoundSpeed2 = Gamma*Pressure/V_domain[nDim+2]; - /*--- Match the pressure, density and energy at the wall. ---*/ + /*--- Compute the acoustic Riemann invariant that is extrapolated + from the domain interior. ---*/ - Pressure = nodes->GetPressure(iPoint); - Density = Pressure / (Gas_Constant * Temperature); - Energy = Pressure / (Density*Gamma_Minus_One) + 0.5 * Velocity2; + Riemann = Two_Gamma_M1*sqrt(SoundSpeed2); + for (iDim = 0; iDim < nDim; iDim++) + Riemann += Velocity[iDim]*UnitNormal[iDim]; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- Speed of sound squared for fictitious inlet state ---*/ - V_inlet[0] = Temperature; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure / Density; - break; - } + SoundSpeed2 = Riemann; + for (iDim = 0; iDim < nDim; iDim++) + SoundSpeed2 -= Vel_Mag*Flow_Dir[iDim]*UnitNormal[iDim]; - case INLET_TYPE::MASS_FLOW: { + SoundSpeed2 = max(0.0,0.5*Gamma_Minus_One*SoundSpeed2); + SoundSpeed2 = SoundSpeed2*SoundSpeed2; - /*--- Impose the wall velocity from the interior wall node. ---*/ + /*--- Pressure for the fictitious inlet state ---*/ - Velocity2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); - Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1]; - } + Pressure = SoundSpeed2*Density/Gamma; - Pressure = nodes->GetPressure(iPoint); - Density = nodes->GetDensity(iPoint); - Temperature = Pressure / ( Gas_Constant * Density); - Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- Energy for the fictitious inlet state ---*/ - V_inlet[0] = Temperature; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - V_inlet[nDim+3] = Energy + Pressure / Density; - - break; - } + Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Vel_Mag*Vel_Mag; + if (tkeNeeded) Energy += GetTke_Inf(); - default: - SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); - break; + /*--- Primitive variables, using the derived quantities ---*/ + Temperature = Pressure / ( Gas_Constant * Density); + V_inlet[0] = Temperature; + for (iDim = 0; iDim < nDim; iDim++) + V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; + V_inlet[nDim+1] = Pressure; + V_inlet[nDim+2] = Density; + V_inlet[nDim+3] = Energy + Pressure/Density; - } + break; } - /*--- Set various quantities in the solver class ---*/ - - conv_numerics->SetPrimitive(V_domain, V_inlet); + default: + SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); + break; + } - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + /*--- Set various quantities in the solver class ---*/ - /*--- Compute the residual using an upwind scheme ---*/ + conv_numerics->SetPrimitive(V_domain, V_inlet); - auto residual = conv_numerics->ComputeResidual(config); + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Update residual value ---*/ + /*--- Compute the residual using an upwind scheme ---*/ - LinSysRes.AddBlock(iPoint, residual); + auto residual = conv_numerics->ComputeResidual(config); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Update residual value ---*/ - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + LinSysRes.AddBlock(iPoint, residual); -// /*--- Viscous contribution, commented out because serious convergence problems ---*/ -// -// if (viscous) { -// -// /*--- Set laminar and eddy viscosity at the infinity ---*/ -// -// V_inlet[nDim+5] = nodes->GetLaminarViscosity(iPoint); -// V_inlet[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_inlet); -// 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.jacobian_i); -// -// } + /*--- Jacobian contribution for implicit integration ---*/ - } + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); } END_SU2_OMP_FOR From 2de834c0fb644f433db04829fc169a4be728595d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 8 Dec 2024 14:13:19 -0800 Subject: [PATCH 20/20] update regressions --- .github/workflows/regression.yml | 6 +++--- TestCases/parallel_regression.py | 4 ++-- TestCases/serial_regression.py | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 1ee1b583e1d..2502204b64c 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -209,7 +209,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: @@ -255,7 +255,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: @@ -300,7 +300,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c fix_cornernode -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 44f11fe177b..5462560119e 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -311,7 +311,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 = [-6.660438, -1.186310, -0.956714, 0.000642, -0.000644, 0.000540, -0.001183, 596.990000, 300.060000, 296.920000, 22.236000, 0.524020, 37.284000, 2.349300] + flatplate_udobj.test_vals = [-6.660355, -1.186227, -0.956763, 0.000642, -0.000643, 0.000540, -0.001184, 596.990000, 300.060000, 296.920000, 22.235000, 0.524150, 37.285000, 2.350100] test_list.append(flatplate_udobj) # Laminar cylinder (steady) @@ -344,7 +344,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.007503, -7.227301, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.007512, -7.227061, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.492864, -7.671632, -0.000000, 2.085796] test_list.append(poiseuille_profile) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 9fc971df334..10312c92cc3 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -198,7 +198,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.009019, -7.262438, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.009021, -7.262796, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.494684, -7.711379, -0.000000, 2.085796] #last 4 columns test_list.append(poiseuille_profile)