From 4f8fcccf6e62f98f28f60fd8a3ce012c77a67b5d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 27 Jun 2021 13:25:37 +0100 Subject: [PATCH 1/8] create a type for "matrix views" (wraps a pointer) to avoid ** --- .../containers/container_decorators.hpp | 87 ++++++++++++------- Common/include/geometry/dual_grid/CPoint.hpp | 2 +- SU2_CFD/include/numerics/CNumerics.hpp | 71 +++++++-------- SU2_CFD/include/numerics/radiation.hpp | 9 +- SU2_CFD/include/output/CFlowOutput.hpp | 34 +++++++- .../include/solvers/CFVMFlowSolverBase.hpp | 34 ++++---- .../include/solvers/CFVMFlowSolverBase.inl | 8 +- .../include/variables/CAdjEulerVariable.hpp | 4 +- SU2_CFD/include/variables/CEulerVariable.hpp | 6 +- SU2_CFD/include/variables/CHeatVariable.hpp | 2 +- .../include/variables/CIncEulerVariable.hpp | 6 +- .../include/variables/CNEMOEulerVariable.hpp | 6 +- SU2_CFD/include/variables/CTurbSAVariable.hpp | 2 +- SU2_CFD/include/variables/CTurbVariable.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 10 +-- .../src/integration/CMultiGridIntegration.cpp | 3 +- .../interfaces/fsi/CFlowTractionInterface.cpp | 2 +- SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp | 4 +- SU2_CFD/src/output/CFlowCompOutput.cpp | 2 +- SU2_CFD/src/output/CFlowIncOutput.cpp | 2 +- SU2_CFD/src/output/CFlowOutput.cpp | 33 ------- SU2_CFD/src/solvers/CAdjEulerSolver.cpp | 10 +-- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 17 ++-- SU2_CFD/src/solvers/CAdjTurbSolver.cpp | 16 ++-- SU2_CFD/src/solvers/CHeatSolver.cpp | 22 ++--- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 12 ++- SU2_CFD/src/solvers/CTurbSASolver.cpp | 4 +- SU2_CFD/src/variables/CTurbSAVariable.cpp | 2 +- 28 files changed, 211 insertions(+), 201 deletions(-) diff --git a/Common/include/containers/container_decorators.hpp b/Common/include/containers/container_decorators.hpp index af46dff61eb..8e4f72a2c1f 100644 --- a/Common/include/containers/container_decorators.hpp +++ b/Common/include/containers/container_decorators.hpp @@ -30,6 +30,43 @@ #include "C2DContainer.hpp" +/*! + * \brief Class to represent a matrix (without owning the data, this just wraps a pointer). + */ +template +class CMatrixView { +public: + using Scalar = typename std::remove_const::type; + using Index = unsigned long; + +private: + T* m_ptr; + Index m_cols; + +public: + CMatrixView(T* ptr = nullptr, Index cols = 0) : m_ptr(ptr), m_cols(cols) {} + + template friend class CMatrixView; + template + CMatrixView(const CMatrixView& other) : m_ptr(other.m_ptr), m_cols(other.m_cols) {} + + explicit CMatrixView(const su2matrix& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {} + + template::value> = 0> + explicit CMatrixView(su2matrix& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {} + + const Scalar* operator[] (Index i) const noexcept { return &m_ptr[i*m_cols]; } + const Scalar& operator() (Index i, Index j) const noexcept { return m_ptr[i*m_cols + j]; } + + template::value> = 0> + Scalar* operator[] (Index i) noexcept { return &m_ptr[i*m_cols]; } + + template::value> = 0> + Scalar& operator() (Index i, Index j) noexcept { return m_ptr[i*m_cols + j]; } + + friend CMatrixView operator+ (CMatrixView mv, Index incr) { return CMatrixView(mv[incr], mv.m_cols); } +}; + /*! * \class C3DContainerDecorator * \brief Decorate a matrix type (Storage) with 3 dimensions. @@ -44,6 +81,9 @@ class C3DContainerDecorator { static constexpr bool IsRowMajor = true; static constexpr bool IsColumnMajor = false; + using Matrix = CMatrixView; + using ConstMatrix = CMatrixView; + using CInnerIter = typename Storage::CInnerIter; template using CInnerIterGather = typename Storage::template CInnerIterGather >; @@ -78,6 +118,18 @@ class C3DContainerDecorator { Scalar& operator() (Index i, Index j, Index k) noexcept { return m_storage(i, j*m_innerSz + k); } const Scalar& operator() (Index i, Index j, Index k) const noexcept { return m_storage(i, j*m_innerSz + k); } + /*! + * \brief Matrix access. + */ + Matrix operator[] (Index i) noexcept { return Matrix(m_storage[i], m_innerSz); } + ConstMatrix operator[] (Index i) const noexcept { return ConstMatrix(m_storage[i], m_innerSz); } + + /*! + * \brief Matrix access with an offset. + */ + Matrix operator() (Index i, Index j) noexcept { return Matrix(m_storage[i]+j*m_innerSz, m_innerSz); } + ConstMatrix operator() (Index i, Index j) const noexcept { return ConstMatrix(m_storage[i]+j*m_innerSz, m_innerSz); } + /*! * \brief Get a scalar iterator to the inner-most dimension of the container. */ @@ -105,42 +157,11 @@ class C3DContainerDecorator { }; /*! - * \brief Some typedefs for the + * \brief Some typedefs for 3D containers */ using C3DIntMatrix = C3DContainerDecorator >; using C3DDoubleMatrix = C3DContainerDecorator; - -/*! - * \class CVectorOfMatrix - * \brief This contrived container is used to store small matrices in a contiguous manner - * but still present the "su2double**" interface to the outside world. - * The "interface" part should be replaced by something more efficient, e.g. a "matrix view". - */ -class CVectorOfMatrix: public C3DDoubleMatrix { -private: - su2matrix interface; - -public: - CVectorOfMatrix() = default; - - CVectorOfMatrix(Index length, Index rows, Index cols, Scalar value = 0) noexcept { - resize(length, rows, cols, value); - } - - void resize(Index length, Index rows, Index cols, Scalar value = 0) noexcept { - C3DDoubleMatrix::resize(length, rows, cols, value); - interface.resize(length,rows); - for(Index i=0; i GetGridVel_Grad(unsigned long iPoint) const { return GridVel_Grad[iPoint]; } /*! * \brief Set the adjoint values of the (geometric) coordinates. diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 7e5fd557df6..2b8f8b0a62f 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -158,28 +158,22 @@ class CNumerics { su2double *TurbPsi_i, /*!< \brief Vector of adjoint turbulent variables at point i. */ *TurbPsi_j; /*!< \brief Vector of adjoint turbulent variables at point j. */ - su2double - **ConsVar_Grad_i, /*!< \brief Gradient of conservative variables at point i. */ - **ConsVar_Grad_j, /*!< \brief Gradient of conservative variables at point j. */ - **ConsVar_Grad; /*!< \brief Gradient of conservative variables which is a scalar. */ - su2double - **PrimVar_Grad_i, /*!< \brief Gradient of primitive variables at point i. */ - **PrimVar_Grad_j; /*!< \brief Gradient of primitive variables at point j. */ - su2double - **PsiVar_Grad_i, /*!< \brief Gradient of adjoint variables at point i. */ - **PsiVar_Grad_j; /*!< \brief Gradient of adjoint variables at point j. */ - su2double - **TurbVar_Grad_i, /*!< \brief Gradient of turbulent variables at point i. */ - **TurbVar_Grad_j; /*!< \brief Gradient of turbulent variables at point j. */ - su2double - **TransVar_Grad_i, /*!< \brief Gradient of turbulent variables at point i. */ - **TransVar_Grad_j; /*!< \brief Gradient of turbulent variables at point j. */ - su2double - **TurbPsi_Grad_i, /*!< \brief Gradient of adjoint turbulent variables at point i. */ - **TurbPsi_Grad_j; /*!< \brief Gradient of adjoint turbulent variables at point j. */ - su2double - **AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ - **AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ + CMatrixView + ConsVar_Grad_i, /*!< \brief Gradient of conservative variables at point i. */ + ConsVar_Grad_j, /*!< \brief Gradient of conservative variables at point j. */ + ConsVar_Grad, /*!< \brief Gradient of conservative variables which is a scalar. */ + PrimVar_Grad_i, /*!< \brief Gradient of primitive variables at point i. */ + PrimVar_Grad_j, /*!< \brief Gradient of primitive variables at point j. */ + PsiVar_Grad_i, /*!< \brief Gradient of adjoint variables at point i. */ + PsiVar_Grad_j, /*!< \brief Gradient of adjoint variables at point j. */ + TurbVar_Grad_i, /*!< \brief Gradient of turbulent variables at point i. */ + TurbVar_Grad_j, /*!< \brief Gradient of turbulent variables at point j. */ + TransVar_Grad_i, /*!< \brief Gradient of turbulent variables at point i. */ + TransVar_Grad_j, /*!< \brief Gradient of turbulent variables at point j. */ + TurbPsi_Grad_i, /*!< \brief Gradient of adjoint turbulent variables at point i. */ + TurbPsi_Grad_j, /*!< \brief Gradient of adjoint turbulent variables at point j. */ + AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ + AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ su2double *Coord_i, /*!< \brief Cartesians coordinates of point i. */ @@ -334,8 +328,8 @@ class CNumerics { * \param[in] val_consvar_grad_i - Gradient of the conservative variable at point i. * \param[in] val_consvar_grad_j - Gradient of the conservative variable at point j. */ - inline void SetConsVarGradient(su2double **val_consvar_grad_i, - su2double **val_consvar_grad_j) { + inline void SetConsVarGradient(CMatrixView val_consvar_grad_i, + CMatrixView val_consvar_grad_j) { ConsVar_Grad_i = val_consvar_grad_i; ConsVar_Grad_j = val_consvar_grad_j; } @@ -344,15 +338,15 @@ class CNumerics { * \brief Set the gradient of the conservative variables. * \param[in] val_consvar_grad - Gradient of the conservative variable which is a scalar. */ - inline void SetConsVarGradient(su2double **val_consvar_grad) { ConsVar_Grad = val_consvar_grad; } + inline void SetConsVarGradient(CMatrixView val_consvar_grad) { ConsVar_Grad = val_consvar_grad; } /*! * \brief Set the gradient of the primitive variables. * \param[in] val_primvar_grad_i - Gradient of the primitive variable at point i. * \param[in] val_primvar_grad_j - Gradient of the primitive variable at point j. */ - void SetPrimVarGradient(su2double **val_primvar_grad_i, - su2double **val_primvar_grad_j) { + void SetPrimVarGradient(CMatrixView val_primvar_grad_i, + CMatrixView val_primvar_grad_j) { PrimVar_Grad_i = val_primvar_grad_i; PrimVar_Grad_j = val_primvar_grad_j; } @@ -372,8 +366,8 @@ class CNumerics { * \param[in] val_psivar_grad_i - Gradient of the adjoint variable at point i. * \param[in] val_psivar_grad_j - Gradient of the adjoint variable at point j. */ - inline void SetAdjointVarGradient(su2double **val_psivar_grad_i, - su2double **val_psivar_grad_j) { + inline void SetAdjointVarGradient(CMatrixView val_psivar_grad_i, + CMatrixView val_psivar_grad_j) { PsiVar_Grad_i = val_psivar_grad_i; PsiVar_Grad_j = val_psivar_grad_j; } @@ -403,8 +397,8 @@ class CNumerics { * \param[in] val_turbvar_grad_i - Gradient of the turbulent variable at point i. * \param[in] val_turbvar_grad_j - Gradient of the turbulent variable at point j. */ - inline void SetTurbVarGradient(su2double **val_turbvar_grad_i, - su2double **val_turbvar_grad_j) { + inline void SetTurbVarGradient(CMatrixView val_turbvar_grad_i, + CMatrixView val_turbvar_grad_j) { TurbVar_Grad_i = val_turbvar_grad_i; TurbVar_Grad_j = val_turbvar_grad_j; } @@ -414,8 +408,8 @@ class CNumerics { * \param[in] val_turbvar_grad_i - Gradient of the turbulent variable at point i. * \param[in] val_turbvar_grad_j - Gradient of the turbulent variable at point j. */ - inline void SetTransVarGradient(su2double **val_transvar_grad_i, - su2double **val_transvar_grad_j) { + inline void SetTransVarGradient(CMatrixView val_transvar_grad_i, + CMatrixView val_transvar_grad_j) { TransVar_Grad_i = val_transvar_grad_i; TransVar_Grad_j = val_transvar_grad_j; } @@ -435,8 +429,8 @@ class CNumerics { * \param[in] val_turbpsivar_grad_i - Gradient of the adjoint turbulent variable at point i. * \param[in] val_turbpsivar_grad_j - Gradient of the adjoint turbulent variable at point j. */ - inline void SetTurbAdjointGradient(su2double **val_turbpsivar_grad_i, - su2double **val_turbpsivar_grad_j) { + inline void SetTurbAdjointGradient(CMatrixView val_turbpsivar_grad_i, + CMatrixView val_turbpsivar_grad_j) { TurbPsi_Grad_i = val_turbpsivar_grad_i; TurbPsi_Grad_j = val_turbpsivar_grad_j; } @@ -685,8 +679,8 @@ class CNumerics { * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i. * \param[in] val_auxvar_grad_j - Gradient of the auxiliary variable at point j. */ - inline void SetAuxVarGrad(su2double **val_auxvar_grad_i, - su2double **val_auxvar_grad_j) { + inline void SetAuxVarGrad(CMatrixView val_auxvar_grad_i, + CMatrixView val_auxvar_grad_j) { AuxVar_Grad_i = val_auxvar_grad_i; AuxVar_Grad_j = val_auxvar_grad_j; } @@ -1549,7 +1543,8 @@ class CNumerics { * \param[in] val_radvar_grad_i - Gradient of the turbulent variable at point i. * \param[in] val_radvar_grad_j - Gradient of the turbulent variable at point j. */ - inline virtual void SetRadVarGradient(const su2double* const* val_radvar_grad_i, const su2double* const* val_radvar_grad_j) { } + inline virtual void SetRadVarGradient(CMatrixView val_radvar_grad_i, + CMatrixView val_radvar_grad_j) { } /*! * \brief Set the gradient of the radiation variables. diff --git a/SU2_CFD/include/numerics/radiation.hpp b/SU2_CFD/include/numerics/radiation.hpp index 4d1cd10ea8f..bbff8fc9191 100644 --- a/SU2_CFD/include/numerics/radiation.hpp +++ b/SU2_CFD/include/numerics/radiation.hpp @@ -35,8 +35,8 @@ class CNumericsRadiation : public CNumerics { bool implicit, incompressible; const su2double *RadVar_i; /*!< \brief Vector of radiation variables at point i. */ const su2double *RadVar_j; /*!< \brief Vector of radiation variables at point j. */ - const su2double* const* RadVar_Grad_i; /*!< \brief Gradient of turbulent variables at point i. */ - const su2double* const* RadVar_Grad_j; /*!< \brief Gradient of turbulent variables at point j. */ + CMatrixView RadVar_Grad_i; /*!< \brief Gradient of turbulent variables at point i. */ + CMatrixView RadVar_Grad_j; /*!< \brief Gradient of turbulent variables at point j. */ su2double Absorption_Coeff; /*!< \brief Absorption coefficient. */ su2double Scattering_Coeff; /*!< \brief Scattering coefficient. */ @@ -61,14 +61,13 @@ class CNumericsRadiation : public CNumerics { RadVar_j = val_radvar_j; } - /*! * \brief Set the gradient of the radiation variables. * \param[in] val_radvar_grad_i - Gradient of the turbulent variable at point i. * \param[in] val_radvar_grad_j - Gradient of the turbulent variable at point j. */ - inline void SetRadVarGradient(const su2double* const* val_radvar_grad_i, - const su2double* const* val_radvar_grad_j) final { + inline void SetRadVarGradient(CMatrixView val_radvar_grad_i, + CMatrixView val_radvar_grad_j) final { RadVar_Grad_i = val_radvar_grad_i; RadVar_Grad_j = val_radvar_grad_j; } diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp index 81168fb872d..e16a694fbae 100644 --- a/SU2_CFD/include/output/CFlowOutput.hpp +++ b/SU2_CFD/include/output/CFlowOutput.hpp @@ -94,7 +94,39 @@ class CFlowOutput : public CFVMOutput{ * \param[in] VelocityGradient - Velocity gradients * \return Value of the Q criteration at the node */ - su2double GetQ_Criterion(su2double** VelocityGradient) const; + template + su2double GetQ_Criterion(const T& VelocityGradient) const { + + /*--- Make a 3D copy of the gradient so we do not have worry about nDim ---*/ + + su2double Grad_Vel[3][3] = {{0.0}}; + + for (unsigned short iDim = 0; iDim < nDim; iDim++) + for (unsigned short jDim = 0 ; jDim < nDim; jDim++) + Grad_Vel[iDim][jDim] = VelocityGradient[iDim][jDim]; + + /*--- Q Criterion Eq 1.2 of HALLER, G. (2005). An objective definition of a vortex. + Journal of Fluid Mechanics, 525, 1-26. doi:10.1017/S0022112004002526 ---*/ + + /*--- Components of the strain rate tensor (symmetric) ---*/ + su2double s11 = Grad_Vel[0][0]; + su2double s12 = 0.5 * (Grad_Vel[0][1] + Grad_Vel[1][0]); + su2double s13 = 0.5 * (Grad_Vel[0][2] + Grad_Vel[2][0]); + su2double s22 = Grad_Vel[1][1]; + su2double s23 = 0.5 * (Grad_Vel[1][2] + Grad_Vel[2][1]); + su2double s33 = Grad_Vel[2][2]; + + /*--- Components of the spin tensor (skew-symmetric) ---*/ + su2double omega12 = 0.5 * (Grad_Vel[0][1] - Grad_Vel[1][0]); + su2double omega13 = 0.5 * (Grad_Vel[0][2] - Grad_Vel[2][0]); + su2double omega23 = 0.5 * (Grad_Vel[1][2] - Grad_Vel[2][1]); + + /*--- Q = ||Omega|| - ||Strain|| ---*/ + su2double Q = 2*(pow(omega12,2) + pow(omega13,2) + pow(omega23,2)) - + (pow(s11,2) + pow(s22,2) + pow(s33,2) + 2*(pow(s12,2) + pow(s13,2) + pow(s23,2))); + + return Q; + } /*! * \brief Write information to meta data file diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 6aecc7b7e6d..089b11fb852 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -997,12 +997,11 @@ class CFVMFlowSolverBase : public CSolver { /*! * \brief Evaluate the vorticity and strain rate magnitude. - * \tparam VelocityOffset: Index in the primitive variables where the velocity starts. + * \tparam VelocityOffset - Index in the primitive variables where the velocity starts. */ - template + template void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) { - const auto& Gradient_Primitive = nodes->GetGradient_Primitive(); auto& StrainMag = nodes->GetStrainMag(); ompMasterAssignBarrier(StrainMag_Max,0.0, Omega_Max,0.0); @@ -1012,31 +1011,28 @@ class CFVMFlowSolverBase : public CSolver { SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { - constexpr size_t u = VelocityOffset; - constexpr size_t v = VelocityOffset+1; - constexpr size_t w = VelocityOffset+2; + const auto VelocityGradient = nodes->GetGradient_Primitive(iPoint, VelocityOffset); + auto Vorticity = nodes->GetVorticity(iPoint); /*--- Vorticity ---*/ - su2double* Vorticity = nodes->GetVorticity(iPoint); - - Vorticity[0] = 0.0; Vorticity[1] = 0.0; - - Vorticity[2] = Gradient_Primitive(iPoint,v,0)-Gradient_Primitive(iPoint,u,1); + Vorticity[0] = 0.0; + Vorticity[1] = 0.0; + Vorticity[2] = VelocityGradient(1,0)-VelocityGradient(0,1); if (nDim == 3) { - Vorticity[0] = Gradient_Primitive(iPoint,w,1)-Gradient_Primitive(iPoint,v,2); - Vorticity[1] = -(Gradient_Primitive(iPoint,w,0)-Gradient_Primitive(iPoint,u,2)); + Vorticity[0] = VelocityGradient(2,1)-VelocityGradient(1,2); + Vorticity[1] = -(VelocityGradient(2,0)-VelocityGradient(0,2)); } /*--- Strain Magnitude ---*/ AD::StartPreacc(); - AD::SetPreaccIn(&Gradient_Primitive[iPoint][VelocityOffset], nDim, nDim); + AD::SetPreaccIn(VelocityGradient, nDim, nDim); su2double Div = 0.0; for (unsigned long iDim = 0; iDim < nDim; iDim++) - Div += Gradient_Primitive(iPoint, iDim+VelocityOffset, iDim); + Div += VelocityGradient(iDim, iDim); Div /= 3.0; StrainMag(iPoint) = 0.0; @@ -1044,7 +1040,7 @@ class CFVMFlowSolverBase : public CSolver { /*--- Add diagonal part ---*/ for (unsigned long iDim = 0; iDim < nDim; iDim++) { - StrainMag(iPoint) += pow(Gradient_Primitive(iPoint, iDim+VelocityOffset, iDim) - Div, 2); + StrainMag(iPoint) += pow(VelocityGradient(iDim, iDim) - Div, 2); } if (nDim == 2) { StrainMag(iPoint) += pow(Div, 2); @@ -1052,11 +1048,11 @@ class CFVMFlowSolverBase : public CSolver { /*--- Add off diagonals ---*/ - StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,u,1) + Gradient_Primitive(iPoint,v,0)), 2); + StrainMag(iPoint) += 2.0*pow(0.5*(VelocityGradient(0,1) + VelocityGradient(1,0)), 2); if (nDim == 3) { - StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,u,2) + Gradient_Primitive(iPoint,w,0)), 2); - StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,v,2) + Gradient_Primitive(iPoint,w,1)), 2); + StrainMag(iPoint) += 2.0*pow(0.5*(VelocityGradient(0,2) + VelocityGradient(2,0)), 2); + StrainMag(iPoint) += 2.0*pow(0.5*(VelocityGradient(1,2) + VelocityGradient(2,1)), 2); } StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint)); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index b4dc537102d..256643fe1d1 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1053,8 +1053,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolver** solve Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, GradTangVel[MAXNDIM] = {0.0}; /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/ - su2double** Grad_Reflected = new su2double*[nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) Grad_Reflected[iVar] = new su2double[nDim]; + su2activematrix Grad_Reflected(nPrimVarGrad, nDim); /*--- Loop over all the vertices on this boundary marker. ---*/ @@ -1260,7 +1259,7 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolver** solve GradNormVel[iDim] * UnitNormal[iVar] + GradTangVel[iDim] * Tangential[iVar]; /*--- Set the primitive gradients of the boundary and reflected state. ---*/ - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), Grad_Reflected); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), CMatrixView(Grad_Reflected)); /*--- Turbulent kinetic energy. ---*/ if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) @@ -1280,9 +1279,6 @@ void CFVMFlowSolverBase::BC_Sym_Plane(CGeometry* geometry, CSolver** solve } // for iVertex END_SU2_OMP_FOR - /*--- Free locally allocated memory ---*/ - for (iVar = 0; iVar < nPrimVarGrad; iVar++) delete[] Grad_Reflected[iVar]; - delete[] Grad_Reflected; } template diff --git a/SU2_CFD/include/variables/CAdjEulerVariable.hpp b/SU2_CFD/include/variables/CAdjEulerVariable.hpp index 61192a61818..c37e74ce899 100644 --- a/SU2_CFD/include/variables/CAdjEulerVariable.hpp +++ b/SU2_CFD/include/variables/CAdjEulerVariable.hpp @@ -169,7 +169,9 @@ class CAdjEulerVariable : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { + return Gradient_Reconstruction[iPoint]; + } /*! * \brief Get the reconstruction gradient for variables at all points. diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 27a09d7615a..2c2af8edef5 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -143,7 +143,9 @@ class CEulerVariable : public CVariable { * \brief Get the value of the primitive variables gradient. * \return Value of the primitive variables gradient. */ - inline su2double **GetGradient_Primitive(unsigned long iPoint) final { return Gradient_Primitive[iPoint]; } + inline CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) final { + return Gradient_Primitive(iPoint,iVar); + } /*! * \brief Get the value of the primitive variables gradient. @@ -156,7 +158,7 @@ class CEulerVariable : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } /*! * \brief A virtual member. diff --git a/SU2_CFD/include/variables/CHeatVariable.hpp b/SU2_CFD/include/variables/CHeatVariable.hpp index 02aa1400b41..fe1c87c2a14 100644 --- a/SU2_CFD/include/variables/CHeatVariable.hpp +++ b/SU2_CFD/include/variables/CHeatVariable.hpp @@ -63,7 +63,7 @@ class CHeatVariable final : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } /*! * \brief Get the reconstruction gradient for primitive variable at all points. diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 418c73da917..eb73839648c 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -120,7 +120,9 @@ class CIncEulerVariable : public CVariable { * \param[in] iPoint - Point index. * \return Value of the primitive variables gradient. */ - inline su2double **GetGradient_Primitive(unsigned long iPoint) final { return Gradient_Primitive[iPoint]; } + inline CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) final { + return Gradient_Primitive(iPoint,iVar); + } /*! * \brief Get the value of the primitive variables gradient. @@ -134,7 +136,7 @@ class CIncEulerVariable : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } /*! * \brief Set the value of the pressure. diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index f01ff2adf18..9c24b9cfed5 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -255,7 +255,7 @@ class CNEMOEulerVariable : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } /*! * \brief Subtract value to the gradient of the primitive variables. @@ -281,7 +281,9 @@ class CNEMOEulerVariable : public CVariable { * \brief Get the value of the primitive variables gradient. * \return Value of the primitive variables gradient. */ - inline su2double **GetGradient_Primitive(unsigned long iPoint) final { return Gradient_Primitive[iPoint]; } + inline CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) final { + return Gradient_Primitive(iPoint,iVar); + } /*! * \brief Get the primitive variable gradients for all points. diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp index bb71ae81f1f..6eb70eb4171 100644 --- a/SU2_CFD/include/variables/CTurbSAVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp @@ -108,7 +108,7 @@ class CTurbSAVariable final : public CTurbVariable { * \brief Set the vortex tilting measure for computation of the EDDES length scale * \param[in] iPoint - Point index. */ - void SetVortex_Tilting(unsigned long iPoint, const su2double* const* PrimGrad_Flow, + void SetVortex_Tilting(unsigned long iPoint, CMatrixView, const su2double* Vorticity, su2double LaminarViscosity) override; /*! diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 3066299c756..010c3414881 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -77,7 +77,7 @@ class CTurbVariable : public CVariable { * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ - inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + inline CMatrixView GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } /*! * \brief Get the reconstruction gradient for primitive variable at all points. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index e83e6c1376a..35d70e3eee4 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -614,7 +614,7 @@ class CVariable { * \param[in] iPoint - Point index. * \return Value of the solution gradient. */ - inline su2double** GetAuxVarGradient(unsigned long iPoint) { + inline CMatrixView GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } @@ -703,7 +703,7 @@ class CVariable { * \param[in] iPoint - Point index. * \return Value of the gradient solution. */ - inline su2double **GetGradient(unsigned long iPoint) { return Gradient[iPoint]; } + inline CMatrixView GetGradient(unsigned long iPoint) { return Gradient[iPoint]; } /*! * \brief Get the value of the solution gradient. @@ -1644,7 +1644,7 @@ class CVariable { * \brief A virtual member. * \return Value of the primitive variables gradient. */ - inline virtual su2double **GetGradient_Primitive(unsigned long iPoint) { return nullptr; } + inline virtual CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) { return nullptr; } /*! * \brief A virtual member. @@ -1656,7 +1656,7 @@ class CVariable { * \brief Get the value of the primitive gradient for MUSCL reconstruction. * \return Value of the primitive gradient for MUSCL reconstruction. */ - inline virtual su2double **GetGradient_Reconstruction(unsigned long iPoint) { return nullptr; } + inline virtual CMatrixView GetGradient_Reconstruction(unsigned long iPoint) { return nullptr; } /*! * \brief Get the reconstruction gradient for primitive variable at all points. @@ -2186,7 +2186,7 @@ class CVariable { inline virtual su2double GetTauWall(unsigned long iPoint) const { return 0.0; } - inline virtual void SetVortex_Tilting(unsigned long iPoint, const su2double* const* PrimGrad_Flow, + inline virtual void SetVortex_Tilting(unsigned long iPoint, CMatrixView PrimGrad_Flow, const su2double* Vorticity, su2double LaminarViscosity) {} inline virtual su2double GetVortex_Tilting(unsigned long iPoint) const { return 0.0; } diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index f3ae3835121..d361e5e4972 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -641,7 +641,6 @@ void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSyst unsigned long Point_Fine, Point_Coarse; unsigned short iVar, iDim, iChildren; su2double Area_Parent, Area_Children; - const su2double* const* Gradient_fine = nullptr; const unsigned short nDim = geo_coarse->GetnDim(); const unsigned short nVar = sol_coarse->GetnVar(); @@ -661,7 +660,7 @@ void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSyst for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); Area_Children = geo_fine->nodes->GetVolume(Point_Fine); - Gradient_fine = sol_fine->GetNodes()->GetGradient(Point_Fine); + auto Gradient_fine = sol_fine->GetNodes()->GetGradient(Point_Fine); for (iVar = 0; iVar < nVar; iVar++) for (iDim = 0; iDim < nDim; iDim++) diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index 139184f58d4..14abc5e6279 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -189,7 +189,7 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); su2double tau[3][3]; - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity); + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->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/NEMO/NEMO_sources.cpp b/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp index 1b6efff9da5..a342081fac8 100644 --- a/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp +++ b/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp @@ -297,7 +297,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi unsigned short iDim, iSpecies, iVar; su2double rho, rhov, vel2, H, yinv, T, Tve, Ru, RuSI; - su2double *Ds, **GV, ktr, kve; + su2double *Ds, ktr, kve; /*--- Rename for convenience ---*/ Ds = Diffusion_Coeff_i; @@ -306,7 +306,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi rho = V_i[RHO_INDEX]; T = V_i[T_INDEX]; Tve = V_i[TVE_INDEX]; - GV = PrimVar_Grad_i; + auto GV = PrimVar_Grad_i; RuSI= UNIVERSAL_GAS_CONSTANT; Ru = 1000.0*RuSI; diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index f4bf29ed537..94a4d5dd68c 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -553,7 +553,7 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv } else { SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } - SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(&(Node_Flow->GetGradient_Primitive(iPoint)[1]))); + SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1))); } LoadCommonFVMOutputs(config, geometry, iPoint); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index d688b97f85b..e81a95a1989 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -650,7 +650,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve } else { SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } - SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(&(Node_Flow->GetGradient_Primitive(iPoint)[1]))); + SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1))); } // Streamwise Periodicity diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index d25a861c8c4..79b5eba9f19 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -816,39 +816,6 @@ void CFlowOutput::Set_CpInverseDesign(CSolver *solver, CGeometry *geometry, CCon } -su2double CFlowOutput::GetQ_Criterion(su2double** VelocityGradient) const { - - /*--- Make a 3D copy of the gradient so we do not have worry about nDim ---*/ - - su2double Grad_Vel[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}}; - - for (unsigned short iDim = 0; iDim < nDim; iDim++) - for (unsigned short jDim = 0 ; jDim < nDim; jDim++) - Grad_Vel[iDim][jDim] = VelocityGradient[iDim][jDim]; - - /*--- Q Criterion Eq 1.2 of HALLER, G. (2005). An objective definition of a vortex. - Journal of Fluid Mechanics, 525, 1-26. doi:10.1017/S0022112004002526 ---*/ - - /*--- Components of the strain rate tensor (symmetric) ---*/ - su2double s11 = Grad_Vel[0][0]; - su2double s12 = 0.5 * (Grad_Vel[0][1] + Grad_Vel[1][0]); - su2double s13 = 0.5 * (Grad_Vel[0][2] + Grad_Vel[2][0]); - su2double s22 = Grad_Vel[1][1]; - su2double s23 = 0.5 * (Grad_Vel[1][2] + Grad_Vel[2][1]); - su2double s33 = Grad_Vel[2][2]; - - /*--- Components of the spin tensor (skew-symmetric) ---*/ - su2double omega12 = 0.5 * (Grad_Vel[0][1] - Grad_Vel[1][0]); - su2double omega13 = 0.5 * (Grad_Vel[0][2] - Grad_Vel[2][0]); - su2double omega23 = 0.5 * (Grad_Vel[1][2] - Grad_Vel[2][1]); - - /*--- Q = ||Omega|| - ||Strain|| ---*/ - su2double Q = 2*(pow(omega12,2) + pow(omega13,2) + pow(omega23,2)) - - (pow(s11,2) + pow(s22,2) + pow(s33,2) + 2*(pow(s12,2) + pow(s13,2) + pow(s23,2))); - - return Q; -} - void CFlowOutput::WriteAdditionalFiles(CConfig *config, CGeometry *geometry, CSolver **solver_container){ if (config->GetFixed_CL_Mode() || config->GetFixed_CM_Mode()){ diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index e7b503dbd1d..d43c1fd83d6 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -1699,7 +1699,7 @@ void CAdjEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont CNumerics* numerics = numerics_container[CONV_TERM]; - su2double **Gradient_i, **Gradient_j, Project_Grad_i, Project_Grad_j, *Limiter_i = nullptr, + su2double Project_Grad_i, Project_Grad_j, *Limiter_i = nullptr, *Limiter_j = nullptr, *Psi_i = nullptr, *Psi_j = nullptr, *V_i, *V_j; unsigned long iEdge, iPoint, jPoint, counter_local = 0, counter_global = 0; unsigned short iDim, iVar; @@ -1746,8 +1746,8 @@ void CAdjEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Adjoint variables using gradient reconstruction and limiters ---*/ - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); + const auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint); + const auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint); if (limiter) { Limiter_i = nodes->GetLimiter(iPoint); @@ -2135,7 +2135,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver unsigned short iPos, jPos; unsigned short iDim, iMarker, iNeigh; su2double *d = nullptr, *Normal = nullptr, *Psi = nullptr, *U = nullptr, Enthalpy, conspsi = 0.0, Mach_Inf, - Area, **PrimVar_Grad = nullptr, *ConsPsi_Grad = nullptr, + Area, *ConsPsi_Grad = nullptr, ConsPsi, d_press, grad_v, v_gradconspsi, UnitNormal[3], *GridVel = nullptr, eps, 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, *USens, D[3][3], Dd[3], scale = 1.0; @@ -2245,7 +2245,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); Area = GeometryToolbox::Norm(nDim, Normal); - PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + const auto PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint)[0]; ConsPsi = nodes->GetAuxVar(iPoint); diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index f6014e666dd..2b8446543f0 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -595,14 +595,13 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con unsigned long iVertex, iPoint; unsigned short iDim, jDim, iMarker, iPos, jPos; - su2double *d = nullptr, **PsiVar_Grad = nullptr, **PrimVar_Grad = nullptr, div_phi, *Normal = nullptr, Area, + su2double *d = nullptr, div_phi, *Normal = nullptr, Area, 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], 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; su2double *USens = new su2double[nVar]; su2double *UnitNormal = new su2double[nDim]; @@ -694,8 +693,8 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con if (geometry->nodes->GetDomain(iPoint)) { - PsiVar_Grad = nodes->GetGradient(iPoint); - PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + const auto PsiVar_Grad = nodes->GetGradient(iPoint); + const auto PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); Laminar_Viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); @@ -782,7 +781,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/ - GridVel_Grad = geometry->nodes->GetGridVel_Grad(iPoint); + const auto GridVel_Grad = geometry->nodes->GetGridVel_Grad(iPoint); for (iDim = 0; iDim < nDim; iDim++) { normal_grad_gridvel[iDim] = 0.0; for (jDim = 0; jDim < nDim; jDim++) @@ -1185,7 +1184,7 @@ void CAdjNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_contai unsigned short iDim, iVar, jVar, jDim; unsigned long iVertex, iPoint, total_index, Point_Normal; - su2double *d, l1psi, vartheta, Sigma_5, **PsiVar_Grad, phi[3] = {}; + su2double *d, l1psi, vartheta, Sigma_5, phi[3] = {}; su2double sq_vel, ProjGridVel, Enthalpy = 0.0, *GridVel; su2double ViscDens, XiDens, Density, Pressure = 0.0, dPhiE_dn; su2double Laminar_Viscosity = 0.0, Eddy_Viscosity = 0.0; @@ -1332,7 +1331,7 @@ void CAdjNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_contai /*--- Store the adjoint velocity and energy gradients for clarity ---*/ - PsiVar_Grad = nodes->GetGradient(iPoint); + const auto PsiVar_Grad = nodes->GetGradient(iPoint); for (iDim = 0; iDim < nDim; iDim++) { GradPsiE[iDim] = PsiVar_Grad[nVar-1][iDim]; for (jDim = 0; jDim < nDim; jDim++) @@ -1547,7 +1546,7 @@ void CAdjNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_cont unsigned long iVertex, iPoint, total_index; unsigned short iDim, iVar, jVar, jDim; su2double *d, q, *U, dVisc_T, rho, pressure, div_phi, - force_stress, Sigma_5, **PsiVar_Grad, phi[3] = {0.0,0.0,0.0}; + force_stress, Sigma_5, phi[3] = {0.0,0.0,0.0}; su2double phis1, phis2, sq_vel, ProjVel, Enthalpy, *GridVel, phi_u, d_n; su2double Energy, ViscDens, XiDens, Density, SoundSpeed, Pressure, dPhiE_dn, Laminar_Viscosity, Eddy_Viscosity, Sigma_xx, Sigma_yy, Sigma_zz, Sigma_xy, Sigma_xz, Sigma_yz, @@ -1699,7 +1698,7 @@ void CAdjNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_cont /*--- Additional contributions to adjoint density (weak imposition) ---*/ /*--- Acquire gradient information ---*/ - PsiVar_Grad = nodes->GetGradient(iPoint); + auto PsiVar_Grad = nodes->GetGradient(iPoint); GradP = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)[nVar-1]; GradDens = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)[nVar]; diff --git a/SU2_CFD/src/solvers/CAdjTurbSolver.cpp b/SU2_CFD/src/solvers/CAdjTurbSolver.cpp index 29afe5c2ab1..94b8a4e4ee9 100644 --- a/SU2_CFD/src/solvers/CAdjTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjTurbSolver.cpp @@ -288,7 +288,7 @@ void CAdjTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_conta CNumerics* numerics = numerics_container[CONV_TERM]; unsigned long iEdge, iPoint, jPoint; - su2double *U_i, *U_j, *TurbPsi_i, *TurbPsi_j, **TurbVar_Grad_i, **TurbVar_Grad_j; + su2double *U_i, *U_j, *TurbPsi_i, *TurbPsi_j; for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { @@ -310,8 +310,8 @@ void CAdjTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetTurbAdjointVar(TurbPsi_i, TurbPsi_j); /*--- Gradient of turbulent variables w/o reconstruction ---*/ - TurbVar_Grad_i = solver_container[TURB_SOL]->GetNodes()->GetGradient(iPoint); - TurbVar_Grad_j = solver_container[TURB_SOL]->GetNodes()->GetGradient(jPoint); + const auto TurbVar_Grad_i = solver_container[TURB_SOL]->GetNodes()->GetGradient(iPoint); + const auto TurbVar_Grad_j = solver_container[TURB_SOL]->GetNodes()->GetGradient(jPoint); numerics->SetTurbVarGradient(TurbVar_Grad_i, TurbVar_Grad_j); /*--- Set normal vectors and length ---*/ @@ -387,8 +387,8 @@ void CAdjTurbSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta //CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM]; unsigned long iPoint; - su2double *U_i, **GradPrimVar_i, *TurbVar_i; - su2double **TurbVar_Grad_i, *TurbPsi_i, **PsiVar_Grad_i; // Gradients + su2double *U_i, *TurbVar_i; + su2double *TurbPsi_i; // Gradients /*--- Piecewise source term ---*/ for (iPoint = 0; iPoint < nPointDomain; iPoint++) { @@ -398,7 +398,7 @@ void CAdjTurbSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetConservative(U_i, nullptr); /*--- Gradient of primitive variables w/o reconstruction ---*/ - GradPrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + auto GradPrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); numerics->SetPrimVarGradient(GradPrimVar_i, nullptr); /*--- Laminar viscosity of the fluid ---*/ @@ -409,7 +409,7 @@ void CAdjTurbSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta numerics->SetTurbVar(TurbVar_i, nullptr); /*--- Gradient of Turbulent Variables w/o reconstruction ---*/ - TurbVar_Grad_i = solver_container[TURB_SOL]->GetNodes()->GetGradient(iPoint); + auto TurbVar_Grad_i = solver_container[TURB_SOL]->GetNodes()->GetGradient(iPoint); numerics->SetTurbVarGradient(TurbVar_Grad_i, nullptr); /*--- Turbulent adjoint variables w/o reconstruction ---*/ @@ -418,7 +418,7 @@ void CAdjTurbSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Gradient of Adjoint flow variables w/o reconstruction (for non-conservative terms depending on gradients of flow adjoint vars.) ---*/ - PsiVar_Grad_i = solver_container[ADJFLOW_SOL]->GetNodes()->GetGradient(iPoint); + auto PsiVar_Grad_i = solver_container[ADJFLOW_SOL]->GetNodes()->GetGradient(iPoint); numerics->SetAdjointVarGradient(PsiVar_Grad_i, nullptr); /*--- Set volume and distances to the surface ---*/ diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index e31233ae6e4..a9579695bf0 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -352,8 +352,8 @@ void CHeatSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe CNumerics* numerics = numerics_container[CONV_TERM]; - su2double *V_i, *V_j, Temp_i, Temp_i_Corrected, Temp_j, Temp_j_Corrected, **Gradient_i, **Gradient_j, Project_Grad_i, Project_Grad_j, - **Temp_i_Grad, **Temp_j_Grad, Project_Temp_i_Grad, Project_Temp_j_Grad; + su2double *V_i, *V_j, Temp_i, Temp_i_Corrected, Temp_j, Temp_j_Corrected, + Project_Grad_i, Project_Grad_j, Project_Temp_i_Grad, Project_Temp_j_Grad; su2double Vector_i[MAXNDIM] = {0}; su2double Vector_j[MAXNDIM] = {0}; @@ -377,8 +377,8 @@ void CHeatSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe V_i = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); V_j = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(jPoint); - Temp_i_Grad = nodes->GetGradient(iPoint); - Temp_j_Grad = nodes->GetGradient(jPoint); + const auto Temp_i_Grad = nodes->GetGradient(iPoint); + const auto Temp_j_Grad = nodes->GetGradient(jPoint); numerics->SetConsVarGradient(Temp_i_Grad, Temp_j_Grad); Temp_i = nodes->GetTemperature(iPoint); @@ -392,10 +392,10 @@ void CHeatSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe Vector_j[iDim] = 0.5*(geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim)); } - Gradient_i = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Reconstruction(iPoint); - Gradient_j = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Reconstruction(jPoint); - Temp_i_Grad = nodes->GetGradient_Reconstruction(iPoint); - Temp_j_Grad = nodes->GetGradient_Reconstruction(jPoint); + const auto Gradient_i = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Reconstruction(iPoint); + const auto Gradient_j = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Reconstruction(jPoint); + const auto Temp_i_Grad = nodes->GetGradient_Reconstruction(iPoint); + const auto Temp_j_Grad = nodes->GetGradient_Reconstruction(jPoint); /*Loop to correct the flow variables*/ for (auto iVar = 0u; iVar < nVarFlow; iVar++) { @@ -450,7 +450,7 @@ void CHeatSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain CNumerics* numerics = numerics_container[VISC_TERM]; su2double laminar_viscosity, Prandtl_Lam, Prandtl_Turb, eddy_viscosity_i, eddy_viscosity_j, - thermal_diffusivity_i, thermal_diffusivity_j, Temp_i, Temp_j, **Temp_i_Grad, **Temp_j_Grad; + thermal_diffusivity_i, thermal_diffusivity_j, Temp_i, Temp_j; const bool turb = ((config->GetKind_Solver() == INC_RANS) || (config->GetKind_Solver() == DISC_ADJ_INC_RANS)); @@ -471,8 +471,8 @@ void CHeatSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain geometry->nodes->GetCoord(jPoint)); numerics->SetNormal(geometry->edges->GetNormal(iEdge)); - Temp_i_Grad = nodes->GetGradient(iPoint); - Temp_j_Grad = nodes->GetGradient(jPoint); + const auto Temp_i_Grad = nodes->GetGradient(iPoint); + const auto Temp_j_Grad = nodes->GetGradient(jPoint); numerics->SetConsVarGradient(Temp_i_Grad, Temp_j_Grad); /*--- Primitive variables w/o reconstruction ---*/ diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 3a10d22f102..b4c9ff25aef 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -290,8 +290,7 @@ void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, unsigned short T_INDEX, TVE_INDEX, RHOCVTR_INDEX; unsigned long iVertex, iPoint, total_index; su2double dTdn, dTvedn, ktr, kve, pcontrol, Wall_HeatFlux; - su2double *Normal, Area; - su2double **GradV, *V; + su2double *Normal, Area, *V; /*--- Assign booleans ---*/ implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); @@ -330,7 +329,7 @@ void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, // Note: Contributions from qtr and qve are used for proportional control // to drive the solution toward the specified heatflux more quickly. V = nodes->GetPrimitive(iPoint); - GradV = nodes->GetGradient_Primitive(iPoint); + const auto GradV = nodes->GetGradient_Primitive(iPoint); dTdn = 0.0; dTvedn = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -442,7 +441,7 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, su2double rho, Ys; su2double *Normal, Area; su2double *Ds, *V, *dYdn, SdYdn; - su2double **GradV, **GradY; + su2double **GradY; /*--- Assign booleans ---*/ implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); @@ -498,7 +497,7 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, /*--- Get temperature gradient information ---*/ V = nodes->GetPrimitive(iPoint); - GradV = nodes->GetGradient_Primitive(iPoint); + const auto GradV = nodes->GetGradient_Primitive(iPoint); dTdn = 0.0; dTvedn = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -964,7 +963,6 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, su2double Viscosity, Eddy_Visc, Lambda; su2double Density, GasConstant; - const su2double* const* Grad_PrimVar; su2double Vector_Tangent_dT[MAXNDIM] = {0.0}, Vector_Tangent_dTve[MAXNDIM] = {0.0}, Vector_Tangent_HF[MAXNDIM] = {0.0}; su2double dTn, dTven; su2double rhoCvtr, rhoCvve; @@ -1058,7 +1056,7 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, kve = kve*(1.0+scl); /*--- Retrieve Primitive Gradients ---*/ - Grad_PrimVar = nodes->GetGradient_Primitive(iPoint); + const auto Grad_PrimVar = nodes->GetGradient_Primitive(iPoint); /*--- Calculate specific gas constant --- */ //TODO: Move to fluidmodel? diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 21d307bdac0..0cef0e5ecb9 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -1672,7 +1672,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC su2double cb1 = 0.1355, ct3 = 1.2, ct4 = 0.5; su2double sigma = 2./3., cb2 = 0.622, f_max=1.0, f_min=0.1, a1=0.15, a2=0.3; su2double cw1 = 0.0, Ji = 0.0, Ji_2 = 0.0, Ji_3 = 0.0, fv1 = 0.0, fv2 = 0.0, ft2 = 0.0, psi_2 = 0.0; - const su2double *coord_i = nullptr, *coord_j = nullptr, *const *primVarGrad = nullptr, *vorticity = nullptr; + const su2double *coord_i = nullptr, *coord_j = nullptr, *vorticity = nullptr; su2double delta[3] = {0.0}, ratioOmega[3] = {0.0}, vortexTiltingMeasure = 0.0; SU2_OMP_FOR_DYN(omp_chunk_size) @@ -1681,7 +1681,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC coord_i = geometry->nodes->GetCoord(iPoint); nNeigh = geometry->nodes->GetnPoint(iPoint); wallDistance = geometry->nodes->GetWall_Distance(iPoint); - primVarGrad = solver[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + const auto primVarGrad = solver[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); vorticity = solver[FLOW_SOL]->GetNodes()->GetVorticity(iPoint); density = solver[FLOW_SOL]->GetNodes()->GetDensity(iPoint); laminarViscosity = solver[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp index 11b5faa92af..84983cf72ba 100644 --- a/SU2_CFD/src/variables/CTurbSAVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp @@ -51,7 +51,7 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi Vortex_Tilting.resize(nPoint); } -void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, const su2double* const* PrimGrad_Flow, +void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView PrimGrad_Flow, const su2double* Vorticity, su2double LaminarViscosity) { su2double Strain[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}, Omega, StrainDotVort[3], numVecVort[3]; From 6971ab0e4c0b234239424b680c05dc07a3d1dfba Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sun, 27 Jun 2021 13:34:16 +0100 Subject: [PATCH 2/8] small fix --- Common/include/containers/container_decorators.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Common/include/containers/container_decorators.hpp b/Common/include/containers/container_decorators.hpp index 8e4f72a2c1f..2d8c5297a1b 100644 --- a/Common/include/containers/container_decorators.hpp +++ b/Common/include/containers/container_decorators.hpp @@ -50,11 +50,11 @@ class CMatrixView { template CMatrixView(const CMatrixView& other) : m_ptr(other.m_ptr), m_cols(other.m_cols) {} - explicit CMatrixView(const su2matrix& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {} - - template::value> = 0> explicit CMatrixView(su2matrix& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {} + template::value> = 0> + explicit CMatrixView(const su2matrix& mat) : m_ptr(mat.data()), m_cols(mat.cols()) {} + const Scalar* operator[] (Index i) const noexcept { return &m_ptr[i*m_cols]; } const Scalar& operator() (Index i, Index j) const noexcept { return m_ptr[i*m_cols + j]; } From 31e94a09910163801be6fa77652be7ad95a854b8 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 27 Jun 2021 18:00:17 +0100 Subject: [PATCH 3/8] const correctness and less duplication in CNumerics --- SU2_CFD/include/numerics/CNumerics.hpp | 161 ++++++++++-------- SU2_CFD/include/numerics/heat.hpp | 56 +----- SU2_CFD/include/numerics/radiation.hpp | 14 +- .../numerics/turbulent/turb_diffusion.hpp | 3 - SU2_CFD/src/drivers/CDriver.cpp | 4 +- SU2_CFD/src/numerics/CNumerics.cpp | 77 +-------- SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp | 6 +- SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp | 12 +- SU2_CFD/src/numerics/NEMO/convection/msw.cpp | 5 +- SU2_CFD/src/numerics/NEMO/convection/roe.cpp | 3 +- SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 4 + SU2_CFD/src/numerics/flow/flow_sources.cpp | 4 +- SU2_CFD/src/numerics/heat.cpp | 132 ++------------ SU2_CFD/src/numerics/radiation.cpp | 53 +----- .../src/numerics/turbulent/turb_diffusion.cpp | 38 +---- 15 files changed, 140 insertions(+), 432 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 2b8f8b0a62f..35822fbde2a 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -53,22 +53,18 @@ class CNumerics { su2double Gamma_Minus_One; /*!< \brief Fluids's Gamma - 1.0 . */ su2double Minf; /*!< \brief Free stream Mach number . */ su2double Gas_Constant; /*!< \brief Gas constant. */ - su2double *Vector; /*!< \brief Auxiliary vector. */ su2double Prandtl_Lam; /*!< \brief Laminar Prandtl's number. */ su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl's number. */ su2double - **Flux_Tensor, /*!< \brief Flux tensor (used for viscous and inviscid purposes. */ *Proj_Flux_Tensor; /*!< \brief Flux tensor projected in a direction. */ su2double **tau; /*!< \brief Viscous stress tensor. */ const su2double delta [3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; /*!< \brief Identity matrix. */ - su2double + const su2double *Diffusion_Coeff_i, /*!< \brief Species diffusion coefficients at point i. */ *Diffusion_Coeff_j; /*!< \brief Species diffusion coefficients at point j. */ su2double Laminar_Viscosity_i, /*!< \brief Laminar viscosity at point i. */ - Laminar_Viscosity_j, /*!< \brief Laminar viscosity at point j. */ - Laminar_Viscosity_id, /*!< \brief Variation of laminar viscosity at point i. */ - Laminar_Viscosity_jd; /*!< \brief Variation of laminar viscosity at point j. */ + Laminar_Viscosity_j; /*!< \brief Laminar viscosity at point j. */ su2double Thermal_Conductivity_i, /*!< \brief Thermal conductivity at point i. */ Thermal_Conductivity_j, /*!< \brief Thermal conductivity at point j. */ @@ -77,9 +73,6 @@ class CNumerics { Thermal_Diffusivity_i, /*!< \brief Thermal diffusivity at point i. */ Thermal_Diffusivity_j; /*!< \brief Thermal diffusivity at point j. */ su2double - SpecificHeat_i, /*!< \brief Specific heat c_p at point j. */ - SpecificHeat_j; /*!< \brief Specific heat c_p at point j. */ - su2double Cp_i, /*!< \brief Cp at point i. */ Cp_j; /*!< \brief Cp at point j. */ su2double @@ -92,9 +85,6 @@ class CNumerics { Pressure_i, /*!< \brief Pressure at point i. */ Pressure_j; /*!< \brief Pressure at point j. */ su2double - GravityForce_i, /*!< \brief Gravity force at point i. */ - GravityForce_j; /*!< \brief Gravity force at point j. */ - su2double Density_i, /*!< \brief Density at point i. */ Density_j; /*!< \brief Density at point j. */ su2double @@ -107,9 +97,6 @@ class CNumerics { Lambda_i, /*!< \brief Spectral radius at point i. */ Lambda_j; /*!< \brief Spectral radius at point j. */ su2double - LambdaComb_i, /*!< \brief Spectral radius at point i. */ - LambdaComb_j; /*!< \brief Spectral radius at point j. */ - su2double SoundSpeed_i, /*!< \brief Sound speed at point i. */ SoundSpeed_j; /*!< \brief Sound speed at point j. */ su2double @@ -121,41 +108,34 @@ class CNumerics { su2double Temp_i, /*!< \brief Temperature at point i. */ Temp_j; /*!< \brief Temperature at point j. */ - su2double + const su2double *Und_Lapl_i, /*!< \brief Undivided laplacians at point i. */ *Und_Lapl_j; /*!< \brief Undivided laplacians at point j. */ su2double Sensor_i, /*!< \brief Pressure sensor at point i. */ Sensor_j; /*!< \brief Pressure sensor at point j. */ - su2double + const su2double *GridVel_i, /*!< \brief Grid velocity at point i. */ *GridVel_j; /*!< \brief Grid velocity at point j. */ - su2double + const su2double *U_i, /*!< \brief Vector of conservative variables at point i. */ - *U_id, /*!< \brief Vector of derivative of conservative variables at point i. */ - *U_j, /*!< \brief Vector of conservative variables at point j. */ - *U_jd; /*!< \brief Vector of derivative of conservative variables at point j. */ - su2double + *U_j; /*!< \brief Vector of conservative variables at point j. */ + const su2double *V_i, /*!< \brief Vector of primitive variables at point i. */ *V_j; /*!< \brief Vector of primitive variables at point j. */ - su2double + const su2double *S_i, /*!< \brief Vector of secondary variables at point i. */ *S_j; /*!< \brief Vector of secondary variables at point j. */ - su2double + const su2double *Psi_i, /*!< \brief Vector of adjoint variables at point i. */ *Psi_j; /*!< \brief Vector of adjoint variables at point j. */ - su2double - *DeltaU_i, /*!< \brief Vector of linearized variables at point i. */ - *DeltaU_j; /*!< \brief Vector of linearized variables at point j. */ - su2double + const su2double *TurbVar_i, /*!< \brief Vector of turbulent variables at point i. */ - *TurbVar_id, /*!< \brief Vector of derivative of turbulent variables at point i. */ - *TurbVar_j, /*!< \brief Vector of turbulent variables at point j. */ - *TurbVar_jd; /*!< \brief Vector of derivative of turbulent variables at point j. */ - su2double + *TurbVar_j; /*!< \brief Vector of turbulent variables at point j. */ + const su2double *TransVar_i, /*!< \brief Vector of turbulent variables at point i. */ *TransVar_j; /*!< \brief Vector of turbulent variables at point j. */ - su2double + const su2double *TurbPsi_i, /*!< \brief Vector of adjoint turbulent variables at point i. */ *TurbPsi_j; /*!< \brief Vector of adjoint turbulent variables at point j. */ CMatrixView @@ -175,37 +155,31 @@ class CNumerics { AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ - su2double + const su2double *Coord_i, /*!< \brief Cartesians coordinates of point i. */ *Coord_j; /*!< \brief Cartesians coordinates of point j. */ unsigned short Neighbor_i, /*!< \brief Number of neighbors of the point i. */ Neighbor_j; /*!< \brief Number of neighbors of the point j. */ - const su2double - *Normal; /*!< \brief Normal vector, its norm is the area of the face. */ - su2double - *UnitNormal, /*!< \brief Unitary normal vector. */ - *UnitNormald; /*!< \brief Derivative of unitary normal vector. */ + const su2double *Normal = nullptr; /*!< \brief Normal vector, its norm is the area of the face. */ + su2double UnitNormal[MAXNDIM] = {0.0}; /*!< \brief Unitary normal vector. */ su2double TimeStep, /*!< \brief Time step useful in dual time method. */ Area, /*!< \brief Area of the face i-j. */ Volume; /*!< \brief Volume of the control volume around point i. */ su2double vel2_inf; /*!< \brief value of the square of freestream speed. */ - su2double + const su2double *WindGust_i, /*!< \brief Wind gust at point i. */ - *WindGust_j; /*!< \brief Wind gust at point j. */ - su2double + *WindGust_j, /*!< \brief Wind gust at point j. */ *WindGustDer_i, /*!< \brief Wind gust derivatives at point i. */ *WindGustDer_j; /*!< \brief Wind gust derivatives at point j. */ - su2double *Vorticity_i, *Vorticity_j; /*!< \brief Vorticity. */ + const su2double *Vorticity_i, *Vorticity_j; /*!< \brief Vorticity. */ su2double StrainMag_i, StrainMag_j; /*!< \brief Strain rate magnitude. */ su2double Dissipation_i, Dissipation_j; /*!< \brief Dissipation. */ su2double Dissipation_ij; su2double roughness_i = 0.0, /*!< \brief Roughness of the wall nearest to point i. */ roughness_j = 0.0; /*!< \brief Roughness of the wall nearest to point j. */ - su2double *l, *m; - su2double MeanPerturbedRSM[3][3];/*!< \brief Perturbed Reynolds stress tensor */ bool using_uq; /*!< \brief Flag for UQ methodology */ su2double PerturbedStrainMag; /*!< \brief Strain magnitude calculated using perturbed stress tensor */ @@ -278,7 +252,7 @@ class CNumerics { * \brief Set the value of the vorticity * \param[in] val_vorticity - Value of the vorticity. */ - void SetVorticity(su2double *val_vorticity_i, su2double *val_vorticity_j) { + void SetVorticity(const su2double *val_vorticity_i, const su2double *val_vorticity_j) { Vorticity_i = val_vorticity_i; Vorticity_j = val_vorticity_j; } @@ -298,7 +272,7 @@ class CNumerics { * \param[in] val_u_i - Value of the conservative variable at point i. * \param[in] val_u_j - Value of the conservative variable at point j. */ - inline void SetConservative(su2double *val_u_i, su2double *val_u_j) { + inline void SetConservative(const su2double *val_u_i, const su2double *val_u_j) { U_i = val_u_i; U_j = val_u_j; } @@ -308,7 +282,7 @@ class CNumerics { * \param[in] val_v_i - Value of the primitive variable at point i. * \param[in] val_v_j - Value of the primitive variable at point j. */ - inline void SetPrimitive(su2double *val_v_i, su2double *val_v_j) { + inline void SetPrimitive(const su2double *val_v_i, const su2double *val_v_j) { V_i = val_v_i; V_j = val_v_j; } @@ -318,7 +292,7 @@ class CNumerics { * \param[in] val_v_i - Value of the primitive variable at point i. * \param[in] val_v_j - Value of the primitive variable at point j. */ - inline void SetSecondary(su2double *val_s_i, su2double *val_s_j) { + inline void SetSecondary(const su2double *val_s_i, const su2double *val_s_j) { S_i = val_s_i; S_j = val_s_j; } @@ -356,7 +330,7 @@ class CNumerics { * \param[in] val_psi_i - Value of the adjoint variable at point i. * \param[in] val_psi_j - Value of the adjoint variable at point j. */ - inline void SetAdjointVar(su2double *val_psi_i, su2double *val_psi_j) { + inline void SetAdjointVar(const su2double *val_psi_i, const su2double *val_psi_j) { Psi_i = val_psi_i; Psi_j = val_psi_j; } @@ -377,7 +351,7 @@ class CNumerics { * \param[in] val_turbvar_i - Value of the turbulent variable at point i. * \param[in] val_turbvar_j - Value of the turbulent variable at point j. */ - inline void SetTurbVar(su2double *val_turbvar_i, su2double *val_turbvar_j) { + inline void SetTurbVar(const su2double *val_turbvar_i, const su2double *val_turbvar_j) { TurbVar_i = val_turbvar_i; TurbVar_j = val_turbvar_j; } @@ -387,7 +361,7 @@ class CNumerics { * \param[in] val_transvar_i - Value of the turbulent variable at point i. * \param[in] val_transvar_j - Value of the turbulent variable at point j. */ - inline void SetTransVar(su2double *val_transvar_i, su2double *val_transvar_j) { + inline void SetTransVar(const su2double *val_transvar_i, const su2double *val_transvar_j) { TransVar_i = val_transvar_i; TransVar_j = val_transvar_j; } @@ -419,7 +393,7 @@ class CNumerics { * \param[in] val_turbpsivar_i - Value of the adjoint turbulent variable at point i. * \param[in] val_turbpsivar_j - Value of the adjoint turbulent variable at point j. */ - inline void SetTurbAdjointVar(su2double *val_turbpsivar_i, su2double *val_turbpsivar_j) { + inline void SetTurbAdjointVar(const su2double *val_turbpsivar_i, const su2double *val_turbpsivar_j) { TurbPsi_i = val_turbpsivar_i; TurbPsi_j = val_turbpsivar_j; } @@ -653,6 +627,58 @@ class CNumerics { } } + /*! + * \brief Perturb the Reynolds stress tensor based on parameters. + * \param[in] nDim - Dimension of the flow problem, 2 or 3. + * \param[in] nVar - Number of variables. + * \param[in] normal - Area vector. + * \param[in] coord_i - Coordinate of point i. + * \param[in] coord_j - Coordinate of point j. + * \param[in] grad_i - Gradient at point i. + * \param[in] grad_j - Gradient at point j. + * \param[in] correct - Correct + * \param[in] var_i - Variable at point i. + * \param[in] var_j - Variable at point j. + * \param[out] projNormal - Average gradient projected on normal. + * \param[out] projCorrected - Corrected gradient projected on normal. + * \return (Edge_Vector DOT normal) / |Edge_Vector|^2. + */ + template + FORCEINLINE static su2double ComputeProjectedGradient(int nDim, int nVar, const Vec1& normal, + const Vec1& coord_i, const Vec1& coord_j, + const Mat& grad_i, const Mat& grad_j, + bool correct, + const Vec2& var_i, const Vec2& var_j, + su2double* projNormal, + su2double* projCorrected) { + nDim = (nDim > 2)? 3 : 2; + su2double edgeVec[MAXNDIM], dist_ij_2 = 0.0, proj_vector_ij = 0.0; + + for (int iDim = 0; iDim < nDim; iDim++) { + edgeVec[iDim] = coord_j[iDim] - coord_i[iDim]; + dist_ij_2 += pow(edgeVec[iDim], 2); + proj_vector_ij += edgeVec[iDim] * normal[iDim]; + } + proj_vector_ij /= max(dist_ij_2,EPS); + + /*--- Mean gradient approximation. ---*/ + for (int iVar = 0; iVar < nVar; iVar++) { + projNormal[iVar] = 0.0; + su2double edgeProj = 0.0; + + for (int iDim = 0; iDim < nDim; iDim++) { + su2double meanGrad = 0.5 * (grad_i[iVar][iDim] + grad_j[iVar][iDim]); + projNormal[iVar] += meanGrad * normal[iDim]; + if (correct) edgeProj += meanGrad * edgeVec[iDim]; + } + + projCorrected[iVar] = projNormal[iVar]; + if (correct) projCorrected[iVar] -= (edgeProj - (var_j[iVar]-var_i[iVar])) * proj_vector_ij; + } + + return proj_vector_ij; + } + /*! * \brief Set the value of the first blending function. * \param[in] val_F1_i - Value of the first Menter blending function at point i. @@ -690,8 +716,8 @@ class CNumerics { * \param[in] val_diffusioncoeff_i - Value of the diffusion coefficients at i. * \param[in] val_diffusioncoeff_j - Value of the diffusion coefficients at j */ - inline void SetDiffusionCoeff(su2double* val_diffusioncoeff_i, - su2double* val_diffusioncoeff_j) { + inline void SetDiffusionCoeff(const su2double* val_diffusioncoeff_i, + const su2double* val_diffusioncoeff_j) { Diffusion_Coeff_i = val_diffusioncoeff_i; Diffusion_Coeff_j = val_diffusioncoeff_j; } @@ -750,8 +776,8 @@ class CNumerics { */ inline void SetSpecificHeat(su2double val_specific_heat_i, su2double val_specific_heat_j) { - SpecificHeat_i = val_specific_heat_i; - SpecificHeat_j = val_specific_heat_j; + Cp_i = val_specific_heat_i; + Cp_j = val_specific_heat_j; } /*! @@ -800,7 +826,7 @@ class CNumerics { * \param[in] val_coord_i - Coordinates of the point i. * \param[in] val_coord_j - Coordinates of the point j. */ - inline void SetCoord(su2double *val_coord_i, su2double *val_coord_j) { + inline void SetCoord(const su2double *val_coord_i, const su2double *val_coord_j) { Coord_i = val_coord_i; Coord_j = val_coord_j; } @@ -810,7 +836,7 @@ class CNumerics { * \param[in] val_gridvel_i - Grid velocity of the point i. * \param[in] val_gridvel_j - Grid velocity of the point j. */ - inline void SetGridVel(su2double *val_gridvel_i, su2double *val_gridvel_j) { + inline void SetGridVel(const su2double *val_gridvel_i, const su2double *val_gridvel_j) { GridVel_i = val_gridvel_i; GridVel_j = val_gridvel_j; } @@ -820,7 +846,7 @@ class CNumerics { * \param[in] val_windgust_i - Wind gust of the point i. * \param[in] val_windgust_j - Wind gust of the point j. */ - inline void SetWindGust(su2double *val_windgust_i, su2double *val_windgust_j) { + inline void SetWindGust(const su2double *val_windgust_i, const su2double *val_windgust_j) { WindGust_i = val_windgust_i; WindGust_j = val_windgust_j; } @@ -830,7 +856,7 @@ class CNumerics { * \param[in] val_windgust_i - Wind gust derivatives of the point i. * \param[in] val_windgust_j - Wind gust derivatives of the point j. */ - inline void SetWindGustDer(su2double *val_windgustder_i, su2double *val_windgustder_j) { + inline void SetWindGustDer(const su2double *val_windgustder_i, const su2double *val_windgustder_j) { WindGustDer_i = val_windgustder_i; WindGustDer_j = val_windgustder_j; } @@ -910,7 +936,7 @@ class CNumerics { * \param[in] val_und_lapl_i Undivided laplacian at point i. * \param[in] val_und_lapl_j Undivided laplacian at point j. */ - inline void SetUndivided_Laplacian(su2double *val_und_lapl_i, su2double *val_und_lapl_j) { + inline void SetUndivided_Laplacian(const su2double *val_und_lapl_i, const su2double *val_und_lapl_j) { Und_Lapl_i = val_und_lapl_i; Und_Lapl_j = val_und_lapl_j; } @@ -962,15 +988,6 @@ class CNumerics { */ inline su2double GetDissipation() const { return Dissipation_ij; } - /*! - * \brief Get the inviscid fluxes. - * \param[in] val_density - Value of the density. - * \param[in] val_velocity - Value of the velocity. - * \param[in] val_pressure - Value of the pressure. - * \param[in] val_enthalpy - Value of the enthalpy. - */ - void GetInviscidFlux(su2double val_density, const su2double *val_velocity, su2double val_pressure, su2double val_enthalpy); - /*! * \brief Compute the projected inviscid flux vector. * \param[in] val_density - Pointer to the density. @@ -1504,7 +1521,7 @@ class CNumerics { * \brief Computes a basis of orthogonal vectors from a supplied vector * \param[in] config - Normal vector */ - void CreateBasis(const su2double *val_Normal); + void CreateBasis(const su2double *val_Normal, su2double* l, su2double* m); /*! * \brief Set the value of the Tauwall diff --git a/SU2_CFD/include/numerics/heat.hpp b/SU2_CFD/include/numerics/heat.hpp index 61d81227e6c..ae7fb9ce1f9 100644 --- a/SU2_CFD/include/numerics/heat.hpp +++ b/SU2_CFD/include/numerics/heat.hpp @@ -123,53 +123,7 @@ class CUpwSca_Heat : public CNumerics { */ class CAvgGrad_Heat : public CNumerics { private: - su2double **Mean_GradHeatVar; - su2double *Proj_Mean_GradHeatVar_Normal, *Proj_Mean_GradHeatVar_Corrected; - su2double *Edge_Vector; - bool implicit; - su2double dist_ij_2, proj_vector_ij, Thermal_Diffusivity_Mean; - unsigned short iVar, iDim; - -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGrad_Heat(void) override; - - /*! - * \brief Compute the viscous heat residual using an average of gradients with correction. - * \param[out] val_residual - Pointer to the total residual. - * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[out] Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - * \param[in] config - Definition of the particular problem. - */ - void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) override; -}; - -/*! - * \class CAvgGradCorrected_Heat - * \brief Class for computing viscous term using average of gradients with correction (heat equation). - * \ingroup ViscDiscr - * \author O. Burghardt. - * \version 7.1.1 "Blackbird" - */ -class CAvgGradCorrected_Heat : public CNumerics { -private: - su2double **Mean_GradHeatVar; - su2double *Proj_Mean_GradHeatVar_Kappa, *Proj_Mean_GradHeatVar_Edge, *Proj_Mean_GradHeatVar_Corrected; - su2double *Edge_Vector; - bool implicit; - su2double dist_ij_2, proj_vector_ij, Thermal_Diffusivity_Mean; - unsigned short iVar, iDim; + bool implicit, correct; public: /*! @@ -177,13 +131,9 @@ class CAvgGradCorrected_Heat : public CNumerics { * \param[in] val_nDim - Number of dimensions of the problem. * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] correct - Correct the gradient. */ - CAvgGradCorrected_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGradCorrected_Heat(void) override; + CAvgGrad_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config, bool correct); /*! * \brief Compute the viscous heat residual using an average of gradients with correction. diff --git a/SU2_CFD/include/numerics/radiation.hpp b/SU2_CFD/include/numerics/radiation.hpp index bbff8fc9191..22a6f9bbd02 100644 --- a/SU2_CFD/include/numerics/radiation.hpp +++ b/SU2_CFD/include/numerics/radiation.hpp @@ -100,14 +100,7 @@ class CSourceP1 final : public CNumericsRadiation { class CAvgGradCorrected_P1 final : public CNumericsRadiation { private: - su2double **Mean_GradP1Var; /*!< \brief Average of gradients at cell face */ - su2double *Edge_Vector, /*!< \brief Vector from node i to node j. */ - *Proj_Mean_GradP1Var_Edge, /*!< \brief Mean_gradTurbVar DOT normal */ - *Proj_Mean_GradP1Var_Kappa, - *Proj_Mean_GradP1Var_Corrected; /*!< \brief Mean_gradTurbVar DOT normal, corrected if required*/ - su2double dist_ij, /*!< \brief Length of the edge and face. */ - proj_vector_ij; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - su2double GammaP1; /*!< \brief P1 parameter */ + su2double GammaP1; /*!< \brief P1 parameter */ public: /*! @@ -118,11 +111,6 @@ class CAvgGradCorrected_P1 final : public CNumericsRadiation { */ CAvgGradCorrected_P1(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config); - /*! - * \brief Destructor of the class. - */ - ~CAvgGradCorrected_P1(void) override; - /*! * \brief Compute the viscous residual of the P1 equation. * \param[out] val_residual - Pointer to the total residual. diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index 090e9212a91..fb07044dc04 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -47,11 +47,8 @@ class CAvgGrad_Scalar : public CNumerics { protected: su2double - Edge_Vector[MAXNDIM] = {0.0}, /*!< \brief Vector from node i to node j. */ *Proj_Mean_GradTurbVar_Normal = nullptr, /*!< \brief Mean_gradTurbVar DOT normal. */ - *Proj_Mean_GradTurbVar_Edge = nullptr, /*!< \brief Mean_gradTurbVar DOT Edge_Vector. */ *Proj_Mean_GradTurbVar = nullptr, /*!< \brief Mean_gradTurbVar DOT normal, corrected if required. */ - dist_ij_2 = 0.0, /*!< \brief |Edge_Vector|^2 */ proj_vector_ij = 0.0, /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ *Flux = nullptr, /*!< \brief Final result, diffusive flux/residual. */ **Jacobian_i = nullptr, /*!< \brief Flux Jacobian w.r.t. node i. */ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 677f19a7bb8..0b40964ac3a 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1932,8 +1932,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - numerics[iMGlevel][HEAT_SOL][visc_term] = new CAvgGradCorrected_Heat(nDim, nVar_Heat, config); - numerics[iMGlevel][HEAT_SOL][visc_bound_term] = new CAvgGrad_Heat(nDim, nVar_Heat, config); + numerics[iMGlevel][HEAT_SOL][visc_term] = new CAvgGrad_Heat(nDim, nVar_Heat, config, true); + numerics[iMGlevel][HEAT_SOL][visc_bound_term] = new CAvgGrad_Heat(nDim, nVar_Heat, config, false); switch (config->GetKind_ConvNumScheme_Heat()) { diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index aec05d975f7..863828c7350 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -33,20 +33,10 @@ CNumerics::CNumerics(void) { - Normal = nullptr; - UnitNormal = nullptr; - UnitNormald = nullptr; - Proj_Flux_Tensor = nullptr; - Flux_Tensor = nullptr; tau = nullptr; - Vector = nullptr; - - l = nullptr; - m = nullptr; - using_uq = false; nemo = false; @@ -56,9 +46,7 @@ CNumerics::CNumerics(void) { CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) { - unsigned short iVar, iDim; - - Normal = nullptr; + unsigned short iDim; nDim = val_nDim; nVar = val_nVar; @@ -69,13 +57,6 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Prandtl_Turb = config->GetPrandtl_Turb(); Gas_Constant = config->GetGas_ConstantND(); - UnitNormal = new su2double [nDim] (); - UnitNormald = new su2double [nDim] (); - - Flux_Tensor = new su2double* [nVar]; - for (iVar = 0; iVar < (nVar); iVar++) - Flux_Tensor[iVar] = new su2double [nDim] (); - tau = new su2double* [nDim]; for (iDim = 0; iDim < nDim; iDim++) tau[iDim] = new su2double [nDim] (); @@ -85,11 +66,6 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, turb_ke_i = 0.0; turb_ke_j = 0.0; - Vector = new su2double[nDim] (); - - l = new su2double [nDim] (); - m = new su2double [nDim] (); - Dissipation_ij = 1.0; /* --- Initializing variables for the UQ methodology --- */ @@ -103,63 +79,14 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, CNumerics::~CNumerics(void) { - delete [] UnitNormal; - delete [] UnitNormald; - // visc delete [] Proj_Flux_Tensor; - if (Flux_Tensor) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) - delete [] Flux_Tensor[iVar]; - delete [] Flux_Tensor; - } - if (tau) { for (unsigned short iDim = 0; iDim < nDim; iDim++) delete [] tau[iDim]; delete [] tau; } - - delete [] Vector; - - delete [] l; - delete [] m; -} - -void CNumerics::GetInviscidFlux(su2double val_density, const su2double *val_velocity, - su2double val_pressure, su2double val_enthalpy) { - if (nDim == 3) { - Flux_Tensor[0][0] = val_density*val_velocity[0]; - Flux_Tensor[1][0] = Flux_Tensor[0][0]*val_velocity[0]+val_pressure; - Flux_Tensor[2][0] = Flux_Tensor[0][0]*val_velocity[1]; - Flux_Tensor[3][0] = Flux_Tensor[0][0]*val_velocity[2]; - Flux_Tensor[4][0] = Flux_Tensor[0][0]*val_enthalpy; - - Flux_Tensor[0][1] = val_density*val_velocity[1]; - Flux_Tensor[1][1] = Flux_Tensor[0][1]*val_velocity[0]; - Flux_Tensor[2][1] = Flux_Tensor[0][1]*val_velocity[1]+val_pressure; - Flux_Tensor[3][1] = Flux_Tensor[0][1]*val_velocity[2]; - Flux_Tensor[4][1] = Flux_Tensor[0][1]*val_enthalpy; - - Flux_Tensor[0][2] = val_density*val_velocity[2]; - Flux_Tensor[1][2] = Flux_Tensor[0][2]*val_velocity[0]; - Flux_Tensor[2][2] = Flux_Tensor[0][2]*val_velocity[1]; - Flux_Tensor[3][2] = Flux_Tensor[0][2]*val_velocity[2]+val_pressure; - Flux_Tensor[4][2] = Flux_Tensor[0][2]*val_enthalpy; - - } - else { - Flux_Tensor[0][0] = val_density*val_velocity[0]; - Flux_Tensor[1][0] = Flux_Tensor[0][0]*val_velocity[0]+val_pressure; - Flux_Tensor[2][0] = Flux_Tensor[0][0]*val_velocity[1]; - Flux_Tensor[3][0] = Flux_Tensor[0][0]*val_enthalpy; - - Flux_Tensor[0][1] = val_density*val_velocity[1]; - Flux_Tensor[1][1] = Flux_Tensor[0][1]*val_velocity[0]; - Flux_Tensor[2][1] = Flux_Tensor[0][1]*val_velocity[1]+val_pressure; - Flux_Tensor[3][1] = Flux_Tensor[0][1]*val_enthalpy; - } } void CNumerics::GetInviscidProjFlux(const su2double *val_density, @@ -1678,7 +1605,7 @@ void CNumerics::GetPrimitive2Conservative (const su2double *val_Mean_PrimVar, } } -void CNumerics::CreateBasis(const su2double *val_Normal) { +void CNumerics::CreateBasis(const su2double *val_Normal, su2double* l, su2double* m) { unsigned short iDim; su2double modm, modl; diff --git a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp index 9ef54677d83..5c7fcdc1a3c 100644 --- a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp +++ b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp @@ -244,6 +244,8 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, su2double rho, T, Tve, RuSI, Ru; auto& Ms = fluidmodel->GetSpeciesMolarMass(); + su2activematrix Flux_Tensor(nVar,nDim); + /*--- Initialize ---*/ for (iVar = 0; iVar < nVar; iVar++) { Proj_Flux_Tensor[iVar] = 0.0; @@ -282,8 +284,10 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar, //kve += Cpve*(val_eddy_viscosity/Prandtl_Turb); /*--- Pre-compute mixture quantities ---*/ + + su2double Vector[MAXNDIM] = {0.0}; + for (iDim = 0; iDim < nDim; iDim++) { - Vector[iDim] = 0.0; for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { Vector[iDim] += rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim]; } diff --git a/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp b/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp index a342081fac8..239df6b3005 100644 --- a/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp +++ b/SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp @@ -297,10 +297,10 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi unsigned short iDim, iSpecies, iVar; su2double rho, rhov, vel2, H, yinv, T, Tve, Ru, RuSI; - su2double *Ds, ktr, kve; + su2double ktr, kve; /*--- Rename for convenience ---*/ - Ds = Diffusion_Coeff_i; + auto Ds = Diffusion_Coeff_i; ktr = Thermal_Conductivity_i; kve = Thermal_Conductivity_ve_i; rho = V_i[RHO_INDEX]; @@ -346,6 +346,10 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi if (viscous) { if (!rans){ turb_ke_i = 0.0; } + su2double Vector = 0.0; + for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) + Vector += rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1]; + su2double sumJhs_y = 0.0; su2double sumJeve_y = 0.0; su2double Mass = 0.0; @@ -362,8 +366,8 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi /*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/ for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { - sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies]; - sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies]; + sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector) * hs[iSpecies]; + sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector) * eve_i[iSpecies]; } for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) diff --git a/SU2_CFD/src/numerics/NEMO/convection/msw.cpp b/SU2_CFD/src/numerics/NEMO/convection/msw.cpp index b0503a30856..57f9ef86d0a 100644 --- a/SU2_CFD/src/numerics/NEMO/convection/msw.cpp +++ b/SU2_CFD/src/numerics/NEMO/convection/msw.cpp @@ -198,7 +198,8 @@ CNumerics::ResidualType<> CUpwMSW_NEMO::ComputeResidual(const CConfig *config) { epsilon*epsilon)); /*--- Compute projected P, invP, and Lambda ---*/ - CreateBasis(UnitNormal); + su2double l[MAXNDIM], m[MAXNDIM]; + CreateBasis(UnitNormal,l,m); GetPMatrix(Ust_i, Vst_i, dPdUst_i, UnitNormal, l, m, P_Tensor); GetPMatrix_inv(Ust_i, Vst_i, dPdUst_i, UnitNormal, l, m, invP_Tensor); @@ -232,7 +233,7 @@ CNumerics::ResidualType<> CUpwMSW_NEMO::ComputeResidual(const CConfig *config) { epsilon*epsilon)); /*--- Compute projected P, invP, and Lambda ---*/ - CreateBasis(UnitNormal); + CreateBasis(UnitNormal,l,m); GetPMatrix(Ust_j, Vst_j, dPdUst_j, UnitNormal, l, m, P_Tensor); GetPMatrix_inv(Ust_j, Vst_j, dPdUst_j, UnitNormal, l, m, invP_Tensor); diff --git a/SU2_CFD/src/numerics/NEMO/convection/roe.cpp b/SU2_CFD/src/numerics/NEMO/convection/roe.cpp index bdf97a17e9c..841b3a0b3b5 100644 --- a/SU2_CFD/src/numerics/NEMO/convection/roe.cpp +++ b/SU2_CFD/src/numerics/NEMO/convection/roe.cpp @@ -109,7 +109,8 @@ CNumerics::ResidualType<> CUpwRoe_NEMO::ComputeResidual(const CConfig *config) { fluidmodel->ComputedPdU(RoeV, roe_eves, RoedPdU); /*--- Calculate dual grid tangent vectors for P & invP ---*/ - CreateBasis(UnitNormal); + su2double l[MAXNDIM], m[MAXNDIM]; + CreateBasis(UnitNormal,l,m); /*--- Compute the inviscid projected fluxes ---*/ GetInviscidProjFlux(U_i, V_i, Normal, ProjFlux_i); diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 9062cf1f04d..6dbcea6cda6 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -212,6 +212,8 @@ void CAvgGrad_Base::GetViscousProjFlux(const su2double *val_primvar, /*--- Primitive variables -> [Temp vel_x vel_y vel_z Pressure] ---*/ + su2double Flux_Tensor[5][3]; + if (nDim == 2) { Flux_Tensor[0][0] = 0.0; Flux_Tensor[1][0] = tau[0][0]; @@ -676,6 +678,8 @@ void CAvgGradInc_Flow::GetViscousIncProjFlux(const su2double* const *val_gradpri /*--- Gradient of primitive variables -> [Pressure vel_x vel_y vel_z Temperature] ---*/ + su2double Flux_Tensor[5][3]; + if (nDim == 2) { Flux_Tensor[0][0] = 0.0; Flux_Tensor[1][0] = tau[0][0]; diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 24f6356f824..887cd8f6079 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -674,7 +674,7 @@ CNumerics::ResidualType<> CSourceWindGust::ComputeResidual(const CConfig* config smx = rho*(du_gust_dt + (u+u_gust)*du_gust_dx + (v+v_gust)*du_gust_dy); smy = rho*(dv_gust_dt + (u+u_gust)*dv_gust_dx + (v+v_gust)*dv_gust_dy); //smz = rho*(dw_gust_dt + (u+u_gust)*dw_gust_dx + (v+v_gust)*dw_gust_dy) + (w+w_gust)*dw_gust_dz; - + se = u*smx + v*smy + p*(du_gust_dx + dv_gust_dy); //se = u*smx + v*smy + w*smz + p*(du_gust_dx + dv_gust_dy + dw_gust_dz); @@ -777,7 +777,7 @@ CNumerics::ResidualType<> CSourceIncStreamwisePeriodic_Outlet::ComputeResidual(c /*--- Force the area avg inlet Temp to match the Inc_Temperature_Init with additional residual contribution ---*/ const su2double delta_T = SPvals.Streamwise_Periodic_InletTemperature - config->GetInc_Temperature_Init()/config->GetTemperature_Ref(); - residual[nDim+1] += 0.5 * abs(local_Massflow) * SpecificHeat_i * delta_T; + residual[nDim+1] += 0.5 * abs(local_Massflow) * Cp_i * delta_T; return ResidualType<>(residual, jacobian, nullptr); } diff --git a/SU2_CFD/src/numerics/heat.cpp b/SU2_CFD/src/numerics/heat.cpp index 05bddd0ff91..389f6c2cec2 100644 --- a/SU2_CFD/src/numerics/heat.cpp +++ b/SU2_CFD/src/numerics/heat.cpp @@ -169,33 +169,17 @@ void CUpwSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jaco } -CAvgGrad_Heat::CAvgGrad_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : +CAvgGrad_Heat::CAvgGrad_Heat(unsigned short val_nDim, unsigned short val_nVar, + const CConfig *config, bool correct_) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Heat() == EULER_IMPLICIT); - - Edge_Vector = new su2double [nDim]; - Proj_Mean_GradHeatVar_Normal = new su2double [nVar]; - Proj_Mean_GradHeatVar_Corrected = new su2double [nVar]; - Mean_GradHeatVar = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Mean_GradHeatVar[iVar] = new su2double [nDim]; - -} - -CAvgGrad_Heat::~CAvgGrad_Heat(void) { - - delete [] Edge_Vector; - delete [] Proj_Mean_GradHeatVar_Normal; - delete [] Proj_Mean_GradHeatVar_Corrected; - for (iVar = 0; iVar < nVar; iVar++) - delete [] Mean_GradHeatVar[iVar]; - delete [] Mean_GradHeatVar; - + correct = correct_; } void CAvgGrad_Heat::ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) { + constexpr int nVar = 1; AD::StartPreacc(); AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); @@ -204,31 +188,15 @@ void CAvgGrad_Heat::ComputeResidual(su2double *val_residual, su2double **Jacobia AD::SetPreaccIn(ConsVar_Grad_i[0],nDim); AD::SetPreaccIn(ConsVar_Grad_j[0],nDim); AD::SetPreaccIn(Thermal_Diffusivity_i); AD::SetPreaccIn(Thermal_Conductivity_j); - Thermal_Diffusivity_Mean = 0.5*(Thermal_Diffusivity_i + Thermal_Diffusivity_j); + su2double NormalGrad[nVar], CorrectedGrad[nVar]; - /*--- Compute vector going from iPoint to jPoint ---*/ + auto proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ConsVar_Grad_i, + ConsVar_Grad_j, correct, &Temp_i, &Temp_j, + NormalGrad, CorrectedGrad); - dist_ij_2 = 0; proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; - - /*--- Mean gradient approximation. Projection of the mean gradient in the direction of the edge ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradHeatVar_Normal[iVar] = 0.0; - Proj_Mean_GradHeatVar_Corrected[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradHeatVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); - Proj_Mean_GradHeatVar_Normal[iVar] += Mean_GradHeatVar[iVar][iDim]*Normal[iDim]; - } - Proj_Mean_GradHeatVar_Corrected[iVar] = Proj_Mean_GradHeatVar_Normal[iVar]; - } + const su2double Thermal_Diffusivity_Mean = 0.5*(Thermal_Diffusivity_i + Thermal_Diffusivity_j); - val_residual[0] = Thermal_Diffusivity_Mean*Proj_Mean_GradHeatVar_Corrected[0]; + val_residual[0] = Thermal_Diffusivity_Mean*CorrectedGrad[0]; /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ if (implicit) { @@ -240,83 +208,3 @@ void CAvgGrad_Heat::ComputeResidual(su2double *val_residual, su2double **Jacobia AD::EndPreacc(); } - -CAvgGradCorrected_Heat::CAvgGradCorrected_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : - CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Heat() == EULER_IMPLICIT); - - Edge_Vector = new su2double [nDim]; - Proj_Mean_GradHeatVar_Edge = new su2double [nVar]; - Proj_Mean_GradHeatVar_Kappa = new su2double [nVar]; - Proj_Mean_GradHeatVar_Corrected = new su2double [nVar]; - Mean_GradHeatVar = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Mean_GradHeatVar[iVar] = new su2double [nDim]; - -} - -CAvgGradCorrected_Heat::~CAvgGradCorrected_Heat(void) { - - delete [] Edge_Vector; - delete [] Proj_Mean_GradHeatVar_Edge; - delete [] Proj_Mean_GradHeatVar_Kappa; - delete [] Proj_Mean_GradHeatVar_Corrected; - for (iVar = 0; iVar < nVar; iVar++) - delete [] Mean_GradHeatVar[iVar]; - delete [] Mean_GradHeatVar; - -} - -void CAvgGradCorrected_Heat::ComputeResidual(su2double *val_residual, su2double **Jacobian_i, - su2double **Jacobian_j, CConfig *config) { - - AD::StartPreacc(); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(Temp_i); AD::SetPreaccIn(Temp_j); - AD::SetPreaccIn(ConsVar_Grad_i[0],nDim); AD::SetPreaccIn(ConsVar_Grad_j[0],nDim); - AD::SetPreaccIn(Thermal_Diffusivity_i); AD::SetPreaccIn(Thermal_Diffusivity_j); - - Thermal_Diffusivity_Mean = 0.5*(Thermal_Diffusivity_i + Thermal_Diffusivity_j); - - /*--- Compute vector going from iPoint to jPoint ---*/ - - dist_ij_2 = 0; proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; - - /*--- Mean gradient approximation. Projection of the mean gradient - in the direction of the edge ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradHeatVar_Edge[iVar] = 0.0; - Proj_Mean_GradHeatVar_Kappa[iVar] = 0.0; - - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradHeatVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); - Proj_Mean_GradHeatVar_Kappa[iVar] += Mean_GradHeatVar[iVar][iDim]*Normal[iDim]; - Proj_Mean_GradHeatVar_Edge[iVar] += Mean_GradHeatVar[iVar][iDim]*Edge_Vector[iDim]; - } - Proj_Mean_GradHeatVar_Corrected[iVar] = Proj_Mean_GradHeatVar_Kappa[iVar]; - Proj_Mean_GradHeatVar_Corrected[iVar] -= Proj_Mean_GradHeatVar_Edge[iVar]*proj_vector_ij - - (Temp_j-Temp_i)*proj_vector_ij; - } - - val_residual[0] = Thermal_Diffusivity_Mean*Proj_Mean_GradHeatVar_Corrected[0]; - - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ - - if (implicit) { - Jacobian_i[0][0] = -Thermal_Diffusivity_Mean*proj_vector_ij; - Jacobian_j[0][0] = Thermal_Diffusivity_Mean*proj_vector_ij; - } - - AD::SetPreaccOut(val_residual, nVar); - AD::EndPreacc(); -} diff --git a/SU2_CFD/src/numerics/radiation.cpp b/SU2_CFD/src/numerics/radiation.cpp index 1a1d098a82d..f1049c2db10 100644 --- a/SU2_CFD/src/numerics/radiation.cpp +++ b/SU2_CFD/src/numerics/radiation.cpp @@ -74,30 +74,11 @@ CAvgGradCorrected_P1::CAvgGradCorrected_P1(unsigned short val_nDim, unsigned sho CNumericsRadiation(val_nDim, val_nVar, config) { GammaP1 = 1.0 / (3.0*(Absorption_Coeff + Scattering_Coeff)); - - Edge_Vector = new su2double [nDim]; - Proj_Mean_GradP1Var_Edge = new su2double [nVar]; - Proj_Mean_GradP1Var_Kappa = new su2double [nVar]; - Proj_Mean_GradP1Var_Corrected = new su2double [nVar]; - Mean_GradP1Var = new su2double* [nVar]; - for (unsigned short iVar = 0; iVar < nVar; iVar++) - Mean_GradP1Var[iVar] = new su2double [nDim]; -} - -CAvgGradCorrected_P1::~CAvgGradCorrected_P1(void) { - - delete [] Edge_Vector; - delete [] Proj_Mean_GradP1Var_Edge; - delete [] Proj_Mean_GradP1Var_Kappa; - delete [] Proj_Mean_GradP1Var_Corrected; - for (unsigned short iVar = 0; iVar < nVar; iVar++) - delete [] Mean_GradP1Var[iVar]; - delete [] Mean_GradP1Var; } void CAvgGradCorrected_P1::ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) { - unsigned short iVar, iDim; + constexpr int nVar = 1; AD::StartPreacc(); AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); @@ -105,35 +86,13 @@ void CAvgGradCorrected_P1::ComputeResidual(su2double *val_residual, su2double ** AD::SetPreaccIn(RadVar_i,nVar); AD::SetPreaccIn(RadVar_j,nVar); AD::SetPreaccIn(RadVar_Grad_i,nVar,nDim); AD::SetPreaccIn(RadVar_Grad_j,nVar,nDim); - /*--- Compute vector going from iPoint to jPoint ---*/ + su2double NormalGrad[nVar], CorrectedGrad[nVar]; - dist_ij = 0; proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij; - - /*--- Mean gradient approximation. Projection of the mean gradient - in the direction of the edge ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradP1Var_Edge[iVar] = 0.0; - Proj_Mean_GradP1Var_Kappa[iVar] = 0.0; - - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradP1Var[iVar][iDim] = 0.5*(RadVar_Grad_i[iVar][iDim] + RadVar_Grad_j[iVar][iDim]); - Proj_Mean_GradP1Var_Kappa[iVar] += Mean_GradP1Var[iVar][iDim]*Normal[iDim]; - Proj_Mean_GradP1Var_Edge[iVar] += Mean_GradP1Var[iVar][iDim]*Edge_Vector[iDim]; - } - Proj_Mean_GradP1Var_Corrected[iVar] = Proj_Mean_GradP1Var_Kappa[iVar]; - Proj_Mean_GradP1Var_Corrected[iVar] -= Proj_Mean_GradP1Var_Edge[iVar]*proj_vector_ij - - (RadVar_j[iVar]-RadVar_i[iVar])*proj_vector_ij; - } + auto proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, RadVar_Grad_i, + RadVar_Grad_j, true, RadVar_i, RadVar_j, + NormalGrad, CorrectedGrad); - val_residual[0] = GammaP1*Proj_Mean_GradP1Var_Corrected[0]; + val_residual[0] = GammaP1*CorrectedGrad[0]; /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ diff --git a/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp b/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp index 98cd23203ce..aa97f45bb14 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_diffusion.cpp @@ -38,7 +38,6 @@ CAvgGrad_Scalar::CAvgGrad_Scalar(unsigned short val_nDim, incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) { Proj_Mean_GradTurbVar_Normal = new su2double [nVar] (); - Proj_Mean_GradTurbVar_Edge = new su2double [nVar] (); Proj_Mean_GradTurbVar = new su2double [nVar] (); Flux = new su2double [nVar] (); @@ -53,7 +52,6 @@ CAvgGrad_Scalar::CAvgGrad_Scalar(unsigned short val_nDim, CAvgGrad_Scalar::~CAvgGrad_Scalar(void) { delete [] Proj_Mean_GradTurbVar_Normal; - delete [] Proj_Mean_GradTurbVar_Edge; delete [] Proj_Mean_GradTurbVar; delete [] Flux; @@ -69,8 +67,6 @@ CAvgGrad_Scalar::~CAvgGrad_Scalar(void) { CNumerics::ResidualType<> CAvgGrad_Scalar::ComputeResidual(const CConfig* config) { - unsigned short iVar, iDim; - AD::StartPreacc(); AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); AD::SetPreaccIn(Normal, nDim); @@ -96,37 +92,9 @@ CNumerics::ResidualType<> CAvgGrad_Scalar::ComputeResidual(const CConfig* config Eddy_Viscosity_i = V_i[nDim+6]; Eddy_Viscosity_j = V_j[nDim+6]; } - /*--- Compute vector going from iPoint to jPoint ---*/ - - dist_ij_2 = 0; proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; - - /*--- Mean gradient approximation ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradTurbVar_Normal[iVar] = 0.0; - Proj_Mean_GradTurbVar_Edge[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - su2double Mean_GradTurbVar = 0.5*(TurbVar_Grad_i[iVar][iDim] + - TurbVar_Grad_j[iVar][iDim]); - - Proj_Mean_GradTurbVar_Normal[iVar] += Mean_GradTurbVar * Normal[iDim]; - - if (correct_gradient) - Proj_Mean_GradTurbVar_Edge[iVar] += Mean_GradTurbVar * Edge_Vector[iDim]; - } - Proj_Mean_GradTurbVar[iVar] = Proj_Mean_GradTurbVar_Normal[iVar]; - if (correct_gradient) { - Proj_Mean_GradTurbVar[iVar] -= Proj_Mean_GradTurbVar_Edge[iVar]*proj_vector_ij - - (TurbVar_j[iVar]-TurbVar_i[iVar])*proj_vector_ij; - } - } - + proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, TurbVar_Grad_i, + TurbVar_Grad_j, correct_gradient, TurbVar_i, TurbVar_j, + Proj_Mean_GradTurbVar_Normal, Proj_Mean_GradTurbVar); FinishResidualCalc(config); AD::SetPreaccOut(Flux, nVar); From 7d1eb8c2fdf063f87bf488cae28b209ce86eba68 Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sun, 27 Jun 2021 18:07:11 +0100 Subject: [PATCH 4/8] correct brief --- SU2_CFD/include/numerics/CNumerics.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 35822fbde2a..5eee0764f15 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -628,7 +628,7 @@ class CNumerics { } /*! - * \brief Perturb the Reynolds stress tensor based on parameters. + * \brief Project average gradient onto normal (with or w/o correction) for viscous fluxes of scalar quantities. * \param[in] nDim - Dimension of the flow problem, 2 or 3. * \param[in] nVar - Number of variables. * \param[in] normal - Area vector. From 72e119a496df0e01dd3bbedb79344cdd2970c579 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 27 Jun 2021 22:22:31 +0100 Subject: [PATCH 5/8] more cleaning + update regressions --- SU2_CFD/include/numerics/heat.hpp | 31 +----- SU2_CFD/src/numerics/heat.cpp | 98 +++++++++---------- .../chtPinArray_2d/of_grad_findiff.csv.ref | 2 +- TestCases/serial_regression.py | 4 +- TestCases/serial_regression_AD.py | 2 +- 5 files changed, 54 insertions(+), 83 deletions(-) diff --git a/SU2_CFD/include/numerics/heat.hpp b/SU2_CFD/include/numerics/heat.hpp index ae7fb9ce1f9..a6e56488f41 100644 --- a/SU2_CFD/include/numerics/heat.hpp +++ b/SU2_CFD/include/numerics/heat.hpp @@ -37,18 +37,10 @@ * \version 7.1.1 "Blackbird" */ class CCentSca_Heat : public CNumerics { - private: - unsigned short iDim; /*!< \brief Iteration on dimension and variables. */ - su2double *Diff_Lapl, /*!< \brief Diference of conservative variables and undivided laplacians. */ - *MeanVelocity, ProjVelocity, - ProjVelocity_i, ProjVelocity_j, /*!< \brief Mean and projected velocities. */ - Param_Kappa_4, /*!< \brief Artificial dissipation parameters. */ - Local_Lambda_i, Local_Lambda_j, - MeanLambda, /*!< \brief Local eingenvalues. */ - cte_0, cte_1; /*!< \brief Artificial dissipation values. */ - bool implicit, /*!< \brief Implicit calculation. */ - dynamic_grid; /*!< \brief Modification for grid movement. */ + su2double Param_Kappa_4; /*!< \brief Artificial dissipation parameters. */ + bool implicit; /*!< \brief Implicit calculation. */ + bool dynamic_grid; /*!< \brief Modification for grid movement. */ public: /*! @@ -57,12 +49,7 @@ class CCentSca_Heat : public CNumerics { * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CCentSca_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CCentSca_Heat(void) override; + CCentSca_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config); /*! * \brief Compute the flow residual using a JST method. @@ -85,10 +72,7 @@ class CCentSca_Heat : public CNumerics { */ class CUpwSca_Heat : public CNumerics { private: - su2double *Velocity_i, *Velocity_j; bool implicit, dynamic_grid; - su2double q_ij, a0, a1; - unsigned short iDim; public: /*! @@ -97,12 +81,7 @@ class CUpwSca_Heat : public CNumerics { * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CUpwSca_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CUpwSca_Heat(void) override; + CUpwSca_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config); /*! * \brief Compute the scalar upwind flux between two nodes i and j. diff --git a/SU2_CFD/src/numerics/heat.cpp b/SU2_CFD/src/numerics/heat.cpp index 389f6c2cec2..31c48037503 100644 --- a/SU2_CFD/src/numerics/heat.cpp +++ b/SU2_CFD/src/numerics/heat.cpp @@ -28,32 +28,17 @@ #include "../../include/numerics/heat.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" -CCentSca_Heat::CCentSca_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : +CCentSca_Heat::CCentSca_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT); - /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - - MeanVelocity = new su2double [nDim]; - - Laminar_Viscosity_i = config->GetViscosity_FreeStreamND(); - Laminar_Viscosity_j = config->GetViscosity_FreeStreamND(); - Param_Kappa_4 = config->GetKappa_4th_Heat(); - Diff_Lapl = new su2double [nVar]; -} - -CCentSca_Heat::~CCentSca_Heat(void) { - - delete [] MeanVelocity; - delete [] Diff_Lapl; } void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { - AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+3); AD::SetPreaccIn(V_j, nDim+3); AD::SetPreaccIn(Temp_i); AD::SetPreaccIn(Temp_j); @@ -71,16 +56,27 @@ void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jac /*--- Projected velocities at the current edge ---*/ - ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjVelocity_i += V_i[iDim+1]*Normal[iDim]; - ProjVelocity_j += V_j[iDim+1]*Normal[iDim]; + su2double ProjVelocity_i = 0.0; + su2double ProjVelocity_j = 0.0; + + if (dynamic_grid) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + su2double Velocity_i = V_i[iDim+1] - GridVel_i[iDim]; + su2double Velocity_j = V_j[iDim+1] - GridVel_j[iDim]; + ProjVelocity_i += Velocity_i*Normal[iDim]; + ProjVelocity_j += Velocity_j*Normal[iDim]; + } + } + else { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + ProjVelocity_i += V_i[iDim+1]*Normal[iDim]; + ProjVelocity_j += V_j[iDim+1]*Normal[iDim]; + } } - Area = GeometryToolbox::Norm(nDim, Normal); /*--- Computing the second order centered scheme part ---*/ - ProjVelocity = 0.5*(ProjVelocity_i+ProjVelocity_j); + su2double ProjVelocity = 0.5*(ProjVelocity_i+ProjVelocity_j); val_residual[0] = 0.5*(Temp_i + Temp_j)*ProjVelocity; @@ -91,20 +87,21 @@ void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jac /*--- Adding artificial dissipation to stabilize the centered scheme ---*/ - Diff_Lapl[0] = Und_Lapl_i[0]-Und_Lapl_j[0]; + su2double Diff_Lapl = Und_Lapl_i[0]-Und_Lapl_j[0]; + su2double Area2 = GeometryToolbox::SquaredNorm(nDim, Normal); - SoundSpeed_i = sqrt(ProjVelocity_i*ProjVelocity_i + (BetaInc2_i/DensityInc_i)*Area*Area); - SoundSpeed_j = sqrt(ProjVelocity_j*ProjVelocity_j + (BetaInc2_j/DensityInc_j)*Area*Area); + su2double SoundSpeed_i = sqrt(pow(ProjVelocity_i,2) + (BetaInc2_i/DensityInc_i)*Area2); + su2double SoundSpeed_j = sqrt(pow(ProjVelocity_j,2) + (BetaInc2_j/DensityInc_j)*Area2); - Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i; - Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j; - MeanLambda = 0.5*(Local_Lambda_i+Local_Lambda_j); + su2double Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i; + su2double Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j; + su2double MeanLambda = 0.5*(Local_Lambda_i+Local_Lambda_j); - val_residual[0] += -Param_Kappa_4*Diff_Lapl[0]*MeanLambda; + val_residual[0] += -Param_Kappa_4*Diff_Lapl*MeanLambda; if (implicit) { - cte_0 = Param_Kappa_4*su2double(Neighbor_i+1)*MeanLambda; - cte_1 = Param_Kappa_4*su2double(Neighbor_j+1)*MeanLambda; + su2double cte_0 = Param_Kappa_4*su2double(Neighbor_i+1)*MeanLambda; + su2double cte_1 = Param_Kappa_4*su2double(Neighbor_j+1)*MeanLambda; val_Jacobian_i[0][0] += cte_0; val_Jacobian_j[0][0] -= cte_1; @@ -115,32 +112,17 @@ void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jac } -CUpwSca_Heat::CUpwSca_Heat(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : +CUpwSca_Heat::CUpwSca_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT); - /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - Velocity_i = new su2double [nDim]; - Velocity_j = new su2double [nDim]; - - Laminar_Viscosity_i = config->GetViscosity_FreeStreamND(); - Laminar_Viscosity_j = config->GetViscosity_FreeStreamND(); -} - -CUpwSca_Heat::~CUpwSca_Heat(void) { - - delete [] Velocity_i; - delete [] Velocity_j; - } void CUpwSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { - q_ij = 0.0; - AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+1); AD::SetPreaccIn(V_j, nDim+1); AD::SetPreaccIn(Temp_i); AD::SetPreaccIn(Temp_j); @@ -149,14 +131,24 @@ void CUpwSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jaco AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); } - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - Velocity_j[iDim] = V_j[iDim+1]; - q_ij += 0.5*(Velocity_i[iDim]+Velocity_j[iDim])*Normal[iDim]; + su2double q_ij = 0.0; + + if (dynamic_grid) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + su2double Velocity_i = V_i[iDim+1] - GridVel_i[iDim]; + su2double Velocity_j = V_j[iDim+1] - GridVel_j[iDim]; + q_ij += 0.5*(Velocity_i+Velocity_j)*Normal[iDim]; + } } + else { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + q_ij += 0.5*(V_i[iDim+1]+V_j[iDim+1])*Normal[iDim]; + } + } + + su2double a0 = 0.5*(q_ij+fabs(q_ij)); + su2double a1 = 0.5*(q_ij-fabs(q_ij)); - a0 = 0.5*(q_ij+fabs(q_ij)); - a1 = 0.5*(q_ij-fabs(q_ij)); val_residual[0] = a0*Temp_i+a1*Temp_j; if (implicit) { diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref index 46ec27a655b..25785d5d61f 100644 --- a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref +++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref @@ -1,2 +1,2 @@ "VARIABLE" , "AVG_DENSITY[0]", "AVG_ENTHALPY[0]", "AVG_NORMALVEL[0]", "DRAG[0]" , "EFFICIENCY[0]" , "FORCE_X[0]" , "FORCE_Y[0]" , "FORCE_Z[0]" , "LIFT[0]" , "MOMENT_X[0]" , "MOMENT_Y[0]" , "MOMENT_Z[0]" , "SIDEFORCE[0]" , "SURFACE_MACH[0]", "SURFACE_MASSFLOW[0]", "SURFACE_MOM_DISTORTION[0]", "SURFACE_PRESSURE_DROP[0]", "SURFACE_SECONDARY[0]", "SURFACE_SECOND_OVER_UNIFORM[0]", "SURFACE_STATIC_PRESSURE[0]", "SURFACE_STATIC_TEMPERATURE[0]", "SURFACE_TOTAL_PRESSURE[0]", "SURFACE_TOTAL_TEMPERATURE[0]", "SURFACE_UNIFORMITY[0]", "AVG_TEMPERATURE[1]", "MAXIMUM_HEATFLUX[1]", "TOTAL_HEATFLUX[1]", "FINDIFF_STEP" -0 , 0.0 , -99999.96982514858, 2.2204999999731917e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -0.04999999997368221, -1.1100000001884e-08, -2.069999999187999 , 0.0 , 2.12000000054946 , 3.7100000016554446 , 330.00000030369847 , -39.999997625272954 , 315.0000011942211 , -40.000008993956726 , -1.400000004814217 , -149.9999996212864, 0.0 , -469.9999976764957, 1e-08 +0 , 0.0 , -99999.96982514858, 3.330700000025998e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -0.04999999997368221, -1.1100000001884e-08, -2.069999999187999 , 0.0 , 2.12000000054946 , 3.7100000016554446 , 330.00000030369847 , -39.999997625272954 , 315.0000011942211 , -40.000008993956726 , -1.400000004814217 , -149.9999996212864, 0.0 , -469.9999976764957, 1e-08 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 040d9fba020..7d955ea5e9d 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -345,7 +345,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-11.451011, -12.798258, -5.863895, 1.049989, 0.019163, -1.925028] + turb_naca0012_sst.test_vals = [-11.451010, -12.798258, -5.863895, 1.049989, 0.019163, -1.925018] turb_naca0012_sst.su2_exec = "SU2_CFD" turb_naca0012_sst.new_output = True turb_naca0012_sst.timeout = 3200 @@ -1821,7 +1821,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-11.451011, -12.798258, -5.863895, 1.049989, 0.019163, -1.925028] + pywrapper_turb_naca0012_sst.test_vals = [-11.451010, -12.798258, -5.863895, 1.049989, 0.019163, -1.925018] pywrapper_turb_naca0012_sst.su2_exec = "SU2_CFD.py -f" pywrapper_turb_naca0012_sst.new_output = True pywrapper_turb_naca0012_sst.timeout = 3200 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 6b24acb58da..10ba9d26bc6 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -241,7 +241,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [-2.271573, 0.671242, -3.172000, -8.231600] #last 4 columns + discadj_heat.test_vals = [-2.271569, 0.671288, -3.172000, -8.231500] #last 4 columns discadj_heat.su2_exec = "SU2_CFD_AD" discadj_heat.timeout = 1600 discadj_heat.tol = 0.00001 From a4885544558c69d68a25754eb50ce11dee3efa49 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 27 Jun 2021 22:26:50 +0100 Subject: [PATCH 6/8] not needed member variable --- SU2_CFD/include/numerics/turbulent/turb_convection.hpp | 1 - SU2_CFD/src/numerics/turbulent/turb_convection.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 9270e12e6b4..9e2cea2032a 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -47,7 +47,6 @@ class CUpwScalar : public CNumerics { protected: su2double - q_ij = 0.0, /*!< \brief Projected velocity at the face. */ a0 = 0.0, /*!< \brief The maximum of the face-normal velocity and 0 */ a1 = 0.0, /*!< \brief The minimum of the face-normal velocity and 0 */ *Flux = nullptr, /*!< \brief Final result, diffusive flux/residual. */ diff --git a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp b/SU2_CFD/src/numerics/turbulent/turb_convection.cpp index a2cb1f2abb9..a73a2dfd8b5 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_convection.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_convection.cpp @@ -74,7 +74,7 @@ CNumerics::ResidualType<> CUpwScalar::ComputeResidual(const CConfig* config) { Density_i = V_i[nDim+2]; Density_j = V_j[nDim+2]; - q_ij = 0.0; + su2double q_ij = 0.0; if (dynamic_grid) { for (iDim = 0; iDim < nDim; iDim++) { su2double Velocity_i = V_i[iDim+1] - GridVel_i[iDim]; From 520b97826a5fdc04d120b3bc68831f24656ae09b Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Thu, 15 Jul 2021 19:26:23 +0100 Subject: [PATCH 7/8] Apply suggestions from code review Co-authored-by: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com> --- SU2_CFD/src/numerics/heat.cpp | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/SU2_CFD/src/numerics/heat.cpp b/SU2_CFD/src/numerics/heat.cpp index 31c48037503..c53da856827 100644 --- a/SU2_CFD/src/numerics/heat.cpp +++ b/SU2_CFD/src/numerics/heat.cpp @@ -76,7 +76,7 @@ void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jac /*--- Computing the second order centered scheme part ---*/ - su2double ProjVelocity = 0.5*(ProjVelocity_i+ProjVelocity_j); + const su2double ProjVelocity = 0.5*(ProjVelocity_i+ProjVelocity_j); val_residual[0] = 0.5*(Temp_i + Temp_j)*ProjVelocity; @@ -87,21 +87,21 @@ void CCentSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jac /*--- Adding artificial dissipation to stabilize the centered scheme ---*/ - su2double Diff_Lapl = Und_Lapl_i[0]-Und_Lapl_j[0]; - su2double Area2 = GeometryToolbox::SquaredNorm(nDim, Normal); + const su2double Diff_Lapl = Und_Lapl_i[0]-Und_Lapl_j[0]; + const su2double Area2 = GeometryToolbox::SquaredNorm(nDim, Normal); - su2double SoundSpeed_i = sqrt(pow(ProjVelocity_i,2) + (BetaInc2_i/DensityInc_i)*Area2); - su2double SoundSpeed_j = sqrt(pow(ProjVelocity_j,2) + (BetaInc2_j/DensityInc_j)*Area2); + const su2double SoundSpeed_i = sqrt(pow(ProjVelocity_i,2) + (BetaInc2_i/DensityInc_i)*Area2); + const su2double SoundSpeed_j = sqrt(pow(ProjVelocity_j,2) + (BetaInc2_j/DensityInc_j)*Area2); - su2double Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i; - su2double Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j; - su2double MeanLambda = 0.5*(Local_Lambda_i+Local_Lambda_j); + const su2double Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i; + const su2double Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j; + const su2double MeanLambda = 0.5*(Local_Lambda_i+Local_Lambda_j); val_residual[0] += -Param_Kappa_4*Diff_Lapl*MeanLambda; if (implicit) { - su2double cte_0 = Param_Kappa_4*su2double(Neighbor_i+1)*MeanLambda; - su2double cte_1 = Param_Kappa_4*su2double(Neighbor_j+1)*MeanLambda; + const su2double cte_0 = Param_Kappa_4*su2double(Neighbor_i+1)*MeanLambda; + const su2double cte_1 = Param_Kappa_4*su2double(Neighbor_j+1)*MeanLambda; val_Jacobian_i[0][0] += cte_0; val_Jacobian_j[0][0] -= cte_1; @@ -135,8 +135,8 @@ void CUpwSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jaco if (dynamic_grid) { for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double Velocity_i = V_i[iDim+1] - GridVel_i[iDim]; - su2double Velocity_j = V_j[iDim+1] - GridVel_j[iDim]; + const su2double Velocity_i = V_i[iDim+1] - GridVel_i[iDim]; + const su2double Velocity_j = V_j[iDim+1] - GridVel_j[iDim]; q_ij += 0.5*(Velocity_i+Velocity_j)*Normal[iDim]; } } @@ -146,8 +146,8 @@ void CUpwSca_Heat::ComputeResidual(su2double *val_residual, su2double **val_Jaco } } - su2double a0 = 0.5*(q_ij+fabs(q_ij)); - su2double a1 = 0.5*(q_ij-fabs(q_ij)); + const su2double a0 = 0.5*(q_ij+fabs(q_ij)); + const su2double a1 = 0.5*(q_ij-fabs(q_ij)); val_residual[0] = a0*Temp_i+a1*Temp_j; From 0ad7a0293d49524d61e259749a93be007b8ebced Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Thu, 15 Jul 2021 20:22:52 +0100 Subject: [PATCH 8/8] Update SU2_CFD/include/numerics/CNumerics.hpp --- SU2_CFD/include/numerics/CNumerics.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 5eee0764f15..1e8b9e03778 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -651,6 +651,7 @@ class CNumerics { const Vec2& var_i, const Vec2& var_j, su2double* projNormal, su2double* projCorrected) { + assert(nDim == 2 || nDim == 3); nDim = (nDim > 2)? 3 : 2; su2double edgeVec[MAXNDIM], dist_ij_2 = 0.0, proj_vector_ij = 0.0;