From 294deb1755ded710e75f22b3344a8dbc3754d05b Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Tue, 5 Oct 2021 22:46:46 -0700 Subject: [PATCH 01/12] initial commit towards generalized solver indices+ nemo turb --- .../numerics/turbulent/turb_convection.hpp | 12 +- SU2_CFD/include/variables/CEulerVariable.hpp | 13 ++ .../include/variables/CIncEulerVariable.hpp | 15 +++ .../include/variables/CNEMOEulerVariable.hpp | 19 ++- SU2_CFD/src/drivers/CDriver.cpp | 122 +++++++++--------- 5 files changed, 120 insertions(+), 61 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 9bb218827c1..95f0f9e160b 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -36,7 +36,8 @@ * \ingroup ConvDiscr * \author A. Bueno. */ -class CUpwSca_TurbSA final : public CUpwScalar { +template +class CUpwSca_TurbSA_generic final : public CUpwScalar { private: /*! * \brief Adds any extra variables to AD @@ -59,6 +60,9 @@ class CUpwSca_TurbSA final : public CUpwScalar { CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); }; +using CUpwSca_TurbSA = CUpwSca_TurbSA_generic; +using CUpwSca_IncTurbSA = CUpwSca_TurbSA_generic; +using CUpwSca_NEMOTurbSA = CUpwSca_TurbSA_generic; /*! * \class CUpwSca_TurbSST @@ -66,7 +70,8 @@ class CUpwSca_TurbSA final : public CUpwScalar { * \ingroup ConvDiscr * \author A. Campos. */ -class CUpwSca_TurbSST final : public CUpwScalar { +template +class CUpwSca_TurbSST_generic final : public CUpwScalar { private: /*! * \brief Adds any extra variables to AD @@ -89,3 +94,6 @@ class CUpwSca_TurbSST final : public CUpwScalar { CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); }; +using CUpwSca_TurbSST = CUpwSca_TurbSST_generic; +using CUpwSca_IncTurbSST = CUpwSca_TurbSST_generic; +using CUpwSca_NEMOTurbSST = CUpwSca_TurbSST_generic; \ No newline at end of file diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 4aab05de791..e3a9c549b3f 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -41,6 +41,19 @@ class CEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 12; + struct CSolverIndices{ + CFlowIndices(nDim){ + inline unsigned short GetTIndex() const {return 0;} + inline unsigned short GetVelIndex() const {return 1;} + inline unsigned short GetPIndex() const {return nDim+1;} + inline unsigned short GetRhoIndex() const {return nDim+2;} + inline unsigned short GetHIndex() const {return nDim+3;} + inline unsigned short GetAIndex() const {return nDim+4;} + inline unsigned short GetLamViscIndex() const {return nDim+5;} + inline unsigned short GetEddyViscIndex() const {return nDim+6;} + } + } + protected: /*!< \brief Secondary variables (dPdrho_e, dPde_rho, dTdrho_e, dTde_rho, dmudrho_T, dmudT_rho, dktdrho_T, dktdT_rho) * in compressible (Euler: 2, NS: 8) flows. */ diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 0d75dbe3686..7aeaeabf743 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -41,6 +41,21 @@ class CIncEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 12; + struct CSolverIndices{ + CIncFlowIndices(nDim){ + inline unsigned short GetPIndex() const {return 0;} + inline unsigned short GetVelIndex() const {return 1;} + inline unsigned short GetTIndex() const {return nDim+1;} + inline unsigned short GetRhoIndex() const {return nDim+2;} + inline unsigned short GetBetaIndex() const {return nDim+3;} + inline unsigned short GetLamViscIndex() const {return nDim+4;} + inline unsigned short GetEddyViscIndex() const {return nDim+5;} + inline unsigned short GetKIndex() const {return nDim+6;} + inline unsigned short GetCpTotalIndex() const {return nDim+7;} + inline unsigned short GetCvTotalIndex() const {return nDim+8;} + } + } + protected: VectorType Streamwise_Periodic_RecoveredPressure, /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */ Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */ diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index 9f5f65a9c17..23a6272aac2 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -34,7 +34,7 @@ /*! * \class CNEMOEulerVariable * \brief Main class for defining the variables of the NEMO Euler's solver. - * \note Primitive variables (rhos_s, T, Tve, vx, vy, vw, P, rho, h, rhoCvtr, rhoCvve) + * \note Primitive variables (rhos_s, T, Tve, vx, vy, vw, P, rho, h, a, rhoCvtr, rhoCvve) * \ingroup Euler_Equations * \author S. R. Copeland, F. Palacios, W. Maier, C. Garbacz */ @@ -42,6 +42,23 @@ class CNEMOEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 25; + struct CSolverIndices{ + CNEMOIndices(nDim, nSpecies){ + inline unsigned short GetRhosIndex() const {return 0;} + inline unsigned short GetTIndex() const {return nSpecies;} + inline unsigned short GetTveIndex() const {return nSpecies+1;} + inline unsigned short GetVelIndex() const {return nSpecies+2} + inline unsigned short GetPIndex() const {return nSpecies+nDim+2;} + inline unsigned short GetRhoIndex() const {return nSpecies+nDim+3;} + inline unsigned short GetHIndex() const {return nSpecies+nDim+4;} + inline unsigned short GetAIndex() const {return nSpecies+nDim+5;} + inline unsigned short GetRhoCvtrIndex() const {return nSpecies+nDim+6;} + inline unsigned short GetRhoCvveIndex() const {return nSpecies+nDim+7;} + inline unsigned short GetLamViscIndex() const {return nSpecies+nDim+8;} + inline unsigned short GetEddyViscIndex() const {return nSpecies+nDim+9;} + } + } + protected: bool ionization; /*!< \brief Presence of charged species in gas mixture. */ bool monoatomic = false; /*!< \brief Presence of single species gas. */ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 85e2edab153..2ab71d29a4f 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1197,6 +1197,67 @@ void CDriver::Integration_Postprocessing(CIntegration ***integration, CGeometry } +template +void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, CConfig *config, CNumerics ****&numerics) { + + /*--- Definition of the convective scheme for each equation and mesh level ---*/ + + switch (config->GetKind_ConvNumScheme_Turb()) { + case NO_UPWIND: + SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); + break; + case SPACE_UPWIND : + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras || neg_spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras ) { + numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); + } + else if (menter_sst) numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); + } + break; + default: + SU2_MPI::Error("Invalid convective scheme for the turbulence equations.", CURRENT_FUNCTION); + break; + } + + /*--- Definition of the viscous scheme for each equation and mesh level ---*/ + + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras){ + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); + } + else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); + else if (menter_sst) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); + } + + /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ + + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); + else if (e_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); + else if (comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); + else if (e_comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); + else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); + else if (menter_sst) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config); + numerics[iMGlevel][TURB_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Turb, config); + } + + /*--- Definition of the boundary condition method ---*/ + + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { + numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); + } + else if (neg_spalart_allmaras) { + numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); + } + else if (menter_sst) { + numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, config); + } + } +} void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const { if (rank == MASTER_NODE) @@ -1833,64 +1894,9 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Solver definition for the turbulent model problem ---*/ if (turbulent) { - - /*--- Definition of the convective scheme for each equation and mesh level ---*/ - - switch (config->GetKind_ConvNumScheme_Turb()) { - case NO_UPWIND: - SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); - break; - case SPACE_UPWIND : - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras || neg_spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras ) { - numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - } - else if (menter_sst) numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - } - break; - default: - SU2_MPI::Error("Invalid convective scheme for the turbulence equations.", CURRENT_FUNCTION); - break; - } - - /*--- Definition of the viscous scheme for each equation and mesh level ---*/ - - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras){ - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); - } - else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); - else if (menter_sst) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); - } - - /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ - - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); - else if (e_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); - else if (comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); - else if (e_comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); - else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); - else if (menter_sst) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config); - numerics[iMGlevel][TURB_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Turb, config); - } - - /*--- Definition of the boundary condition method ---*/ - - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { - numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); - } - else if (neg_spalart_allmaras) { - numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); - } - else if (menter_sst) { - numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, config); - } - } + if (incompressible) InstantiateTurbulentNumerics(nVar_Turb, config, numerics); + else if (nemo) InstantiateTurbulentNumerics(nVar_Turb, config, numerics); + else InstantiateTurbulentNumerics(nVar_Turb, config, numerics); } /*--- Solver definition for the transition model problem ---*/ From 77dcc80eadd1135cc5b20bfba3469e0323733332 Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sat, 9 Oct 2021 11:38:21 +0100 Subject: [PATCH 02/12] Update SU2_CFD/include/variables/CEulerVariable.hpp --- SU2_CFD/include/variables/CEulerVariable.hpp | 25 ++++++++++---------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index e3a9c549b3f..10dde042d3b 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -41,18 +41,19 @@ class CEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 12; - struct CSolverIndices{ - CFlowIndices(nDim){ - inline unsigned short GetTIndex() const {return 0;} - inline unsigned short GetVelIndex() const {return 1;} - inline unsigned short GetPIndex() const {return nDim+1;} - inline unsigned short GetRhoIndex() const {return nDim+2;} - inline unsigned short GetHIndex() const {return nDim+3;} - inline unsigned short GetAIndex() const {return nDim+4;} - inline unsigned short GetLamViscIndex() const {return nDim+5;} - inline unsigned short GetEddyViscIndex() const {return nDim+6;} - } - } + template + struct CIndices { + const IndexType nDim; + CIndices(IndexType ndim, IndexType) : nDim(ndim) {} + inline IndexType Temperature() const { return 0; } + inline IndexType Velocity() const { return 1; } + inline IndexType Pressure() const { return nDim+1; } + inline IndexType Density() const { return nDim+2; } + inline IndexType Enthalpy() const { return nDim+3; } + inline IndexType SoundSpeed() const { return nDim+4; } + inline IndexType LamViscosity() const { return nDim+5; } + inline IndexType EddyViscosity() const { return nDim+6; } + }; protected: /*!< \brief Secondary variables (dPdrho_e, dPde_rho, dTdrho_e, dTde_rho, dmudrho_T, dmudT_rho, dktdrho_T, dktdT_rho) From bb356910cd21c459c4fe93eea78bc03ce7d23560 Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Mon, 1 Nov 2021 18:37:48 -0700 Subject: [PATCH 03/12] Update incomp and nemo index structure --- .../include/variables/CIncEulerVariable.hpp | 29 ++++++++-------- .../include/variables/CNEMOEulerVariable.hpp | 33 ++++++++++--------- 2 files changed, 32 insertions(+), 30 deletions(-) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 913289007e2..a3738e6631e 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -41,20 +41,21 @@ class CIncEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 12; - struct CSolverIndices{ - CIncFlowIndices(nDim){ - inline unsigned short GetPIndex() const {return 0;} - inline unsigned short GetVelIndex() const {return 1;} - inline unsigned short GetTIndex() const {return nDim+1;} - inline unsigned short GetRhoIndex() const {return nDim+2;} - inline unsigned short GetBetaIndex() const {return nDim+3;} - inline unsigned short GetLamViscIndex() const {return nDim+4;} - inline unsigned short GetEddyViscIndex() const {return nDim+5;} - inline unsigned short GetKIndex() const {return nDim+6;} - inline unsigned short GetCpTotalIndex() const {return nDim+7;} - inline unsigned short GetCvTotalIndex() const {return nDim+8;} - } - } + template + struct CIndices { + const IndexType nDim; + CIndices(IndexType ndim, IndexType) : nDim(ndim) {} + inline IndexType Pressure() const {return 0;} + inline IndexType Velocity() const {return 1;} + inline IndexType Temperature() const {return nDim+1;} + inline IndexType Density() const {return nDim+2;} + inline IndexType Beta() const {return nDim+3;} + inline IndexType LamViscosity() const {return nDim+4;} + inline IndexType EddyViscosity() const {return nDim+5;} + inline IndexType ThemalConduct() const {return nDim+6;} + inline IndexType CpTotal() const {return nDim+7;} + inline IndexType CvTotal() const {return nDim+8;} + }; protected: VectorType Streamwise_Periodic_RecoveredPressure, /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */ diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index 17d379fc272..0d95a74df7c 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -42,22 +42,23 @@ class CNEMOEulerVariable : public CFlowVariable { public: static constexpr size_t MAXNVAR = 25; - struct CSolverIndices{ - CNEMOIndices(nDim, nSpecies){ - inline unsigned short GetRhosIndex() const {return 0;} - inline unsigned short GetTIndex() const {return nSpecies;} - inline unsigned short GetTveIndex() const {return nSpecies+1;} - inline unsigned short GetVelIndex() const {return nSpecies+2} - inline unsigned short GetPIndex() const {return nSpecies+nDim+2;} - inline unsigned short GetRhoIndex() const {return nSpecies+nDim+3;} - inline unsigned short GetHIndex() const {return nSpecies+nDim+4;} - inline unsigned short GetAIndex() const {return nSpecies+nDim+5;} - inline unsigned short GetRhoCvtrIndex() const {return nSpecies+nDim+6;} - inline unsigned short GetRhoCvveIndex() const {return nSpecies+nDim+7;} - inline unsigned short GetLamViscIndex() const {return nSpecies+nDim+8;} - inline unsigned short GetEddyViscIndex() const {return nSpecies+nDim+9;} - } - } + template + struct CIndices { + const IndexType nDim, nSpecies; + CIndices(IndexType ndim, IndexType nSpecies, IndexType) : nDim(ndim), nSpecies(nspecies) {} + inline IndexType SpeciesDensities() const {return 0; } + inline IndexType Temperature() const { return nSpecies; } + inline IndexType Temperature_ve() const { return nSpecies+1; } + inline IndexType Velocity() const { return nSpecies+2; } + inline IndexType Pressure() const {return nSpecies+nDim+2;} + inline IndexType Density() const {return nSpecies+nDim+3;} + inline IndexType Enthalpy() const {return nSpecies+nDim+4;} + inline IndexType SoundSpeed() const {return nSpecies+nDim+5;} + inline IndexType RhoCvtr() const {return nSpecies+nDim+6;} + inline IndexType RhoCvve() const {return nSpecies+nDim+7;} + inline IndexType LamViscosity() const {return nSpecies+nDim+8;} + inline IndexType EddyViscosity() const {return nSpecies+nDim+9;} + }; protected: bool ionization; /*!< \brief Presence of charged species in gas mixture. */ From 3213f855b8c6fed57f91de63c99248a1973e0008 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 6 Nov 2021 23:35:57 +0000 Subject: [PATCH 04/12] put the code in a compilable state --- SU2_CFD/include/drivers/CDriver.hpp | 7 + .../numerics/scalar/scalar_convection.hpp | 70 ++++++++- .../numerics/turbulent/turb_convection.hpp | 78 ++++++++-- .../include/variables/CNEMOEulerVariable.hpp | 10 +- SU2_CFD/src/drivers/CDriver.cpp | 146 +++++++++++------- SU2_CFD/src/meson.build | 2 - .../src/numerics/scalar/scalar_convection.cpp | 95 ------------ .../numerics/turbulent/turb_convection.cpp | 72 --------- 8 files changed, 225 insertions(+), 255 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 52aeb70fb63..089d321cb8b 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -213,6 +213,13 @@ class CDriver { */ void Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const; + /*! + * \brief Helper to instantiate turbulence numerics specialized for different flow solvers. + */ + template + void InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig *config, + const CSolver* turb_solver, CNumerics ****&numerics) const; + /*! * \brief Definition and allocation of all solver classes. * \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved). diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index 0b2095ac40f..f7e70877b57 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -44,6 +44,7 @@ * \ingroup ConvDiscr * \author C. Pederson, A. Bueno., and A. Campos. */ +template class CUpwScalar : public CNumerics { protected: su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0 */ @@ -69,21 +70,80 @@ class CUpwScalar : public CNumerics { 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] ndim - Number of dimensions of the problem. + * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CUpwScalar(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); + CUpwScalar(unsigned short ndim, unsigned short nvar, const CConfig* config) + : CNumerics(ndim, nvar, config), + implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), + incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), + dynamic_grid(config->GetDynamic_Grid()) { + Flux = new su2double[nVar]; + Jacobian_i = new su2double*[nVar]; + Jacobian_j = new su2double*[nVar]; + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double[nVar]; + Jacobian_j[iVar] = new su2double[nVar]; + } + } /*! * \brief Destructor of the class. */ - ~CUpwScalar(void) override; + ~CUpwScalar() override { + delete[] Flux; + if (Jacobian_i != nullptr) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + delete[] Jacobian_i[iVar]; + delete[] Jacobian_j[iVar]; + } + delete[] Jacobian_i; + delete[] Jacobian_j; + } + } /*! * \brief Compute the scalar upwind flux between two nodes i and j. * \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; + CNumerics::ResidualType<> ComputeResidual(const CConfig* config) final { + AD::StartPreacc(); + AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_j, nVar); + if (dynamic_grid) { + AD::SetPreaccIn(GridVel_i, nDim); + AD::SetPreaccIn(GridVel_j, nDim); + } + + ExtraADPreaccIn(); + + Density_i = V_i[nDim + 2]; + Density_j = V_j[nDim + 2]; + + su2double q_ij = 0.0; + if (dynamic_grid) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + su2double Velocity_i = V_i[iDim + 1] - GridVel_i[iDim]; + su2double Velocity_j = V_j[iDim + 1] - GridVel_j[iDim]; + q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; + } + } else { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + q_ij += 0.5 * (V_i[iDim + 1] + V_j[iDim + 1]) * Normal[iDim]; + } + } + + a0 = 0.5 * (q_ij + fabs(q_ij)); + a1 = 0.5 * (q_ij - fabs(q_ij)); + + FinishResidualCalc(config); + + AD::SetPreaccOut(Flux, nVar); + AD::EndPreacc(); + + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); + } }; diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 9bc1ba6f585..3a80fbe11b6 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -37,18 +37,41 @@ * \author A. Bueno. */ template -class CUpwSca_TurbSA_generic final : public CUpwScalar { +class CUpwSca_TurbSA final : public CUpwScalar { private: + using Base = CUpwScalar; + using Base::nDim; + using Base::V_i; + using Base::V_j; + using Base::a0; + using Base::a1; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::implicit; + /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn() override; + void ExtraADPreaccIn() override { + AD::SetPreaccIn(V_i, nDim+1); + AD::SetPreaccIn(V_j, nDim+1); + } /*! * \brief SA specific steps in the ComputeResidual method * \param[in] config - Definition of the particular problem. */ - void FinishResidualCalc(const CConfig* config) override; + void FinishResidualCalc(const CConfig* config) override { + Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0]; + + if (implicit) { + Jacobian_i[0][0] = a0; + Jacobian_j[0][0] = a1; + } + } public: /*! @@ -57,12 +80,9 @@ class CUpwSca_TurbSA_generic final : public CUpwScalar { * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - + CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CUpwScalar(val_nDim, val_nVar, config) {} }; -using CUpwSca_TurbSA = CUpwSca_TurbSA_generic; -using CUpwSca_IncTurbSA = CUpwSca_TurbSA_generic; -using CUpwSca_NEMOTurbSA = CUpwSca_TurbSA_generic; /*! * \class CUpwSca_TurbSST @@ -71,18 +91,47 @@ using CUpwSca_NEMOTurbSA = CUpwSca_TurbSA_generic; * \author A. Campos. */ template -class CUpwSca_TurbSST_generic final : public CUpwScalar { +class CUpwSca_TurbSST final : public CUpwScalar { private: + using Base = CUpwScalar; + using Base::nDim; + using Base::V_i; + using Base::V_j; + using Base::a0; + using Base::a1; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::Density_i; + using Base::Density_j; + using Base::implicit; + /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn() override; + void ExtraADPreaccIn() override { + AD::SetPreaccIn(V_i, nDim+3); + AD::SetPreaccIn(V_j, nDim+3); + } /*! * \brief SST specific steps in the ComputeResidual method * \param[in] config - Definition of the particular problem. */ - void FinishResidualCalc(const CConfig* config) override; + void FinishResidualCalc(const CConfig* config) override { + Flux[0] = a0*Density_i*ScalarVar_i[0] + a1*Density_j*ScalarVar_j[0]; + Flux[1] = a0*Density_i*ScalarVar_i[1] + a1*Density_j*ScalarVar_j[1]; + + if (implicit) { + Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0; + + Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1; + } + } public: /*! @@ -91,9 +140,6 @@ class CUpwSca_TurbSST_generic final : public CUpwScalar { * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - + CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CUpwScalar(val_nDim, val_nVar, config) {} }; -using CUpwSca_TurbSST = CUpwSca_TurbSST_generic; -using CUpwSca_IncTurbSST = CUpwSca_TurbSST_generic; -using CUpwSca_NEMOTurbSST = CUpwSca_TurbSST_generic; \ No newline at end of file diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index 0d95a74df7c..fd0f68fee09 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -45,11 +45,11 @@ class CNEMOEulerVariable : public CFlowVariable { template struct CIndices { const IndexType nDim, nSpecies; - CIndices(IndexType ndim, IndexType nSpecies, IndexType) : nDim(ndim), nSpecies(nspecies) {} - inline IndexType SpeciesDensities() const {return 0; } - inline IndexType Temperature() const { return nSpecies; } - inline IndexType Temperature_ve() const { return nSpecies+1; } - inline IndexType Velocity() const { return nSpecies+2; } + CIndices(IndexType ndim, IndexType nspecies) : nDim(ndim), nSpecies(nspecies) {} + inline IndexType SpeciesDensities() const {return 0;} + inline IndexType Temperature() const {return nSpecies;} + inline IndexType Temperature_ve() const {return nSpecies+1;} + inline IndexType Velocity() const {return nSpecies+2;} inline IndexType Pressure() const {return nSpecies+nDim+2;} inline IndexType Density() const {return nSpecies+nDim+3;} inline IndexType Enthalpy() const {return nSpecies+nDim+4;} diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 25b63fd7d29..b0fcf7786b9 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1198,7 +1198,46 @@ void CDriver::Integration_Postprocessing(CIntegration ***integration, CGeometry } template -void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, CConfig *config, CNumerics ****&numerics) { +void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, const CConfig *config, + const CSolver* turb_solver, CNumerics ****&numerics) const { + const int conv_term = CONV_TERM + offset; + const int visc_term = VISC_TERM + offset; + + const int source_first_term = SOURCE_FIRST_TERM + offset; + const int source_second_term = SOURCE_SECOND_TERM + offset; + + const int conv_bound_term = CONV_BOUND_TERM + offset; + const int visc_bound_term = VISC_BOUND_TERM + offset; + + bool spalart_allmaras, neg_spalart_allmaras, e_spalart_allmaras, comp_spalart_allmaras, e_comp_spalart_allmaras, menter_sst; + spalart_allmaras = neg_spalart_allmaras = e_spalart_allmaras = comp_spalart_allmaras = e_comp_spalart_allmaras = menter_sst = false; + + /*--- Assign turbulence model booleans ---*/ + + switch (config->GetKind_Turb_Model()) { + case TURB_MODEL::NONE: + SU2_MPI::Error("No turbulence model selected.", CURRENT_FUNCTION); + break; + case TURB_MODEL::SA: spalart_allmaras = true; break; + case TURB_MODEL::SA_NEG: neg_spalart_allmaras = true; break; + case TURB_MODEL::SA_E: e_spalart_allmaras = true; break; + case TURB_MODEL::SA_COMP: comp_spalart_allmaras = true; break; + case TURB_MODEL::SA_E_COMP: e_comp_spalart_allmaras = true; break; + case TURB_MODEL::SST: menter_sst = true; break; + case TURB_MODEL::SST_SUST: menter_sst = true; break; + } + + /*--- If the Menter SST model is used, store the constants of the model and determine the + free stream values of the turbulent kinetic energy and dissipation rate. ---*/ + + const su2double *constants = nullptr; + su2double kine_Inf = 0.0, omega_Inf = 0.0; + + if (menter_sst) { + constants = turb_solver->GetConstants(); + kine_Inf = turb_solver->GetTke_Inf(); + omega_Inf = turb_solver->GetOmega_Inf(); + } /*--- Definition of the convective scheme for each equation and mesh level ---*/ @@ -1207,9 +1246,9 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, CConfig *co SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); break; case SPACE_UPWIND : - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras || neg_spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras ) { - numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); + for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras || neg_spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { + numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); } else if (menter_sst) numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); } @@ -1221,43 +1260,58 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, CConfig *co /*--- Definition of the viscous scheme for each equation and mesh level ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras){ - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); + for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); } - else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); - else if (menter_sst) numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); + else if (neg_spalart_allmaras) + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); + else if (menter_sst) + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); - else if (e_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); - else if (comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); - else if (e_comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); - else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); - else if (menter_sst) numerics[iMGlevel][TURB_SOL][source_first_term] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config); + for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + auto& turb_source_first_term = numerics[iMGlevel][TURB_SOL][source_first_term]; + + if (spalart_allmaras) + turb_source_first_term = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); + else if (e_spalart_allmaras) + turb_source_first_term = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); + else if (comp_spalart_allmaras) + turb_source_first_term = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); + else if (e_comp_spalart_allmaras) + turb_source_first_term = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); + else if (neg_spalart_allmaras) + turb_source_first_term = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); + else if (menter_sst) + turb_source_first_term = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config); + numerics[iMGlevel][TURB_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Turb, config); } /*--- Definition of the boundary condition method ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); } else if (neg_spalart_allmaras) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); } else if (menter_sst) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, config); } } } +/*--- Explicit instantiation of the template above, needed because it is defined in a cpp file, instead of hpp. ---*/ +template void CDriver::InstantiateTurbulentNumerics(unsigned short, int, const CConfig*, + const CSolver*, CNumerics****&) const; + void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const { if (rank == MASTER_NODE) @@ -1280,21 +1334,16 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol numerics = new CNumerics***[config->GetnMGLevels()+1] (); - const su2double *constants = nullptr; - su2double kine_Inf = 0.0, omega_Inf = 0.0; - bool compressible = false; bool incompressible = false; bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS); bool roe_low_dissipation = (config->GetKind_RoeLowDiss() != NO_ROELOWDISS); /*--- Initialize some useful booleans ---*/ - bool euler, ns, NEMO_euler, NEMO_ns, turbulent, adj_euler, adj_ns, adj_turb, fem_euler, fem_ns, fem_turbulent; - bool spalart_allmaras, neg_spalart_allmaras, e_spalart_allmaras, comp_spalart_allmaras, e_comp_spalart_allmaras, menter_sst; + bool euler, ns, NEMO_euler, NEMO_ns, turbulent, adj_euler, adj_ns, adj_turb, fem_euler, fem_ns; bool fem, heat, transition, template_solver; - euler = ns = NEMO_euler = NEMO_ns = turbulent = adj_euler = adj_ns = adj_turb = fem_euler = fem_ns = fem_turbulent = false; - spalart_allmaras = neg_spalart_allmaras = e_spalart_allmaras = comp_spalart_allmaras = e_comp_spalart_allmaras = menter_sst = false; + euler = ns = NEMO_euler = NEMO_ns = turbulent = adj_euler = adj_ns = adj_turb = fem_euler = fem_ns = false; fem = heat = transition = template_solver = false; /*--- Assign booleans ---*/ @@ -1346,7 +1395,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol case FEM_RANS: case DISC_ADJ_FEM_RANS: - fem_ns = compressible = fem_turbulent = true; break; + fem_ns = compressible = true; break; case FEM_LES: fem_ns = compressible = true; break; @@ -1372,31 +1421,6 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } - /*--- Assign turbulence model booleans ---*/ - - if (turbulent || fem_turbulent) - switch (config->GetKind_Turb_Model()) { - case TURB_MODEL::NONE: - SU2_MPI::Error("No turbulence model selected.", CURRENT_FUNCTION); - break; - case TURB_MODEL::SA: spalart_allmaras = true; break; - case TURB_MODEL::SA_NEG: neg_spalart_allmaras = true; break; - case TURB_MODEL::SA_E: e_spalart_allmaras = true; break; - case TURB_MODEL::SA_COMP: comp_spalart_allmaras = true; break; - case TURB_MODEL::SA_E_COMP: e_comp_spalart_allmaras = true; break; - case TURB_MODEL::SST: menter_sst = true; break; - case TURB_MODEL::SST_SUST: menter_sst = true; break; - } - - /*--- If the Menter SST model is used, store the constants of the model and determine the - free stream values of the turbulent kinetic energy and dissipation rate. ---*/ - - if (menter_sst) { - constants = solver[MESH_0][TURB_SOL]->GetConstants(); - kine_Inf = solver[MESH_0][TURB_SOL]->GetTke_Inf(); - omega_Inf = solver[MESH_0][TURB_SOL]->GetOmega_Inf(); - } - /*--- Number of variables for the template ---*/ if (template_solver) nVar_Flow = solver[MESH_0][FLOW_SOL]->GetnVar(); @@ -1412,12 +1436,11 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (fem_euler) nVar_Flow = solver[MESH_0][FLOW_SOL]->GetnVar(); if (fem_ns) nVar_Flow = solver[MESH_0][FLOW_SOL]->GetnVar(); - //if (fem_turbulent) nVar_Turb = solver_container[MESH_0][FEM_TURB_SOL]->GetnVar(); if (fem) nVar_FEM = solver[MESH_0][FEA_SOL]->GetnVar(); if (heat) nVar_Heat = solver[MESH_0][HEAT_SOL]->GetnVar(); - if (config->AddRadiation()) nVar_Rad = solver[MESH_0][RAD_SOL]->GetnVar(); + if (config->AddRadiation()) nVar_Rad = solver[MESH_0][RAD_SOL]->GetnVar(); /*--- Number of variables for adjoint problem ---*/ @@ -1430,7 +1453,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (NEMO_euler || NEMO_ns) nPrimVar_NEMO = solver[MESH_0][FLOW_SOL]->GetnPrimVar(); if (NEMO_euler || NEMO_ns) nPrimVarGrad_NEMO = solver[MESH_0][FLOW_SOL]->GetnPrimVarGrad(); - /*--- Definition of the Class for the numerical method: numerics_container[INSTANCE_LEVEL][MESH_LEVEL][EQUATION][EQ_TERM] ---*/ + /*--- Definition of the Class for the numerical method: + numerics_container[INSTANCE_LEVEL][MESH_LEVEL][EQUATION][EQ_TERM] ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { numerics[iMGlevel] = new CNumerics** [MAX_SOLS]; @@ -1894,9 +1918,11 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Solver definition for the turbulent model problem ---*/ if (turbulent) { - if (incompressible) InstantiateTurbulentNumerics(nVar_Turb, config, numerics); - else if (nemo) InstantiateTurbulentNumerics(nVar_Turb, config, numerics); - else InstantiateTurbulentNumerics(nVar_Turb, config, numerics); + /// TODO: Replace the dummy int type here for the real thing commented below. + InstantiateTurbulentNumerics(nVar_Turb, offset, config, solver[MESH_0][TURB_SOL], numerics); +// if (incompressible) InstantiateTurbulentNumerics(...); +// else if (nemo) InstantiateTurbulentNumerics(...); +// else InstantiateTurbulentNumerics(...); } /*--- Solver definition for the transition model problem ---*/ @@ -2094,7 +2120,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Solver definition for the turbulent adjoint problem ---*/ if (adj_turb) { - if (!spalart_allmaras) + if (config->GetKind_Turb_Model() != TURB_MODEL::SA) SU2_MPI::Error("Only the SA turbulence model can be used with the continuous adjoint solver.", CURRENT_FUNCTION); /*--- Definition of the convective scheme for each equation and mesh level ---*/ diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index c9ace58e3d6..8d444ab8054 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -134,10 +134,8 @@ su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/continuous_adjoint/adj_convection.cpp', 'numerics/continuous_adjoint/adj_diffusion.cpp', 'numerics/continuous_adjoint/adj_sources.cpp', - 'numerics/scalar/scalar_convection.cpp', 'numerics/scalar/scalar_diffusion.cpp', 'numerics/scalar/scalar_sources.cpp', - 'numerics/turbulent/turb_convection.cpp', 'numerics/turbulent/turb_diffusion.cpp', 'numerics/turbulent/turb_sources.cpp', 'numerics/elasticity/CFEAElasticity.cpp', diff --git a/SU2_CFD/src/numerics/scalar/scalar_convection.cpp b/SU2_CFD/src/numerics/scalar/scalar_convection.cpp index 57ce254f086..8b137891791 100644 --- a/SU2_CFD/src/numerics/scalar/scalar_convection.cpp +++ b/SU2_CFD/src/numerics/scalar/scalar_convection.cpp @@ -1,96 +1 @@ -/*! - * \file scalar_convection.cpp - * \brief Implementation of numerics classes to compute convective - * fluxes in scalar problems. - * \author F. Palacios, T. Economon - * \version 7.2.1 "Blackbird" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2021, 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/numerics/scalar/scalar_convection.hpp" - -CUpwScalar::CUpwScalar(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CNumerics(val_nDim, val_nVar, config), - implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), - incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), - dynamic_grid(config->GetDynamic_Grid()) { - Flux = new su2double[nVar]; - Jacobian_i = new su2double*[nVar]; - Jacobian_j = new su2double*[nVar]; - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double[nVar]; - Jacobian_j[iVar] = new su2double[nVar]; - } -} - -CUpwScalar::~CUpwScalar(void) { - delete[] Flux; - if (Jacobian_i != nullptr) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - delete[] Jacobian_i[iVar]; - delete[] Jacobian_j[iVar]; - } - delete[] Jacobian_i; - delete[] Jacobian_j; - } -} - -CNumerics::ResidualType<> CUpwScalar::ComputeResidual(const CConfig* config) { - unsigned short iDim; - - AD::StartPreacc(); - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_j, nVar); - if (dynamic_grid) { - AD::SetPreaccIn(GridVel_i, nDim); - AD::SetPreaccIn(GridVel_j, nDim); - } - - ExtraADPreaccIn(); - - Density_i = V_i[nDim + 2]; - Density_j = V_j[nDim + 2]; - - su2double q_ij = 0.0; - if (dynamic_grid) { - for (iDim = 0; iDim < nDim; iDim++) { - su2double Velocity_i = V_i[iDim + 1] - GridVel_i[iDim]; - su2double Velocity_j = V_j[iDim + 1] - GridVel_j[iDim]; - q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; - } - } else { - for (iDim = 0; iDim < nDim; iDim++) { - q_ij += 0.5 * (V_i[iDim + 1] + V_j[iDim + 1]) * Normal[iDim]; - } - } - - a0 = 0.5 * (q_ij + fabs(q_ij)); - a1 = 0.5 * (q_ij - fabs(q_ij)); - - FinishResidualCalc(config); - - AD::SetPreaccOut(Flux, nVar); - AD::EndPreacc(); - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); -} diff --git a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp b/SU2_CFD/src/numerics/turbulent/turb_convection.cpp index c3e523b1606..8b137891791 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_convection.cpp @@ -1,73 +1 @@ -/*! - * \file turb_convection.cpp - * \brief Implementation of numerics classes to compute convective - * fluxes in turbulence problems. - * \author F. Palacios, T. Economon - * \version 7.2.1 "Blackbird" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2021, 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/numerics/turbulent/turb_convection.hpp" - -CUpwSca_TurbSA::CUpwSca_TurbSA(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CUpwScalar(val_nDim, val_nVar, config) { } - -void CUpwSca_TurbSA::ExtraADPreaccIn() { - AD::SetPreaccIn(V_i, nDim+1); - AD::SetPreaccIn(V_j, nDim+1); -} - -void CUpwSca_TurbSA::FinishResidualCalc(const CConfig* config) { - - Flux[0] = a0*ScalarVar_i[0]+a1*ScalarVar_j[0]; - - if (implicit) { - Jacobian_i[0][0] = a0; - Jacobian_j[0][0] = a1; - } -} - -CUpwSca_TurbSST::CUpwSca_TurbSST(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CUpwScalar(val_nDim, val_nVar, config) { } - -void CUpwSca_TurbSST::ExtraADPreaccIn() { - AD::SetPreaccIn(V_i, nDim+3); - AD::SetPreaccIn(V_j, nDim+3); -} - -void CUpwSca_TurbSST::FinishResidualCalc(const CConfig* config) { - - Flux[0] = a0*Density_i*ScalarVar_i[0]+a1*Density_j*ScalarVar_j[0]; - Flux[1] = a0*Density_i*ScalarVar_i[1]+a1*Density_j*ScalarVar_j[1]; - - if (implicit) { - Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0; - - Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1; - } -} From 15a6915fb93cda29f94a3b4bc3735a50d3d96226 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 12:18:16 +0000 Subject: [PATCH 05/12] use the index classes in the variables --- SU2_CFD/include/variables/CEulerVariable.hpp | 33 +++++---- .../include/variables/CIncEulerVariable.hpp | 68 +++++++++++-------- SU2_CFD/include/variables/CIncNSVariable.hpp | 18 +++-- .../include/variables/CNEMOEulerVariable.hpp | 4 +- SU2_CFD/include/variables/CNEMONSVariable.hpp | 7 -- SU2_CFD/include/variables/CNSVariable.hpp | 33 +++++---- SU2_CFD/include/variables/CVariable.hpp | 12 ---- SU2_CFD/src/variables/CEulerVariable.cpp | 3 +- SU2_CFD/src/variables/CIncEulerVariable.cpp | 14 ++-- SU2_CFD/src/variables/CIncNSVariable.cpp | 13 ++-- SU2_CFD/src/variables/CNEMOEulerVariable.cpp | 1 + SU2_CFD/src/variables/CNSVariable.cpp | 10 +-- 12 files changed, 113 insertions(+), 103 deletions(-) diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 733d3a0e8c0..31f5a68bb6f 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -51,11 +51,15 @@ class CEulerVariable : public CFlowVariable { inline IndexType Density() const { return nDim+2; } inline IndexType Enthalpy() const { return nDim+3; } inline IndexType SoundSpeed() const { return nDim+4; } - inline IndexType LamViscosity() const { return nDim+5; } + inline IndexType LaminarViscosity() const { return nDim+5; } inline IndexType EddyViscosity() const { return nDim+6; } + inline IndexType ThermalConductivity() const { return nDim+7; } + inline IndexType CpTotal() const { return nDim+8; } }; protected: + const CIndices indices; + /*!< \brief Secondary variables (dPdrho_e, dPde_rho, dTdrho_e, dTde_rho, dmudrho_T, dmudT_rho, dktdrho_T, dktdT_rho) * in compressible (Euler: 2, NS: 8) flows. */ MatrixType Secondary; @@ -111,7 +115,8 @@ class CEulerVariable : public CFlowVariable { * \brief Set the value of the enthalpy. */ inline void SetEnthalpy(unsigned long iPoint) final { - Primitive(iPoint,nDim+3) = (Solution(iPoint,nVar-1) + Primitive(iPoint,nDim+1)) / Solution(iPoint,0); + Primitive(iPoint, indices.Enthalpy()) = + (Solution(iPoint,nVar-1) + Primitive(iPoint, indices.Pressure())) / Solution(iPoint,0); } /*! @@ -166,8 +171,8 @@ class CEulerVariable : public CFlowVariable { * \brief Set the value of the density for the incompressible flows. */ inline bool SetDensity(unsigned long iPoint) final { - Primitive(iPoint,nDim+2) = Solution(iPoint,0); - return Primitive(iPoint,nDim+2) <= 0.0; + Primitive(iPoint, indices.Density()) = Solution(iPoint,0); + return Primitive(iPoint, indices.Density()) <= 0.0; } /*! @@ -175,7 +180,7 @@ class CEulerVariable : public CFlowVariable { * \param[in] temperature - how agitated the particles are :) */ inline bool SetTemperature(unsigned long iPoint, su2double temperature) final { - Primitive(iPoint,0) = temperature; + Primitive(iPoint, indices.Temperature()) = temperature; return temperature <= 0.0; } @@ -183,19 +188,19 @@ class CEulerVariable : public CFlowVariable { * \brief Get the flow pressure. * \return Value of the flow pressure. */ - inline su2double GetPressure(unsigned long iPoint) const final { return Primitive(iPoint,nDim+1); } + inline su2double GetPressure(unsigned long iPoint) const final { return Primitive(iPoint, indices.Pressure()); } /*! * \brief Get the speed of the sound. * \return Value of speed of the sound. */ - inline su2double GetSoundSpeed(unsigned long iPoint) const final { return Primitive(iPoint,nDim+4); } + inline su2double GetSoundSpeed(unsigned long iPoint) const final { return Primitive(iPoint, indices.SoundSpeed()); } /*! * \brief Get the enthalpy of the flow. * \return Value of the enthalpy of the flow. */ - inline su2double GetEnthalpy(unsigned long iPoint) const final { return Primitive(iPoint,nDim+3); } + inline su2double GetEnthalpy(unsigned long iPoint) const final { return Primitive(iPoint, indices.Enthalpy()); } /*! * \brief Get the density of the flow. @@ -213,14 +218,16 @@ class CEulerVariable : public CFlowVariable { * \brief Get the temperature of the flow. * \return Value of the temperature of the flow. */ - inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint,0); } + inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint,indices.Temperature()); } /*! * \brief Get the velocity of the flow. * \param[in] iDim - Index of the dimension. * \return Value of the velocity for the dimension iDim. */ - inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { return Primitive(iPoint,iDim+1); } + inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { + return Primitive(iPoint,iDim+indices.Velocity()); + } /*! * \brief Get the projected velocity in a unitary vector direction (compressible solver). @@ -230,7 +237,7 @@ class CEulerVariable : public CFlowVariable { inline su2double GetProjVel(unsigned long iPoint, const su2double *val_vector) const final { su2double ProjVel = 0.0; for (unsigned long iDim = 0; iDim < nDim; iDim++) - ProjVel += Primitive(iPoint,iDim+1)*val_vector[iDim]; + ProjVel += Primitive(iPoint,iDim+indices.Velocity())*val_vector[iDim]; return ProjVel; } @@ -241,8 +248,8 @@ class CEulerVariable : public CFlowVariable { inline void SetVelocity(unsigned long iPoint) final { Velocity2(iPoint) = 0.0; for (unsigned long iDim = 0; iDim < nDim; iDim++) { - Primitive(iPoint,iDim+1) = Solution(iPoint,iDim+1) / Solution(iPoint,0); - Velocity2(iPoint) += pow(Primitive(iPoint,iDim+1),2); + Primitive(iPoint,iDim+indices.Velocity()) = Solution(iPoint,iDim+1) / Solution(iPoint,0); + Velocity2(iPoint) += pow(Primitive(iPoint,iDim+indices.Velocity()),2); } } diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index a3738e6631e..855aa09694c 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -45,19 +45,21 @@ class CIncEulerVariable : public CFlowVariable { struct CIndices { const IndexType nDim; CIndices(IndexType ndim, IndexType) : nDim(ndim) {} - inline IndexType Pressure() const {return 0;} - inline IndexType Velocity() const {return 1;} - inline IndexType Temperature() const {return nDim+1;} - inline IndexType Density() const {return nDim+2;} - inline IndexType Beta() const {return nDim+3;} - inline IndexType LamViscosity() const {return nDim+4;} - inline IndexType EddyViscosity() const {return nDim+5;} - inline IndexType ThemalConduct() const {return nDim+6;} - inline IndexType CpTotal() const {return nDim+7;} - inline IndexType CvTotal() const {return nDim+8;} + inline IndexType Pressure() const { return 0; } + inline IndexType Velocity() const { return 1; } + inline IndexType Temperature() const { return nDim+1; } + inline IndexType Density() const { return nDim+2; } + inline IndexType Beta() const { return nDim+3; } + inline IndexType LaminarViscosity() const { return nDim+4; } + inline IndexType EddyViscosity() const { return nDim+5; } + inline IndexType ThermalConductivity() const { return nDim+6; } + inline IndexType CpTotal() const { return nDim+7; } + inline IndexType CvTotal() const { return nDim+8; } }; protected: + const CIndices indices; + VectorType Streamwise_Periodic_RecoveredPressure, /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */ Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */ public: @@ -78,16 +80,15 @@ class CIncEulerVariable : public CFlowVariable { * \brief Set the value of the pressure. * \param[in] iPoint - Point index. */ - inline void SetPressure(unsigned long iPoint) final { Primitive(iPoint,0) = Solution(iPoint,0); } + inline void SetPressure(unsigned long iPoint) final { Primitive(iPoint, indices.Pressure()) = Solution(iPoint,0); } /*! * \brief Set the value of the density for the incompressible flows. * \param[in] iPoint - Point index. */ inline bool SetDensity(unsigned long iPoint, su2double val_density) final { - Primitive(iPoint,nDim+2) = val_density; - if (Primitive(iPoint,nDim+2) > 0.0) return false; - else return true; + Primitive(iPoint, indices.Density()) = val_density; + return val_density <= 0.0; } /*! @@ -97,8 +98,8 @@ class CIncEulerVariable : public CFlowVariable { inline void SetVelocity(unsigned long iPoint) final { Velocity2(iPoint) = 0.0; for (unsigned long iDim = 0; iDim < nDim; iDim++) { - Primitive(iPoint,iDim+1) = Solution(iPoint,iDim+1); - Velocity2(iPoint) += pow(Primitive(iPoint,iDim+1),2); + Primitive(iPoint, iDim+indices.Velocity()) = Solution(iPoint,iDim+1); + Velocity2(iPoint) += pow(Primitive(iPoint, iDim+indices.Velocity()), 2); } } @@ -107,47 +108,50 @@ class CIncEulerVariable : public CFlowVariable { * \param[in] iPoint - Point index. */ inline bool SetTemperature(unsigned long iPoint, su2double val_temperature) final { - Primitive(iPoint,nDim+1) = val_temperature; - if (Primitive(iPoint,nDim+1) > 0.0) return false; - else return true; + Primitive(iPoint, indices.Temperature()) = val_temperature; + return val_temperature <= 0.0; } /*! * \brief Set the value of the beta coeffient for incompressible flows. * \param[in] iPoint - Point index. */ - inline void SetBetaInc2(unsigned long iPoint, su2double val_betainc2) final { Primitive(iPoint,nDim+3) = val_betainc2; } + inline void SetBetaInc2(unsigned long iPoint, su2double betainc2) final { + Primitive(iPoint, indices.Beta()) = betainc2; + } /*! * \brief Get the flow pressure. * \return Value of the flow pressure. */ - inline su2double GetPressure(unsigned long iPoint) const final { return Primitive(iPoint,0); } + inline su2double GetPressure(unsigned long iPoint) const final { return Primitive(iPoint, indices.Pressure()); } /*! * \brief Get the value of beta squared for the incompressible flow * \return Value of beta squared. */ - inline su2double GetBetaInc2(unsigned long iPoint) const final { return Primitive(iPoint,nDim+3); } + inline su2double GetBetaInc2(unsigned long iPoint) const final { return Primitive(iPoint, indices.Beta()); } /*! * \brief Get the density of the flow. * \return Value of the density of the flow. */ - inline su2double GetDensity(unsigned long iPoint) const final { return Primitive(iPoint,nDim+2); } + inline su2double GetDensity(unsigned long iPoint) const final { return Primitive(iPoint, indices.Density()); } /*! * \brief Get the temperature of the flow. * \return Value of the temperature of the flow. */ - inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint,nDim+1); } + inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint, indices.Temperature()); } /*! * \brief Get the velocity of the flow. * \param[in] iDim - Index of the dimension. * \return Value of the velocity for the dimension iDim. */ - inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { return Primitive(iPoint,iDim+1); } + inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { + return Primitive(iPoint, iDim+indices.Velocity()); + } /*! * \brief Get the projected velocity in a unitary vector direction (compressible solver). @@ -157,7 +161,7 @@ class CIncEulerVariable : public CFlowVariable { inline su2double GetProjVel(unsigned long iPoint, const su2double *val_vector) const final { su2double ProjVel = 0.0; for (unsigned long iDim = 0; iDim < nDim; iDim++) - ProjVel += Primitive(iPoint,iDim+1)*val_vector[iDim]; + ProjVel += Primitive(iPoint, iDim+indices.Velocity())*val_vector[iDim]; return ProjVel; } @@ -186,24 +190,28 @@ class CIncEulerVariable : public CFlowVariable { /*! * \brief Set the specific heat Cp. */ - inline void SetSpecificHeatCp(unsigned long iPoint, su2double val_Cp) final { Primitive(iPoint, nDim+7) = val_Cp; } + inline void SetSpecificHeatCp(unsigned long iPoint, su2double val_Cp) final { + Primitive(iPoint, indices.CpTotal()) = val_Cp; + } /*! * \brief Set the specific heat Cv. */ - inline void SetSpecificHeatCv(unsigned long iPoint, su2double val_Cv) final { Primitive(iPoint, nDim+8) = val_Cv; } + inline void SetSpecificHeatCv(unsigned long iPoint, su2double val_Cv) final { + Primitive(iPoint, indices.CvTotal()) = val_Cv; + } /*! * \brief Get the specific heat at constant P of the flow. * \return Value of the specific heat at constant P of the flow. */ - inline su2double GetSpecificHeatCp(unsigned long iPoint) const final { return Primitive(iPoint, nDim+7); } + inline su2double GetSpecificHeatCp(unsigned long iPoint) const final { return Primitive(iPoint, indices.CpTotal()); } /*! * \brief Get the specific heat at constant V of the flow. * \return Value of the specific heat at constant V of the flow. */ - inline su2double GetSpecificHeatCv(unsigned long iPoint) const final { return Primitive(iPoint, nDim+8); } + inline su2double GetSpecificHeatCv(unsigned long iPoint) const final { return Primitive(iPoint, indices.CvTotal()); } /*! * \brief Set the recovered pressure for streamwise periodic flow. diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index d38294f7214..10b52b750f6 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/SU2_CFD/include/variables/CIncNSVariable.hpp @@ -59,7 +59,7 @@ class CIncNSVariable final : public CIncEulerVariable { * \brief Set the laminar viscosity. */ inline void SetLaminarViscosity(unsigned long iPoint, su2double laminarViscosity) override { - Primitive(iPoint,nDim+4) = laminarViscosity; + Primitive(iPoint, indices.LaminarViscosity()) = laminarViscosity; } /*! @@ -67,33 +67,39 @@ class CIncNSVariable final : public CIncEulerVariable { * \param[in] eddy_visc - Value of the eddy viscosity. */ inline void SetEddyViscosity(unsigned long iPoint, su2double eddy_visc) override { - Primitive(iPoint,nDim+5) = eddy_visc; + Primitive(iPoint, indices.EddyViscosity()) = eddy_visc; } /*! * \brief Get the laminar viscosity of the flow. * \return Value of the laminar viscosity of the flow. */ - inline su2double GetLaminarViscosity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+4); } + inline su2double GetLaminarViscosity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.LaminarViscosity()); + } /*! * \brief Get the eddy viscosity of the flow. * \return The eddy viscosity of the flow. */ - inline su2double GetEddyViscosity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+5); } + inline su2double GetEddyViscosity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.EddyViscosity()); + } /*! * \brief Set the thermal conductivity. */ inline void SetThermalConductivity(unsigned long iPoint, su2double thermalConductivity) override { - Primitive(iPoint,nDim+6) = thermalConductivity; + Primitive(iPoint, indices.ThermalConductivity()) = thermalConductivity; } /*! * \brief Get the thermal conductivity of the flow. * \return Value of the laminar viscosity of the flow. */ - inline su2double GetThermalConductivity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+6); } + inline su2double GetThermalConductivity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.ThermalConductivity()); + } /*! * \brief Set all the primitive variables for incompressible flows diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index fd0f68fee09..9037011884d 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -56,11 +56,13 @@ class CNEMOEulerVariable : public CFlowVariable { inline IndexType SoundSpeed() const {return nSpecies+nDim+5;} inline IndexType RhoCvtr() const {return nSpecies+nDim+6;} inline IndexType RhoCvve() const {return nSpecies+nDim+7;} - inline IndexType LamViscosity() const {return nSpecies+nDim+8;} + inline IndexType LaminarViscosity() const {return nSpecies+nDim+8;} inline IndexType EddyViscosity() const {return nSpecies+nDim+9;} }; protected: + const CIndices indices; + bool ionization; /*!< \brief Presence of charged species in gas mixture. */ bool monoatomic = false; /*!< \brief Presence of single species gas. */ diff --git a/SU2_CFD/include/variables/CNEMONSVariable.hpp b/SU2_CFD/include/variables/CNEMONSVariable.hpp index 00a0bfd5286..2b157e3a815 100644 --- a/SU2_CFD/include/variables/CNEMONSVariable.hpp +++ b/SU2_CFD/include/variables/CNEMONSVariable.hpp @@ -128,11 +128,4 @@ class CNEMONSVariable final : public CNEMOEulerVariable { * \return Value of the laminar viscosity of the flow. */ inline su2double GetThermalConductivity_ve(unsigned long iPoint) const override { return ThermalCond_ve(iPoint); } - - /*! - * \brief Set the temperature at the wall - */ - inline void SetWallTemperature(unsigned long iPoint, su2double temperature_wall) override { - Primitive(iPoint,T_INDEX) = temperature_wall; - } }; diff --git a/SU2_CFD/include/variables/CNSVariable.hpp b/SU2_CFD/include/variables/CNSVariable.hpp index 8db072702bd..a7d2e102492 100644 --- a/SU2_CFD/include/variables/CNSVariable.hpp +++ b/SU2_CFD/include/variables/CNSVariable.hpp @@ -62,56 +62,61 @@ class CNSVariable final : public CEulerVariable { * \brief Set the laminar viscosity. */ inline void SetLaminarViscosity(unsigned long iPoint, su2double laminarViscosity) override { - Primitive(iPoint,nDim+5) = laminarViscosity; + Primitive(iPoint, indices.LaminarViscosity()) = laminarViscosity; } /*! * \brief Set the laminar viscosity. */ inline void SetThermalConductivity(unsigned long iPoint, su2double thermalConductivity) override { - Primitive(iPoint,nDim+7) = thermalConductivity; + Primitive(iPoint, indices.ThermalConductivity()) = thermalConductivity; } /*! * \brief Set the specific heat Cp. */ - inline void SetSpecificHeatCp(unsigned long iPoint, su2double val_Cp) override { Primitive(iPoint,nDim+8) = val_Cp; } + inline void SetSpecificHeatCp(unsigned long iPoint, su2double val_Cp) override { + Primitive(iPoint, indices.CpTotal()) = val_Cp; + } /*! * \overload * \param[in] eddy_visc - Value of the eddy viscosity. */ - inline void SetEddyViscosity(unsigned long iPoint, su2double eddy_visc) override { Primitive(iPoint,nDim+6) = eddy_visc; } + inline void SetEddyViscosity(unsigned long iPoint, su2double eddy_visc) override { + Primitive(iPoint, indices.EddyViscosity()) = eddy_visc; + } /*! * \brief Get the laminar viscosity of the flow. * \return Value of the laminar viscosity of the flow. */ - inline su2double GetLaminarViscosity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+5); } + inline su2double GetLaminarViscosity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.LaminarViscosity()); + } /*! * \brief Get the thermal conductivity of the flow. * \return Value of the laminar viscosity of the flow. */ - inline su2double GetThermalConductivity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+7); } + inline su2double GetThermalConductivity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.ThermalConductivity()); + } /*! * \brief Get the eddy viscosity of the flow. * \return The eddy viscosity of the flow. */ - inline su2double GetEddyViscosity(unsigned long iPoint) const override { return Primitive(iPoint,nDim+6); } + inline su2double GetEddyViscosity(unsigned long iPoint) const override { + return Primitive(iPoint, indices.EddyViscosity()); + } /*! * \brief Get the specific heat at constant P of the flow. * \return Value of the specific heat at constant P of the flow. */ - inline su2double GetSpecificHeatCp(unsigned long iPoint) const override { return Primitive(iPoint,nDim+8); } - - /*! - * \brief Set the temperature at the wall - */ - inline void SetWallTemperature(unsigned long iPoint, su2double temperature_wall) override { - Primitive(iPoint,0) = temperature_wall; + inline su2double GetSpecificHeatCp(unsigned long iPoint) const override { + return Primitive(iPoint, indices.CpTotal()); } /*! diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 1f1dd851cdf..5898eebf6a9 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1413,18 +1413,6 @@ class CVariable { */ inline virtual void SetPrimitive(unsigned long iPoint, CConfig *config) {} - /*! - * \brief A virtual member. - * \param[in] Temperature_Wall - Value of the Temperature at the wall - */ - inline virtual void SetWallTemperature(unsigned long iPoint, su2double Temperature_Wall) {} - - /*! - * \brief A virtual member. - * \param[in] Temperature_Wall - Value of the Temperature at the wall - */ - inline virtual void SetWallTemperature(unsigned long iPoint, su2double* Temperature_Wall) {} - /*! * \brief Set the thermal coefficient. * \param[in] config - Configuration parameters. diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index d27405ae44a..75ea6dbcdb5 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -30,7 +30,8 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2double energy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) - : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config) { + : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config), + indices(ndim, 0) { const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 92805a7d11b..84d801f333d 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -30,7 +30,8 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double temperature, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) - : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config) { + : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config), + indices(ndim, 0) { const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); @@ -65,8 +66,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { - unsigned long iVar; - bool check_dens = false, check_temp = false, physical = true; + bool physical = true; /*--- Set the value of the pressure ---*/ @@ -74,8 +74,8 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel /*--- Set the value of the temperature directly ---*/ - su2double Temperature = Solution(iPoint,nDim+1); - check_temp = SetTemperature(iPoint,Temperature); + su2double Temperature = Solution(iPoint, nDim+1); + const auto check_temp = SetTemperature(iPoint, Temperature); /*--- Use the fluid model to compute the new value of density. Note that the thermodynamic pressure is constant and decoupled @@ -87,7 +87,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel /*--- Set the value of the density ---*/ - check_dens = SetDensity(iPoint, FluidModel->GetDensity()); + const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity()); /*--- Non-physical solution found. Revert to old values. ---*/ @@ -95,7 +95,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel /*--- Copy the old solution ---*/ - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0ul; iVar < nVar; iVar++) Solution(iPoint, iVar) = Solution_Old(iPoint, iVar); /*--- Recompute the primitive variables ---*/ diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index fc606d952bb..33c7c039cef 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -51,8 +51,7 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel) { - unsigned short iVar; - bool check_dens = false, check_temp = false, physical = true; + bool physical = true; /*--- Set the value of the pressure ---*/ @@ -60,8 +59,8 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Set the value of the temperature directly ---*/ - su2double Temperature = Solution(iPoint,nDim+1); - check_temp = SetTemperature(iPoint,Temperature); + su2double Temperature = Solution(iPoint, nDim+1); + const auto check_temp = SetTemperature(iPoint,Temperature); /*--- Use the fluid model to compute the new value of density. Note that the thermodynamic pressure is constant and decoupled @@ -73,7 +72,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Set the value of the density ---*/ - check_dens = SetDensity(iPoint, FluidModel->GetDensity()); + const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity()); /*--- Non-physical solution found. Revert to old values. ---*/ @@ -81,12 +80,12 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Copy the old solution ---*/ - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0ul; iVar < nVar; iVar++) Solution(iPoint,iVar) = Solution_Old(iPoint,iVar); /*--- Recompute the primitive variables ---*/ - Temperature = Solution(iPoint,nDim+1); + Temperature = Solution(iPoint, nDim+1); SetTemperature(iPoint, Temperature); FluidModel->SetTDState_T(Temperature); SetDensity(iPoint, FluidModel->GetDensity()); diff --git a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp index 19c24d2e1c0..aea7cebdfa6 100644 --- a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp +++ b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp @@ -41,6 +41,7 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure, const CConfig *config, CNEMOGas *fluidmodel) : CFlowVariable(npoint, ndim, nvar, nvarprim, nvarprimgrad, config), + indices(ndim, config->GetnSpecies()), implicit(config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT) { unsigned short iDim, iSpecies; diff --git a/SU2_CFD/src/variables/CNSVariable.cpp b/SU2_CFD/src/variables/CNSVariable.cpp index 8b9f8cec0bb..01a94284eda 100644 --- a/SU2_CFD/src/variables/CNSVariable.cpp +++ b/SU2_CFD/src/variables/CNSVariable.cpp @@ -66,9 +66,9 @@ void CNSVariable::SetRoe_Dissipation_NTS(unsigned long iPoint, /*--- Density ---*/ AD::SetPreaccIn(Solution(iPoint,0)); /*--- Laminar viscosity --- */ - AD::SetPreaccIn(Primitive(iPoint,nDim+5)); + AD::SetPreaccIn(Primitive(iPoint, indices.LaminarViscosity())); /*--- Eddy viscosity ---*/ - AD::SetPreaccIn(Primitive(iPoint,nDim+6)); + AD::SetPreaccIn(Primitive(iPoint, indices.EddyViscosity())); /*--- Central/upwind blending based on: * Zhixiang Xiao, Jian Liu, Jingbo Huang, and Song Fu. "Numerical @@ -111,15 +111,15 @@ void CNSVariable::SetRoe_Dissipation_FD(unsigned long iPoint, su2double val_wall AD::SetPreaccIn(Gradient_Primitive[iPoint], nVar, nDim); AD::SetPreaccIn(val_wall_dist); /*--- Eddy viscosity ---*/ - AD::SetPreaccIn(Primitive(iPoint,nDim+5)); + AD::SetPreaccIn(Primitive(iPoint, indices.EddyViscosity())); /*--- Laminar viscosity --- */ - AD::SetPreaccIn(Primitive(iPoint,nDim+6)); + AD::SetPreaccIn(Primitive(iPoint, indices.LaminarViscosity())); su2double uijuij = 0.0; for(unsigned long iDim = 0; iDim < nDim; ++iDim) for(unsigned long jDim = 0; jDim < nDim; ++jDim) - uijuij += pow(Gradient_Primitive(iPoint,1+iDim,jDim),2); + uijuij += pow(Gradient_Primitive(iPoint, indices.Velocity()+iDim, jDim),2); uijuij = max(sqrt(uijuij),1e-10); From ac3f474bfd5ffd25c1c608f790ae428b490ef2c4 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 12:43:23 +0000 Subject: [PATCH 06/12] using the indices in scalar convection --- .../numerics/scalar/scalar_convection.hpp | 20 +++++++------ .../numerics/turbulent/turb_convection.hpp | 21 +++++-------- SU2_CFD/src/drivers/CDriver.cpp | 30 ++++++++++++++----- 3 files changed, 40 insertions(+), 31 deletions(-) diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index f7e70877b57..8da8296f98c 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -47,8 +47,9 @@ template class CUpwScalar : public CNumerics { protected: - su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0 */ - su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0 */ + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ + su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */ su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */ su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */ su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */ @@ -56,7 +57,8 @@ class CUpwScalar : public CNumerics { const bool implicit = false, incompressible = false, dynamic_grid = false; /*! - * \brief A pure virtual function; Adds any extra variables to AD + * \brief A pure virtual function. Derived classes must use it to register the additional + * variables they use as preaccumulation inputs, e.g. the density for SST. */ virtual void ExtraADPreaccIn() = 0; @@ -76,6 +78,7 @@ class CUpwScalar : public CNumerics { */ CUpwScalar(unsigned short ndim, unsigned short nvar, const CConfig* config) : CNumerics(ndim, nvar, config), + idx(ndim, config->GetnSpecies()), implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), dynamic_grid(config->GetDynamic_Grid()) { @@ -117,22 +120,21 @@ class CUpwScalar : public CNumerics { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); } + AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); + AD::SetPreaccIn(&V_j[idx.Velocity()], nDim); ExtraADPreaccIn(); - Density_i = V_i[nDim + 2]; - Density_j = V_j[nDim + 2]; - su2double q_ij = 0.0; if (dynamic_grid) { for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double Velocity_i = V_i[iDim + 1] - GridVel_i[iDim]; - su2double Velocity_j = V_j[iDim + 1] - GridVel_j[iDim]; + su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim]; + su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim]; q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; } } else { for (unsigned short iDim = 0; iDim < nDim; iDim++) { - q_ij += 0.5 * (V_i[iDim + 1] + V_j[iDim + 1]) * Normal[iDim]; + q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim]; } } diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 3a80fbe11b6..4f24885cf44 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -40,9 +40,6 @@ template class CUpwSca_TurbSA final : public CUpwScalar { private: using Base = CUpwScalar; - using Base::nDim; - using Base::V_i; - using Base::V_j; using Base::a0; using Base::a1; using Base::Flux; @@ -53,12 +50,9 @@ class CUpwSca_TurbSA final : public CUpwScalar { using Base::implicit; /*! - * \brief Adds any extra variables to AD + * \brief Adds any extra variables to AD. */ - void ExtraADPreaccIn() override { - AD::SetPreaccIn(V_i, nDim+1); - AD::SetPreaccIn(V_j, nDim+1); - } + void ExtraADPreaccIn() override {} /*! * \brief SA specific steps in the ComputeResidual method @@ -104,16 +98,15 @@ class CUpwSca_TurbSST final : public CUpwScalar { using Base::Jacobian_j; using Base::ScalarVar_i; using Base::ScalarVar_j; - using Base::Density_i; - using Base::Density_j; using Base::implicit; + using Base::idx; /*! * \brief Adds any extra variables to AD */ void ExtraADPreaccIn() override { - AD::SetPreaccIn(V_i, nDim+3); - AD::SetPreaccIn(V_j, nDim+3); + AD::SetPreaccIn(V_i[idx.Density()]); + AD::SetPreaccIn(V_j[idx.Density()]); } /*! @@ -121,8 +114,8 @@ class CUpwSca_TurbSST final : public CUpwScalar { * \param[in] config - Definition of the particular problem. */ void FinishResidualCalc(const CConfig* config) override { - Flux[0] = a0*Density_i*ScalarVar_i[0] + a1*Density_j*ScalarVar_j[0]; - Flux[1] = a0*Density_i*ScalarVar_i[1] + a1*Density_j*ScalarVar_j[1]; + Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0]; + Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1]; if (implicit) { Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index b0fcf7786b9..2a483f5735c 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -51,6 +51,10 @@ #include "../../include/interfaces/fsi/CFlowTractionInterface.hpp" #include "../../include/interfaces/fsi/CDiscAdjFlowTractionInterface.hpp" +#include "../../include/variables/CEulerVariable.hpp" +#include "../../include/variables/CIncEulerVariable.hpp" +#include "../../include/variables/CNEMOEulerVariable.hpp" + #include "../../include/numerics/template.hpp" #include "../../include/numerics/transition.hpp" #include "../../include/numerics/radiation.hpp" @@ -1309,8 +1313,14 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, } } /*--- Explicit instantiation of the template above, needed because it is defined in a cpp file, instead of hpp. ---*/ -template void CDriver::InstantiateTurbulentNumerics(unsigned short, int, const CConfig*, - const CSolver*, CNumerics****&) const; +template void CDriver::InstantiateTurbulentNumerics>( + unsigned short, int, const CConfig*, const CSolver*, CNumerics****&) const; + +template void CDriver::InstantiateTurbulentNumerics>( + unsigned short, int, const CConfig*, const CSolver*, CNumerics****&) const; + +template void CDriver::InstantiateTurbulentNumerics>( + unsigned short, int, const CConfig*, const CSolver*, CNumerics****&) const; void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSolver ***solver, CNumerics ****&numerics) const { @@ -1765,7 +1775,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Solver definition for the Potential, Euler, Navier-Stokes NEMO problems ---*/ - if ((NEMO_euler) || (NEMO_ns)) { + if (NEMO_euler || NEMO_ns) { /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Flow()) { @@ -1918,11 +1928,15 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Solver definition for the turbulent model problem ---*/ if (turbulent) { - /// TODO: Replace the dummy int type here for the real thing commented below. - InstantiateTurbulentNumerics(nVar_Turb, offset, config, solver[MESH_0][TURB_SOL], numerics); -// if (incompressible) InstantiateTurbulentNumerics(...); -// else if (nemo) InstantiateTurbulentNumerics(...); -// else InstantiateTurbulentNumerics(...); + if (incompressible) + InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + solver[MESH_0][TURB_SOL], numerics); + else if (NEMO_ns) + InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + solver[MESH_0][TURB_SOL], numerics); + else + InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + solver[MESH_0][TURB_SOL], numerics); } /*--- Solver definition for the transition model problem ---*/ From 9790d119149fbe557f16a5d1fdc6b6352677dc83 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 15:32:53 +0000 Subject: [PATCH 07/12] viscous scalar and turbulent numerics --- .../numerics/scalar/scalar_diffusion.hpp | 73 +++++++- .../numerics/turbulent/turb_diffusion.hpp | 169 +++++++++++++++--- SU2_CFD/obj/Makefile.am | 4 - SU2_CFD/src/drivers/CDriver.cpp | 13 +- SU2_CFD/src/meson.build | 2 - .../src/numerics/scalar/scalar_convection.cpp | 1 - .../src/numerics/scalar/scalar_diffusion.cpp | 108 ----------- .../numerics/turbulent/turb_convection.cpp | 1 - .../src/numerics/turbulent/turb_diffusion.cpp | 150 ---------------- 9 files changed, 219 insertions(+), 302 deletions(-) delete mode 100644 SU2_CFD/src/numerics/scalar/scalar_convection.cpp delete mode 100644 SU2_CFD/src/numerics/scalar/scalar_diffusion.cpp delete mode 100644 SU2_CFD/src/numerics/turbulent/turb_convection.cpp delete mode 100644 SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 6f7800e3989..04c96d6de42 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -44,14 +44,18 @@ * \ingroup ViscDiscr * \author C. Pederson, A. Bueno, and F. Palacios */ +template class CAvgGrad_Scalar : public CNumerics { protected: - su2double* Proj_Mean_GradScalarVar_Normal = nullptr; /*!< \brief Mean_gradScalarVar DOT normal. */ - su2double* Proj_Mean_GradScalarVar = nullptr; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ - su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */ - su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */ - su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */ + enum : unsigned short {MAXNVAR = 8}; + + const FlowIndices idx; + su2double Proj_Mean_GradScalarVar_Normal[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal. */ + su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ + su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ + su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ + su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ + su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ const bool correct_gradient = false, implicit = false, incompressible = false; @@ -75,17 +79,68 @@ class CAvgGrad_Scalar : public CNumerics { * \param[in] correct_gradient - Whether to correct gradient for skewness. * \param[in] config - Definition of the particular problem. */ - CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_gradient, const CConfig* config); + CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, + const CConfig* config) + : CNumerics(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()), + correct_gradient(correct_grad), + implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), + incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) { + if (nVar > MAXNVAR) { + SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION); + } + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double[nVar](); + Jacobian_j[iVar] = new su2double[nVar](); + } + } /*! * \brief Destructor of the class. */ - ~CAvgGrad_Scalar(void) override; + ~CAvgGrad_Scalar() override { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + delete[] Jacobian_i[iVar]; + delete[] Jacobian_j[iVar]; + } + } /*! * \brief Compute the viscous residual using an average of gradients without correction. * \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; + ResidualType<> ComputeResidual(const CConfig* config) final { + AD::StartPreacc(); + AD::SetPreaccIn(Coord_i, nDim); + AD::SetPreaccIn(Coord_j, nDim); + AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim); + if (correct_gradient) { + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_j, nVar); + } + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); + AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]); + + Density_i = V_i[idx.Density()]; + Density_j = V_j[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; + Laminar_Viscosity_j = V_j[idx.LaminarViscosity()]; + Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; + Eddy_Viscosity_j = V_j[idx.EddyViscosity()]; + + ExtraADPreaccIn(); + + proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j, + correct_gradient, ScalarVar_i, ScalarVar_j, + Proj_Mean_GradScalarVar_Normal, Proj_Mean_GradScalarVar); + FinishResidualCalc(config); + + AD::SetPreaccOut(Flux, nVar); + AD::EndPreacc(); + + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); + } }; diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index 15ad07944eb..143dcdf26ea 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -25,32 +25,61 @@ * You should have received a copy of the GNU Lesser General Public * License along with SU2. If not, see . */ - #pragma once #include "../scalar/scalar_diffusion.hpp" - /*! * \class CAvgGrad_TurbSA * \brief Class for computing viscous term using average of gradients (Spalart-Allmaras Turbulence model). * \ingroup ViscDiscr * \author A. Bueno. */ -class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { +template +class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { private: + using Base = CAvgGrad_Scalar; + using Base::Laminar_Viscosity_i; + using Base::Laminar_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::implicit; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + const su2double sigma = 2.0/3.0; /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn(void) override; + 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; + void FinishResidualCalc(const CConfig* config) override { + /*--- Compute mean effective viscosity ---*/ + + su2double nu_i = Laminar_Viscosity_i/Density_i; + su2double nu_j = Laminar_Viscosity_j/Density_j; + su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); + + /// TODO: I smell bug... + Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; + + /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ + + if (implicit) { + Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; + Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; + } + } public: /*! @@ -61,7 +90,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config); + bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) {} }; /*! @@ -70,21 +100,65 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { * \ingroup ViscDiscr * \author F. Palacios */ -class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { +template +class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { private: + using Base = CAvgGrad_Scalar; + using Base::Laminar_Viscosity_i; + using Base::Laminar_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_Mean_GradScalarVar_Normal; + using Base::proj_vector_ij; + using Base::implicit; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + const su2double sigma = 2.0/3.0; const su2double cn1 = 16.0; /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn(void) override; + void ExtraADPreaccIn() override {} /*! - * \brief SA specific steps in the ComputeResidual method + * \brief SA-neg specific steps in the ComputeResidual method * \param[in] config - Definition of the particular problem. */ - void FinishResidualCalc(const CConfig* config) override; + void FinishResidualCalc(const CConfig* config) override { + /*--- Compute mean effective viscosity ---*/ + + su2double nu_i = Laminar_Viscosity_i/Density_i; + su2double nu_j = Laminar_Viscosity_j/Density_j; + + su2double nu_ij = 0.5*(nu_i+nu_j); + su2double nu_tilde_ij = 0.5*(ScalarVar_i[0] + ScalarVar_j[0]); + + su2double nu_e; + + if (nu_tilde_ij > 0.0) { + nu_e = nu_ij + nu_tilde_ij; + } + else { + su2double Xi = nu_tilde_ij/nu_ij; + su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi); + nu_e = nu_ij + fn*nu_tilde_ij; + } + + Flux[0] = nu_e*Proj_Mean_GradScalarVar_Normal[0]/sigma; + + /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ + + if (implicit) { + Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; + Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; + } + } public: /*! @@ -95,7 +169,8 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSA_Neg(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config); + bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) {} }; /*! @@ -104,26 +179,73 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { * \ingroup ViscDiscr * \author A. Bueno. */ -class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { +template +class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { private: - const su2double - sigma_k1 = 0.0, /*!< \brief Constants for the viscous terms, k-w (1), k-eps (2)*/ - sigma_k2 = 0.0, - sigma_om1 = 0.0, - sigma_om2 = 0.0; + 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::implicit; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + + const su2double sigma_k1; /*!< \brief Constants for the viscous terms, k-w (1), k-eps (2)*/ + const su2double sigma_k2; + const su2double sigma_om1; + const su2double sigma_om2; su2double F1_i, F1_j; /*!< \brief Menter's first blending function */ /*! * \brief Adds any extra variables to AD */ - void ExtraADPreaccIn(void) override; + void ExtraADPreaccIn() override { + AD::SetPreaccIn(F1_i, F1_j); + } /*! * \brief SST specific steps in the ComputeResidual method * \param[in] config - Definition of the particular problem. */ - void FinishResidualCalc(const CConfig* config) override; + void FinishResidualCalc(const CConfig* config) override { + /*--- Compute the blended constant for the viscous terms ---*/ + const su2double sigma_kine_i = F1_i*sigma_k1 + (1.0 - F1_i)*sigma_k2; + const su2double sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2; + const su2double sigma_omega_i = F1_i*sigma_om1 + (1.0 - F1_i)*sigma_om2; + const su2double sigma_omega_j = F1_j*sigma_om1 + (1.0 - F1_j)*sigma_om2; + + /*--- Compute mean effective dynamic viscosity ---*/ + const su2double diff_i_kine = Laminar_Viscosity_i + sigma_kine_i*Eddy_Viscosity_i; + const su2double diff_j_kine = Laminar_Viscosity_j + sigma_kine_j*Eddy_Viscosity_j; + const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i; + const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j; + + const su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine); + const su2double diff_omega = 0.5*(diff_i_omega + diff_j_omega); + + Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; + Flux[1] = diff_omega*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_kine*proj_on_rho_i; Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_on_rho_i; + + const su2double proj_on_rho_j = proj_vector_ij/Density_j; + Jacobian_j[0][0] = diff_kine*proj_on_rho_j; Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_on_rho_j; + } + } public: /*! @@ -135,7 +257,13 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { * \param[in] config - Definition of the particular problem. */ CAvgGrad_TurbSST(unsigned short val_nDim, unsigned short val_nVar, - const su2double* constants, bool correct_grad, const CConfig* config); + const su2double* constants, bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), + sigma_k1(constants[0]), + sigma_k2(constants[1]), + sigma_om1(constants[2]), + sigma_om2(constants[3]) { + } /*! * \brief Sets value of first blending function. @@ -143,5 +271,4 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { void SetF1blending(su2double val_F1_i, su2double val_F1_j) override { F1_i = val_F1_i; F1_j = val_F1_j; } - }; diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 9f8c2b1eecb..4400cb6f3be 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -106,11 +106,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/numerics/continuous_adjoint/adj_convection.cpp \ ../src/numerics/continuous_adjoint/adj_diffusion.cpp \ ../src/numerics/continuous_adjoint/adj_sources.cpp \ - ../src/numerics/scalar/scalar_convection.cpp \ - ../src/numerics/scalar/scalar_diffusion.cpp \ ../src/numerics/scalar/scalar_sources.cpp \ - ../src/numerics/turbulent/turb_convection.cpp \ - ../src/numerics/turbulent/turb_diffusion.cpp \ ../src/numerics/turbulent/turb_sources.cpp \ ../src/numerics/elasticity/CFEAElasticity.cpp \ ../src/numerics/elasticity/CFEALinearElasticity.cpp \ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 2a483f5735c..992a5a84822 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1266,12 +1266,12 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, true, config); } else if (neg_spalart_allmaras) - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, true, config); else if (menter_sst) - numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); + numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, true, config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1300,15 +1300,16 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA(nDim, nVar_Turb, false, config); } else if (neg_spalart_allmaras) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg(nDim, nVar_Turb, false, config); } else if (menter_sst) { numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST(nDim, nVar_Turb, config); - numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, config); + numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSST(nDim, nVar_Turb, constants, false, + config); } } } diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 8d444ab8054..ba1ce1e9e9d 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -134,9 +134,7 @@ su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/continuous_adjoint/adj_convection.cpp', 'numerics/continuous_adjoint/adj_diffusion.cpp', 'numerics/continuous_adjoint/adj_sources.cpp', - 'numerics/scalar/scalar_diffusion.cpp', 'numerics/scalar/scalar_sources.cpp', - 'numerics/turbulent/turb_diffusion.cpp', 'numerics/turbulent/turb_sources.cpp', 'numerics/elasticity/CFEAElasticity.cpp', 'numerics/elasticity/CFEALinearElasticity.cpp', diff --git a/SU2_CFD/src/numerics/scalar/scalar_convection.cpp b/SU2_CFD/src/numerics/scalar/scalar_convection.cpp deleted file mode 100644 index 8b137891791..00000000000 --- a/SU2_CFD/src/numerics/scalar/scalar_convection.cpp +++ /dev/null @@ -1 +0,0 @@ - diff --git a/SU2_CFD/src/numerics/scalar/scalar_diffusion.cpp b/SU2_CFD/src/numerics/scalar/scalar_diffusion.cpp deleted file mode 100644 index 8d62e778bc0..00000000000 --- a/SU2_CFD/src/numerics/scalar/scalar_diffusion.cpp +++ /dev/null @@ -1,108 +0,0 @@ -/*! - * \file scalar_diffusion.cpp - * \brief Implementation of numerics classes to compute viscous - * fluxes in scalar problems. - * \author F. Palacios, T. Economon - * \version 7.2.1 "Blackbird" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2021, 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/numerics/scalar/scalar_diffusion.hpp" - -CAvgGrad_Scalar::CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, - const CConfig* config) - : CNumerics(val_nDim, val_nVar, config), - correct_gradient(correct_grad), - implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), - incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) { - Proj_Mean_GradScalarVar_Normal = new su2double[nVar](); - Proj_Mean_GradScalarVar = new su2double[nVar](); - - Flux = new su2double[nVar](); - Jacobian_i = new su2double*[nVar]; - Jacobian_j = new su2double*[nVar]; - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double[nVar](); - Jacobian_j[iVar] = new su2double[nVar](); - } -} - -CAvgGrad_Scalar::~CAvgGrad_Scalar(void) { - delete[] Proj_Mean_GradScalarVar_Normal; - delete[] Proj_Mean_GradScalarVar; - - delete[] Flux; - if (Jacobian_i != nullptr) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - delete[] Jacobian_i[iVar]; - delete[] Jacobian_j[iVar]; - } - delete[] Jacobian_i; - delete[] Jacobian_j; - } -} - -CNumerics::ResidualType<> CAvgGrad_Scalar::ComputeResidual(const CConfig* config) { - AD::StartPreacc(); - AD::SetPreaccIn(Coord_i, nDim); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); - AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim); - if (correct_gradient) { - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_j, nVar); - } - ExtraADPreaccIn(); - - if (incompressible) { - AD::SetPreaccIn(V_i, nDim + 6); - AD::SetPreaccIn(V_j, nDim + 6); - - Density_i = V_i[nDim + 2]; - Density_j = V_j[nDim + 2]; - Laminar_Viscosity_i = V_i[nDim + 4]; - Laminar_Viscosity_j = V_j[nDim + 4]; - Eddy_Viscosity_i = V_i[nDim + 5]; - Eddy_Viscosity_j = V_j[nDim + 5]; - } else { - AD::SetPreaccIn(V_i, nDim + 7); - AD::SetPreaccIn(V_j, nDim + 7); - - Density_i = V_i[nDim + 2]; - Density_j = V_j[nDim + 2]; - Laminar_Viscosity_i = V_i[nDim + 5]; - Laminar_Viscosity_j = V_j[nDim + 5]; - Eddy_Viscosity_i = V_i[nDim + 6]; - Eddy_Viscosity_j = V_j[nDim + 6]; - } - - proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j, - correct_gradient, ScalarVar_i, ScalarVar_j, Proj_Mean_GradScalarVar_Normal, - Proj_Mean_GradScalarVar); - FinishResidualCalc(config); - - AD::SetPreaccOut(Flux, nVar); - AD::EndPreacc(); - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); -} diff --git a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp b/SU2_CFD/src/numerics/turbulent/turb_convection.cpp deleted file mode 100644 index 8b137891791..00000000000 --- a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp +++ /dev/null @@ -1 +0,0 @@ - diff --git a/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp b/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp deleted file mode 100644 index ad1c8459914..00000000000 --- a/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp +++ /dev/null @@ -1,150 +0,0 @@ -/*! - * \file turb_diffusion.cpp - * \brief Implementation of numerics classes to compute viscous - * fluxes in turbulence problems. - * \author F. Palacios, T. Economon - * \version 7.2.1 "Blackbird" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2021, 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/numerics/turbulent/turb_diffusion.hpp" - - -CAvgGrad_TurbSA::CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config) : - CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) { } - -void CAvgGrad_TurbSA::ExtraADPreaccIn() { } - -void CAvgGrad_TurbSA::FinishResidualCalc(const CConfig* config) { - - /*--- Compute mean effective viscosity ---*/ - - su2double nu_i = Laminar_Viscosity_i/Density_i; - su2double nu_j = Laminar_Viscosity_j/Density_j; - su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); - - Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; - - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ - - if (implicit) { - Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; - Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; - } - -} - -CAvgGrad_TurbSA_Neg::CAvgGrad_TurbSA_Neg(unsigned short val_nDim, - unsigned short val_nVar, - bool correct_grad, - const CConfig* config) : - CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) { } - -void CAvgGrad_TurbSA_Neg::ExtraADPreaccIn() { } - -void CAvgGrad_TurbSA_Neg::FinishResidualCalc(const CConfig* config) { - - /*--- Compute mean effective viscosity ---*/ - - su2double nu_i = Laminar_Viscosity_i/Density_i; - su2double nu_j = Laminar_Viscosity_j/Density_j; - - su2double nu_ij = 0.5*(nu_i+nu_j); - su2double nu_tilde_ij = 0.5*(ScalarVar_i[0]+ScalarVar_j[0]); - - su2double nu_e; - - if (nu_tilde_ij > 0.0) { - nu_e = nu_ij + nu_tilde_ij; - } - else { - su2double Xi = nu_tilde_ij/nu_ij; - su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi); - nu_e = nu_ij + fn*nu_tilde_ij; - } - - Flux[0] = nu_e*Proj_Mean_GradScalarVar_Normal[0]/sigma; - - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ - - if (implicit) { - Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; - Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; - } - -} - -CAvgGrad_TurbSST::CAvgGrad_TurbSST(unsigned short val_nDim, - unsigned short val_nVar, - const su2double *constants, - bool correct_grad, - const CConfig* config) : - CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), - sigma_k1(constants[0]), - sigma_k2(constants[1]), - sigma_om1(constants[2]), - sigma_om2(constants[3]) { - -} - -void CAvgGrad_TurbSST::ExtraADPreaccIn() { - AD::SetPreaccIn(F1_i); AD::SetPreaccIn(F1_j); -} - -void CAvgGrad_TurbSST::FinishResidualCalc(const CConfig* config) { - - su2double sigma_kine_i, sigma_kine_j, sigma_omega_i, sigma_omega_j; - su2double diff_i_kine, diff_i_omega, diff_j_kine, diff_j_omega; - - /*--- Compute the blended constant for the viscous terms ---*/ - sigma_kine_i = F1_i*sigma_k1 + (1.0 - F1_i)*sigma_k2; - sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2; - sigma_omega_i = F1_i*sigma_om1 + (1.0 - F1_i)*sigma_om2; - sigma_omega_j = F1_j*sigma_om1 + (1.0 - F1_j)*sigma_om2; - - /*--- Compute mean effective dynamic viscosity ---*/ - diff_i_kine = Laminar_Viscosity_i + sigma_kine_i*Eddy_Viscosity_i; - diff_j_kine = Laminar_Viscosity_j + sigma_kine_j*Eddy_Viscosity_j; - diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i; - diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j; - - su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine); - su2double diff_omega = 0.5*(diff_i_omega + diff_j_omega); - - Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; - Flux[1] = diff_omega*Proj_Mean_GradScalarVar[1]; - - /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ - if (implicit) { - su2double proj_on_rho = proj_vector_ij/Density_i; - - Jacobian_i[0][0] = -diff_kine*proj_on_rho; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_on_rho; - - proj_on_rho = proj_vector_ij/Density_j; - - Jacobian_j[0][0] = diff_kine*proj_on_rho; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_on_rho; - } - -} From 60eed7c317eda3aaba55e2cd4c87690bbd1b9bab Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 16:48:32 +0000 Subject: [PATCH 08/12] fix old bug in SA-neg diffusion and abstract "velocity gradient" --- SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp | 4 +--- SU2_CFD/include/variables/CEulerVariable.hpp | 8 ++++++++ SU2_CFD/include/variables/CIncEulerVariable.hpp | 8 ++++++++ SU2_CFD/include/variables/CNEMOEulerVariable.hpp | 8 ++++++++ SU2_CFD/include/variables/CVariable.hpp | 9 +++++++++ SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp | 2 +- SU2_CFD/src/output/CFlowCompOutput.cpp | 2 +- SU2_CFD/src/output/CFlowIncOutput.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 2 +- SU2_CFD/src/solvers/CTurbSASolver.cpp | 6 +++--- 10 files changed, 41 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index 143dcdf26ea..40eac034996 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -70,7 +70,6 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { su2double nu_j = Laminar_Viscosity_j/Density_j; su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); - /// TODO: I smell bug... Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ @@ -111,7 +110,6 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { using Base::ScalarVar_i; using Base::ScalarVar_j; using Base::Proj_Mean_GradScalarVar; - using Base::Proj_Mean_GradScalarVar_Normal; using Base::proj_vector_ij; using Base::implicit; using Base::Flux; @@ -150,7 +148,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { nu_e = nu_ij + fn*nu_tilde_ij; } - Flux[0] = nu_e*Proj_Mean_GradScalarVar_Normal[0]/sigma; + Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 31f5a68bb6f..7cc3deedf33 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -229,6 +229,14 @@ class CEulerVariable : public CFlowVariable { return Primitive(iPoint,iDim+indices.Velocity()); } + /*! + * \brief Get the velocity gradient. + * \return Value of the velocity gradient. + */ + inline CMatrixView GetVelocityGradient(unsigned long iPoint) const final { + return Gradient_Primitive(iPoint, indices.Velocity()); + } + /*! * \brief Get the projected velocity in a unitary vector direction (compressible solver). * \param[in] val_vector - Direction of projection. diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 855aa09694c..b750d5fc846 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -153,6 +153,14 @@ class CIncEulerVariable : public CFlowVariable { return Primitive(iPoint, iDim+indices.Velocity()); } + /*! + * \brief Get the velocity gradient. + * \return Value of the velocity gradient. + */ + inline CMatrixView GetVelocityGradient(unsigned long iPoint) const final { + return Gradient_Primitive(iPoint, indices.Velocity()); + } + /*! * \brief Get the projected velocity in a unitary vector direction (compressible solver). * \param[in] val_vector - Direction of projection. diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index 9037011884d..2c52c6bb5b9 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -235,6 +235,14 @@ class CNEMOEulerVariable : public CFlowVariable { */ inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { return Primitive(iPoint,VEL_INDEX+iDim); } + /*! + * \brief Get the velocity gradient. + * \return Value of the velocity gradient. + */ + inline CMatrixView GetVelocityGradient(unsigned long iPoint) const final { + return Gradient_Primitive(iPoint, indices.Velocity()); + } + /*! * \brief Get the projected velocity in a unitary vector direction (compressible solver). * \param[in] val_vector - Direction of projection. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 5898eebf6a9..e2a003710bb 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1042,6 +1042,15 @@ class CVariable { */ inline virtual su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \return Value of the velocity gradient. + */ + inline virtual CMatrixView GetVelocityGradient(unsigned long iPoint) const { + return CMatrixView(); + } + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index da75e7397e3..8494cc02512 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -189,7 +189,7 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); su2double tau[3][3]; - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow,1), Viscosity); + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { // Viscous component in the tn vector --> Units of force (non-dimensional). diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index d9effc33b46..79a5a3228e1 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -557,7 +557,7 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv } else { SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } - SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1))); + SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetVelocityGradient(iPoint))); } LoadCommonFVMOutputs(config, geometry, iPoint); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 5e460daf22b..df1693ee2ad 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -659,7 +659,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve } else { SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } - SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1))); + SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetVelocityGradient(iPoint))); } // Streamwise Periodicity diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 97d69a0fc6a..624210748e9 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3984,7 +3984,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) if (viscous_flow) { su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); su2double Tau[3][3]; - CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity); + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity); for (iDim = 0; iDim < nDim; iDim++) { auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], iNormal); } diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 1c8fbec34c4..737b5c2ce27 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -1248,7 +1248,7 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde); visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), - nodes->GetGradient(iPoint)); + nodes->GetGradient(iPoint)); /*--- Compute residual, and Jacobians ---*/ @@ -1364,7 +1364,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC const auto coord_i = geometry->nodes->GetCoord(iPoint); const auto nNeigh = geometry->nodes->GetnPoint(iPoint); const auto wallDistance = geometry->nodes->GetWall_Distance(iPoint); - const auto primVarGrad = flowNodes->GetGradient_Primitive(iPoint); + const auto velocityGrad = flowNodes->GetVelocityGradient(iPoint); const auto vorticity = flowNodes->GetVorticity(iPoint); const auto density = flowNodes->GetDensity(iPoint); const auto laminarViscosity = flowNodes->GetLaminarViscosity(iPoint); @@ -1375,7 +1375,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC su2double uijuij = 0.0; for(auto iDim = 0u; iDim < nDim; iDim++){ for(auto jDim = 0u; jDim < nDim; jDim++){ - uijuij += primVarGrad[1+iDim][jDim]*primVarGrad[1+iDim][jDim]; + uijuij += pow(velocityGrad[iDim][jDim], 2); } } uijuij = sqrt(fabs(uijuij)); From 36282d70597d9cb85046f62c5cb182a8be2a6f86 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 22:03:55 +0000 Subject: [PATCH 09/12] update regressions after bug fix --- TestCases/hybrid_regression.py | 2 +- TestCases/hybrid_regression_AD.py | 2 +- TestCases/parallel_regression.py | 4 ++-- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 2 +- TestCases/serial_regression_AD.py | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index e164cd26953..ee54aa341e4 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -183,7 +183,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.531271, -14.899968, 1.064330, 0.019756, 20, -1.825644, 0, -12.818150] + turb_naca0012_sa.test_vals = [-8.627052, -10.377936, 1.064491, 0.019710, 20.000000, -1.763095, 20.000000, -4.794176] test_list.append(turb_naca0012_sa) # NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253) diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index fd74fd0226e..e987e96fc6e 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -75,7 +75,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.230632, 0.696530, 0.177890, -0.000016, 5, -3.007652, 5, -7.666352] + discadj_rans_naca0012_sa.test_vals = [-2.230631, 0.644953, 0.177890, -0.000016, 5.000000, -3.007652, 5.000000, -7.631910] test_list.append(discadj_rans_naca0012_sa) # Adjoint turbulent NACA0012 SST diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 057c10dd359..8f6bcd72aa4 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -330,7 +330,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.147929, -14.466781, 1.064330, 0.0197560, 20, -1.741326, 20, -4.864272] + turb_naca0012_sa.test_vals = [-8.621456, -10.378269, 1.064502, 0.019710, 20.000000, -1.811700, 20.000000, -5.171326] turb_naca0012_sa.su2_exec = "parallel_computation.py -f" turb_naca0012_sa.timeout = 3200 turb_naca0012_sa.tol = 0.00001 @@ -641,7 +641,7 @@ def main(): turbmod_sa_neg_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_neg_rae2822.cfg_file = "turb_SA_NEG_RAE2822.cfg" turbmod_sa_neg_rae2822.test_iter = 20 - turbmod_sa_neg_rae2822.test_vals = [-2.004688, 0.742320, 0.497307, -5.265640, 0.809478, 0.061986] + turbmod_sa_neg_rae2822.test_vals = [-2.004689, 0.742306, 0.497308, -5.265793, 0.809463, 0.062016] turbmod_sa_neg_rae2822.su2_exec = "mpirun -n 2 SU2_CFD" turbmod_sa_neg_rae2822.timeout = 1600 turbmod_sa_neg_rae2822.new_output = True diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index ef74841a3b9..fbabc4483df 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -95,7 +95,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.230578, 0.696567, 0.181590, -0.000018, 5, -3.421214, 5, -6.798877] + discadj_rans_naca0012_sa.test_vals = [-2.230578, 0.645001, 0.181590, -0.000018, 5.000000, -3.421214, 5.000000, -6.769609] discadj_rans_naca0012_sa.su2_exec = "parallel_computation.py -f" discadj_rans_naca0012_sa.timeout = 1600 discadj_rans_naca0012_sa.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 3a2e775ca79..f6ce7708bfb 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -355,7 +355,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.133933, -14.498178, 1.064330, 0.019756, 20, -1.898609, 20, -2.719825] + turb_naca0012_sa.test_vals = [-8.629583, -10.377793, 1.064488, 0.019711, 20.000000, -2.173971, 20.000000, -5.213344] turb_naca0012_sa.su2_exec = "SU2_CFD" turb_naca0012_sa.new_output = True turb_naca0012_sa.timeout = 3200 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 82869238f45..2ecd3e84c62 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -84,7 +84,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.230556, 0.696586, 0.180740, -0.000018, 5, -4.275186, 5, -8.926013] #last 8 columns + discadj_rans_naca0012_sa.test_vals = [-2.230555, 0.645023, 0.180740, -0.000018, 5.000000, -4.275184, 5.000000, -8.892454] #last 8 columns discadj_rans_naca0012_sa.su2_exec = "SU2_CFD_AD" discadj_rans_naca0012_sa.timeout = 1600 discadj_rans_naca0012_sa.tol = 0.00001 From 40c558647ac2563e2b2a6c423ff85fe42cf4175e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 7 Nov 2021 23:00:08 +0000 Subject: [PATCH 10/12] turbulence sources --- .../numerics/turbulent/turb_sources.hpp | 23 +- SU2_CFD/src/drivers/CDriver.cpp | 13 +- .../src/numerics/turbulent/turb_sources.cpp | 222 ++++++++++-------- 3 files changed, 149 insertions(+), 109 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 9f3ebbfea16..d4a0517f057 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -125,8 +125,11 @@ class CSourceBase_TurbSA : public CNumerics { * \ingroup SourceDiscr * \author A. Bueno. */ +template class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { private: + const FlowIndices idx; + su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; su2double norm2_Grad; @@ -135,7 +138,7 @@ class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { unsigned short iDim; bool transition; bool axisymmetric; - + public: /*! * \brief Constructor of the class. @@ -161,8 +164,11 @@ class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { * \author E.Molina, A. Bueno. * \version 7.2.1 "Blackbird" */ +template class CSourcePieceWise_TurbSA_COMP final : public CSourceBase_TurbSA { private: + const FlowIndices idx; + su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; su2double norm2_Grad; @@ -196,8 +202,11 @@ class CSourcePieceWise_TurbSA_COMP final : public CSourceBase_TurbSA { * \author E.Molina, A. Bueno. * \version 7.2.1 "Blackbird" */ +template class CSourcePieceWise_TurbSA_E final : public CSourceBase_TurbSA { private: + const FlowIndices idx; + su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; su2double norm2_Grad; @@ -229,8 +238,11 @@ class CSourcePieceWise_TurbSA_E final : public CSourceBase_TurbSA { * \author E.Molina, A. Bueno. * \version 7.2.1 "Blackbird" */ +template class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { private: + const FlowIndices idx; + su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; su2double norm2_Grad; @@ -238,7 +250,6 @@ class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { su2double dr, dg, dfw; su2double Sbar; su2double aux_cc, CompCorrection, c5; - unsigned short jDim; public: /*! @@ -263,8 +274,11 @@ class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { * \ingroup SourceDiscr * \author F. Palacios */ +template class CSourcePieceWise_TurbSA_Neg : public CSourceBase_TurbSA { private: + const FlowIndices idx; + su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; su2double norm2_Grad; @@ -295,8 +309,11 @@ class CSourcePieceWise_TurbSA_Neg : public CSourceBase_TurbSA { * \ingroup SourceDiscr * \author A. Campos. */ +template class CSourcePieceWise_TurbSST final : public CNumerics { private: + const FlowIndices idx; + su2double F1_i, F1_j, F2_i, @@ -334,7 +351,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*! * \brief Add contribution due to axisymmetric formulation to 2D residual */ - inline void ResidualAxisymmetric(su2double alfa_blended, su2double zeta){ + inline void ResidualAxisymmetric(su2double alfa_blended, su2double zeta) { if (Coord_i[1] < EPS) return; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 992a5a84822..3eb9d78a69e 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1280,17 +1280,18 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset, auto& turb_source_first_term = numerics[iMGlevel][TURB_SOL][source_first_term]; if (spalart_allmaras) - turb_source_first_term = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); + turb_source_first_term = new CSourcePieceWise_TurbSA(nDim, nVar_Turb, config); else if (e_spalart_allmaras) - turb_source_first_term = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); + turb_source_first_term = new CSourcePieceWise_TurbSA_E(nDim, nVar_Turb, config); else if (comp_spalart_allmaras) - turb_source_first_term = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); + turb_source_first_term = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config); else if (e_comp_spalart_allmaras) - turb_source_first_term = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); + turb_source_first_term = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config); else if (neg_spalart_allmaras) - turb_source_first_term = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); + turb_source_first_term = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config); else if (menter_sst) - turb_source_first_term = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config); + turb_source_first_term = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, + config); numerics[iMGlevel][TURB_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Turb, config); } diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index 37965ef34fa..5295b86c6e2 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -28,6 +28,10 @@ #include "../../../include/numerics/turbulent/turb_sources.hpp" +#include "../../../include/variables/CEulerVariable.hpp" +#include "../../../include/variables/CIncEulerVariable.hpp" +#include "../../../include/variables/CNEMOEulerVariable.hpp" + CSourceBase_TurbSA::CSourceBase_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : @@ -58,15 +62,18 @@ CSourceBase_TurbSA::CSourceBase_TurbSA(unsigned short val_nDim, } -CSourcePieceWise_TurbSA::CSourcePieceWise_TurbSA(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config) { +template +CSourcePieceWise_TurbSA::CSourcePieceWise_TurbSA(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CSourceBase_TurbSA(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()) { transition = (config->GetKind_Trans_Model() == BC); } -CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig* config) { // AD::StartPreacc(); // AD::SetPreaccIn(V_i, nDim+6); @@ -79,14 +86,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig // Set the boolean here depending on whether the point is closest to a rough wall or not. roughwall = (roughness_i > 0.0); - if (incompressible) { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - } - else { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; Residual = 0.0; Production = 0.0; @@ -215,12 +216,16 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig } -CSourcePieceWise_TurbSA_COMP::CSourcePieceWise_TurbSA_COMP(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), c5(3.5) { } +template +CSourcePieceWise_TurbSA_COMP::CSourcePieceWise_TurbSA_COMP(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CSourceBase_TurbSA(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()), + c5(3.5) { } -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CConfig* config) { // AD::StartPreacc(); // AD::SetPreaccIn(V_i, nDim+6); @@ -230,14 +235,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CC // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - if (incompressible) { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - } - else { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; Residual = 0.0; Production = 0.0; @@ -297,12 +296,14 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CC Residual = Production - Destruction + CrossProduction; /*--- Compressibility Correction term ---*/ - Pressure_i = V_i[nDim+1]; + Pressure_i = V_i[idx.Pressure()]; SoundSpeed_i = sqrt(Pressure_i*Gamma/Density_i); aux_cc=0; for(iDim=0;iDim CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CC } -CSourcePieceWise_TurbSA_E::CSourcePieceWise_TurbSA_E(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config) { } +template +CSourcePieceWise_TurbSA_E::CSourcePieceWise_TurbSA_E(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CSourceBase_TurbSA(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()) { } -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E::ComputeResidual(const CConfig* config) { unsigned short iDim, jDim; @@ -352,14 +356,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E::ComputeResidual(const CConf // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - if (incompressible) { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - } - else { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; Residual = 0.0; Production = 0.0; @@ -379,9 +377,13 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E::ComputeResidual(const CConf Sbar = 0.0; for(iDim=0;iDim CSourcePieceWise_TurbSA_E::ComputeResidual(const CConf } -CSourcePieceWise_TurbSA_E_COMP::CSourcePieceWise_TurbSA_E_COMP(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config) { } +template +CSourcePieceWise_TurbSA_E_COMP::CSourcePieceWise_TurbSA_E_COMP(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CSourceBase_TurbSA(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()) { } -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const CConfig* config) { - unsigned short iDim; + unsigned short iDim, jDim; // AD::StartPreacc(); // AD::SetPreaccIn(V_i, nDim+6); @@ -479,14 +484,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - if (incompressible) { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - } - else { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; Residual = 0.0; Production = 0.0; @@ -506,9 +505,13 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const Sbar = 0.0; for(iDim=0;iDim CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const Residual = Production - Destruction + CrossProduction; /*--- Compressibility Correction term ---*/ - Pressure_i = V_i[nDim+1]; + Pressure_i = V_i[idx.Pressure()]; SoundSpeed_i = sqrt(Pressure_i*Gamma/Density_i); aux_cc=0; for(iDim=0;iDim CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const } -CSourcePieceWise_TurbSA_Neg::CSourcePieceWise_TurbSA_Neg(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config) { } +template +CSourcePieceWise_TurbSA_Neg::CSourcePieceWise_TurbSA_Neg(unsigned short val_nDim, + unsigned short val_nVar, + const CConfig* config) : + CSourceBase_TurbSA(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()) { } -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_Neg::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSA_Neg::ComputeResidual(const CConfig* config) { unsigned short iDim; @@ -619,14 +627,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_Neg::ComputeResidual(const CCo // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - if (incompressible) { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - } - else { - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; Residual = 0.0; Production = 0.0; @@ -751,13 +753,15 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA_Neg::ComputeResidual(const CCo } -CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, - unsigned short val_nVar, - const su2double *constants, - su2double val_kine_Inf, - su2double val_omega_Inf, - const CConfig* config) : - CNumerics(val_nDim, val_nVar, config) { +template +CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, + unsigned short val_nVar, + const su2double *constants, + su2double val_kine_Inf, + su2double val_omega_Inf, + const CConfig* config) : + CNumerics(val_nDim, val_nVar, config), + idx(val_nDim, config->GetnSpecies()) { incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); sustaining_terms = (config->GetKind_Turb_Model() == TURB_MODEL::SST_SUST); @@ -785,7 +789,8 @@ CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, } -CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfig* config) { +template +CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfig* config) { AD::StartPreacc(); AD::SetPreaccIn(StrainMag_i); @@ -793,7 +798,7 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); AD::SetPreaccIn(F1_i); AD::SetPreaccIn(F2_i); AD::SetPreaccIn(CDkw_i); - AD::SetPreaccIn(PrimVar_Grad_i, nDim+1, nDim); + AD::SetPreaccIn(PrimVar_Grad_i, nDim+idx.Velocity(), nDim); AD::SetPreaccIn(Vorticity_i, 3); unsigned short iDim; @@ -803,20 +808,11 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi Vorticity_i[1]*Vorticity_i[1] + Vorticity_i[2]*Vorticity_i[2]); - if (incompressible) { - AD::SetPreaccIn(V_i, nDim+6); + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+4]; - Eddy_Viscosity_i = V_i[nDim+5]; - } - else { - AD::SetPreaccIn(V_i, nDim+7); - - Density_i = V_i[nDim+2]; - Laminar_Viscosity_i = V_i[nDim+5]; - Eddy_Viscosity_i = V_i[nDim+6]; - } + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; + Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; Residual[0] = 0.0; Residual[1] = 0.0; Jacobian_i[0][0] = 0.0; Jacobian_i[0][1] = 0.0; @@ -833,13 +829,13 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi diverg = 0.0; for (iDim = 0; iDim < nDim; iDim++) - diverg += PrimVar_Grad_i[iDim+1][iDim]; + diverg += PrimVar_Grad_i[iDim+idx.Velocity()][iDim]; /* if using UQ methodolgy, calculate production using perturbed Reynolds stress matrix */ if (using_uq){ ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, - PrimVar_Grad_i+1, Density_i, Eddy_Viscosity_i, + PrimVar_Grad_i+idx.Velocity(), Density_i, Eddy_Viscosity_i, ScalarVar_i[0], MeanPerturbedRSM); SetPerturbedStrainMag(ScalarVar_i[0]); pk = Eddy_Viscosity_i*PerturbedStrainMag*PerturbedStrainMag @@ -912,7 +908,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } -void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke){ +template +void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke) { /*--- Compute norm of perturbed strain rate tensor. ---*/ @@ -928,3 +925,28 @@ void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke){ PerturbedStrainMag = sqrt(2.0*PerturbedStrainMag); } + +/*--- Explicit instantiations until we don't move this to the hpp. ---*/ +template class CSourcePieceWise_TurbSA >; +template class CSourcePieceWise_TurbSA >; +template class CSourcePieceWise_TurbSA >; + +template class CSourcePieceWise_TurbSA_COMP >; +template class CSourcePieceWise_TurbSA_COMP >; +template class CSourcePieceWise_TurbSA_COMP >; + +template class CSourcePieceWise_TurbSA_E >; +template class CSourcePieceWise_TurbSA_E >; +template class CSourcePieceWise_TurbSA_E >; + +template class CSourcePieceWise_TurbSA_E_COMP >; +template class CSourcePieceWise_TurbSA_E_COMP >; +template class CSourcePieceWise_TurbSA_E_COMP >; + +template class CSourcePieceWise_TurbSA_Neg >; +template class CSourcePieceWise_TurbSA_Neg >; +template class CSourcePieceWise_TurbSA_Neg >; + +template class CSourcePieceWise_TurbSST >; +template class CSourcePieceWise_TurbSST >; +template class CSourcePieceWise_TurbSST >; From f7d054e21add37d4e23c0c38a63d32834691dc45 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 8 Nov 2021 12:35:19 +0000 Subject: [PATCH 11/12] doxygen --- .../numerics/scalar/scalar_convection.hpp | 34 ++++++------------- .../numerics/scalar/scalar_diffusion.hpp | 25 +++++--------- .../numerics/turbulent/turb_sources.hpp | 12 +++---- 3 files changed, 25 insertions(+), 46 deletions(-) diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index 8da8296f98c..5e6945af59d 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -47,12 +47,15 @@ template class CUpwScalar : public CNumerics { protected: + enum : unsigned short {MAXNVAR = 8}; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */ - su2double* Flux = nullptr; /*!< \brief Final result, diffusive flux/residual. */ - su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */ - su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */ + su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ + su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ + su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ + su2double JacobianBuffer[2*MAXNVAR*MAXNVAR]; /*!< \brief Static storage for the two Jacobians. */ const bool implicit = false, incompressible = false, dynamic_grid = false; @@ -82,27 +85,12 @@ class CUpwScalar : public CNumerics { implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT), incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), dynamic_grid(config->GetDynamic_Grid()) { - Flux = new su2double[nVar]; - Jacobian_i = new su2double*[nVar]; - Jacobian_j = new su2double*[nVar]; - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double[nVar]; - Jacobian_j[iVar] = new su2double[nVar]; + if (nVar > MAXNVAR) { + SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION); } - } - - /*! - * \brief Destructor of the class. - */ - ~CUpwScalar() override { - delete[] Flux; - if (Jacobian_i != nullptr) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - delete[] Jacobian_i[iVar]; - delete[] Jacobian_j[iVar]; - } - delete[] Jacobian_i; - delete[] Jacobian_j; + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar]; + Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR]; } } diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 04c96d6de42..e3eef508d1b 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -49,13 +49,13 @@ class CAvgGrad_Scalar : public CNumerics { protected: enum : unsigned short {MAXNVAR = 8}; - const FlowIndices idx; - su2double Proj_Mean_GradScalarVar_Normal[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal. */ - su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ + su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */ const bool correct_gradient = false, implicit = false, incompressible = false; @@ -90,18 +90,8 @@ class CAvgGrad_Scalar : public CNumerics { SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION); } for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double[nVar](); - Jacobian_j[iVar] = new su2double[nVar](); - } - } - - /*! - * \brief Destructor of the class. - */ - ~CAvgGrad_Scalar() override { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - delete[] Jacobian_i[iVar]; - delete[] Jacobian_j[iVar]; + Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar]; + Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR]; } } @@ -133,9 +123,10 @@ class CAvgGrad_Scalar : public CNumerics { ExtraADPreaccIn(); + su2double ProjGradScalarVarNoCorr[MAXNVAR]; proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j, - correct_gradient, ScalarVar_i, ScalarVar_j, - Proj_Mean_GradScalarVar_Normal, Proj_Mean_GradScalarVar); + correct_gradient, ScalarVar_i, ScalarVar_j, ProjGradScalarVarNoCorr, + Proj_Mean_GradScalarVar); FinishResidualCalc(config); AD::SetPreaccOut(Flux, nVar); diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index d4a0517f057..2929966e284 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -128,7 +128,7 @@ class CSourceBase_TurbSA : public CNumerics { template class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; @@ -167,7 +167,7 @@ class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { template class CSourcePieceWise_TurbSA_COMP final : public CSourceBase_TurbSA { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; @@ -205,7 +205,7 @@ class CSourcePieceWise_TurbSA_COMP final : public CSourceBase_TurbSA { template class CSourcePieceWise_TurbSA_E final : public CSourceBase_TurbSA { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; @@ -241,7 +241,7 @@ class CSourcePieceWise_TurbSA_E final : public CSourceBase_TurbSA { template class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; @@ -277,7 +277,7 @@ class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { template class CSourcePieceWise_TurbSA_Neg : public CSourceBase_TurbSA { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; su2double r, g, g_6, glim, fw; @@ -312,7 +312,7 @@ class CSourcePieceWise_TurbSA_Neg : public CSourceBase_TurbSA { template class CSourcePieceWise_TurbSST final : public CNumerics { private: - const FlowIndices idx; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double F1_i, F1_j, From 57a736ab66894d50baee204ab177c64d72355aa2 Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Mon, 8 Nov 2021 16:44:33 +0000 Subject: [PATCH 12/12] Apply suggestions from code review Co-authored-by: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com> --- .../numerics/turbulent/turb_diffusion.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index 40eac034996..2b563ca5b9d 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -66,9 +66,9 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { void FinishResidualCalc(const CConfig* config) override { /*--- Compute mean effective viscosity ---*/ - su2double nu_i = Laminar_Viscosity_i/Density_i; - su2double nu_j = Laminar_Viscosity_j/Density_j; - su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); + const su2double nu_i = Laminar_Viscosity_i/Density_i; + const su2double nu_j = Laminar_Viscosity_j/Density_j; + const su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; @@ -131,11 +131,11 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { void FinishResidualCalc(const CConfig* config) override { /*--- Compute mean effective viscosity ---*/ - su2double nu_i = Laminar_Viscosity_i/Density_i; - su2double nu_j = Laminar_Viscosity_j/Density_j; + const su2double nu_i = Laminar_Viscosity_i/Density_i; + const su2double nu_j = Laminar_Viscosity_j/Density_j; - su2double nu_ij = 0.5*(nu_i+nu_j); - su2double nu_tilde_ij = 0.5*(ScalarVar_i[0] + ScalarVar_j[0]); + const su2double nu_ij = 0.5*(nu_i+nu_j); + const su2double nu_tilde_ij = 0.5*(ScalarVar_i[0] + ScalarVar_j[0]); su2double nu_e; @@ -143,8 +143,8 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { nu_e = nu_ij + nu_tilde_ij; } else { - su2double Xi = nu_tilde_ij/nu_ij; - su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi); + const su2double Xi = nu_tilde_ij/nu_ij; + const su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi); nu_e = nu_ij + fn*nu_tilde_ij; }