From f2cc18a6cfe98168311449584f9f1628bd1edfac Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Mon, 30 Nov 2020 12:10:43 +0100 Subject: [PATCH 01/18] CNumerics::MeanRateOfStrain computed by CNumerics::CompMROfSMat and removed CAvgGrad_Base::GetMeanRateOfStrainMatrix, CSourcePieceWise_TurbSST::GetMeanRateOfStrainMatrix The new CNumerics member is initialized only if using_uq, in analogy to the member MeanReynoldsStress. Generalize this later. --- SU2_CFD/include/numerics/CNumerics.hpp | 8 ++++ .../include/numerics/flow/flow_diffusion.hpp | 6 --- .../numerics/turbulent/turb_sources.hpp | 6 --- SU2_CFD/src/numerics/CNumerics.cpp | 32 +++++++++++++ SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 43 ++---------------- .../src/numerics/turbulent/turb_sources.cpp | 45 ++----------------- 6 files changed, 46 insertions(+), 94 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 37eff768daa2..f6b18885bf81 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -215,6 +215,7 @@ class CNumerics { su2double *l, *m; + su2double **MeanRateOfStrain; /*!< Mean rate of strain tensor. */ su2double **MeanReynoldsStress; /*!< \brief Mean Reynolds stress tensor */ su2double **MeanPerturbedRSM; /*!< \brief Perturbed Reynolds stress tensor */ bool using_uq, /*!< \brief Flag for UQ methodology */ @@ -459,6 +460,13 @@ class CNumerics { TurbPsi_Grad_j = val_turbpsivar_grad_j; } + /*! + * \brief Set the mean rate of strain matrix. + * \details The parameter primvargrad can be either PrimVar_Grad_i or Mean_GradPrimVar. + * \param[in] primvargrad - A primitive variable gradient matrix. + */ + void ComputeMeanRateOfStrainMatrix(su2double** primvargrad); + /*! * \brief Set the value of the first blending function. * \param[in] val_F1_i - Value of the first Menter blending function at point i. diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index e5e356c62522..08a2141a0966 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -190,12 +190,6 @@ class CAvgGrad_Base : public CNumerics { */ void SetPerturbedRSM(su2double turb_ke, const CConfig* config); - /*! - * \brief Get the mean rate of strain matrix based on velocity gradients - * \param[in] S_ij - */ - void GetMeanRateOfStrainMatrix(su2double **S_ij) const; - public: /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index fe8886cdb98f..97f55d37c411 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -339,12 +339,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics { */ void SetPerturbedStrainMag(su2double turb_ke); - /*! - * \brief Get the mean rate of strain matrix based on velocity gradients - * \param[in] S_ij - */ - void GetMeanRateOfStrainMatrix(su2double **S_ij); - public: /*! * \brief Constructor of the class. diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 75955dcbb3df..7115a849d5ad 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -109,6 +109,7 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, /* --- Initializing variables for the UQ methodology --- */ using_uq = config->GetUsing_UQ(); if (using_uq){ + MeanRateOfStrain = new su2double* [3]; MeanReynoldsStress = new su2double* [3]; MeanPerturbedRSM = new su2double* [3]; A_ij = new su2double* [3]; @@ -120,6 +121,7 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Barycentric_Coord = new su2double [2]; New_Coord = new su2double [2]; for (iDim = 0; iDim < 3; iDim++){ + MeanRateOfStrain[iDim] = new su2double [3]; MeanReynoldsStress[iDim] = new su2double [3]; MeanPerturbedRSM[iDim] = new su2double [3]; A_ij[iDim] = new su2double [3]; @@ -184,6 +186,7 @@ CNumerics::~CNumerics(void) { if (using_uq) { for (unsigned short iDim = 0; iDim < 3; iDim++){ + delete [] MeanRateOfStrain[iDim]; delete [] MeanReynoldsStress[iDim]; delete [] MeanPerturbedRSM[iDim]; delete [] A_ij[iDim]; @@ -192,6 +195,7 @@ CNumerics::~CNumerics(void) { delete [] New_Eig_Vec[iDim]; delete [] Corners[iDim]; } + delete [] MeanRateOfStrain; delete [] MeanReynoldsStress; delete [] MeanPerturbedRSM; delete [] A_ij; @@ -472,6 +476,34 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou AD::EndPassive(wasActive); } +void CNumerics::ComputeMeanRateOfStrainMatrix(su2double **primvargrad){ + + /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ + + if (nDim == 3){ + MeanRateOfStrain[0][0] = primvargrad[1][0]; + MeanRateOfStrain[1][1] = primvargrad[2][1]; + MeanRateOfStrain[2][2] = primvargrad[3][2]; + MeanRateOfStrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); + MeanRateOfStrain[0][2] = 0.5 * (primvargrad[1][2] + primvargrad[3][0]); + MeanRateOfStrain[1][2] = 0.5 * (primvargrad[2][2] + primvargrad[3][1]); + MeanRateOfStrain[1][0] = MeanRateOfStrain[0][1]; + MeanRateOfStrain[2][1] = MeanRateOfStrain[1][2]; + MeanRateOfStrain[2][0] = MeanRateOfStrain[0][2]; + } + else { // nDim==2 + MeanRateOfStrain[0][0] = primvargrad[1][0]; + MeanRateOfStrain[1][1] = primvargrad[2][1]; + MeanRateOfStrain[2][2] = 0.0; + MeanRateOfStrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); + MeanRateOfStrain[0][2] = 0.0; + MeanRateOfStrain[1][2] = 0.0; + MeanRateOfStrain[1][0] = MeanRateOfStrain[0][1]; + MeanRateOfStrain[2][1] = MeanRateOfStrain[1][2]; + MeanRateOfStrain[2][0] = MeanRateOfStrain[0][2]; + } +} + void CNumerics::GetPreconditioner(const su2double *val_density, const su2double *val_velocity, const su2double *val_betainc2, const su2double *val_cp, const su2double *val_temperature, const su2double *val_drhodt, diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index d1a9ac4f1628..13b296b10db2 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -217,67 +217,30 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, tau[iDim][jDim] = tau[iDim][jDim]*(val_tau_wall/WallShearStress); } -void CAvgGrad_Base::GetMeanRateOfStrainMatrix(su2double **S_ij) const -{ - /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ - - if (nDim == 3){ - S_ij[0][0] = Mean_GradPrimVar[1][0]; - S_ij[1][1] = Mean_GradPrimVar[2][1]; - S_ij[2][2] = Mean_GradPrimVar[3][2]; - S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]); - S_ij[0][2] = 0.5 * (Mean_GradPrimVar[1][2] + Mean_GradPrimVar[3][0]); - S_ij[1][2] = 0.5 * (Mean_GradPrimVar[2][2] + Mean_GradPrimVar[3][1]); - S_ij[1][0] = S_ij[0][1]; - S_ij[2][1] = S_ij[1][2]; - S_ij[2][0] = S_ij[0][2]; - } - else { - S_ij[0][0] = Mean_GradPrimVar[1][0]; - S_ij[1][1] = Mean_GradPrimVar[2][1]; - S_ij[2][2] = 0.0; - S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]); - S_ij[0][2] = 0.0; - S_ij[1][2] = 0.0; - S_ij[1][0] = S_ij[0][1]; - S_ij[2][1] = S_ij[1][2]; - S_ij[2][0] = S_ij[0][2]; - - } -} - void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ unsigned short iDim, jDim; - su2double **S_ij = new su2double* [3]; su2double muT = Mean_Eddy_Viscosity; su2double divVel = 0; su2double density; su2double TWO3 = 2.0/3.0; density = Mean_PrimVar[nDim+2]; - for (iDim = 0; iDim < 3; iDim++){ - S_ij[iDim] = new su2double [3]; - } - - GetMeanRateOfStrainMatrix(S_ij); + ComputeMeanRateOfStrainMatrix(Mean_GradPrimVar); /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */ for (iDim = 0; iDim < 3; iDim++){ - divVel += S_ij[iDim][iDim]; + divVel += MeanRateOfStrain[iDim][iDim]; } for (iDim = 0; iDim < 3; iDim++){ for (jDim = 0; jDim < 3; jDim++){ MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] - - muT / density * (2 * S_ij[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); + - muT / density * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); } } - for (iDim = 0; iDim < 3; iDim++) - delete [] S_ij[iDim]; - delete [] S_ij; } void CAvgGrad_Base::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index 2ab0f0aab780..ee88bee7fa6a 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -919,65 +919,26 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } -void CSourcePieceWise_TurbSST::GetMeanRateOfStrainMatrix(su2double **S_ij) -{ - /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ - - if (nDim == 3){ - S_ij[0][0] = PrimVar_Grad_i[1][0]; - S_ij[1][1] = PrimVar_Grad_i[2][1]; - S_ij[2][2] = PrimVar_Grad_i[3][2]; - S_ij[0][1] = 0.5 * (PrimVar_Grad_i[1][1] + PrimVar_Grad_i[2][0]); - S_ij[0][2] = 0.5 * (PrimVar_Grad_i[1][2] + PrimVar_Grad_i[3][0]); - S_ij[1][2] = 0.5 * (PrimVar_Grad_i[2][2] + PrimVar_Grad_i[3][1]); - S_ij[1][0] = S_ij[0][1]; - S_ij[2][1] = S_ij[1][2]; - S_ij[2][0] = S_ij[0][2]; - } - else { - S_ij[0][0] = PrimVar_Grad_i[1][0]; - S_ij[1][1] = PrimVar_Grad_i[2][1]; - S_ij[2][2] = 0.0; - S_ij[0][1] = 0.5 * (PrimVar_Grad_i[1][1] + PrimVar_Grad_i[2][0]); - S_ij[0][2] = 0.0; - S_ij[1][2] = 0.0; - S_ij[1][0] = S_ij[0][1]; - S_ij[2][1] = S_ij[1][2]; - S_ij[2][0] = S_ij[0][2]; - - } -} - void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ unsigned short iDim, jDim; - su2double **S_ij = new su2double* [3]; su2double divVel = 0; su2double TWO3 = 2.0/3.0; - - - for (iDim = 0; iDim < 3; iDim++){ - S_ij[iDim] = new su2double [3]; - } - - GetMeanRateOfStrainMatrix(S_ij); + ComputeMeanRateOfStrainMatrix(PrimVar_Grad_i); /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */ for (iDim = 0; iDim < 3; iDim++){ - divVel += S_ij[iDim][iDim]; + divVel += MeanRateOfStrain[iDim][iDim]; } for (iDim = 0; iDim < 3; iDim++){ for (jDim = 0; jDim < 3; jDim++){ MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] - - Eddy_Viscosity_i / Density_i * (2 * S_ij[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); + - Eddy_Viscosity_i / Density_i * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); } } - for (iDim = 0; iDim < 3; iDim++) - delete [] S_ij[iDim]; - delete [] S_ij; } void CSourcePieceWise_TurbSST::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ From 8b893d339978333cbda52d022500da53f6af628f Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Mon, 30 Nov 2020 13:29:16 +0100 Subject: [PATCH 02/18] CNumerics::ComputeReynoldsStressMatrix implemented and used in CAvgGrad_Base::SetReynoldsStressMatrix, CSourcePieceWise_TurbSST::SetReynoldsStressMatrix --- SU2_CFD/include/numerics/CNumerics.hpp | 9 ++++++++ SU2_CFD/src/numerics/CNumerics.cpp | 19 +++++++++++++++ SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 23 +------------------ .../src/numerics/turbulent/turb_sources.cpp | 19 +-------------- 4 files changed, 30 insertions(+), 40 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index f6b18885bf81..5a7da88e246e 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -467,6 +467,15 @@ class CNumerics { */ void ComputeMeanRateOfStrainMatrix(su2double** primvargrad); + /*! + * \brief Set the mean Reynolds stress matrix +(u_i' u_j')~. + * \details The mean rate of strain matrix must be already set. + * \param[in] turb_ke - Turbulent kinetic energy + * \param[in] eddy_visc - Eddy viscosity + * \param[in] density - Density + */ + void ComputeReynoldsStressMatrix(su2double turb_ke, su2double eddy_visc, su2double density); + /*! * \brief Set the value of the first blending function. * \param[in] val_F1_i - Value of the first Menter blending function at point i. diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 7115a849d5ad..7ed54f64aea4 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -504,6 +504,25 @@ void CNumerics::ComputeMeanRateOfStrainMatrix(su2double **primvargrad){ } } +void CNumerics::ComputeReynoldsStressMatrix(su2double turb_ke, su2double eddy_visc, su2double density){ + su2double divVel = 0; + su2double TWO3 = 2.0/3.0; + + /* --- Using the rate of strain matrix already set, calculate Reynolds stress tensor --- */ + + for (unsigned short iDim = 0; iDim < 3; iDim++){ + divVel += MeanRateOfStrain[iDim][iDim]; + } + + for (unsigned short iDim = 0; iDim < 3; iDim++){ + for (unsigned short jDim = 0; jDim < 3; jDim++){ + MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] + - eddy_visc / density * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); + } + } + +} + void CNumerics::GetPreconditioner(const su2double *val_density, const su2double *val_velocity, const su2double *val_betainc2, const su2double *val_cp, const su2double *val_temperature, const su2double *val_drhodt, diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 13b296b10db2..cfd73bf39311 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -218,29 +218,8 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, } void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ - - unsigned short iDim, jDim; - su2double muT = Mean_Eddy_Viscosity; - su2double divVel = 0; - su2double density; - su2double TWO3 = 2.0/3.0; - density = Mean_PrimVar[nDim+2]; - ComputeMeanRateOfStrainMatrix(Mean_GradPrimVar); - - /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */ - - for (iDim = 0; iDim < 3; iDim++){ - divVel += MeanRateOfStrain[iDim][iDim]; - } - - for (iDim = 0; iDim < 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] - - muT / density * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); - } - } - + ComputeReynoldsStressMatrix(turb_ke, Mean_Eddy_Viscosity, Mean_PrimVar[nDim+2]); } void CAvgGrad_Base::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index ee88bee7fa6a..0342c184c0a6 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -920,25 +920,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - unsigned short iDim, jDim; - su2double divVel = 0; - su2double TWO3 = 2.0/3.0; - ComputeMeanRateOfStrainMatrix(PrimVar_Grad_i); - - /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */ - - for (iDim = 0; iDim < 3; iDim++){ - divVel += MeanRateOfStrain[iDim][iDim]; - } - - for (iDim = 0; iDim < 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] - - Eddy_Viscosity_i / Density_i * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); - } - } - + ComputeReynoldsStressMatrix(turb_ke, Eddy_Viscosity_i, Density_i); } void CSourcePieceWise_TurbSST::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ From c2ded5afeb500c9c9f96f2824f0a826fca50ae20 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Tue, 1 Dec 2020 10:02:52 +0100 Subject: [PATCH 03/18] Allocation of Reynolds matrix if CConfig::GetUsing_ReyStress this applies to the mean rate of strain matrix as well Currently they are only allocated if UQ methodology is used --- Common/include/CConfig.hpp | 8 ++++++ SU2_CFD/include/numerics/CNumerics.hpp | 3 +- SU2_CFD/src/numerics/CNumerics.cpp | 29 ++++++++++++++------ SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 4 ++- 4 files changed, 34 insertions(+), 10 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 265f7dc91949..511d0f94b5f8 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -8828,6 +8828,14 @@ class CConfig { */ bool GetPrintInlet_InterpolatedData(void) const { return PrintInlet_InterpolatedData; } + /*! + * \brief Get information about using the Reynolds stress tensor and mean rate of strain matrix. + * \return TRUE means that they will be used. + */ + bool GetUsing_ReynoldsStress(void) const { + return (using_uq); + } + /*! * \brief Get information about using UQ methodology * \return TRUE means that UQ methodology of eigenspace perturbation will be used diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 5a7da88e246e..f86a841844cf 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -215,7 +215,8 @@ class CNumerics { su2double *l, *m; - su2double **MeanRateOfStrain; /*!< Mean rate of strain tensor. */ + bool using_reynoldsstress; /*!< \brief Flag for usage of mean Reynolds stress and strain rate matrices */ + su2double **MeanRateOfStrain; /*!< \brief Mean rate of strain tensor. */ su2double **MeanReynoldsStress; /*!< \brief Mean Reynolds stress tensor */ su2double **MeanPerturbedRSM; /*!< \brief Perturbed Reynolds stress tensor */ bool using_uq, /*!< \brief Flag for UQ methodology */ diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 7ed54f64aea4..e96c81418f8a 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -49,6 +49,7 @@ CNumerics::CNumerics(void) { l = nullptr; m = nullptr; + using_reynoldsstress = false; using_uq = false; nemo = false; @@ -106,11 +107,20 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Dissipation_ij = 1.0; - /* --- Initializing variables for the UQ methodology --- */ - using_uq = config->GetUsing_UQ(); - if (using_uq){ + /* --- Initializing Reynolds stress matrix and strain rate matrix --- */ + using_reynoldsstress= config->GetUsing_ReynoldsStress(); + if (using_reynoldsstress){ MeanRateOfStrain = new su2double* [3]; MeanReynoldsStress = new su2double* [3]; + for (iDim = 0; iDim < 3; iDim++){ + MeanRateOfStrain[iDim] = new su2double [3]; + MeanReynoldsStress[iDim] = new su2double [3]; + } + } + + /* --- Initializing additional variables for the UQ methodology --- */ + using_uq = config->GetUsing_UQ(); + if (using_uq){ MeanPerturbedRSM = new su2double* [3]; A_ij = new su2double* [3]; newA_ij = new su2double* [3]; @@ -121,8 +131,6 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Barycentric_Coord = new su2double [2]; New_Coord = new su2double [2]; for (iDim = 0; iDim < 3; iDim++){ - MeanRateOfStrain[iDim] = new su2double [3]; - MeanReynoldsStress[iDim] = new su2double [3]; MeanPerturbedRSM[iDim] = new su2double [3]; A_ij[iDim] = new su2double [3]; newA_ij[iDim] = new su2double [3]; @@ -184,10 +192,17 @@ CNumerics::~CNumerics(void) { delete [] l; delete [] m; - if (using_uq) { + if (using_reynoldsstress){ for (unsigned short iDim = 0; iDim < 3; iDim++){ delete [] MeanRateOfStrain[iDim]; delete [] MeanReynoldsStress[iDim]; + } + delete [] MeanRateOfStrain; + delete [] MeanReynoldsStress; + } + + if (using_uq) { + for (unsigned short iDim = 0; iDim < 3; iDim++){ delete [] MeanPerturbedRSM[iDim]; delete [] A_ij[iDim]; delete [] newA_ij[iDim]; @@ -195,8 +210,6 @@ CNumerics::~CNumerics(void) { delete [] New_Eig_Vec[iDim]; delete [] Corners[iDim]; } - delete [] MeanRateOfStrain; - delete [] MeanReynoldsStress; delete [] MeanPerturbedRSM; delete [] A_ij; delete [] newA_ij; diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index cfd73bf39311..d1837edb263c 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -131,7 +131,9 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, for (iDim = 0 ; iDim < nDim; iDim++) div_vel += val_gradprimvar[iDim+1][iDim]; - /* --- If UQ methodology is used, calculate tau using the perturbed reynolds stress tensor --- */ + /* --- If UQ methodology is used, use the perturbed Reynolds stress tensor + * for the turbulent part of tau. Otherwise both the laminar and turbulent + * parts of tau can be computed with the total viscosity. --- */ if (using_uq){ for (iDim = 0 ; iDim < nDim; iDim++) From a182c89567c826c533ad5000ad23bd1d3fb6120f Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Tue, 1 Dec 2020 11:02:11 +0100 Subject: [PATCH 04/18] generalized CNumerics methods for Reynolds stress and strain rate The methods CNumerics::ComputeStressTensor and CNumerics::ComputeMeanRateOfStrainMatrix are now static, and their input (rate of strain or primitive variable gradients, respectively) and output (stress tensor or rate of strain, resp) are given as parameters. We try to use this for both laminar and turbulent stresses. --- SU2_CFD/include/numerics/CNumerics.hpp | 23 +++++--- SU2_CFD/src/numerics/CNumerics.cpp | 59 ++++++++++--------- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 35 ++++++----- .../src/numerics/turbulent/turb_sources.cpp | 9 ++- 4 files changed, 69 insertions(+), 57 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index f86a841844cf..cc6d71938acc 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -462,20 +462,27 @@ class CNumerics { } /*! - * \brief Set the mean rate of strain matrix. - * \details The parameter primvargrad can be either PrimVar_Grad_i or Mean_GradPrimVar. + * \brief Compute the mean rate of strain matrix. + * \details The parameter primvargrad can be e.g. PrimVar_Grad_i or Mean_GradPrimVar. + * \param[in] nDim - 2 or 3 + * \param[out] rateofstrain - Rate of strain matrix * \param[in] primvargrad - A primitive variable gradient matrix. */ - void ComputeMeanRateOfStrainMatrix(su2double** primvargrad); + static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double* const* primvargrad); /*! - * \brief Set the mean Reynolds stress matrix +(u_i' u_j')~. - * \details The mean rate of strain matrix must be already set. - * \param[in] turb_ke - Turbulent kinetic energy - * \param[in] eddy_visc - Eddy viscosity + * \brief Compute the stress tensor from the rate of strain tensor. + * \details If the Reynolds stress tensor is defined as +(u_i' u_j')~, divide the result + * of this function by (-rho). + * \param[in] nDim - 2 or 3 + * \param[out] stress - Stress tensor + * \param[in] rateofstrain - Rate of strain tensor + * \param[in] viscosity - Viscosity * \param[in] density - Density + * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor */ - void ComputeReynoldsStressMatrix(su2double turb_ke, su2double eddy_visc, su2double density); + static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* rateofstrain, + su2double viscosity, su2double density, su2double turb_ke=0.0); /*! * \brief Set the value of the first blending function. diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index e96c81418f8a..b5b483785a94 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -489,48 +489,49 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou AD::EndPassive(wasActive); } -void CNumerics::ComputeMeanRateOfStrainMatrix(su2double **primvargrad){ +void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double * const* primvargrad){ /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ if (nDim == 3){ - MeanRateOfStrain[0][0] = primvargrad[1][0]; - MeanRateOfStrain[1][1] = primvargrad[2][1]; - MeanRateOfStrain[2][2] = primvargrad[3][2]; - MeanRateOfStrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); - MeanRateOfStrain[0][2] = 0.5 * (primvargrad[1][2] + primvargrad[3][0]); - MeanRateOfStrain[1][2] = 0.5 * (primvargrad[2][2] + primvargrad[3][1]); - MeanRateOfStrain[1][0] = MeanRateOfStrain[0][1]; - MeanRateOfStrain[2][1] = MeanRateOfStrain[1][2]; - MeanRateOfStrain[2][0] = MeanRateOfStrain[0][2]; + rateofstrain[0][0] = primvargrad[1][0]; + rateofstrain[1][1] = primvargrad[2][1]; + rateofstrain[2][2] = primvargrad[3][2]; + rateofstrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); + rateofstrain[0][2] = 0.5 * (primvargrad[1][2] + primvargrad[3][0]); + rateofstrain[1][2] = 0.5 * (primvargrad[2][2] + primvargrad[3][1]); + rateofstrain[1][0] = rateofstrain[0][1]; + rateofstrain[2][1] = rateofstrain[1][2]; + rateofstrain[2][0] = rateofstrain[0][2]; } else { // nDim==2 - MeanRateOfStrain[0][0] = primvargrad[1][0]; - MeanRateOfStrain[1][1] = primvargrad[2][1]; - MeanRateOfStrain[2][2] = 0.0; - MeanRateOfStrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); - MeanRateOfStrain[0][2] = 0.0; - MeanRateOfStrain[1][2] = 0.0; - MeanRateOfStrain[1][0] = MeanRateOfStrain[0][1]; - MeanRateOfStrain[2][1] = MeanRateOfStrain[1][2]; - MeanRateOfStrain[2][0] = MeanRateOfStrain[0][2]; + rateofstrain[0][0] = primvargrad[1][0]; + rateofstrain[1][1] = primvargrad[2][1]; + rateofstrain[2][2] = 0.0; + rateofstrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); + rateofstrain[0][2] = 0.0; + rateofstrain[1][2] = 0.0; + rateofstrain[1][0] = rateofstrain[0][1]; + rateofstrain[2][1] = rateofstrain[1][2]; + rateofstrain[2][0] = rateofstrain[0][2]; } } -void CNumerics::ComputeReynoldsStressMatrix(su2double turb_ke, su2double eddy_visc, su2double density){ - su2double divVel = 0; +void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* rateofstrain, + su2double viscosity, su2double density, su2double turb_ke){ su2double TWO3 = 2.0/3.0; - /* --- Using the rate of strain matrix already set, calculate Reynolds stress tensor --- */ - - for (unsigned short iDim = 0; iDim < 3; iDim++){ - divVel += MeanRateOfStrain[iDim][iDim]; + su2double divVel = 0; + for (unsigned short iDim = 0; iDim < nDim; iDim++){ + divVel += rateofstrain[iDim][iDim]; } - for (unsigned short iDim = 0; iDim < 3; iDim++){ - for (unsigned short jDim = 0; jDim < 3; jDim++){ - MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim] - - eddy_visc / density * (2 * MeanRateOfStrain[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]); + for (unsigned short iDim = 0; iDim < nDim; iDim++){ + for (unsigned short jDim = 0; jDim < nDim; jDim++){ + stress[iDim][jDim] = + viscosity * 2 * rateofstrain[iDim][jDim] + - TWO3 * viscosity * divVel * (iDim==jDim) + - TWO3 * density * turb_ke * (iDim==jDim); } } diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index d1837edb263c..c0a47aadd4ae 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -123,30 +123,23 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, const su2double val_laminar_viscosity, const su2double val_eddy_viscosity) { - unsigned short iDim, jDim; const su2double Density = val_primvar[nDim+2]; - const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; - - su2double div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - div_vel += val_gradprimvar[iDim+1][iDim]; /* --- If UQ methodology is used, use the perturbed Reynolds stress tensor * for the turbulent part of tau. Otherwise both the laminar and turbulent * parts of tau can be computed with the total viscosity. --- */ + ComputeMeanRateOfStrainMatrix(nDim,MeanRateOfStrain,val_gradprimvar); if (using_uq){ - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = val_laminar_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] ) - - TWO3*val_laminar_viscosity*div_vel*delta[iDim][jDim] - Density * MeanPerturbedRSM[iDim][jDim]; - + ComputeStressTensor(nDim, tau, MeanRateOfStrain, val_laminar_viscosity, Density, 0.0); // laminar part + // add turbulent part which was perturbed + for (unsigned short iDim = 0 ; iDim < nDim; iDim++) + for (unsigned short jDim = 0 ; jDim < nDim; jDim++) + tau[iDim][jDim] += (-Density) * MeanPerturbedRSM[iDim][jDim]; } else { - - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = total_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] ) - - TWO3*total_viscosity*div_vel*delta[iDim][jDim]; + // compute both parts in one step + const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; + ComputeStressTensor(nDim, tau, MeanRateOfStrain, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? } } @@ -220,8 +213,14 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, } void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeMeanRateOfStrainMatrix(Mean_GradPrimVar); - ComputeReynoldsStressMatrix(turb_ke, Mean_Eddy_Viscosity, Mean_PrimVar[nDim+2]); + ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain, Mean_GradPrimVar); + su2double meandensity = Mean_PrimVar[nDim+2]; + ComputeStressTensor(nDim, MeanReynoldsStress, MeanRateOfStrain, Mean_Eddy_Viscosity, meandensity, turb_ke); + for(unsigned short iDim=0; iDim CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeMeanRateOfStrainMatrix(PrimVar_Grad_i); - ComputeReynoldsStressMatrix(turb_ke, Eddy_Viscosity_i, Density_i); + ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain, PrimVar_Grad_i); + ComputeStressTensor(nDim, MeanReynoldsStress, MeanRateOfStrain, Eddy_Viscosity_i, Density_i, turb_ke); + for(unsigned short iDim=0; iDim Date: Tue, 1 Dec 2020 13:34:23 +0100 Subject: [PATCH 05/18] Using new Reynolds function at another place --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c871c5a841a4..3bb25b78a210 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -202,16 +202,8 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - - div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - div_vel += PrimVar_Grad_i[iDim+1][iDim]; - - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = (total_viscosity*(PrimVar_Grad_i[jDim+1][iDim] + - PrimVar_Grad_i[iDim+1][jDim] ) - -TWO3*total_viscosity*div_vel*delta[iDim][jDim]); + ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain,PrimVar_Grad_i); + ComputeStressTensor(nDim, tau, MeanRateOfStrain, total_viscosity, 0.0, 0.0); /*--- Viscous terms. ---*/ From abbcdb0acee3b4d385a0d29a5fe841e26f0533bb Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Tue, 1 Dec 2020 13:42:55 +0100 Subject: [PATCH 06/18] Got rid of CNumerics::MeanRateOfStrain again CNumerics::ComputeStressTensor takes primitive variable gradient instead of rate of strain matrix now --- SU2_CFD/include/numerics/CNumerics.hpp | 11 +++++------ SU2_CFD/src/numerics/CNumerics.cpp | 12 ++++-------- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 8 +++----- SU2_CFD/src/numerics/flow/flow_sources.cpp | 3 +-- SU2_CFD/src/numerics/turbulent/turb_sources.cpp | 3 +-- 5 files changed, 14 insertions(+), 23 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index cc6d71938acc..e34df2c98485 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -215,8 +215,7 @@ class CNumerics { su2double *l, *m; - bool using_reynoldsstress; /*!< \brief Flag for usage of mean Reynolds stress and strain rate matrices */ - su2double **MeanRateOfStrain; /*!< \brief Mean rate of strain tensor. */ + bool using_reynoldsstress; /*!< \brief Flag for usage of mean Reynolds stress matrix */ su2double **MeanReynoldsStress; /*!< \brief Mean Reynolds stress tensor */ su2double **MeanPerturbedRSM; /*!< \brief Perturbed Reynolds stress tensor */ bool using_uq, /*!< \brief Flag for UQ methodology */ @@ -471,17 +470,17 @@ class CNumerics { static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double* const* primvargrad); /*! - * \brief Compute the stress tensor from the rate of strain tensor. - * \details If the Reynolds stress tensor is defined as +(u_i' u_j')~, divide the result + * \brief Compute the stress tensor from the velocity gradients. + * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result * of this function by (-rho). * \param[in] nDim - 2 or 3 * \param[out] stress - Stress tensor - * \param[in] rateofstrain - Rate of strain tensor + * \param[in] primvargrad - A primitive variable gradient matrix. * \param[in] viscosity - Viscosity * \param[in] density - Density * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor */ - static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* rateofstrain, + static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, su2double viscosity, su2double density, su2double turb_ke=0.0); /*! diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index b5b483785a94..431cdb4a6068 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -107,13 +107,11 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Dissipation_ij = 1.0; - /* --- Initializing Reynolds stress matrix and strain rate matrix --- */ + /* --- Initializing Reynolds stress matrix --- */ using_reynoldsstress= config->GetUsing_ReynoldsStress(); if (using_reynoldsstress){ - MeanRateOfStrain = new su2double* [3]; MeanReynoldsStress = new su2double* [3]; for (iDim = 0; iDim < 3; iDim++){ - MeanRateOfStrain[iDim] = new su2double [3]; MeanReynoldsStress[iDim] = new su2double [3]; } } @@ -194,10 +192,8 @@ CNumerics::~CNumerics(void) { if (using_reynoldsstress){ for (unsigned short iDim = 0; iDim < 3; iDim++){ - delete [] MeanRateOfStrain[iDim]; delete [] MeanReynoldsStress[iDim]; } - delete [] MeanRateOfStrain; delete [] MeanReynoldsStress; } @@ -517,19 +513,19 @@ void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** r } } -void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* rateofstrain, +void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, su2double viscosity, su2double density, su2double turb_ke){ su2double TWO3 = 2.0/3.0; su2double divVel = 0; for (unsigned short iDim = 0; iDim < nDim; iDim++){ - divVel += rateofstrain[iDim][iDim]; + divVel += primvargrad[iDim+1][iDim]; } for (unsigned short iDim = 0; iDim < nDim; iDim++){ for (unsigned short jDim = 0; jDim < nDim; jDim++){ stress[iDim][jDim] = - viscosity * 2 * rateofstrain[iDim][jDim] + viscosity * (primvargrad[iDim+1][jDim]+primvargrad[jDim+1][iDim]) - TWO3 * viscosity * divVel * (iDim==jDim) - TWO3 * density * turb_ke * (iDim==jDim); } diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index c0a47aadd4ae..eb5ec42b569e 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -129,9 +129,8 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, * for the turbulent part of tau. Otherwise both the laminar and turbulent * parts of tau can be computed with the total viscosity. --- */ - ComputeMeanRateOfStrainMatrix(nDim,MeanRateOfStrain,val_gradprimvar); if (using_uq){ - ComputeStressTensor(nDim, tau, MeanRateOfStrain, val_laminar_viscosity, Density, 0.0); // laminar part + ComputeStressTensor(nDim, tau, val_gradprimvar, val_laminar_viscosity, Density, 0.0); // laminar part // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) @@ -139,7 +138,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } else { // compute both parts in one step const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; - ComputeStressTensor(nDim, tau, MeanRateOfStrain, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? + ComputeStressTensor(nDim, tau, val_gradprimvar, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? } } @@ -213,9 +212,8 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, } void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain, Mean_GradPrimVar); su2double meandensity = Mean_PrimVar[nDim+2]; - ComputeStressTensor(nDim, MeanReynoldsStress, MeanRateOfStrain, Mean_Eddy_Viscosity, meandensity, turb_ke); + ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar, Mean_Eddy_Viscosity, meandensity, turb_ke); for(unsigned short iDim=0; iDim CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain,PrimVar_Grad_i); - ComputeStressTensor(nDim, tau, MeanRateOfStrain, total_viscosity, 0.0, 0.0); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i, total_viscosity, 0.0, 0.0); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index a27880c13d2e..8221466ceb8b 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -920,8 +920,7 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeMeanRateOfStrainMatrix(nDim, MeanRateOfStrain, PrimVar_Grad_i); - ComputeStressTensor(nDim, MeanReynoldsStress, MeanRateOfStrain, Eddy_Viscosity_i, Density_i, turb_ke); + ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i, Eddy_Viscosity_i, Density_i, turb_ke); for(unsigned short iDim=0; iDim Date: Tue, 1 Dec 2020 15:08:55 +0100 Subject: [PATCH 07/18] Used ComputeStressTensor to replace some explicit computations --- SU2_CFD/include/numerics/CNumerics.hpp | 4 +-- .../include/solvers/CFVMFlowSolverBase.inl | 19 ++++------ .../interfaces/fsi/CFlowTractionInterface.cpp | 17 ++++----- SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp | 20 ++--------- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 2 +- SU2_CFD/src/numerics/flow/flow_sources.cpp | 6 ++-- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 18 ++++------ SU2_CFD/src/solvers/CNEMONSSolver.cpp | 13 ++----- SU2_CFD/src/solvers/CNSSolver.cpp | 36 ++++++------------- 9 files changed, 40 insertions(+), 95 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e34df2c98485..ea211bc0072d 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -472,7 +472,7 @@ class CNumerics { /*! * \brief Compute the stress tensor from the velocity gradients. * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result - * of this function by (-rho). + * of this function by (-rho). The argument density is only used if turb_ke is not 0. * \param[in] nDim - 2 or 3 * \param[out] stress - Stress tensor * \param[in] primvargrad - A primitive variable gradient matrix. @@ -481,7 +481,7 @@ class CNumerics { * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor */ static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, - su2double viscosity, su2double density, su2double turb_ke=0.0); + su2double viscosity, su2double density=0.0, su2double turb_ke=0.0); /*! * \brief Set the value of the first blending function. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index c6190ac980f0..db89f3061bda 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2030,12 +2030,11 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr unsigned long iVertex, iPoint, iPointNormal; unsigned short iMarker, iMarker_Monitoring, iDim, jDim; unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0; - su2double Viscosity = 0.0, div_vel, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0, + su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0, RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion, UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp, thermal_conductivity, thermal_conductivity_tr, thermal_conductivity_ve = 0.0, - MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor, - delta[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; + MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor; const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr; su2double **Grad_PrimVar = nullptr, dTn, dTven; @@ -2190,16 +2189,10 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - - div_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim]; - - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0; jDim < nDim; jDim++) { - Tau[iDim][jDim] = Viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) - - TWO3 * Viscosity * div_vel * delta[iDim][jDim]; - } - } + su2double *Tau_pointer[3] = {&(Tau[0][0]), &(Tau[1][0]), &(Tau[2][0])}; + su2double *Grad_Vel_pointer[4] = {nullptr, &(Grad_Vel[0][0]), &(Grad_Vel[1][0]), &(Grad_Vel[2][0])}; + CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_Vel_pointer, Viscosity); + // Grad_Vel is not a primitive variable gradient, so we have to shift the index. /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index 14f687a65bad..e229c04ed7ec 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -188,22 +188,17 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); - const su2double* const* GradVel = &flow_nodes->GetGradient_Primitive(Point_Flow)[1]; - - // Divergence of the velocity - su2double DivVel = 0.0; - for (auto iVar = 0u; iVar < nVar; iVar++) DivVel += GradVel[iVar][iVar]; - + su2double tau_data[9]; // avoid dynamic allocation + su2double *tau[3]; + tau[0] = tau_data; tau[1] = tau_data+3; tau[2] = tau_data+6; + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow),Viscosity); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { - // Viscous stress - su2double delta_ij = (iVar == jVar); - su2double tau_ij = Viscosity*(GradVel[jVar][iVar] + GradVel[iVar][jVar] - TWO3*DivVel*delta_ij); - // Viscous component in the tn vector --> Units of force (non-dimensional). - Donor_Variable[iVar] += tau_ij * Normal_Flow[jVar]; + Donor_Variable[iVar] += tau[iVar][jVar] * Normal_Flow[jVar]; } } + } // Redimensionalize and take into account ramp transfer of the loads diff --git a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp index 9a7cf845061d..e027a7dea994 100644 --- a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp +++ b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp @@ -237,7 +237,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, // rather than the standard V = [r1, ... , rn, T, Tve, ... ] unsigned short iSpecies, iVar, iDim, jDim; - su2double *Ds, *V, **GV, mu, ktr, kve, div_vel; + su2double *Ds, *V, **GV, mu, ktr, kve; su2double rho, T, Tve, RuSI, Ru; auto& Ms = fluidmodel->GetSpeciesMolarMass(); @@ -279,12 +279,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, //Cpve = V[RHOCVVE_INDEX]+Ru/Mass; //kve += Cpve*(val_eddy_viscosity/Prandtl_Turb); - /*--- Calculate the velocity divergence ---*/ - div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - div_vel += GV[VEL_INDEX+iDim][iDim]; - - + /*--- Pre-compute mixture quantities ---*/ for (iDim = 0; iDim < nDim; iDim++) { Vector[iDim] = 0.0; @@ -294,16 +289,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, } /*--- Compute the viscous stress tensor ---*/ - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - tau[iDim][jDim] = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - tau[iDim][jDim] += mu * (val_gradprimvar[VEL_INDEX+jDim][iDim] + - val_gradprimvar[VEL_INDEX+iDim][jDim]); - } - tau[iDim][iDim] -= TWO3*mu*div_vel; - } + ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX-1, mu); /*--- Populate entries in the viscous flux vector ---*/ for (iDim = 0; iDim < nDim; iDim++) { diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index eb5ec42b569e..9af8758a3d50 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -130,7 +130,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, * parts of tau can be computed with the total viscosity. --- */ if (using_uq){ - ComputeStressTensor(nDim, tau, val_gradprimvar, val_laminar_viscosity, Density, 0.0); // laminar part + ComputeStressTensor(nDim, tau, val_gradprimvar, val_laminar_viscosity); // laminar part // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index dcfedd4cefb2..15d81da84f20 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -136,7 +136,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) { su2double yinv, Velocity_i[3]; - unsigned short iDim, jDim, iVar, jVar; + unsigned short iDim, iVar, jVar; if (Coord_i[1] > EPS) { @@ -197,12 +197,12 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo Eddy_Viscosity_i = V_i[nDim+5]; Thermal_Conductivity_i = V_i[nDim+6]; - su2double total_viscosity, div_vel; + su2double total_viscosity; total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeStressTensor(nDim, tau, PrimVar_Grad_i, total_viscosity, 0.0, 0.0); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i, total_viscosity); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index e090482cec00..b7aa40670988 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -783,17 +783,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con else val_turb_ke = 0.0; - div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) { - Velocity[iDim] = U[iDim+1]/Density; - div_vel += PrimVar_Grad[iDim+1][iDim]; - } - - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = Laminar_Viscosity*(PrimVar_Grad[jDim+1][iDim] + PrimVar_Grad[iDim+1][jDim]) - - TWO3*Laminar_Viscosity*div_vel*delta[iDim][jDim] - - TWO3*Density*val_turb_ke*delta[iDim][jDim]; + CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad, Laminar_Viscosity, Density, val_turb_ke); /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ @@ -811,6 +801,12 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con /*--- Form Sigma_Psi5v ---*/ + div_vel = 0.0; + for (iDim = 0 ; iDim < nDim; iDim++) { + Velocity[iDim] = U[iDim+1]/Density; + div_vel += PrimVar_Grad[iDim+1][iDim]; + } + gradPsi5_v = 0.0; for (iDim = 0; iDim < nDim; iDim++) { gradPsi5_v += PsiVar_Grad[nDim+1][iDim]*Velocity[iDim]; diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 3c7ac0ef4567..92a424eb01ef 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -1019,7 +1019,6 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, su2double TauElem[3], TauTangent[3]; su2double Tau[3][3]; su2double TauNormal; - su2double div_vel=0, Delta; bool ionization = config->GetIonization(); @@ -1143,17 +1142,9 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar ++) Res_Visc[iVar] = 0.0; - div_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - div_vel += Grad_PrimVar[VEL_INDEX+iDim][iDim]; - + su2double *Tau_pointer[3] = {&(Tau[0][0]), &(Tau[1][0]), &(Tau[2][0])}; + CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_PrimVar+VEL_INDEX-1, Viscosity); for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Delta = 0.0; if (iDim == jDim) Delta = 1.0; - Tau[iDim][jDim] = Viscosity*(Grad_PrimVar[VEL_INDEX+jDim][iDim] + - Grad_PrimVar[VEL_INDEX+iDim][jDim] ) - - TWO3*Viscosity*div_vel*Delta; - } TauElem[iDim] = 0.0; for (jDim = 0; jDim < nDim; jDim++) TauElem[iDim] += Tau[iDim][jDim]*UnitNormal[jDim]; diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 1a501e7628e3..3cc7cc31b384 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -471,23 +471,13 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint); su2double total_viscosity = laminar_viscosity + eddy_viscosity; - const auto Grad_Vel = &nodes->GetGradient_Primitive(iPoint)[1]; - - /*--- Divergence of the velocity ---*/ - - su2double div_vel = 0.0; - for (auto iDim = 0u; iDim < nDim; iDim++) - div_vel += Grad_Vel[iDim][iDim]; - /*--- Compute the viscous stress tensor ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; - for (auto iDim = 0u; iDim < nDim; iDim++) { - for (auto jDim = 0u; jDim < nDim; jDim++) { - tau[iDim][jDim] = total_viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]); - } - tau[iDim][iDim] -= TWO3*total_viscosity*div_vel; - } + su2double *tau_pointer[MAXNDIM]; // avoid dynamic allocation + for(unsigned short iDim=0; iDimGetGradient_Primitive(iPoint), total_viscosity); /*--- Dot product of the stress tensor with the grid velocity ---*/ @@ -988,21 +978,15 @@ void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, C /*--- Compute the shear stress at the wall in the regular fashion by using the stress tensor on the surface ---*/ - su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); - - const auto GradVel = &nodes->GetGradient_Primitive(iPoint)[1]; + su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0}; + su2double *tau_pointer[MAXNDIM]; // avoid dynamic allocation + for(unsigned short iDim=0; iDimGetLaminarViscosity(iPoint); + CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint), Lam_Visc_Wall); - su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0}; for (auto iDim = 0u; iDim < nDim; iDim++) { - for (auto jDim = 0u; jDim < nDim; jDim++) { - tau[iDim][jDim] = Lam_Visc_Wall * (GradVel[jDim][iDim] + GradVel[iDim][jDim]); - } - tau[iDim][iDim] -= TWO3*Lam_Visc_Wall*div_vel; - TauElem[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], UnitNormal); } From ca63718514aced95afcb9beacf5372272d05b0f0 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Wed, 2 Dec 2020 11:26:19 +0100 Subject: [PATCH 08/18] Something is wrong with UQ methodology, trying to fix... --- SU2_CFD/include/numerics/CNumerics.hpp | 5 +++-- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 9 +++++++-- SU2_CFD/src/numerics/CNumerics.cpp | 7 ++++++- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 6 +++--- SU2_CFD/src/numerics/turbulent/turb_sources.cpp | 6 +++--- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 4 +++- SU2_CFD/src/solvers/CNSSolver.cpp | 4 ++-- 7 files changed, 27 insertions(+), 14 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index ea211bc0072d..094f8c2b2649 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -473,15 +473,16 @@ class CNumerics { * \brief Compute the stress tensor from the velocity gradients. * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result * of this function by (-rho). The argument density is only used if turb_ke is not 0. - * \param[in] nDim - 2 or 3 + * \param[in] nDim - Dimension of the flow problem, 2 or 3 * \param[out] stress - Stress tensor * \param[in] primvargrad - A primitive variable gradient matrix. * \param[in] viscosity - Viscosity * \param[in] density - Density * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor + * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2. */ static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, - su2double viscosity, su2double density=0.0, su2double turb_ke=0.0); + su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false); /*! * \brief Set the value of the first blending function. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index db89f3061bda..6226fa1595ef 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2189,8 +2189,13 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - su2double *Tau_pointer[3] = {&(Tau[0][0]), &(Tau[1][0]), &(Tau[2][0])}; - su2double *Grad_Vel_pointer[4] = {nullptr, &(Grad_Vel[0][0]), &(Grad_Vel[1][0]), &(Grad_Vel[2][0])}; + su2double *Tau_pointer[3]; + su2double *Grad_Vel_pointer[4]; + Grad_Vel_pointer[0] = nullptr; + for(iDim=0;iDim<3;iDim++){ + Tau_pointer[iDim] = Tau[iDim]; + Grad_Vel_pointer[iDim+1] = Grad_Vel[iDim]; + } CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_Vel_pointer, Viscosity); // Grad_Vel is not a primitive variable gradient, so we have to shift the index. diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 431cdb4a6068..5b70417150ce 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -514,7 +514,7 @@ void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** r } void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, - su2double viscosity, su2double density, su2double turb_ke){ + su2double viscosity, su2double density, su2double turb_ke, bool reynolds3x3){ su2double TWO3 = 2.0/3.0; su2double divVel = 0; @@ -531,6 +531,11 @@ void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, con } } + if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix + stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0; + stress[2][2] = -TWO3 * ( viscosity * divVel + density * turb_ke ); + } + } void CNumerics::GetPreconditioner(const su2double *val_density, const su2double *val_velocity, diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 9af8758a3d50..beaea9567305 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -213,9 +213,9 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ su2double meandensity = Mean_PrimVar[nDim+2]; - ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar, Mean_Eddy_Viscosity, meandensity, turb_ke); - for(unsigned short iDim=0; iDim CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i, Eddy_Viscosity_i, Density_i, turb_ke); - for(unsigned short iDim=0; iDimGetGradient_Primitive(iPoint), total_viscosity); /*--- Dot product of the stress tensor with the grid velocity ---*/ @@ -981,7 +981,7 @@ void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, C su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0}; su2double *tau_pointer[MAXNDIM]; // avoid dynamic allocation for(unsigned short iDim=0; iDimGetLaminarViscosity(iPoint); CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint), Lam_Visc_Wall); From 219f7e12c4389fa59c48b326108bf064099046be Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Wed, 9 Dec 2020 08:52:28 +0100 Subject: [PATCH 09/18] ComputeStressTensor takes gradient of velocity, not primvar now and conversions of multidim arrays into pointers were modified --- SU2_CFD/include/numerics/CNumerics.hpp | 13 +++++++--- .../include/solvers/CFVMFlowSolverBase.inl | 10 ++----- .../interfaces/fsi/CFlowTractionInterface.cpp | 5 ++-- SU2_CFD/src/numerics/CNumerics.cpp | 26 +++++++++---------- SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp | 2 +- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 6 ++--- SU2_CFD/src/numerics/flow/flow_sources.cpp | 2 +- .../src/numerics/turbulent/turb_sources.cpp | 2 +- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 2 +- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 6 ++--- SU2_CFD/src/solvers/CNSSolver.cpp | 4 +-- 11 files changed, 37 insertions(+), 41 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 094f8c2b2649..57f2bee9ecf3 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -465,23 +465,28 @@ class CNumerics { * \details The parameter primvargrad can be e.g. PrimVar_Grad_i or Mean_GradPrimVar. * \param[in] nDim - 2 or 3 * \param[out] rateofstrain - Rate of strain matrix - * \param[in] primvargrad - A primitive variable gradient matrix. + * \param[in] velgrad - A velocity gradient matrix. */ - static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double* const* primvargrad); + static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double* const* velgrad); /*! * \brief Compute the stress tensor from the velocity gradients. * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result * of this function by (-rho). The argument density is only used if turb_ke is not 0. + * To select the velocity gradient components from a primitive variable gradient PrimVar_Grad_i, + * write PrimVar_Grad_i+1. + * If nDim==2, we use the same formula but only only access the entries [0][0]..[1][1] of + * stress and velgrad. If reynolds3x3 is true, the other non-diagonal entries of stress + * set to zero, and stress[2][2] to some value. * \param[in] nDim - Dimension of the flow problem, 2 or 3 * \param[out] stress - Stress tensor - * \param[in] primvargrad - A primitive variable gradient matrix. + * \param[in] velgrad - A velocity gradient matrix. * \param[in] viscosity - Viscosity * \param[in] density - Density * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2. */ - static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, + static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* velgrad, su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false); /*! diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 6226fa1595ef..82f8031f75f7 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2189,15 +2189,9 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - su2double *Tau_pointer[3]; - su2double *Grad_Vel_pointer[4]; - Grad_Vel_pointer[0] = nullptr; - for(iDim=0;iDim<3;iDim++){ - Tau_pointer[iDim] = Tau[iDim]; - Grad_Vel_pointer[iDim+1] = Grad_Vel[iDim]; - } + su2double *Tau_pointer[3] = {Tau[0], Tau[1], (nDim==3)?Tau[2]:nullptr}; + su2double *Grad_Vel_pointer[3] = {Grad_Vel[0], Grad_Vel[1], (nDim==3)?Grad_Vel[2]:nullptr}; CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_Vel_pointer, Viscosity); - // Grad_Vel is not a primitive variable gradient, so we have to shift the index. /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index e229c04ed7ec..66d179786383 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -189,9 +189,8 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); su2double tau_data[9]; // avoid dynamic allocation - su2double *tau[3]; - tau[0] = tau_data; tau[1] = tau_data+3; tau[2] = tau_data+6; - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow),Viscosity); + su2double *tau[3] = {tau_data, tau_data+3, tau_data+6}; + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,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/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 5b70417150ce..77d160729de9 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -485,26 +485,26 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou AD::EndPassive(wasActive); } -void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double * const* primvargrad){ +void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double * const* velgrad){ /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ if (nDim == 3){ - rateofstrain[0][0] = primvargrad[1][0]; - rateofstrain[1][1] = primvargrad[2][1]; - rateofstrain[2][2] = primvargrad[3][2]; - rateofstrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); - rateofstrain[0][2] = 0.5 * (primvargrad[1][2] + primvargrad[3][0]); - rateofstrain[1][2] = 0.5 * (primvargrad[2][2] + primvargrad[3][1]); + rateofstrain[0][0] = velgrad[0][0]; + rateofstrain[1][1] = velgrad[1][1]; + rateofstrain[2][2] = velgrad[2][2]; + rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); + rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]); + rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]); rateofstrain[1][0] = rateofstrain[0][1]; rateofstrain[2][1] = rateofstrain[1][2]; rateofstrain[2][0] = rateofstrain[0][2]; } else { // nDim==2 - rateofstrain[0][0] = primvargrad[1][0]; - rateofstrain[1][1] = primvargrad[2][1]; + rateofstrain[0][0] = velgrad[0][0]; + rateofstrain[1][1] = velgrad[1][1]; rateofstrain[2][2] = 0.0; - rateofstrain[0][1] = 0.5 * (primvargrad[1][1] + primvargrad[2][0]); + rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); rateofstrain[0][2] = 0.0; rateofstrain[1][2] = 0.0; rateofstrain[1][0] = rateofstrain[0][1]; @@ -513,19 +513,19 @@ void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** r } } -void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* primvargrad, +void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* velgrad, su2double viscosity, su2double density, su2double turb_ke, bool reynolds3x3){ su2double TWO3 = 2.0/3.0; su2double divVel = 0; for (unsigned short iDim = 0; iDim < nDim; iDim++){ - divVel += primvargrad[iDim+1][iDim]; + divVel += velgrad[iDim][iDim]; } for (unsigned short iDim = 0; iDim < nDim; iDim++){ for (unsigned short jDim = 0; jDim < nDim; jDim++){ stress[iDim][jDim] = - viscosity * (primvargrad[iDim+1][jDim]+primvargrad[jDim+1][iDim]) + viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]) - TWO3 * viscosity * divVel * (iDim==jDim) - TWO3 * density * turb_ke * (iDim==jDim); } diff --git a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp index e027a7dea994..ad18af5b5865 100644 --- a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp +++ b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp @@ -289,7 +289,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, } /*--- Compute the viscous stress tensor ---*/ - ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX-1, mu); + ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX, mu); /*--- Populate entries in the viscous flux vector ---*/ for (iDim = 0; iDim < nDim; iDim++) { diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index beaea9567305..fd0630b08ec2 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -130,7 +130,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, * parts of tau can be computed with the total viscosity. --- */ if (using_uq){ - ComputeStressTensor(nDim, tau, val_gradprimvar, val_laminar_viscosity); // laminar part + ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity); // laminar part // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) @@ -138,7 +138,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } else { // compute both parts in one step const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; - ComputeStressTensor(nDim, tau, val_gradprimvar, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? } } @@ -213,7 +213,7 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ su2double meandensity = Mean_PrimVar[nDim+2]; - ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar, Mean_Eddy_Viscosity, meandensity, turb_ke, true); + ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar+1, Mean_Eddy_Viscosity, meandensity, turb_ke, true); for(unsigned short iDim=0; iDim<3; iDim++){ for(unsigned short jDim=0; jDim<3; jDim++){ MeanReynoldsStress[iDim][jDim] /= (-meandensity); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 15d81da84f20..6800da6f91f7 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -202,7 +202,7 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); /*--- The full stress tensor is needed for variable density ---*/ - ComputeStressTensor(nDim, tau, PrimVar_Grad_i, total_viscosity); + ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity); /*--- Viscous terms. ---*/ diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index 3fed0ea69582..3316e4d94e39 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -920,7 +920,7 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i, Eddy_Viscosity_i, Density_i, turb_ke, true); + ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i+1, Eddy_Viscosity_i, Density_i, turb_ke, true); for(unsigned short iDim=0; iDim<3; iDim++){ for(unsigned short jDim=0; jDim<3; jDim++){ MeanReynoldsStress[iDim][jDim] /= (-Density_i); diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index b7aa40670988..b2cea1adc96c 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -783,7 +783,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con else val_turb_ke = 0.0; - CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad, Laminar_Viscosity, Density, val_turb_ke); + CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke); /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 49dcd26f374c..09fa1191ac7a 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -1142,10 +1142,8 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar ++) Res_Visc[iVar] = 0.0; - su2double *Tau_pointer[3]; - for(iDim=0; iDimGetGradient_Primitive(iPoint), total_viscosity); + CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint)+1, total_viscosity); /*--- Dot product of the stress tensor with the grid velocity ---*/ @@ -984,7 +984,7 @@ void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, C tau_pointer[iDim] = tau[iDim]; su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint), Lam_Visc_Wall); + CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint)+1, Lam_Visc_Wall); for (auto iDim = 0u; iDim < nDim; iDim++) { TauElem[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], UnitNormal); From db059dc5b4b4dbcaf57eb6f1dd168523ff8d6dfd Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Thu, 10 Dec 2020 09:43:47 +0100 Subject: [PATCH 10/18] Removed unused variable delta --- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index e9b073bcf381..cea3942bc1fa 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -604,7 +604,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con normal_grad_psi5, normal_grad_T, sigma_partial, Laminar_Viscosity = 0.0, heat_flux_factor, temp_sens = 0.0, *Psi = nullptr, *U = nullptr, Enthalpy, gradPsi5_v, psi5_tau_partial, psi5_tau_grad_vel, source_v_1, Density, Pressure = 0.0, div_vel, val_turb_ke, vartheta, vartheta_partial, psi5_p_div_vel, Omega[3], rho_v[3] = {0.0,0.0,0.0}, - CrossProduct[3], delta[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}, r, ru, rv, rw, rE, p, T, dp_dr, dp_dru, + CrossProduct[3], r, ru, rv, rw, rE, p, T, dp_dr, dp_dru, dp_drv, dp_drw, dp_drE, dH_dr, dH_dru, dH_drv, dH_drw, dH_drE, H, D[3][3], Dd[3], Mach_Inf, eps, scale = 1.0, RefVel2, RefDensity, Mach2Vel, *Velocity_Inf, factor; const su2double* const* GridVel_Grad; From 2ba0c73fd340581d7a6cd2560fc840b8b3641509 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Thu, 10 Dec 2020 10:10:35 +0100 Subject: [PATCH 11/18] Corrected spacing --- Common/include/CConfig.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index dc6f8878b50c..4bcf3d98eb8d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -8831,7 +8831,7 @@ class CConfig { * \return TRUE means that they will be used. */ bool GetUsing_ReynoldsStress(void) const { - return (using_uq); + return (using_uq); } /*! From 62444fde4470036d43b31167a84c941786f64d53 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Thu, 10 Dec 2020 14:03:11 +0100 Subject: [PATCH 12/18] Allocation of CNum::MeanReyStress for UQ not treated separately from the other UQ quantities any more. This partly reverts commit c2ded5af. --- Common/include/CConfig.hpp | 8 -------- SU2_CFD/include/numerics/CNumerics.hpp | 1 - SU2_CFD/src/numerics/CNumerics.cpp | 23 +++++------------------ 3 files changed, 5 insertions(+), 27 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 4bcf3d98eb8d..f7c1580057fa 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -8826,14 +8826,6 @@ class CConfig { */ bool GetPrintInlet_InterpolatedData(void) const { return PrintInlet_InterpolatedData; } - /*! - * \brief Get information about using the Reynolds stress tensor and mean rate of strain matrix. - * \return TRUE means that they will be used. - */ - bool GetUsing_ReynoldsStress(void) const { - return (using_uq); - } - /*! * \brief Get information about using UQ methodology * \return TRUE means that UQ methodology of eigenspace perturbation will be used diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index d8cf33c705a1..f57f0ec3eac4 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -215,7 +215,6 @@ class CNumerics { su2double *l, *m; - bool using_reynoldsstress; /*!< \brief Flag for usage of mean Reynolds stress matrix */ su2double **MeanReynoldsStress; /*!< \brief Mean Reynolds stress tensor */ su2double **MeanPerturbedRSM; /*!< \brief Perturbed Reynolds stress tensor */ bool using_uq, /*!< \brief Flag for UQ methodology */ diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 25aeebc1dba4..0edfa2c0564b 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -49,7 +49,6 @@ CNumerics::CNumerics(void) { l = nullptr; m = nullptr; - using_reynoldsstress = false; using_uq = false; nemo = false; @@ -107,18 +106,10 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Dissipation_ij = 1.0; - /* --- Initializing Reynolds stress matrix --- */ - using_reynoldsstress= config->GetUsing_ReynoldsStress(); - if (using_reynoldsstress){ - MeanReynoldsStress = new su2double* [3]; - for (iDim = 0; iDim < 3; iDim++){ - MeanReynoldsStress[iDim] = new su2double [3]; - } - } - - /* --- Initializing additional variables for the UQ methodology --- */ + /* --- Initializing variables for the UQ methodology --- */ using_uq = config->GetUsing_UQ(); if (using_uq){ + MeanReynoldsStress = new su2double* [3]; MeanPerturbedRSM = new su2double* [3]; A_ij = new su2double* [3]; newA_ij = new su2double* [3]; @@ -129,6 +120,7 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Barycentric_Coord = new su2double [2]; New_Coord = new su2double [2]; for (iDim = 0; iDim < 3; iDim++){ + MeanReynoldsStress[iDim] = new su2double [3]; MeanPerturbedRSM[iDim] = new su2double [3]; A_ij[iDim] = new su2double [3]; newA_ij[iDim] = new su2double [3]; @@ -190,15 +182,9 @@ CNumerics::~CNumerics(void) { delete [] l; delete [] m; - if (using_reynoldsstress){ - for (unsigned short iDim = 0; iDim < 3; iDim++){ - delete [] MeanReynoldsStress[iDim]; - } - delete [] MeanReynoldsStress; - } - if (using_uq) { for (unsigned short iDim = 0; iDim < 3; iDim++){ + delete [] MeanReynoldsStress[iDim]; delete [] MeanPerturbedRSM[iDim]; delete [] A_ij[iDim]; delete [] newA_ij[iDim]; @@ -206,6 +192,7 @@ CNumerics::~CNumerics(void) { delete [] New_Eig_Vec[iDim]; delete [] Corners[iDim]; } + delete [] MeanReynoldsStress; delete [] MeanPerturbedRSM; delete [] A_ij; delete [] newA_ij; From 078427d9e2e7c3b5462b0f9cb233c771f33ef146 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Thu, 10 Dec 2020 15:30:56 +0100 Subject: [PATCH 13/18] ComputeStressTensor, CompMeanRateOfStrMat templated and inline so that they accept arguments of type su2double**, su2double[3][3], etc. --- SU2_CFD/include/numerics/CNumerics.hpp | 62 +++++++++++++++++-- .../include/solvers/CFVMFlowSolverBase.inl | 4 +- .../interfaces/fsi/CFlowTractionInterface.cpp | 3 +- SU2_CFD/src/numerics/CNumerics.cpp | 53 ---------------- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 3 +- SU2_CFD/src/solvers/CNSSolver.cpp | 11 +--- 6 files changed, 62 insertions(+), 74 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index f57f0ec3eac4..1f485d46acd4 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -465,8 +465,37 @@ class CNumerics { * \param[in] nDim - 2 or 3 * \param[out] rateofstrain - Rate of strain matrix * \param[in] velgrad - A velocity gradient matrix. - */ - static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double* const* velgrad); + * \tparam TWOINDICES_1 - any type that supports the [][] interface + * \tparam TWOINDICES_2 - any type that supports the [][] interface + */ + template + inline static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, TWOINDICES_1& rateofstrain, const TWOINDICES_2& velgrad){ + + /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ + + if (nDim == 3){ + rateofstrain[0][0] = velgrad[0][0]; + rateofstrain[1][1] = velgrad[1][1]; + rateofstrain[2][2] = velgrad[2][2]; + rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); + rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]); + rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]); + rateofstrain[1][0] = rateofstrain[0][1]; + rateofstrain[2][1] = rateofstrain[1][2]; + rateofstrain[2][0] = rateofstrain[0][2]; + } + else { // nDim==2 + rateofstrain[0][0] = velgrad[0][0]; + rateofstrain[1][1] = velgrad[1][1]; + rateofstrain[2][2] = 0.0; + rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); + rateofstrain[0][2] = 0.0; + rateofstrain[1][2] = 0.0; + rateofstrain[1][0] = rateofstrain[0][1]; + rateofstrain[2][1] = rateofstrain[1][2]; + rateofstrain[2][0] = rateofstrain[0][2]; + } + } /*! * \brief Compute the stress tensor from the velocity gradients. @@ -484,9 +513,32 @@ class CNumerics { * \param[in] density - Density * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2. - */ - static void ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* velgrad, - su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false); + * \tparam TWOINDICES_1 - any type that supports the [][] interface + * \tparam TWOINDICES_2 - any type that supports the [][] interface + */ + template + inline static void ComputeStressTensor(unsigned short nDim, TWOINDICES_1& stress, const TWOINDICES_2& velgrad, + su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false){ + su2double divVel = 0; + for (unsigned short iDim = 0; iDim < nDim; iDim++){ + divVel += velgrad[iDim][iDim]; + } + + for (unsigned short iDim = 0; iDim < nDim; iDim++){ + for (unsigned short jDim = 0; jDim < nDim; jDim++){ + stress[iDim][jDim] = + viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]) + - 2./3. * viscosity * divVel * (iDim==jDim) + - 2./3. * density * turb_ke * (iDim==jDim); + } + } + + if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix + stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0; + stress[2][2] = -2./3. * ( viscosity * divVel + density * turb_ke ); + } + + } /*! * \brief Set the value of the first blending function. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 9763a0e21547..94b8780ca034 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2189,9 +2189,7 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - su2double *Tau_pointer[3] = {Tau[0], Tau[1], (nDim==3)?Tau[2]:nullptr}; - su2double *Grad_Vel_pointer[3] = {Grad_Vel[0], Grad_Vel[1], (nDim==3)?Grad_Vel[2]:nullptr}; - CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_Vel_pointer, Viscosity); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity); /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index f748a68c4443..8aa5415b4408 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -188,8 +188,7 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); - su2double tau_data[9]; // avoid dynamic allocation - su2double *tau[3] = {tau_data, tau_data+3, tau_data+6}; + su2double tau[3][3]; CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 0edfa2c0564b..e337969b5cd6 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -472,59 +472,6 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou AD::EndPassive(wasActive); } -void CNumerics::ComputeMeanRateOfStrainMatrix(unsigned short nDim, su2double** rateofstrain, const su2double * const* velgrad){ - - /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ - - if (nDim == 3){ - rateofstrain[0][0] = velgrad[0][0]; - rateofstrain[1][1] = velgrad[1][1]; - rateofstrain[2][2] = velgrad[2][2]; - rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); - rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]); - rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]); - rateofstrain[1][0] = rateofstrain[0][1]; - rateofstrain[2][1] = rateofstrain[1][2]; - rateofstrain[2][0] = rateofstrain[0][2]; - } - else { // nDim==2 - rateofstrain[0][0] = velgrad[0][0]; - rateofstrain[1][1] = velgrad[1][1]; - rateofstrain[2][2] = 0.0; - rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]); - rateofstrain[0][2] = 0.0; - rateofstrain[1][2] = 0.0; - rateofstrain[1][0] = rateofstrain[0][1]; - rateofstrain[2][1] = rateofstrain[1][2]; - rateofstrain[2][0] = rateofstrain[0][2]; - } -} - -void CNumerics::ComputeStressTensor(unsigned short nDim, su2double** stress, const su2double* const* velgrad, - su2double viscosity, su2double density, su2double turb_ke, bool reynolds3x3){ - su2double TWO3 = 2.0/3.0; - - su2double divVel = 0; - for (unsigned short iDim = 0; iDim < nDim; iDim++){ - divVel += velgrad[iDim][iDim]; - } - - for (unsigned short iDim = 0; iDim < nDim; iDim++){ - for (unsigned short jDim = 0; jDim < nDim; jDim++){ - stress[iDim][jDim] = - viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]) - - TWO3 * viscosity * divVel * (iDim==jDim) - - TWO3 * density * turb_ke * (iDim==jDim); - } - } - - if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix - stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0; - stress[2][2] = -TWO3 * ( viscosity * divVel + density * turb_ke ); - } - -} - void CNumerics::GetPreconditioner(const su2double *val_density, const su2double *val_velocity, const su2double *val_betainc2, const su2double *val_cp, const su2double *val_temperature, const su2double *val_drhodt, diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 6835040074ae..c8add4efc55d 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -1142,8 +1142,7 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar ++) Res_Visc[iVar] = 0.0; - su2double *Tau_pointer[3] = {Tau[0], Tau[1], (nDim==3)?Tau[2]:nullptr}; - CNumerics::ComputeStressTensor(nDim, Tau_pointer, Grad_PrimVar+VEL_INDEX, Viscosity); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity); for (iDim = 0; iDim < nDim; iDim++) { TauElem[iDim] = 0.0; for (jDim = 0; jDim < nDim; jDim++) diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 7f0b38831b06..f67e22a19ea6 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -474,10 +474,7 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne /*--- Compute the viscous stress tensor ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; - su2double *tau_pointer[MAXNDIM]; // avoid dynamic allocation - for(unsigned short iDim=0; iDimGetGradient_Primitive(iPoint)+1, total_viscosity); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetGradient_Primitive(iPoint)+1, total_viscosity); /*--- Dot product of the stress tensor with the grid velocity ---*/ @@ -979,12 +976,8 @@ void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, C by using the stress tensor on the surface ---*/ su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0}; - su2double *tau_pointer[MAXNDIM]; // avoid dynamic allocation - for(unsigned short iDim=0; iDimGetLaminarViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau_pointer, nodes->GetGradient_Primitive(iPoint)+1, Lam_Visc_Wall); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetGradient_Primitive(iPoint)+1, Lam_Visc_Wall); for (auto iDim = 0u; iDim < nDim; iDim++) { TauElem[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], UnitNormal); From 2e90bae2ddaba9c708b8a8dc27fdca797df6a9c1 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Thu, 10 Dec 2020 15:58:42 +0100 Subject: [PATCH 14/18] Subtraction from the stress tensor diagonal modified --- SU2_CFD/include/numerics/CNumerics.hpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 1f485d46acd4..368fb4cd384b 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -523,19 +523,18 @@ class CNumerics { for (unsigned short iDim = 0; iDim < nDim; iDim++){ divVel += velgrad[iDim][iDim]; } + su2double pTerm = 2./3. * (divVel * viscosity + density * turb_ke); for (unsigned short iDim = 0; iDim < nDim; iDim++){ for (unsigned short jDim = 0; jDim < nDim; jDim++){ - stress[iDim][jDim] = - viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]) - - 2./3. * viscosity * divVel * (iDim==jDim) - - 2./3. * density * turb_ke * (iDim==jDim); + stress[iDim][jDim] = viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]); } + stress[iDim][iDim] -= pTerm; } if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0; - stress[2][2] = -2./3. * ( viscosity * divVel + density * turb_ke ); + stress[2][2] = -pTerm; } } From 64bd077aa320efcd74bace14dfa14cdd57b60302 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Fri, 11 Dec 2020 08:28:45 +0100 Subject: [PATCH 15/18] Using CNum::CompStressT in CSolver::CompVertexTractions --- SU2_CFD/src/solvers/CSolver.cpp | 27 ++++----------------------- 1 file changed, 4 insertions(+), 23 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 8b4f08df584d..6743086dca2f 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4142,11 +4142,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ (config->GetKind_Solver() == DISC_ADJ_RANS)); // Parameters for the calculations - su2double Pn = 0.0, div_vel = 0.0; - su2double Viscosity = 0.0; - su2double Tau[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}}; - su2double Grad_Vel[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}}; - su2double delta[3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{0.0, 0.0, 1.0}}; + su2double Pn = 0.0; su2double auxForce[3] = {1.0, 0.0, 0.0}; unsigned short iMarker; @@ -4195,26 +4191,11 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ // Calculate tn in the fluid nodes for the viscous term if (viscous_flow) { - - Viscosity = base_nodes->GetLaminarViscosity(iPoint); - - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Grad_Vel[iDim][jDim] = base_nodes->GetGradient_Primitive(iPoint, iDim+1, jDim); - } - } - - // Divergence of the velocity - div_vel = 0.0; for (iDim = 0; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim]; - + su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); + su2double Tau[3][3]; + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity); for (iDim = 0; iDim < nDim; iDim++) { for (jDim = 0 ; jDim < nDim; jDim++) { - - // Viscous stress - Tau[iDim][jDim] = Viscosity*(Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) - - TWO3*Viscosity*div_vel*delta[iDim][jDim]; - - // Viscous component in the tn vector --> Units of force (non-dimensional). auxForce[iDim] += Tau[iDim][jDim]*iNormal[jDim]; } } From 8305f6606029da8077fcdb1d4612189949e6ec04 Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Fri, 11 Dec 2020 09:32:44 +0100 Subject: [PATCH 16/18] Updated regression tests (rans_uq and two others) There were minor differences in the residuals, probably because round-off errors in the computation of the stress tensor accumulate over the solver iterations. turb_naca0012_1c, _2c, _p1c1, _p1c2 in serial, parallel, hybrid regression poiseuille_profile in serial regression stat_fsi in hybrid regression --- TestCases/hybrid_regression.py | 10 +++++----- TestCases/parallel_regression.py | 8 ++++---- TestCases/serial_regression.py | 10 +++++----- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 5f0ed7eed6be..054a7ba1eb48 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -232,7 +232,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.979389, 1.140070, 1.211965, 0.194237] + turb_naca0012_1c.test_vals = [-4.978913, 1.140283, 1.211887, 0.194208] test_list.append(turb_naca0012_1c) # NACA0012 2c @@ -240,7 +240,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.484195, 0.969789, 1.310525, 0.231240] + turb_naca0012_2c.test_vals = [-5.484204, 0.969783, 1.311085, 0.231449] test_list.append(turb_naca0012_2c) # NACA0012 3c @@ -256,7 +256,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.132081, 1.076462, 1.178093, 0.181595] + turb_naca0012_p1c1.test_vals = [-5.132991, 1.076082, 1.177974, 0.181556] test_list.append(turb_naca0012_p1c1) # NACA0012 p1c2 @@ -264,7 +264,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.556648, 0.945129, 1.240986, 0.205071] + turb_naca0012_p1c2.test_vals = [-5.556581, 0.945167, 1.240884, 0.205031] test_list.append(turb_naca0012_p1c2) ###################################### @@ -527,7 +527,7 @@ def main(): stat_fsi.cfg_dir = "fea_fsi/stat_fsi" stat_fsi.cfg_file = "config.cfg" stat_fsi.test_iter = 7 - stat_fsi.test_vals = [-3.242834, -4.866608, 0.000000, 11.000000] + stat_fsi.test_vals = [-3.242851, -4.866383, 0.000000, 11.000000] stat_fsi.multizone = True test_list.append(stat_fsi) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index d7063c7ef4be..fd04078b8f2d 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -683,7 +683,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.973124, 1.141759, 0.861168, 0.014208] + turb_naca0012_1c.test_vals = [-4.973120, 1.141760, 0.861180, 0.014211] turb_naca0012_1c.su2_exec = "parallel_computation.py -f" turb_naca0012_1c.timeout = 1600 turb_naca0012_1c.tol = 0.00001 @@ -694,7 +694,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.484227, 0.967174, 0.901039, 0.029201] + turb_naca0012_2c.test_vals = [-5.484231, 0.967172, 0.901005, 0.029186] turb_naca0012_2c.su2_exec = "parallel_computation.py -f" turb_naca0012_2c.timeout = 1600 turb_naca0012_2c.tol = 0.00001 @@ -716,7 +716,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.128931, 1.076207, 0.954587, 0.040176] + turb_naca0012_p1c1.test_vals = [-5.128763, 1.076245, 0.954666, 0.040210] turb_naca0012_p1c1.su2_exec = "parallel_computation.py -f" turb_naca0012_p1c1.timeout = 1600 turb_naca0012_p1c1.tol = 0.00001 @@ -727,7 +727,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.556651, 0.941872, 0.950627, 0.041946] + turb_naca0012_p1c2.test_vals = [-5.556632, 0.941881, 0.950638, 0.041949] turb_naca0012_p1c2.su2_exec = "parallel_computation.py -f" turb_naca0012_p1c2.timeout = 1600 turb_naca0012_p1c2.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index f5e5d62f60ea..7a15a21d328e 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -250,7 +250,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.494705, -7.711759, -0.000000, 2.085796] #last 4 columns + poiseuille_profile.test_vals = [-12.494720, -7.711373, -0.000000, 2.085796] #last 4 columns poiseuille_profile.su2_exec = "SU2_CFD" poiseuille_profile.new_output = True poiseuille_profile.timeout = 1600 @@ -806,7 +806,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.978069, 1.139128, 0.807062, 0.064726] #last 4 columns + turb_naca0012_1c.test_vals = [-4.978068, 1.139123, 0.806952, 0.064685] #last 4 columns turb_naca0012_1c.su2_exec = "SU2_CFD" turb_naca0012_1c.new_output = True turb_naca0012_1c.timeout = 1600 @@ -818,7 +818,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.484279, 0.967025, 0.823201, 0.070866] #last 4 columns + turb_naca0012_2c.test_vals = [-5.484282, 0.967023, 0.823129, 0.070840] #last 4 columns turb_naca0012_2c.su2_exec = "SU2_CFD" turb_naca0012_2c.new_output = True turb_naca0012_2c.timeout = 1600 @@ -842,7 +842,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.126796, 1.076577, 0.783116, 0.055988] #last 4 columns + turb_naca0012_p1c1.test_vals = [ -5.126540, 1.076620, 0.783153, 0.056001] #last 4 columns turb_naca0012_p1c1.su2_exec = "SU2_CFD" turb_naca0012_p1c1.new_output = True turb_naca0012_p1c1.timeout = 1600 @@ -854,7 +854,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.556585, 0.941677, 0.796006, 0.060817] #last 4 columns + turb_naca0012_p1c2.test_vals = [-5.556554, 0.941694, 0.795964, 0.060801] #last 4 columns turb_naca0012_p1c2.su2_exec = "SU2_CFD" turb_naca0012_p1c2.new_output = True turb_naca0012_p1c2.timeout = 1600 From 09313dc210441ce7c4a60012abcd9777e6c2f4ee Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Fri, 11 Dec 2020 10:16:36 +0100 Subject: [PATCH 17/18] Forgot one regression test in 8305f66 poiseuille_profile in parallel_regression.py --- TestCases/parallel_regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index fd04078b8f2d..08bcb40c6b62 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -225,7 +225,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.492876, -7.672445, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.492859, -7.672756, -0.000000, 2.085796] poiseuille_profile.su2_exec = "parallel_computation.py -f" poiseuille_profile.timeout = 1600 poiseuille_profile.tol = 0.00001 From 8879c13b7e298ae1698fc083d02d66bb521940fc Mon Sep 17 00:00:00 2001 From: Max Aehle Date: Fri, 11 Dec 2020 11:57:31 +0100 Subject: [PATCH 18/18] Using CNum::CompStressT in python_wrapper_structure.cpp this changes the arithmetics --- SU2_CFD/src/python_wrapper_structure.cpp | 33 ++++++------------------ 1 file changed, 8 insertions(+), 25 deletions(-) diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 672073b26873..ec823ff9d19b 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -356,16 +356,9 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) /*--- Parameters for the calculations ---*/ // Pn: Pressure // Pinf: Pressure_infinite - // div_vel: Velocity divergence - // Dij: Dirac delta - su2double Pn = 0.0, div_vel = 0.0, Dij = 0.0; + su2double Pn = 0.0; su2double Viscosity = 0.0; - su2double Grad_Vel[3][3] = { {0.0, 0.0, 0.0} , - {0.0, 0.0, 0.0} , - {0.0, 0.0, 0.0} } ; - su2double Tau[3][3] = { {0.0, 0.0, 0.0} , - {0.0, 0.0, 0.0} , - {0.0, 0.0, 0.0} } ; + su2double Tau[3][3] = {{0.0}}; su2double Pinf = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetPressure_Inf(); @@ -384,11 +377,6 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) /*--- Get the values of pressure and viscosity ---*/ Pn = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetPressure(iPoint); if (viscous_flow) { - for(iDim=0; iDimGetNodes()->GetGradient_Primitive(iPoint, iDim+1, jDim); - } - } Viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); } @@ -399,23 +387,18 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) /*--- Calculate the viscous (shear stress) part of tn in the fluid nodes (force units) ---*/ if ((incompressible || compressible) && viscous_flow) { - div_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - div_vel += Grad_Vel[iDim][iDim]; - if (incompressible) div_vel = 0.0; - + CNumerics::ComputeStressTensor(nDim, Tau, + solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)+1, Viscosity); for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Dij = 0.0; if (iDim == jDim) Dij = 1.0; - Tau[iDim][jDim] = Viscosity*(Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) - TWO3*Viscosity*div_vel*Dij; - PyWrapNodalForce[iDim] += Tau[iDim][jDim]*Normal[jDim]; + for (jDim = 0 ; jDim < nDim; jDim++) { + PyWrapNodalForce[iDim] += Tau[iDim][jDim]*Normal[jDim]; } } } //Divide by local are in case of force density communication. - for(iDim = 0; iDim < nDim; iDim++) { - PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area; + for(iDim = 0; iDim < nDim; iDim++) { + PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area; } halo = false;