diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 39ebcf28c011..9a85586cb336 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -870,6 +870,9 @@ class CConfig { Tke_FreeStream, /*!< \brief Total turbulent kinetic energy of the fluid. */ Intermittency_FreeStream, /*!< \brief Freestream intermittency (for sagt transition model) of the fluid. */ ReThetaT_FreeStream, /*!< \brief Freestream Transition Momentum Thickness Reynolds Number (for LM transition model) of the fluid. */ + TurbulenceIntensity_FreeStream, /*!< \brief Freestream turbulent intensity (for sagt transition model) of the fluid. */ + AmplificationFactor_FreeStream, /*!< \brief Freestream amplifictation factor for the eN 1 equation transition model. */ + Turb2LamViscRatio_FreeStream, /*!< \brief Ratio of turbulent to laminar viscosity. */ NuFactor_FreeStream, /*!< \brief Ratio of turbulent to laminar viscosity. */ NuFactor_Engine, /*!< \brief Ratio of turbulent to laminar viscosity at the engine. */ SecondaryFlow_ActDisk, /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */ @@ -1979,6 +1982,12 @@ class CConfig { * \return Freestream momentum thickness Reynolds number. */ su2double GetReThetaT_FreeStream() const { return ReThetaT_FreeStream; } + + /*! + * \brief Get the value of the non-dimensionalized freestream amplification factor. + * \return Non-dimensionalized freestream amplification factor. + */ + su2double GetAmplificationFactor_FreeStream(void) const { return AmplificationFactor_FreeStream; } /*! * \brief Get the value of the non-dimensionalized freestream turbulence intensity. diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 4e78f588bfa9..691055a6768b 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -107,6 +107,7 @@ class CPoint { ClosestWall_Marker; /*!< \brief Marker index of closest wall element, for given rank and zone index. */ su2vector ClosestWall_Elem; /*!< \brief Element index of closest wall element, for givenrank, zone and marker index. */ + su2activematrix Normals; /*!< \brief Normal of the nearest wall element. */ su2activevector SharpEdge_Distance; /*!< \brief Distance to a sharp edge. */ su2activevector Curvature; /*!< \brief Value of the surface curvature (SU2_GEO). */ @@ -508,6 +509,24 @@ class CPoint { */ inline su2double GetRoughnessHeight(unsigned long iPoint) const { return RoughnessHeight(iPoint); } + /*! + * \brief Set the value of the normal of the nearest wall element. + * \param[in] iPoint - Index of the point. + * \param[in] normal - Value of the normal. + */ + template + inline void SetNormal(unsigned long iPoint, Normals_type const&normal) { + for (unsigned long iDim = 0; iDim < nDim; iDim++) + Normals(iPoint,iDim) = normal[iDim]; + } + + /*! + * \brief Get the value of the normal of the nearest wall element. + * \param[in] iPoint - Index of the point. + * \return normal to the normal of the nearest wall element. + */ + inline su2double *GetNormal(unsigned long iPoint) { return Normals[iPoint]; } + /*! * \brief Set the value of the distance to a sharp edge. * \param[in] iPoint - Index of the point. @@ -908,4 +927,20 @@ class CPoint { } } } + + /*! + * \brief Set wall normal according to stored closest wall information. + * \param[in] normals - Mapping [rank][zone][marker][element] -> normal + */ + template + void SetWallNormals(Normals_type const&normals){ + for (unsigned long iPoint=0; iPoint= 0) + SetNormal(iPoint, normals[rankID][zoneID][markerID][elementID]); + } + } }; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 7e5aa1727def..84c482d7be9b 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1177,10 +1177,12 @@ inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned sh enum class TURB_TRANS_MODEL { NONE, /*!< \brief No transition model. */ LM, /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */ + EN, /*!< \brief Kind of transition model using Amplification factor for Spalart-Allmaras*/ }; static const MapType Trans_Model_Map = { MakePair("NONE", TURB_TRANS_MODEL::NONE) MakePair("LM", TURB_TRANS_MODEL::LM) + MakePair("EN", TURB_TRANS_MODEL::EN) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 33ef724f162f..c26d74036f37 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4823,7 +4823,17 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i for (int i=0; i<7; ++i) eng_cyl[i] /= 12.0; } - if(Turb_Fixed_Values && !OptionIsSet("TURB_FIXED_VALUES_DOMAIN")){ + if (Kind_Trans_Model == TURB_TRANS_MODEL::EN) { + if (saParsedOptions.bc){ + SU2_MPI::Error("Please select only one transition model for the SA turbulence model (Choose either eN or BCM).", CURRENT_FUNCTION); + } else if (!saParsedOptions.ft2) { + SU2_MPI::Error("Please select SA_OPTIONS = WITHFT2 when using SA-Ft2-eN transtion model.", CURRENT_FUNCTION); + } else if ((Kind_Regime == ENUM_REGIME::COMPRESSIBLE) && (Ref_NonDim == 0)) { + SU2_MPI::Error("Please select a non-dimensionalization option other than 'REF_DIMENSIONALIZATION = DIMENSIONAL' when using SA-Ft2-eN transition.", CURRENT_FUNCTION); + } + } + + if (Turb_Fixed_Values && !OptionIsSet("TURB_FIXED_VALUES_DOMAIN")){ SU2_MPI::Error("TURB_FIXED_VALUES activated, but no domain set with TURB_FIXED_VALUES_DOMAIN.", CURRENT_FUNCTION); } @@ -6115,6 +6125,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } switch (Kind_Trans_Model) { case TURB_TRANS_MODEL::NONE: break; + case TURB_TRANS_MODEL::EN: cout << "Low-turbulence Transition model: eN 1 equation model (2014)" << endl; break; case TURB_TRANS_MODEL::LM: { cout << "Transition model: Langtry and Menter's 4 equation model"; if (lmParsedOptions.LM2015) { @@ -6145,6 +6156,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { break; } } + cout << "Hybrid RANS/LES: "; switch (Kind_HybridRANSLES) { case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break; diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 9d526f45ebd3..d2e9120e1032 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -3972,5 +3972,79 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo } } } + + /*--- Function copied from Transition SLM branch ---*/ + su2vector>> WallNormal_container; + WallNormal_container.resize(nZone) = su2vector>(); + for (int iZone = 0; iZone < nZone; iZone++){ + const CConfig* config = config_container[iZone]; + const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0]; + WallNormal_container[iZone].resize(geometry->GetnMarker()); + for (auto iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++){ + if (config->GetViscous_Wall(iMarker)) { + WallNormal_container[iZone][iMarker].resize(geometry->GetnElem_Bound(iMarker), 3); + + for (auto iElem = 0u; iElem < geometry->GetnElem_Bound(iMarker); iElem++) { + su2vector NormalHere; + NormalHere.resize(3) = su2double(0.0); + + for (unsigned short iNode = 0; iNode < geometry->bound[iMarker][iElem]->GetnNodes(); iNode++) { + // Extract global coordinate of the node + unsigned long iPointHere = geometry->bound[iMarker][iElem]->GetNode(iNode); + long iVertexHere = geometry->nodes->GetVertex(iPointHere, iMarker); + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] += geometry->vertex[iMarker][iVertexHere]->GetNormal(iDim); + } + + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] /= geometry->bound[iMarker][iElem]->GetnNodes(); + + su2double NormalMag = 0.0; + for (auto iDim = 0u; iDim < 3; iDim++) + NormalMag += NormalHere[iDim]*NormalHere[iDim]; + NormalMag = sqrt(NormalMag); + + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] /= NormalMag; + + for (auto iDim = 0u; iDim < 3; iDim++) + WallNormal_container[iZone][iMarker](iElem, iDim) = NormalHere[iDim]; + + } + } else { + WallNormal_container[iZone][iMarker].resize(1, 3) = su2double(0.0); + } + } + } + + auto normal_i = + make_pair(nZone, [config_container,geometry_container,iInst,WallNormal_container](unsigned long iZone){ + const CConfig* config = config_container[iZone]; + const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0]; + const auto nMarker = geometry->GetnMarker(); + const auto WallNormal = WallNormal_container[iZone]; + + return make_pair( nMarker, [config,geometry,WallNormal](unsigned long iMarker){ + auto nElem_Bou = geometry->GetnElem_Bound(iMarker); + if (!config->GetViscous_Wall(iMarker)) nElem_Bou = 1; + + return make_pair(nElem_Bou, [WallNormal,iMarker](unsigned long iElem){ + const auto dimensions = 3; + + return make_pair(dimensions, [WallNormal,iMarker,iElem](unsigned short iDim){ + + return WallNormal[iMarker](iElem, iDim); + }); + }); + }); + }); + + NdFlattener<4>Normals_Local(normal_i); + NdFlattener<5> Normals_global(Nd_MPI_Environment(), Normals_Local); + + // use it to update roughnesses + for(int jZone=0; jZonenodes->SetWallNormals(Normals_global); + } } } diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 190bf7ad69a4..9d7feddd8f04 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -132,6 +132,8 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) { RoughnessHeight.resize(npoint) = su2double(0.0); SharpEdge_Distance.resize(npoint) = su2double(0.0); + + Normals.resize(npoint, 3) = su2double(0.0); } void CPoint::SetElems(const vector >& elemsMatrix) { Elem = CCompressedSparsePatternL(elemsMatrix); } diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 236c932b6306..b5d49d9d4b56 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -228,6 +228,7 @@ class CDriver : public CDriverBase { template void InstantiateTransitionNumerics(unsigned short nVar_Trans, int offset, const CConfig* config, const CSolver* trans_solver, CNumerics****& numerics) const; + /*! * \brief Helper to instantiate species transport numerics specialized for different flow solvers. */ diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index c2068ef311b5..664d03ead6da 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -80,6 +80,9 @@ class CNumerics { turb_ke_i, /*!< \brief Turbulent kinetic energy at point i. */ turb_ke_j; /*!< \brief Turbulent kinetic energy at point j. */ su2double + amplification_factor_i, /*!< \brief amplification factor at point i. */ + modified_intermittency_i; /*!< \brief modified intermittency at point i. */ + su2double intermittency_eff_i, /*!< \brief effective intermittency at point i. */ intermittency_i; /*!< \brief intermittency at point i. */ su2double @@ -707,6 +710,39 @@ class CNumerics { */ virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */}; + /*! + * \brief Get the Amplification factor for the e^N transition model. + */ + inline su2double GetAmplificationFactor() const { return amplification_factor_i; } + + /*! + * \brief Set the value of the amplification factor for the e^N transition model. + * \param[in] amplification_factor_i - Value of the amplification factor at point i. + */ + void SetAmplificationFactor(su2double val_amplification_factor_i) { + amplification_factor_i = val_amplification_factor_i; + }; + + /*! + * \brief Get the Modified Intermittency for the e^N transition model. + */ + inline su2double GetModifiedIntermittency() const { return modified_intermittency_i; } + + /*! + * \brief Set the value of the modified intermittency for the e^N transition model. + * \param[in] modified_intermittency_i - Value of the modified intermittency at point i. + */ + void SetModifiedIntermittency(su2double val_modified_intermittency_i) { + modified_intermittency_i = val_modified_intermittency_i; + }; + + /*! + * \brief Set the gradient of the auxiliary variables. + * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i. + * \param[in] val_auxvar_grad_j - Gradient of the auxiliary variable at point j. + */ + inline virtual void SetAuxVar(su2double val_AuxVar) {} + /*! * \brief Set the value of the effective intermittency for the LM model. * \param[in] intermittency_eff_i - Value of the effective intermittency at point i. diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp index b9328d43f639..c314e3bb940e 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp @@ -2,7 +2,7 @@ * \file trans_convection.hpp * \brief Delarations of numerics classes for discretization of * convective fluxes in transition problems. - * \author S. Kang + * \author S. Kang, R. Roos * \version 7.5.1 "Blackbird" * * SU2 Project Website: https://su2code.github.io @@ -27,8 +27,8 @@ */ #pragma once - #include "../turb_convection.hpp" +#include "../../scalar/scalar_convection.hpp" /*! * \class CUpwSca_TransLM @@ -38,3 +38,14 @@ template using CUpwSca_TransLM = CUpwSca_TurbSST; +/*! + * \class CUpwSca_TransLM + * \brief Re-use the SA convective fluxes for the scalar upwind discretization of eN transition model equations. + * \ingroup ConvDiscr + */ +template +using CUpwSca_TransEN = CUpwSca_TurbSST; + +//template +//using CUpwSca_TransEN = CUpwSca_TurbSA; + diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp index 253e35c4dc4f..4d16ebf06a2e 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp @@ -109,7 +109,7 @@ class TransLMCorrelations { break; } case TURB_TRANS_CORRELATION::DEFAULT: - SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has been set in the code.", CURRENT_FUNCTION); break; } @@ -182,7 +182,7 @@ class TransLMCorrelations { break; } case TURB_TRANS_CORRELATION::DEFAULT: - SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has been set in the code.", CURRENT_FUNCTION); break; } diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp index d2155a82e937..c0c50c5442f0 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp @@ -2,7 +2,7 @@ * \file trans_diffusion.hpp * \brief Declarations of numerics classes for discretization of * viscous fluxes in transition problems. - * \author S. Kang + * \author S. Kang, R. Roos * \version 7.5.1 "Blackbird" * * SU2 Project Website: https://su2code.github.io @@ -103,3 +103,101 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar { } }; + +/*! + * \class CAvgGrad_TransEN + * \brief Class for computing viscous term using average of gradient with correction (e^N transition model). + * \ingroup ViscDiscr + * \author R. Roos + */ +template +class CAvgGrad_TransEN final : public CAvgGrad_Scalar { +private: + using Base = CAvgGrad_Scalar; + using Base::Laminar_Viscosity_i; + using Base::Laminar_Viscosity_j; + using Base::Eddy_Viscosity_i; + using Base::Eddy_Viscosity_j; + using Base::Density_i; + using Base::Density_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::Proj_Mean_GradScalarVar; + using Base::proj_vector_ij; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + + const su2double sigma_n = 1.0; + const su2double sigma_y = 1.0; + + /*! + * \brief Adds any extra variables to AD + */ + void ExtraADPreaccIn() override {} + + /*! + * \brief SA specific steps in the ComputeResidual method + * \param[in] config - Definition of the particular problem. + */ + + void FinishResidualCalc(const CConfig* config) override { + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + + /*--- Compute mean effective dynamic viscosity ---*/ + const su2double diff_i_n = sigma_n*(Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double diff_j_n = sigma_n*(Laminar_Viscosity_j + Eddy_Viscosity_j); + const su2double diff_i_gamma = Laminar_Viscosity_i + (Eddy_Viscosity_i/sigma_y); + const su2double diff_j_gamma = Laminar_Viscosity_j + (Eddy_Viscosity_j/sigma_y); + + const su2double diff_n = 0.5*(diff_i_n + diff_j_n); + const su2double diff_gamma = 0.5*(diff_i_gamma + diff_j_gamma); + + Flux[0] = diff_n*Proj_Mean_GradScalarVar[0]; + Flux[1] = diff_gamma*Proj_Mean_GradScalarVar[1]; + + /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ + if (implicit) { + const su2double proj_on_rho_i = proj_vector_ij/Density_i; + Jacobian_i[0][0] = -diff_n*proj_on_rho_i; Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_gamma*proj_on_rho_i; + + const su2double proj_on_rho_j = proj_vector_ij/Density_j; + Jacobian_j[0][0] = diff_n*proj_on_rho_j; Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_gamma*proj_on_rho_j; + } + } + +// void FinishResidualCalc(const CConfig* config) override { +// const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; +// +// /*--- Compute mean effective dynamic viscosity ---*/ +// const su2double diff_i_amplification = (Laminar_Viscosity_i + Eddy_Viscosity_i)/sigma_n; +// const su2double diff_j_amplification = (Laminar_Viscosity_j + Eddy_Viscosity_j)/sigma_n; +// +// const su2double diff_amplification = 0.5*(diff_i_amplification + diff_j_amplification); +// +// Flux[0] = diff_amplification*Proj_Mean_GradScalarVar[0]; +// +// /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ +// +// if (implicit) { +// Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-diff_amplification*proj_vector_ij); +// Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+diff_amplification*proj_vector_ij); +// } +// } + +public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] correct_grad - Whether to correct gradient for skewness. + * \param[in] config - Definition of the particular problem. + */ + CAvgGrad_TransEN(unsigned short val_nDim, unsigned short val_nVar, + bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) { + } + +}; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index 3e47e8172e12..28d44d527f76 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -323,3 +323,156 @@ class CSourcePieceWise_TransLM final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } }; + +/*! + * \class CSourcePieceWise_TranEN + * \brief Class for integrating the source terms of the e^N transition model equations. + * \ingroup SourceDiscr + * \author R. Roos + */ +template +class CSourcePieceWise_TransEN final : public CNumerics { + private: + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + + /*--- eN Closure constants ---*/ + su2double c_1 = 100.0; + su2double c_2 = 0.06; + su2double c_3 = 50.0; + + su2double AuxVar; + + su2double Residual[2]; + su2double* Jacobian_i[2]; + su2double Jacobian_Buffer[4]; // Static storage for the Jacobian (which needs to be pointer for return type). + + public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourcePieceWise_TransEN(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CNumerics(val_nDim, 1, config), + idx(val_nDim, config->GetnSpecies()) { + + /*--- "Allocate" the Jacobian using the static buffer. ---*/ + Jacobian_i[0] = Jacobian_Buffer; + Jacobian_i[1] = Jacobian_Buffer + 2; + } + + /*! + * \brief Residual for source term integration. + * \param[in] config - Definition of the particular problem. + * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. + */ + ResidualType<> ComputeResidual(const CConfig* config) override { + + AD::StartPreacc(); + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.Pressure()], V_i[idx.LaminarViscosity()], StrainMag_i, ScalarVar_i[0], Volume, dist_i); + AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); + AD::SetPreaccIn(Vorticity_i, 3); + AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim); + AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); + + su2double rho = V_i[idx.Density()]; + su2double p = V_i[idx.Pressure()]; + su2double muLam = V_i[idx.LaminarViscosity()]; + su2double muEddy = V_i[idx.EddyViscosity()]; + + const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity_i); + + const su2double vel_u = V_i[idx.Velocity()]; + const su2double vel_v = V_i[1+idx.Velocity()]; + const su2double vel_w = (nDim ==3) ? V_i[2+idx.Velocity()] : 0.0; + + const su2double vel_mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20); + + su2double rhoInf = config->GetDensity_FreeStreamND(); + su2double pInf = config->GetPressure_FreeStreamND(); + const su2double *velInf = config->GetVelocity_FreeStreamND(); + + su2double velInf2 = 0.0; + for(unsigned short iDim = 0; iDim < nDim; ++iDim) { + velInf2 += velInf[iDim]*velInf[iDim]; + } + + Residual[0] = 0.0; + Residual[1] = 0.0; + Jacobian_i[0][0] = 0.0; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = 0.0; + + if (dist_i > 1e-10) { + + /*--- Local pressure-gradient parameter for the boundary layer shape factor. Minimum value of 0.328 for stability ---*/ + const su2double H_L = (pow(dist_i,2)/(muLam/rho))*AuxVar; + + /*--- Integral shape factor ---*/ + const su2double H_12 = min(max(0.26*H_L + 2.4, 2.2), 20.0); + + /*--- F growth parameters ---*/ + const su2double DH_12 = 2.4*H_12 / (H_12 - 1.0); + const su2double lH_12 = (6.54*H_12 - 14.07)/pow(H_12,2); + const su2double mH_12 = (1/lH_12) * (0.058*( pow((H_12 - 4.0),2)/(H_12 - 1) ) - 0.068); + + const su2double F_growth = DH_12*( (1 + mH_12)/2 )* lH_12; + + /*--- F crit parameters ---*/ + const su2double Re_v = (rho*StrainMag_i*pow(dist_i,2))/(muLam + muEddy); + const su2double k_v = 1/(0.4036*pow(H_12,2)-2.5394*H_12+4.3273); + const su2double Re_d2_0 = pow(10,(0.7*tanh((14/(H_12 - 1)) - 9.24) + 2.492/(pow((H_12 - 1),0.43)) + 0.62)); + const su2double Re_v_0 = k_v * Re_d2_0; + + short int F_crit; + if (Re_v < Re_v_0){ + F_crit = 0; + } else { + F_crit = 1; + } + + /*--- Source term expresses stream wise growth of Tollmien_schlichting instabilities ---*/ + const su2double dn_over_dRe_d2 = 0.028*(H_12 - 1) - 0.0345*exp( -(pow( (3.87/(H_12 - 1) - 2.52),2) ) ); + + /*--- F onset parameters ---*/ + const su2double Ncrit = -8.43 - 2.4*log( 2.5*tanh(config->GetTurbulenceIntensity_FreeStream()/2.5) /100); + const su2double F_onset1 = TransVar_i[0]/(Ncrit); + + const su2double F_onset2 = min(F_onset1,2.0); + + const su2double Rt = muEddy/muLam; + const su2double F_onset3 = max(1-( pow((Rt/3.5),3) ),0.0); + + const su2double Fonset = max((F_onset2 - F_onset3), 0.0); + + /*--- F turb parameters ---*/ + const su2double Fturb = exp(- (pow((Rt/2),4))); + + /*--- Production terms ---*/ + const su2double P_amplification = rho*VorticityMag*F_crit*F_growth*dn_over_dRe_d2; + const su2double P_gamma = c_1*rho*StrainMag_i*Fonset*( 1-exp(TransVar_i[1]) ); + const su2double D_gamma = c_2*rho*VorticityMag*Fturb*( c_3*exp(TransVar_i[1]) -1 ); + + /*--- Add Production to residual ---*/ + Residual[0] += P_amplification * Volume; + Residual[1] += (P_gamma - D_gamma) * Volume; + + /*--- Implicit part ---*/ + Jacobian_i[0][0] = (rho*VorticityMag*F_crit*F_growth) * Volume; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = ( -c_1*rho*StrainMag_i*Fonset*exp(TransVar_i[1]) - c_2*rho*VorticityMag*Fturb* c_3*exp(TransVar_i[1]) ) * Volume; + + } + + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(Residual, Jacobian_i, nullptr); + } + + inline void SetAuxVar(su2double val_AuxVar) override { AuxVar = val_AuxVar;} + +}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 541383e6dbd4..d4b070842cc9 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -41,7 +41,6 @@ struct CSAVariables { const su2double cb1 = 0.1355; const su2double cw2 = 0.3; const su2double ct3 = 1.2; - const su2double ct4 = 0.5; const su2double cw3_6 = pow(2, 6); const su2double sigma = 2.0 / 3.0; const su2double cb2 = 0.622; @@ -51,13 +50,20 @@ struct CSAVariables { const su2double CRot = 1.0; const su2double c2 = 0.7, c3 = 0.9; + /*--- List of non-const constants ---*/ + su2double ct4 = 0.5; + su2double amplification, modifiedintermittency; + /*--- List of auxiliary functions ---*/ - su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2; + su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2, Ncrit; /*--- List of helpers ---*/ su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad; su2double intermittency, interDestrFactor; + + /*--- List of booleans ---*/ + bool transEN = false; }; /*! @@ -155,7 +161,18 @@ class CSourceBase_TurbSA : public CNumerics { var.fv2 = 1 - ScalarVar_i[0] / (nu + ScalarVar_i[0] * var.fv1); var.d_fv2 = -(1 / nu - Ji_2 * var.d_fv1) / pow(1 + var.Ji * var.fv1, 2); - /*--- Compute ft2 term ---*/ + /*--- Compute ft2 term. Also includes boolean for e^N transition model that modifies the ft2 term ---*/ + if(TURB_TRANS_MODEL::EN == config->GetKind_Trans_Model()) { + var.transEN = true; +// var.Ncrit = -8.43 - 2.4*log(config->GetTurbulenceIntensity_FreeStream()/100); +// var.amplification = min(amplification_factor_i, var.Ncrit); + var.modifiedintermittency = modified_intermittency_i; + + /*--- Slight deviation from theory to obtain better results. Coder et al: Ct4 = 0.05 ---*/ +// if (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) var.ct4 = 0.025; +// else var.ct4 = 0.05; + } + ft2::get(var); /*--- Compute modified vorticity ---*/ @@ -218,6 +235,7 @@ class CSourceBase_TurbSA : public CNumerics { Residual = (Production - Destruction + CrossProduction) * Volume; Jacobian_i[0] *= Volume; + } AD::SetPreaccOut(Residual); @@ -287,9 +305,17 @@ struct Zero { /*! \brief Non-zero ft2 term according to the literature. */ struct Nonzero { static void get(CSAVariables& var) { - const su2double xsi2 = pow(var.Ji, 2); - var.ft2 = var.ct3 * exp(-var.ct4 * xsi2); - var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; + const su2double xsi2 = pow(var.Ji, 2); + + if (var.transEN == true){ +// var.ft2 = var.ct3 * (1 - exp(2*(var.amplification - var.Ncrit)) ) * exp(-var.ct4 * xsi2); + var.ft2 = var.ct3 * (1 - exp(var.modifiedintermittency) ); + var.d_ft2 = 0.0; +// var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; + } else { + var.ft2 = var.ct3 * exp(-var.ct4 * xsi2); + var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; + } } }; }; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 06cd950aa138..dbcfe1c88a1f 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3157,6 +3157,12 @@ class CSolver { * \return Value of the turbulent frequency. */ inline virtual su2double GetOmega_Inf(void) const { return 0; } + + /*! + * \brief Get value of the Amplification Factor. + * \return Value of the Amplification Factor. + */ + inline virtual su2double GetAmplificationFactor_Inf() const { return 0; } /*! * \brief Get value of the Intermittency. diff --git a/SU2_CFD/include/solvers/CSolverFactory.hpp b/SU2_CFD/include/solvers/CSolverFactory.hpp index bdf8300619f1..26e66acd4764 100644 --- a/SU2_CFD/include/solvers/CSolverFactory.hpp +++ b/SU2_CFD/include/solvers/CSolverFactory.hpp @@ -105,15 +105,15 @@ class CSolverFactory { /*! * \brief Create a transition solver - * \param[in] kindTransModel - Kind of transition solver + * \param[in] kindTransModel- Kind of transition solver * \param[in] solver - The solver container (used to call preprocessing of the flow solver) * \param[in] geometry - The geometry definition * \param[in] config - The configuration * \param[in] iMGLevel - The multigrid level * \param[in] adjoint - Boolean indicating whether a primal or adjoint solver should be allocated - * \return - A pointer to the allocated transition solver + * \return - A pointer to the allocated turbulent solver */ - static CSolver* CreateTransSolver(TURB_TRANS_MODEL kindTransModel , CSolver **solver, CGeometry *geometry, CConfig *config, int iMGLevel, int adjoint); + static CSolver* CreateTransSolver(TURB_TRANS_MODEL kindTransModel, CSolver **solver, CGeometry *geometry, CConfig *config, int iMGLevel, int adjoint); /*! * \brief Create a species solver diff --git a/SU2_CFD/include/solvers/CTransENSolver.hpp b/SU2_CFD/include/solvers/CTransENSolver.hpp new file mode 100644 index 000000000000..12968d246574 --- /dev/null +++ b/SU2_CFD/include/solvers/CTransENSolver.hpp @@ -0,0 +1,201 @@ +/*! + * \file CTransENSolver.hpp + * \brief Headers of the CTransENSolver class + * \author R. Roos + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CTurbSolver.hpp" + +/*! + * \class CTransENSolver + * \brief Main class for defining the e^N transition model solver. + * \ingroup Turbulence_Model + * \author Roos + */ + +class CTransENSolver final : public CTurbSolver { +private: + +public: + /*! + * \overload + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + CTransENSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh); + + /*! + * \brief Restart residual and compute gradients. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void Preprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned short iRKStep, + unsigned short RunTime_EqSystem, + bool Output) override; + + /*! + * \brief Computes the effective intermtittency. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Postprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Compute the viscous flux for the LM equation at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \note Calls a generic implementation after defining a SolverSpecificNumerics object. + */ + void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, CConfig* config) override; + + /*! + * \brief Source term computation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics_container - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Residual(CGeometry *geometry, + CSolver **solver_container, + CNumerics **numerics_container, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Source term computation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Template(CGeometry *geometry, + CSolver **solver_container, + CNumerics *numerics, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Impose the Langtry Menter transition 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_HeatFlux_Wall(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the Navier-Stokes 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 Impose the inlet 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_Inlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the outlet 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_Outlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Get the value of the amplification factor. + * \return Value of the amplification factor. + */ + inline su2double GetAmplificationFactor_Inf(void) const override { return Solution_Inf[0]; } + + /*! + * \brief Load a solution from a restart file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all of the solvers. + * \param[in] config - Definition of the particular problem. + * \param[in] val_iter - Current external iteration number. + * \param[in] val_update_geo - Flag for updating coords and grid velocity. + */ + void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) final; + +}; diff --git a/SU2_CFD/include/variables/CTransENVariable.hpp b/SU2_CFD/include/variables/CTransENVariable.hpp new file mode 100644 index 000000000000..82e492b51912 --- /dev/null +++ b/SU2_CFD/include/variables/CTransENVariable.hpp @@ -0,0 +1,91 @@ +/*! + * \file CTransENVariable.hpp + * \brief Declaration of the variables of the e^N transition model. + * \author F. Palacios, T. Economon + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CTurbVariable.hpp" + +/*! + * \class CTransENVariable + * \brief e^N Transition model variables. + * \ingroup Turbulence_Model + * \author R. Roos + */ + +class CTransENVariable final : public CTurbVariable { +protected: + VectorType AmplificationFactor; + VectorType ModifiedIntermittency; + + VectorType normal_x; + VectorType normal_y; + VectorType normal_z; + +public: + /*! + * \brief Constructor of the class. + * \param[in] AmplificationFactor - Amplification factor(n) (initialization value). + * \param[in] npoint - Number of points/nodes/vertices in the domain. + * \param[in] ndim - Number of dimensions of the problem. + * \param[in] nvar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CTransENVariable(su2double AmplificationFactor,su2double ModifiedIntermittency, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CTransENVariable() override = default; + + /*! + * \brief Set Amplification Factor. + */ + void SetAmplificationFactor(unsigned long iPoint, su2double val_AmplificationFactor) override ; + + /*! + * \brief Set Amplification Factor. + */ + void SetModifiedIntermittency(unsigned long iPoint, su2double val_Gamma) override ; + + void SetNormal(unsigned long iPoint, su2double val_normal_x, su2double val_normal_y, su2double val_normal_z) override; + + /*! + * \brief Value of AmplificationFactor. + */ + inline su2double GetAmplificationFactor(unsigned long iPoint) const override { return AmplificationFactor(iPoint); } + + /*! + * \brief Value of AmplificationFactor. + */ + inline su2double GetModifiedIntermittency(unsigned long iPoint) const override { return ModifiedIntermittency(iPoint); } + + + inline su2double GetNormal_x(unsigned long iPoint) const override {return normal_x(iPoint);}; + inline su2double GetNormal_y(unsigned long iPoint) const override {return normal_y(iPoint);}; + inline su2double GetNormal_z(unsigned long iPoint) const override {return normal_z(iPoint);}; + +}; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index ec2557c8088a..0a10a26a5233 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1129,6 +1129,55 @@ class CVariable { */ inline virtual void SetGammaEff(unsigned long iPoint) {} + /*! + * \brief Set the amplification factor for the eN transition model. + * \param[in] val_dist - Value of the amplification factor. + */ + inline virtual void SetAmplificationFactor(unsigned long iPoint, su2double val_AmplificationFactor) {} + + /*! + * \brief Set the modified intermittency for the eN transition model. + * \param[in] val_dist - Value of the amplification factor. + */ + inline virtual void SetModifiedIntermittency(unsigned long iPoint, su2double val_Gamma) {} + + inline virtual void SetNormal(unsigned long iPoint, su2double val_normal_x, su2double val_normal_y, su2double val_normal_z) {}; + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \return Returns Amplification factor for the eN transition model + */ + inline virtual su2double GetAmplificationFactor(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \return Returns Modified Intermittency for the eN transition model + */ + inline virtual su2double GetModifiedIntermittency(unsigned long iPoint) const { return 0.0; } + + inline virtual su2double GetNormal_x(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetNormal_y(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetNormal_z(unsigned long iPoint) const { return 0.0; } + + /*--- Debug terms ---*/ + inline virtual su2double GetProdN(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetProdG(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestG(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetGammaN(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetHL(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetH12(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetFG(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetFC(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetREV(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetREV0(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDist(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetStrain(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetFonset1(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetFonset(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetFturb(unsigned long iPoint) const { return 0.0; } + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index c8ead1837303..b56e9e8dc5e3 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -165,7 +165,8 @@ CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false), if (rank == MASTER_NODE) cout << "Computing wall distances." << endl; - CGeometry::ComputeWallDistance(config_container, geometry_container); + if (!dry_run) + CGeometry::ComputeWallDistance(config_container, geometry_container); for (iZone = 0; iZone < nZone; iZone++) { @@ -1330,7 +1331,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse const int conv_bound_term = CONV_BOUND_TERM + offset; const int visc_bound_term = VISC_BOUND_TERM + offset; + /*--- Assign transition model booleans ---*/ + const bool LM = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM; + const bool EN = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::EN; /*--- Definition of the convective scheme for each equation and mesh level ---*/ @@ -1340,7 +1344,12 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse break; case SPACE_UPWIND : for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + if (LM) { + numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + } + else if (EN) { + numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransEN(nDim, nVar_Trans, config); + } } break; default: @@ -1351,15 +1360,25 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (LM) { + numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + } + else if (EN) { + numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransEN(nDim, nVar_Trans, true, config); + } } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { auto& trans_source_first_term = numerics[iMGlevel][TRANS_SOL][source_first_term]; - - if (LM) trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); + + if (LM) { + trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); + } + else if (EN) { + trans_source_first_term = new CSourcePieceWise_TransEN(nDim, nVar_Trans, config); + } numerics[iMGlevel][TRANS_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Trans, config); } @@ -1370,9 +1389,14 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse if (LM) { numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, false, config); + } + else if (EN) { + numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransEN(nDim, nVar_Trans, config); + numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransEN(nDim, nVar_Trans, false, config); } } } + /*--- Explicit instantiation of the template above, needed because it is defined in a cpp file, instead of hpp. ---*/ template void CDriver::InstantiateTransitionNumerics>( unsigned short, int, const CConfig*, const CSolver*, CNumerics****&) const; @@ -1498,7 +1522,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver case MAIN_SOLVER::RANS: case MAIN_SOLVER::DISC_ADJ_RANS: ns = compressible = turbulent = true; - transition = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM); + transition = (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE); species = config->GetKind_Species_Model() != SPECIES_MODEL::NONE; break; case MAIN_SOLVER::INC_EULER: @@ -1515,7 +1539,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver case MAIN_SOLVER::DISC_ADJ_INC_RANS: ns = incompressible = turbulent = true; heat = config->GetWeakly_Coupled_Heat(); - transition = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM); + transition = (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE); species = (config->GetKind_Species_Model() != SPECIES_MODEL::NONE); break; case MAIN_SOLVER::FEM_EULER: @@ -2067,7 +2091,6 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } /*--- Solver definition for the species transport problem ---*/ - if (species) { if (incompressible) InstantiateSpeciesNumerics >(nVar_Species, offset, config, diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index cde93a085bbb..ba9e3925f359 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -85,14 +85,14 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe /*--- Solve transition model ---*/ - if (config[val_iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { + if (config[val_iZone]->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TRANS_SYS); integration[val_iZone][val_iInst][TRANS_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TRANS_SYS, val_iZone, val_iInst); } - /*--- Solve the turbulence model ---*/ - + /*--- Solve the turbulence model ---*/ + config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TURB_SYS); integration[val_iZone][val_iInst][TURB_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TURB_SYS, val_iZone, val_iInst); @@ -187,7 +187,7 @@ void CFluidIteration::Update(COutput* output, CIntegration**** integration, CGeo /*--- Update dual time solver for the transition model ---*/ - if (config[val_iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { + if (config[val_iZone]->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { integration[val_iZone][val_iInst][TRANS_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0][TRANS_SOL], config[val_iZone], MESH_0); @@ -283,8 +283,12 @@ void CFluidIteration::Solve(COutput* output, CIntegration**** integration, CGeom Output(output, geometry, solver, config, Inner_Iter, StopCalc, val_iZone, val_iInst); } +// if (Inner_Iter == 9) +// SU2_MPI::Error("10 Iterations done", CURRENT_FUNCTION); + /*--- If the iteration has converged, break the loop ---*/ if (StopCalc) break; + } if (multizone && steady) { diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 946227cd3052..ebcd0df1feef 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -75,6 +75,7 @@ su2_cfd_src += files(['variables/CIncNSVariable.cpp', 'variables/CRadVariable.cpp', 'variables/CRadP1Variable.cpp', 'variables/CTurbSSTVariable.cpp', + 'variables/CTransENVariable.cpp', 'variables/CVariable.cpp', 'variables/CNSVariable.cpp', 'variables/CAdjTurbVariable.cpp', @@ -112,6 +113,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CSpeciesSolver.cpp', 'solvers/CSpeciesFlameletSolver.cpp', 'solvers/CTransLMSolver.cpp', + 'solvers/CTransENSolver.cpp', 'solvers/CTurbSolver.cpp', 'solvers/CTurbSASolver.cpp', 'solvers/CTurbSSTSolver.cpp', diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 20b41e7fc79c..57633588861a 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -981,6 +981,13 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); break; + + case TURB_TRANS_MODEL::EN: + /// DESCRIPTION: Root-mean square residual of the Amplification factor (e^N model). + AddHistoryOutput("RMS_AMPLIFICATION_FACTOR", "rms[n]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the amplification factor (EN model)).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Root-mean square residual of the Modified Intermittency (e^N model). + AddHistoryOutput("RMS_MODIFIED_INTERMITTENCY", "rms[g]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the modified intermittency (EN model)).", HistoryFieldType::RESIDUAL); + break; case TURB_TRANS_MODEL::NONE: break; } @@ -1025,12 +1032,18 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { } switch (config->GetKind_Trans_Model()) { - case TURB_TRANS_MODEL::LM: /// DESCRIPTION: Maximum residual of the intermittency (LM model). AddHistoryOutput("MAX_INTERMITTENCY", "max[LM_1]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the intermittency (LM model).", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + break; + + case TURB_TRANS_MODEL::EN: + /// DESCRIPTION: Maximum residual of the amplification factor (EN model). + AddHistoryOutput("MAX_AMPLIFICATION_FACTOR", "max[n]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the amplification factor (EN model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the modified intermittency (EN model). + AddHistoryOutput("MAX_MODIFIED_INTERMITTENCY", "max[g]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the modified intermittency (EN model).", HistoryFieldType::RESIDUAL); break; case TURB_TRANS_MODEL::NONE: @@ -1084,6 +1097,13 @@ void CFlowOutput::AddHistoryOutputFields_ScalarBGS_RES(const CConfig* config) { /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); break; + + case TURB_TRANS_MODEL::EN: + /// DESCRIPTION: Maximum residual of the amplification factor (EN model). + AddHistoryOutput("BGS_AMPLIFICATION_FACTOR", "bgs[n]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the amplification factor (EN model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the modified intermittency (EN model). + AddHistoryOutput("BGS_MODIFIED_INTERMITTENCY", "bgs[g]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the modified intermittency (EN model).", HistoryFieldType::RESIDUAL); + break; case TURB_TRANS_MODEL::NONE: break; } @@ -1178,6 +1198,17 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co SetHistoryOutputValue("LINSOL_ITER_TRANS", solver[TRANS_SOL]->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL_TRANS", log10(solver[TRANS_SOL]->GetResLinSolver())); break; + + case TURB_TRANS_MODEL::EN: + SetHistoryOutputValue("RMS_AMPLIFICATION_FACTOR", log10(solver[TRANS_SOL]->GetRes_RMS(0))); + SetHistoryOutputValue("MAX_AMPLIFICATION_FACTOR", log10(solver[TRANS_SOL]->GetRes_Max(0))); + SetHistoryOutputValue("RMS_MODIFIED_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("MAX_MODIFIED_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_Max(1))); + if (multiZone) { + SetHistoryOutputValue("BGS_AMPLIFICATION_FACTOR", log10(solver[TRANS_SOL]->GetRes_BGS(0))); + SetHistoryOutputValue("BGS_MODIFIED_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + } + break; case TURB_TRANS_MODEL::NONE: break; } @@ -1242,6 +1273,12 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("INTERMITTENCY", "LM_gamma", "SOLUTION", "LM intermittency"); AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); break; + + case TURB_TRANS_MODEL::EN: + AddVolumeOutput("AMPLIFICATION_FACTOR", "n", "SOLUTION", "e^N Amplification Factor"); + AddVolumeOutput("MODIFIED_INTERMITTENCY", "eN_gamma", "SOLUTION", "e^N Modified Intermittency"); + AddVolumeOutput("TURB_INDEX", "Turb_index", "PRIMITIVE", "Turbulence index"); + break; case TURB_TRANS_MODEL::NONE: break; @@ -1309,6 +1346,11 @@ void CFlowOutput::SetVolumeOutputFieldsScalarResidual(const CConfig* config) { AddVolumeOutput("RES_INTERMITTENCY", "Residual_LM_intermittency", "RESIDUAL", "Residual of LM intermittency"); AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); break; + + case TURB_TRANS_MODEL::EN: + AddVolumeOutput("RES_AMPLIFICATION_FACTOR", "Residual_EN_AMLIFICATION", "RESIDUAL", "Residual of the amplification factor (EN model)"); + AddVolumeOutput("RES_MODIFIED_INTERMITTENCY", "Residual_EN_INTERMITTENCY", "RESIDUAL", "Residual of the modified intermittency (EN model)"); + break; case TURB_TRANS_MODEL::NONE: break; @@ -1531,6 +1573,14 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("RES_INTERMITTENCY", iPoint, trans_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); break; + + case TURB_TRANS_MODEL::EN: + SetVolumeOutputValue("AMPLIFICATION_FACTOR", iPoint, Node_Trans->GetSolution(iPoint, 0)); + SetVolumeOutputValue("RES_AMPLIFICATION_FACTOR", iPoint, trans_solver->LinSysRes(iPoint, 0)); + SetVolumeOutputValue("MODIFIED_INTERMITTENCY", iPoint, Node_Trans->GetSolution(iPoint, 1)); + SetVolumeOutputValue("RES_MODIFIED_INTERMITTENCY", iPoint, trans_solver->LinSysRes(iPoint, 1)); + SetVolumeOutputValue("TURB_INDEX", iPoint, Node_Turb->GetTurbIndex(iPoint)); + break; case TURB_TRANS_MODEL::NONE: break; } @@ -2628,15 +2678,18 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo if (transition) { file << "Transition model: "; switch (Kind_Trans_Model) { - case TURB_TRANS_MODEL::NONE: break; - case TURB_TRANS_MODEL::LM: - file << "Langtry and Menter's transition"; - if (config->GetLMParsedOptions().LM2015) { - file << " w/ cross-flow corrections (2015)\n"; - } else { - file << " (2009)\n"; - } - break; + case TURB_TRANS_MODEL::NONE: break; + case TURB_TRANS_MODEL::LM: + file << "Langtry and Menter's transition"; + if (config->GetLMParsedOptions().LM2015) { + file << " w/ cross-flow corrections (2015)\n"; + } else { + file << " (2009)\n"; + } + break; + case TURB_TRANS_MODEL::EN: + file << "e^N transition\n"; + break; } } break; diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index af9669fa3920..6a67122c8eab 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -36,6 +36,7 @@ #include "../../include/solvers/CTurbSASolver.hpp" #include "../../include/solvers/CTurbSSTSolver.hpp" #include "../../include/solvers/CTransLMSolver.hpp" +#include "../../include/solvers/CTransENSolver.hpp" #include "../../include/solvers/CAdjEulerSolver.hpp" #include "../../include/solvers/CAdjNSSolver.hpp" #include "../../include/solvers/CAdjTurbSolver.hpp" @@ -64,48 +65,48 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon switch (kindMainSolver) { case MAIN_SOLVER::TEMPLATE_SOLVER: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TEMPLATE, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TEMPLATE, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::INC_EULER: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_EULER, solver, geometry, config, iMGLevel); - solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_EULER, solver, geometry, config, iMGLevel); + solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::EULER: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::EULER, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::EULER, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::NEMO_EULER: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NEMO_EULER, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NEMO_EULER, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::INC_NAVIER_STOKES: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); - solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); - solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::NAVIER_STOKES: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::NEMO_NAVIER_STOKES: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NEMO_NAVIER_STOKES, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NEMO_NAVIER_STOKES, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::RANS: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel); - solver[TURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel); - solver[TRANS_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TRANSITION, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel); + solver[TURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel); + solver[TRANS_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TRANSITION, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::INC_RANS: - solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); - solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); - solver[TURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel); - solver[TRANS_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TRANSITION, solver, geometry, config, iMGLevel); - solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); + solver[TURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel); + solver[TRANS_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TRANSITION, solver, geometry, config, iMGLevel); + solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::HEAT_EQUATION: - solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::ADJ_EULER: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::EULER, solver, geometry, config, iMGLevel); @@ -382,11 +383,19 @@ CSolver* CSolverFactory::CreateTransSolver(TURB_TRANS_MODEL kindTransModel, CSol solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); transSolver->Postprocessing(geometry, solver, config, iMGLevel); solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + break; + + case TURB_TRANS_MODEL::EN : + transSolver = new CTransENSolver(geometry, config, iMGLevel); + solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + transSolver->Postprocessing(geometry, solver, config, iMGLevel); + solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); break; - case TURB_TRANS_MODEL::NONE: + + case TURB_TRANS_MODEL::NONE: break; } - } + } return transSolver; } diff --git a/SU2_CFD/src/solvers/CTransENSolver.cpp b/SU2_CFD/src/solvers/CTransENSolver.cpp new file mode 100644 index 000000000000..c1c1758d2556 --- /dev/null +++ b/SU2_CFD/src/solvers/CTransENSolver.cpp @@ -0,0 +1,549 @@ +/*! + * \file CTransENSolver.cpp + * \brief Main subroutines for e^N Transition model solver. + * \author R. Roos + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/solvers/CTransENSolver.hpp" +#include "../../include/variables/CTransENVariable.hpp" +#include "../../include/variables/CFlowVariable.hpp" +#include "../../include/variables/CTurbSAVariable.hpp" +#include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +/*--- This is the implementation of the eN transition model. + The main reference for this model is: Coder, J.G., https://doi.org/10.2514/6.2019-0039. ---*/ + +CTransENSolver::CTransENSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) + : CTurbSolver(geometry, config, true) { + unsigned long iPoint; + ifstream restart_file; + string text_line; + + bool multizone = config->GetMultizone_Problem(); + + /*--- Dimension of the problem --> 1 Transport equation (Amplification factor) ---*/ + nVar = 2; + nPrimVar = 2; + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + /*--- Initialize nVarGrad for deallocation ---*/ + + nVarGrad = nVar; + + /*--- Define geometry constants in the solver structure ---*/ + + nDim = geometry->GetnDim(); + + /*--- Single grid simulation ---*/ + + if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) { + + /*--- Define some auxiliar vector related with the residual ---*/ + + Residual_RMS.resize(nVar,0.0); + Residual_Max.resize(nVar,0.0); + Point_Max.resize(nVar,0); + Point_Max_Coord.resize(nVar,nDim) = su2double(0.0); + + /*--- Initialization of the structure of the whole Jacobian ---*/ + + if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (eN transition model)." << endl; + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); + LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + System.SetxIsZero(true); + + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + + if (config->GetExtraOutput()) { + if (nDim == 2) { nOutputVariables = 13; } + else if (nDim == 3) { nOutputVariables = 19; } + OutputVariables.Initialize(nPoint, nPointDomain, nOutputVariables, 0.0); + OutputHeadingNames = new string[nOutputVariables]; + } + + /*--- Initialize the BGS residuals in multizone problems. ---*/ + if (multizone){ + Residual_BGS.resize(nVar,0.0); + Residual_Max_BGS.resize(nVar,0.0); + Point_Max_BGS.resize(nVar,0); + Point_Max_Coord_BGS.resize(nVar,nDim) = su2double(0.0); + } + + } + + /*--- Initialize lower and upper limits ---*/ + lowerlimit[0] = 1.0e-30; + upperlimit[0] = (-8.43 - 2.4*log(2.5*tanh(config->GetTurbulenceIntensity_FreeStream()/2.5)/100))*2; + + lowerlimit[1] = -5; + upperlimit[1] = 1.0e-30; + + /*--- Far-field flow state quantities and initialization. ---*/ + const su2double AmplificationFactor_Inf = 0; + const su2double ModifiedIntermittency_Inf = 0; + + Solution_Inf[0] = AmplificationFactor_Inf; + Solution_Inf[1] = ModifiedIntermittency_Inf; + + /*--- Initialize the solution to the far-field state everywhere. ---*/ + nodes = new CTransENVariable(AmplificationFactor_Inf,ModifiedIntermittency_Inf, nPoint, nDim, nVar, config); + SetBaseClassPointerToNodes(); + + /*--- MPI solution ---*/ + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); + + /*--- Initializate quantities for SlidingMesh Interface ---*/ + SlidingState.resize(nMarker); + SlidingStateNodes.resize(nMarker); + + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == FLUID_INTERFACE) { + SlidingState[iMarker].resize(nVertex[iMarker], nPrimVar+1) = nullptr; + SlidingStateNodes[iMarker].resize(nVertex[iMarker],0); + } + } + + /*-- Allocation of inlets has to happen in derived classes (not CTurbSolver), + due to arbitrary number of turbulence variables ---*/ + Inlet_TurbVars.resize(nMarker); + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) + Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar) = AmplificationFactor_Inf; + + const su2double CFL = config->GetCFL(MGLevel)*config->GetCFLRedCoeff_Turb(); + for (iPoint = 0; iPoint < nPoint; iPoint++) { + nodes->SetLocalCFL(iPoint, CFL); + } + Min_CFL_Local = CFL; + Max_CFL_Local = CFL; + Avg_CFL_Local = CFL; + + /*--- Add the solver name (max 8 characters) ---*/ + SolverName = "eN model"; + +} + +void CTransENSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, + unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { + config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem); + + /*--- Upwind second order reconstruction and gradients ---*/ + CommonPreprocessing(geometry, config, Output); + + /*--- Setting normals for HL parameter. ---*/ + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + SU2_OMP_FOR_STAT(omp_chunk_size) + + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + auto Normal = geometry->nodes->GetNormal(iPoint); + nodes->SetNormal(iPoint, Normal[0], Normal[1], Normal[2]); + nodes->SetAuxVar(iPoint, 0, flowNodes->GetProjVel(iPoint, Normal)); + } + END_SU2_OMP_FOR + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetAuxVar_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetAuxVar_Gradient_LS(geometry, config); + } + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + su2double AuxVarHere = 0.0; + auto Normal = geometry->nodes->GetNormal(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) + AuxVarHere += Normal[iDim] * nodes->GetAuxVarGradient(iPoint, 0, iDim); + nodes->SetAuxVar(iPoint, 0, AuxVarHere); + } + END_SU2_OMP_FOR + +} + +void CTransENSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { + + /*--- Compute e^N model gradients. ---*/ + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetSolution_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetSolution_Gradient_LS(geometry, config); + } +} + +void CTransENSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, CConfig* config) { + + /*--- Define an object to set solver specific numerics contribution. ---*/ + + auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {}; + + /*--- Now instantiate the generic implementation with the functor above. ---*/ + + Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); +} + + +void CTransENSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, + CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); + + /*--- Pick one numerics object per thread. ---*/ + auto* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; + + /*--- Loop over all points. ---*/ + + AD::StartNoSharedReading(); + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Conservative variables w/o reconstruction ---*/ + + numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), nullptr); + + /*--- Gradient of the primitive and conservative variables ---*/ + + numerics->SetPrimVarGradient(flowNodes->GetGradient_Primitive(iPoint), nullptr); + + /*--- Turbulent variables w/o reconstruction, and its gradient ---*/ + /*--- ScalarVar & ScalarVarGradient : Turbulence model solution(k&w) ---*/ + + numerics->SetScalarVar(turbNodes->GetSolution(iPoint), nullptr); + numerics->SetScalarVarGradient(turbNodes->GetGradient(iPoint), nullptr); + + /*--- Transition variables w/o reconstruction, and its gradient ---*/ + + numerics->SetTransVar(nodes->GetSolution(iPoint), nullptr); + numerics->SetTransVarGradient(nodes->GetGradient(iPoint), nullptr); + + /*--- Set volume ---*/ + + numerics->SetVolume(geometry->nodes->GetVolume(iPoint)); + + /*--- Set distance to the surface ---*/ + + numerics->SetDistance(geometry->nodes->GetWall_Distance(iPoint), 0.0); + + /*--- Set vorticity and strain rate magnitude ---*/ + + numerics->SetVorticity(flowNodes->GetVorticity(iPoint), nullptr); + + numerics->SetStrainMag(flowNodes->GetStrainMag(iPoint), 0.0); + + /*--- Set coordinates (for debugging) ---*/ + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); + + /*--- Set wall normal derivative of the wall normal momentum ---*/ + numerics->SetAuxVar(nodes->GetAuxVar(iPoint, 0)); + + /*--- Compute the source term ---*/ + auto residual = numerics->ComputeResidual(config); + + /*--- Subtract residual and the Jacobian ---*/ + + LinSysRes.SubtractBlock(iPoint, residual); + if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + + } + END_SU2_OMP_FOR + + AD::EndNoSharedReading(); + +} + +void CTransENSolver::Source_Template(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, + CConfig *config, unsigned short iMesh) { +} + +void CTransENSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + 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)) { + + /*--- Allocate the value at the infinity ---*/ + + auto V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Retrieve solution at the farfield boundary node ---*/ + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + + conv_numerics->SetPrimitive(V_domain, V_infty); + + /*--- Set turbulent variable at the wall, and at infinity ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); + + /*--- Set Normal (it is necessary to change the sign) ---*/ + /*--- It's mean wall normal zero flux. */ + + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); + + /*--- Grid Movement ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); + + /*--- Compute residuals and Jacobians ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + + /*--- Add residuals and Jacobians ---*/ + + LinSysRes.AddBlock(iPoint, residual); + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR + +} + +void CTransENSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + +} + +void CTransENSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, + unsigned short val_marker) { + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + 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) ---*/ + + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); + + /*--- Allocate the value at the inlet ---*/ + + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Retrieve solution at the farfield boundary node ---*/ + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + + /*--- Set various quantities in the solver class ---*/ + + conv_numerics->SetPrimitive(V_domain, V_inlet); + + /*--- Non-dimensionalize Inlet_TurbVars if Inlet-Files are used. ---*/ + su2double Inlet_Vars[MAXNVAR]; + Inlet_Vars[0] = Inlet_TurbVars[val_marker][iVertex][0]; + Inlet_Vars[1] = Inlet_TurbVars[val_marker][iVertex][1]; + if (config->GetInlet_Profile_From_File()) { + Inlet_Vars[0] /= pow(config->GetVelocity_Ref(), 2); + Inlet_Vars[1] *= config->GetViscosity_Ref() / (config->GetDensity_Ref() * pow(config->GetVelocity_Ref(), 2)); + } + + /*--- Set the eN variable states. ---*/ + /*--- Load the inlet transition eN model variables (uniform by default). ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Inlet_Vars); + + /*--- Set various other quantities in the solver class ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); + + /*--- Compute the residual using an upwind scheme ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR + +} + +void CTransENSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, + CConfig *config, unsigned short val_marker) { + +BC_Far_Field(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + +} + + +void CTransENSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, + bool val_update_geo) { + + const string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", val_iter); + + /*--- To make this routine safe to call in parallel most of it can only be executed by one thread. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ + + if (config->GetRead_Binary_Restart()) { + Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename); + } else { + Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); + } + + /*--- Skip flow variables ---*/ + + unsigned short skipVars = nDim + solver[MESH_0][FLOW_SOL]->GetnVar() + solver[MESH_0][TURB_SOL] ->GetnVar(); + + /*--- Adjust the number of solution variables in the incompressible + restart. We always carry a space in nVar for the energy equation in the + mean flow solver, but we only write it to the restart if it is active. + Therefore, we must reduce skipVars here if energy is inactive so that + the turbulent variables are read correctly. ---*/ + + const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); + const bool energy = config->GetEnergy_Equation(); + const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + + if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + + /*--- Load data from the restart into correct containers. ---*/ + + unsigned long counter = 0; + for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) { + /*--- Retrieve local index. If this node from the restart file lives + on the current processor, we will load and instantiate the vars. ---*/ + + const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global); + + if (iPoint_Local > -1) { + /*--- We need to store this point's data, so jump to the correct + offset in the buffer of data from the restart file and load it. ---*/ + + const auto index = counter * Restart_Vars[1] + skipVars; + for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); + nodes ->SetAmplificationFactor(iPoint_Local, Restart_Data[index + 2]); +// nodes ->SetModifiedIntermittency(iPoint_Local, Restart_Data[index + 3]); + + /*--- Increment the overall counter for how many points have been loaded. ---*/ + counter++; + } + } + + /*--- Detect a wrong solution file ---*/ + + if (counter != nPointDomain) { + SU2_MPI::Error(string("The solution file ") + restart_filename + string(" does not match with the mesh file!\n") + + string("This can be caused by empty lines at the end of the file."), + CURRENT_FUNCTION); + } + + } // end SU2_OMP_MASTER, pre and postprocessing are thread-safe. + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + /*--- MPI solution and compute the eddy viscosity ---*/ + + solver[MESH_0][TRANS_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][TRANS_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + + /*--- For turbulent+species simulations the solver Pre-/Postprocessing is done by the species solver. ---*/ + if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, + RUNTIME_FLOW_SYS, false); + solver[MESH_0][TURB_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0); + solver[MESH_0][TRANS_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0); + } + + /*--- Interpolate the solution down to the coarse multigrid levels ---*/ + + for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) { + const su2double Area_Parent = geometry[iMesh]->nodes->GetVolume(iPoint); + su2double Solution_Coarse[MAXNVAR] = {0.0}; + + for (auto iChildren = 0ul; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) { + const auto Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren); + const su2double Area_Children = geometry[iMesh - 1]->nodes->GetVolume(Point_Fine); + const su2double* Solution_Fine = solver[iMesh - 1][TURB_SOL]->GetNodes()->GetSolution(Point_Fine); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution_Coarse[iVar] += Solution_Fine[iVar] * Area_Children / Area_Parent; + } + } + solver[iMesh][TURB_SOL]->GetNodes()->SetSolution(iPoint, Solution_Coarse); + } + END_SU2_OMP_FOR + + solver[iMesh][TURB_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][TURB_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + + if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, + false); + solver[iMesh][TURB_SOL]->Postprocessing(geometry[iMesh], solver[iMesh], config, iMesh); + } + } + + /*--- Go back to single threaded execution. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Delete the class memory that is used to load the restart. ---*/ + + delete[] Restart_Vars; + Restart_Vars = nullptr; + delete[] Restart_Data; + Restart_Data = nullptr; + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + +} diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index e1f60cebbb13..997fcb927664 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -289,6 +289,7 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain } AD::EndNoSharedReading(); + } void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, @@ -303,6 +304,7 @@ void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, C /*--- Now instantiate the generic implementation with the functor above. ---*/ Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + } void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, @@ -376,9 +378,16 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai } - /*--- Effective Intermittency ---*/ + /*--- Amplification factor ---*/ - if (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { + if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::EN) { + numerics-> SetAmplificationFactor(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint,0)); + numerics-> SetModifiedIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint,1)); + } + + /*--- Effective Intermittency ---*/ + + if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { numerics->SetIntermittencyEff(solver_container[TRANS_SOL]->GetNodes()->GetIntermittencyEff(iPoint)); numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint, 0)); } @@ -389,7 +398,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai /*--- Store the intermittency ---*/ - if (transition_BC || config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { + if (transition_BC || config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { nodes->SetIntermittency(iPoint,numerics->GetIntermittencyEff()); } diff --git a/SU2_CFD/src/variables/CTransENVariable.cpp b/SU2_CFD/src/variables/CTransENVariable.cpp new file mode 100644 index 000000000000..e0dead934dfa --- /dev/null +++ b/SU2_CFD/src/variables/CTransENVariable.cpp @@ -0,0 +1,65 @@ +/*! + * \file CTransENVariable.cpp + * \brief Definition of the solution fields. + * \author R. Roos + * \version 7.4.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + + +#include "../../include/variables/CTransENVariable.hpp" + +CTransENVariable::CTransENVariable(su2double AmplificationFactor,su2double ModifiedIntermittency, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) + : CTurbVariable(npoint, ndim, nvar, config) { + + for(unsigned long iPoint=0; iPoint