diff --git a/Common/include/grid_movement/CFreeFormDefBox.hpp b/Common/include/grid_movement/CFreeFormDefBox.hpp index 5025d6ec02e0..0a9b59d81e04 100644 --- a/Common/include/grid_movement/CFreeFormDefBox.hpp +++ b/Common/include/grid_movement/CFreeFormDefBox.hpp @@ -756,13 +756,4 @@ class CFreeFormDefBox : public CGridMovement { */ inline unsigned short GetLevel() const { return Level; } - /*! - * \brief Compute the determinant of a 3 by 3 matrix. - * \param[in] val_matrix 3 by 3 matrix. - * \result Determinant of the matrix - */ - inline su2double Determinant_3x3(su2double A00, su2double A01, su2double A02, su2double A10, su2double A11, su2double A12, su2double A20, su2double A21, su2double A22) { - return A00*(A11*A22-A12*A21) - A01*(A10*A22-A12*A20) + A02*(A10*A21-A11*A20); - } - }; diff --git a/Common/include/grid_movement/CVolumetricMovement.hpp b/Common/include/grid_movement/CVolumetricMovement.hpp index 79293aabc109..5f64a6ecf543 100644 --- a/Common/include/grid_movement/CVolumetricMovement.hpp +++ b/Common/include/grid_movement/CVolumetricMovement.hpp @@ -361,26 +361,6 @@ class CVolumetricMovement : public CGridMovement { */ void UpdateGridCoord_Derivatives(CGeometry *geometry, CConfig *config); - /*! - * \brief Compute the determinant of a 3 by 3 matrix. - * 3 by 3 matrix elements - * \param[in] A00 - * \param[in] A01 - * \param[in] A02 - * \param[in] A10 - * \param[in] A11 - * \param[in] A12 - * \param[in] A20 - * \param[in] A21 - * \param[in] A22 - * \result Determinant of the matrix - */ - inline su2double Determinant_3x3(su2double A00, su2double A01, su2double A02, su2double A10, su2double A11, su2double A12, su2double A20, su2double A21, su2double A22) { - return A00*(A11*A22-A12*A21) - A01*(A10*A22-A12*A20) + A02*(A10*A21-A11*A20); - } - - - /*! * \brief Store the number of iterations when moving the mesh. * \param[in] val_nIterMesh - Number of iterations. diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 24d2d5e50c71..36798619b3b3 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -67,8 +67,6 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> unsigned long nElm = 0; /*!< \brief Total number of elements (or number elements on this processor). */ unsigned long nElmDomain = 0; /*!< \brief Total number of elements without Ghost cells. */ unsigned long nVar = 0; /*!< \brief Number of elements in a block. */ - mutable ScalarType dotRes = 0.0; /*!< \brief Result of dot product. to perform a reduction with OpenMP the - variable needs to be declared outside the parallel region. */ /*! * \brief Generic initialization from a scalar or array. @@ -101,7 +99,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> /*! * \brief Default constructor of the class. */ - CSysVector() {} + CSysVector() = default; /*! * \brief Destructor @@ -291,8 +289,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ template ScalarType dot(const VecExpr::CVecExpr& expr) const { + static ScalarType dotRes; /*--- All threads get the same "view" of the vectors and shared variable. ---*/ SU2_OMP_BARRIER + SU2_OMP_MASTER dotRes = 0.0; SU2_OMP_BARRIER diff --git a/Common/include/linear_algebra/blas_structure.hpp b/Common/include/linear_algebra/blas_structure.hpp index 39951ecf33b6..49864652cd36 100644 --- a/Common/include/linear_algebra/blas_structure.hpp +++ b/Common/include/linear_algebra/blas_structure.hpp @@ -143,6 +143,345 @@ class CBlasStructure { } } + /*! + * \brief tred2 + * Author: + * + * Original FORTRAN77 version by Smith, Boyle, Dongarra, Garbow, Ikebe, + * Klema, Moler. + * C++ version by Aashwin Mishra and Jayant Mukhopadhaya. + * + * Reference: + * + * Martin, Reinsch, Wilkinson, + * TRED2, + * Numerische Mathematik, + * Volume 11, pages 181-195, 1968. + * + * James Wilkinson, Christian Reinsch, + * Handbook for Automatic Computation, + * Volume II, Linear Algebra, Part 2, + * Springer, 1971, + * ISBN: 0387054146, + * LC: QA251.W67. + * + * Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, + * Yasuhiko Ikebe, Virginia Klema, Cleve Moler, + * Matrix Eigensystem Routines, EISPACK Guide, + * Lecture Notes in Computer Science, Volume 6, + * Springer Verlag, 1976, + * ISBN13: 978-3540075462, + * LC: QA193.M37 + * + * \param[in,out] V: matrix that needs to be decomposed + * \param[in,out] d: holds eigenvalues + * \param[in,out] e: work vector + * \param[in] n: order of matrix V + */ + template + static void tred2(Mat& V, Vec& d, W& e, int n) { + using Scalar = typename std::decay::type; + + int i,j,k; + + for (j = 0; j < n; j++) { + d[j] = V[n-1][j]; + } + + /* Householder reduction to tridiagonal form. */ + + for (i = n-1; i > 0; i--) { + + /* Scale to avoid under/overflow. */ + + Scalar scale = 0.0; + Scalar h = 0.0; + for (k = 0; k < i; k++) { + scale = scale + fabs(d[k]); + } + if (scale == 0.0) { + e[i] = d[i-1]; + for (j = 0; j < i; j++) { + d[j] = V[i-1][j]; + V[i][j] = 0.0; + V[j][i] = 0.0; + } + } + else { + + /* Generate Householder vector. */ + + for (k = 0; k < i; k++) { + d[k] /= scale; + h += d[k] * d[k]; + } + Scalar f = d[i-1]; + Scalar g = sqrt(h); + if (f > 0) { + g = -g; + } + e[i] = scale * g; + h = h - f * g; + d[i-1] = f - g; + for (j = 0; j < i; j++) { + e[j] = 0.0; + } + + /* Apply similarity transformation to remaining columns. */ + + for (j = 0; j < i; j++) { + f = d[j]; + V[j][i] = f; + g = e[j] + V[j][j] * f; + for (k = j+1; k <= i-1; k++) { + g += V[k][j] * d[k]; + e[k] += V[k][j] * f; + } + e[j] = g; + } + f = 0.0; + for (j = 0; j < i; j++) { + e[j] /= h; + f += e[j] * d[j]; + } + Scalar hh = f / (h + h); + for (j = 0; j < i; j++) { + e[j] -= hh * d[j]; + } + for (j = 0; j < i; j++) { + f = d[j]; + g = e[j]; + for (k = j; k <= i-1; k++) { + V[k][j] -= (f * e[k] + g * d[k]); + } + d[j] = V[i-1][j]; + V[i][j] = 0.0; + } + } + d[i] = h; + } + + /* Accumulate transformations. */ + + for (i = 0; i < n-1; i++) { + V[n-1][i] = V[i][i]; + V[i][i] = 1.0; + Scalar h = d[i+1]; + if (h != 0.0) { + for (k = 0; k <= i; k++) { + d[k] = V[k][i+1] / h; + } + for (j = 0; j <= i; j++) { + Scalar g = 0.0; + for (k = 0; k <= i; k++) { + g += V[k][i+1] * V[k][j]; + } + for (k = 0; k <= i; k++) { + V[k][j] -= g * d[k]; + } + } + } + for (k = 0; k <= i; k++) { + V[k][i+1] = 0.0; + } + } + for (j = 0; j < n; j++) { + d[j] = V[n-1][j]; + V[n-1][j] = 0.0; + } + V[n-1][n-1] = 1.0; + e[0] = 0.0; + } + + /*! + * \brief tql2 + * Author: + * + * Original FORTRAN77 version by Smith, Boyle, Dongarra, Garbow, Ikebe, + * Klema, Moler. + * C++ version by Aashwin Mishra and Jayant Mukhopadhaya. + * + * Reference: + * + * Bowdler, Martin, Reinsch, Wilkinson, + * TQL2, + * Numerische Mathematik, + * Volume 11, pages 293-306, 1968. + * + * James Wilkinson, Christian Reinsch, + * Handbook for Automatic Computation, + * Volume II, Linear Algebra, Part 2, + * Springer, 1971, + * ISBN: 0387054146, + * LC: QA251.W67. + * + * Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, + * Yasuhiko Ikebe, Virginia Klema, Cleve Moler, + * Matrix Eigensystem Routines, EISPACK Guide, + * Lecture Notes in Computer Science, Volume 6, + * Springer Verlag, 1976, + * ISBN13: 978-3540075462, + * LC: QA193.M37 + * + * \param[in,out] V: matrix that will hold the eigenvectors + * \param[in,out] d: array that will hold the ordered eigenvalues + * \param[in,out] e: work vector + * \param[in] n: order of matrix V + */ + template + static void tql2(Mat& V, Vec& d, W& e, int n) { + using Scalar = typename std::decay::type; + + int i,j,k,l; + for (i = 1; i < n; i++) { + e[i-1] = e[i]; + } + e[n-1] = 0.0; + + Scalar f = 0.0; + Scalar tst1 = 0.0; + Scalar eps = pow(2.0,-52.0); + for (l = 0; l < n; l++) { + + /* Find small subdiagonal element */ + + tst1 = max(tst1,(fabs(d[l]) + fabs(e[l]))); + int m = l; + while (m < n) { + if (fabs(e[m]) <= eps*tst1) { + break; + } + m++; + } + + /* If m == l, d[l] is an eigenvalue, */ + /* otherwise, iterate. */ + + if (m > l) { + int iter = 0; + do { + iter = iter + 1; /* (Could check iteration count here.) */ + + /* Compute implicit shift */ + + Scalar g = d[l]; + Scalar p = (d[l+1] - g) / (2.0 * e[l]); + Scalar r = sqrt(p*p+1.0); + if (p < 0) { + r = -r; + } + d[l] = e[l] / (p + r); + d[l+1] = e[l] * (p + r); + Scalar dl1 = d[l+1]; + Scalar h = g - d[l]; + for (i = l+2; i < n; i++) { + d[i] -= h; + } + f = f + h; + + /* Implicit QL transformation. */ + + p = d[m]; + Scalar c = 1.0; + Scalar c2 = c; + Scalar c3 = c; + Scalar el1 = e[l+1]; + Scalar s = 0.0; + Scalar s2 = 0.0; + for (i = m-1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = sqrt(p*p+e[i]*e[i]); + e[i+1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i+1] = h + s * (c * g + s * d[i]); + + /* Accumulate transformation. */ + + for (k = 0; k < n; k++) { + h = V[k][i+1]; + V[k][i+1] = s * V[k][i] + c * h; + V[k][i] = c * V[k][i] - s * h; + } + } + p = -s * s2 * c3 * el1 * e[l] / dl1; + e[l] = s * p; + d[l] = c * p; + + /* Check for convergence. */ + + } while (fabs(e[l]) > eps*tst1); + } + d[l] = d[l] + f; + e[l] = 0.0; + } + + /* Sort eigenvalues and corresponding vectors. */ + + for (i = 0; i < n-1; i++) { + k = i; + Scalar p = d[i]; + for (j = i+1; j < n; j++) { + if (d[j] < p) { + k = j; + p = d[j]; + } + } + if (k != i) { + d[k] = d[i]; + d[i] = p; + for (j = 0; j < n; j++) { + p = V[j][i]; + V[j][i] = V[j][k]; + V[j][k] = p; + } + } + } + } + + /*! + * \brief Decomposes the symmetric matrix A_ij, into eigenvectors and eigenvalues + * \param[in] A_i: symmetric matrix to be decomposed + * \param[in,out] Eig_Vec: stores the eigenvectors + * \param[in,out] Eig_Val: stores the eigenvalues + * \param[in] n: order of matrix A_ij + * \param[in,out] e: work vector + */ + template + static void EigenDecomposition(const Mat& A_ij, Mat& Eig_Vec, Vec& Eig_Val, int n, W& e) { + for (int iDim = 0; iDim < n; iDim++){ + e[iDim] = 0.0; + for (int jDim = 0; jDim < n; jDim++){ + Eig_Vec[iDim][jDim] = A_ij[iDim][jDim]; + } + } + tred2(Eig_Vec, Eig_Val, e, n); + tql2(Eig_Vec, Eig_Val, e, n); + } + + /*! + * \brief Recomposes the eigenvectors and eigenvalues into a matrix + * \param[out] A_ij: recomposed matrix + * \param[in] Eig_Vec: eigenvectors + * \param[in] Eig_Val: eigenvalues + * \param[in] n: order of matrix A_ij + */ + template + static void EigenRecomposition(Mat& A_ij, const Mat& Eig_Vec, const Vec& Eig_Val, int n) { + for (int i = 0; i < n; i++){ + for (int j = 0; j < n; j++){ + A_ij[i][j] = 0.0; + for (int k = 0; k < n; k++) + A_ij[i][j] += Eig_Vec[i][k] * Eig_Val[k] * Eig_Vec[j][k]; + } + } + } + private: #if !(defined(HAVE_LIBXSMM) || defined(HAVE_BLAS) || defined(HAVE_MKL)) || (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 15c91d7f9b00..d227b2019484 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -37,6 +37,7 @@ #include "../fluid/CNEMOGas.hpp" #include "../../include/fluid/CMutationTCLib.hpp" #include "../../include/fluid/CUserDefinedTCLib.hpp" +#include "../../../Common/include/linear_algebra/blas_structure.hpp" using namespace std; @@ -63,10 +64,8 @@ class CNumerics { 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. */ - **delta, /*!< \brief Identity matrix. */ - **delta3; /*!< \brief 3 row Identity matrix. */ + 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 *Diffusion_Coeff_i, /*!< \brief Species diffusion coefficients at point i. */ *Diffusion_Coeff_j; /*!< \brief Species diffusion coefficients at point j. */ @@ -215,19 +214,15 @@ class CNumerics { su2double *l, *m; - su2double **MeanReynoldsStress; /*!< \brief Mean Reynolds stress tensor */ - su2double **MeanPerturbedRSM; /*!< \brief Perturbed Reynolds stress tensor */ - bool using_uq, /*!< \brief Flag for UQ methodology */ - nemo; /*!< \brief Flag for NEMO problems */ + 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 */ unsigned short Eig_Val_Comp; /*!< \brief Component towards which perturbation is perfromed */ su2double uq_delta_b; /*!< \brief Magnitude of perturbation */ su2double uq_urlx; /*!< \brief Under-relaxation factor for numerical stability */ bool uq_permute; /*!< \brief Flag for eigenvector permutation */ - /* Supporting data structures for the eigenspace perturbation for UQ methodology */ - su2double **A_ij, **newA_ij, **Eig_Vec, **New_Eig_Vec, **Corners; - su2double *Eig_Val, *Barycentric_Coord, *New_Coord; + bool nemo; /*!< \brief Flag for NEMO problems */ public: /*! @@ -275,17 +270,6 @@ class CNumerics { */ virtual ~CNumerics(void); - /*! - * \brief Compute the determinant of a 3 by 3 matrix. - * \param[in] val_matrix 3 by 3 matrix. - * \return Determinant of the matrix - */ - inline static su2double Determinant_3x3(su2double A00, su2double A01, su2double A02, - su2double A10, su2double A11, su2double A12, - su2double A20, su2double A21, su2double A22) { - return A00*(A11*A22-A12*A21) - A01*(A10*A22-A12*A20) + A02*(A10*A21-A11*A20); - } - /*! * \brief Set the time step. * \param[in] val_timestep - Value of the time step. @@ -465,15 +449,15 @@ class CNumerics { * \param[in] nDim - 2 or 3 * \param[out] rateofstrain - Rate of strain matrix * \param[in] velgrad - A velocity gradient matrix. - * \tparam TWOINDICES_1 - any type that supports the [][] interface - * \tparam TWOINDICES_2 - any type that supports the [][] interface + * \tparam Mat1 - any type that supports the [][] interface + * \tparam Mat2 - any type that supports the [][] interface */ - template - inline static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, TWOINDICES_1& rateofstrain, const TWOINDICES_2& velgrad){ + template + FORCEINLINE static void ComputeMeanRateOfStrainMatrix(size_t nDim, Mat1& rateofstrain, const Mat2& velgrad){ /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */ - if (nDim == 3){ + if (nDim == 3) { rateofstrain[0][0] = velgrad[0][0]; rateofstrain[1][1] = velgrad[1][1]; rateofstrain[2][2] = velgrad[2][2]; @@ -513,77 +497,168 @@ class CNumerics { * \param[in] density - Density * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2. - * \tparam TWOINDICES_1 - any type that supports the [][] interface - * \tparam TWOINDICES_2 - any type that supports the [][] interface - */ - template - inline static void ComputeStressTensor(unsigned short nDim, TWOINDICES_1& stress, const TWOINDICES_2& velgrad, - su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false){ - su2double divVel = 0; - for (unsigned short iDim = 0; iDim < nDim; iDim++){ + * \tparam Mat1 - any type that supports the [][] interface + * \tparam Mat2 - any type that supports the [][] interface + */ + template + FORCEINLINE static void ComputeStressTensor(size_t nDim, Mat1& stress, const Mat2& velgrad, + Scalar viscosity, Scalar density=0.0, + Scalar turb_ke=0.0, bool reynolds3x3=false){ + Scalar divVel = 0.0; + for (size_t iDim = 0; iDim < nDim; iDim++) { divVel += velgrad[iDim][iDim]; } - su2double pTerm = 2./3. * (divVel * viscosity + density * turb_ke); + Scalar pTerm = 2./3. * (divVel * viscosity + density * turb_ke); - for (unsigned short iDim = 0; iDim < nDim; iDim++){ - for (unsigned short jDim = 0; jDim < nDim; jDim++){ + for (size_t iDim = 0; iDim < nDim; iDim++){ + for (size_t jDim = 0; jDim < nDim; jDim++){ stress[iDim][jDim] = viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]); } stress[iDim][iDim] -= pTerm; } - if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix + if(reynolds3x3 && nDim==2) { // fill the third row and column of Reynolds stress matrix stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0; stress[2][2] = -pTerm; } - } /*! - * \brief Add a correction using a Quadratic Constitutive Relation + * \brief Add a correction using a Quadratic Constitutive Relation to the stress tensor. * - * This function requires that the stress tensor already be - * computed using \ref GetStressTensor - * - * See: Spalart, P. R., "Strategies for Turbulence Modelling and - * Simulation," International Journal of Heat and Fluid Flow, Vol. 21, - * 2000, pp. 252-263 + * See: Spalart, P. R., "Strategies for Turbulence Modelling and Simulation", + * International Journal of Heat and Fluid Flow, Vol. 21, 2000, pp. 252-263 * * \param[in] nDim: 2D or 3D. * \param[in] gradvel: Velocity gradients. * \param[in,out] tau: Shear stress tensor. */ - template - inline static void AddQCR(unsigned short nDim, const U& gradvel, V& tau) { + template + FORCEINLINE static void AddQCR(size_t nDim, const Mat1& gradvel, Mat2& tau) { + using Scalar = typename std::decay::type; - const su2double c_cr1= 0.3; - unsigned short iDim, jDim, kDim; + const Scalar c_cr1 = 0.3; /*--- Denominator Antisymmetric normalized rotation tensor ---*/ - su2double den_aux = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - den_aux += gradvel[iDim][jDim] * gradvel[iDim][jDim]; - den_aux = sqrt(max(den_aux,1E-10)); + Scalar factor = 0.0; + for (size_t iDim = 0; iDim < nDim; iDim++) + for (size_t jDim = 0; jDim < nDim; jDim++) + factor += gradvel[iDim][jDim] * gradvel[iDim][jDim]; + factor = 1.0 / sqrt(max(factor,1E-10)); /*--- Adding the QCR contribution ---*/ - su2double tauQCR[MAXNDIM][MAXNDIM] = {{0.0}}; + Scalar tauQCR[MAXNDIM][MAXNDIM] = {{0.0}}; - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim = 0; jDim < nDim; jDim++){ - for (kDim = 0; kDim < nDim; kDim++){ - su2double O_ik = (gradvel[iDim][kDim] - gradvel[kDim][iDim])/ den_aux; - su2double O_jk = (gradvel[jDim][kDim] - gradvel[kDim][jDim])/ den_aux; - tauQCR[iDim][jDim] = c_cr1 * ((O_ik * tau[jDim][kDim]) + (O_jk * tau[iDim][kDim])); + for (size_t iDim = 0; iDim < nDim; iDim++){ + for (size_t jDim = 0; jDim < nDim; jDim++){ + for (size_t kDim = 0; kDim < nDim; kDim++){ + Scalar O_ik = (gradvel[iDim][kDim] - gradvel[kDim][iDim]) * factor; + Scalar O_jk = (gradvel[jDim][kDim] - gradvel[kDim][jDim]) * factor; + tauQCR[iDim][jDim] += O_ik * tau[jDim][kDim] + O_jk * tau[iDim][kDim]; } } } - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - tau[iDim][jDim] -= tauQCR[iDim][jDim]; + for (size_t iDim = 0; iDim < nDim; iDim++) + for (size_t jDim = 0; jDim < nDim; jDim++) + tau[iDim][jDim] -= c_cr1 * tauQCR[iDim][jDim]; + } + + /*! + * \brief Perturb the Reynolds stress tensor based on parameters. + * \param[in] nDim - Dimension of the flow problem, 2 or 3. + * \param[in] uq_eigval_comp - Component 1C 2C 3C. + * \param[in] uq_permute - Whether to swap order of eigen vectors. + * \param[in] uq_delta_b - Delta_b parameter. + * \param[in] uq_urlx - Relaxation factor. + * \param[in] velgrad - A velocity gradient matrix. + * \param[in] density - Density. + * \param[in] viscosity - Eddy viscosity. + * \param[in] turb_ke: Turbulent kinetic energy. + * \param[out] MeanPerturbedRSM - Perturbed stress tensor. + */ + template + NEVERINLINE static void ComputePerturbedRSM(size_t nDim, size_t uq_eigval_comp, bool uq_permute, su2double uq_delta_b, + su2double uq_urlx, const Mat1& velgrad, Scalar density, + Scalar viscosity, Scalar turb_ke, Mat2& MeanPerturbedRSM) { + Scalar MeanReynoldsStress[3][3]; + ComputeStressTensor(nDim, MeanReynoldsStress, velgrad, viscosity, density, turb_ke, true); + for (size_t iDim = 0; iDim < 3; iDim++) + for (size_t jDim = 0; jDim < 3; jDim++) + MeanReynoldsStress[iDim][jDim] /= -density; + + /* --- Calculate anisotropic part of Reynolds Stress tensor --- */ + + Scalar A_ij[3][3]; + for (size_t iDim = 0; iDim < 3; iDim++) { + for (size_t jDim = 0; jDim < 3; jDim++) { + A_ij[iDim][jDim] = .5 * MeanReynoldsStress[iDim][jDim] / turb_ke; + } + A_ij[iDim][iDim] -= 1.0/3.0; + } + + /* --- Get ordered eigenvectors and eigenvalues of A_ij --- */ + + Scalar work[3], Eig_Vec[3][3], Eig_Val[3]; + CBlasStructure::EigenDecomposition(A_ij, Eig_Vec, Eig_Val, 3, work); + + /* compute convex combination coefficients */ + Scalar c1c = Eig_Val[2] - Eig_Val[1]; + Scalar c2c = 2.0 * (Eig_Val[1] - Eig_Val[0]); + Scalar c3c = 3.0 * Eig_Val[0] + 1.0; + + /* define barycentric traingle corner points */ + Scalar Corners[3][2]; + Corners[0][0] = 1.0; + Corners[0][1] = 0.0; + Corners[1][0] = 0.0; + Corners[1][1] = 0.0; + Corners[2][0] = 0.5; + Corners[2][1] = 0.866025; + + /* define barycentric coordinates */ + Scalar Barycentric_Coord[2]; + Barycentric_Coord[0] = Corners[0][0] * c1c + Corners[1][0] * c2c + Corners[2][0] * c3c; + Barycentric_Coord[1] = Corners[0][1] * c1c + Corners[1][1] * c2c + Corners[2][1] * c3c; + + /* component 1C, 2C, 3C, converted to index of the "corners" */ + Scalar New_Coord[2]; + New_Coord[0] = Corners[uq_eigval_comp-1][0]; + New_Coord[1] = Corners[uq_eigval_comp-1][1]; + + /* calculate perturbed barycentric coordinates */ + Barycentric_Coord[0] = Barycentric_Coord[0] + (uq_delta_b) * (New_Coord[0] - Barycentric_Coord[0]); + Barycentric_Coord[1] = Barycentric_Coord[1] + (uq_delta_b) * (New_Coord[1] - Barycentric_Coord[1]); + + /* rebuild c1c,c2c,c3c based on perturbed barycentric coordinates */ + c3c = Barycentric_Coord[1] / Corners[2][1]; + c1c = Barycentric_Coord[0] - Corners[2][0] * c3c; + c2c = 1 - c1c - c3c; + + /* build new anisotropy eigenvalues */ + Eig_Val[0] = (c3c - 1) / 3.0; + Eig_Val[1] = 0.5 *c2c + Eig_Val[0]; + Eig_Val[2] = c1c + Eig_Val[1]; + + /* permute eigenvectors if required */ + if (uq_permute) { + for (size_t jDim = 0; jDim < 3; jDim++) + swap(Eig_Vec[0][jDim], Eig_Vec[2][jDim]); + } + + CBlasStructure::EigenRecomposition(A_ij, Eig_Vec, Eig_Val, 3); + + /* compute perturbed Reynolds stress matrix; using under-relaxation factor (uq_urlx)*/ + for (size_t iDim = 0; iDim < 3; iDim++) { + for (size_t jDim = 0; jDim < 3; jDim++) { + auto delta_ij = (jDim==iDim)? 1.0 : 0.0; + MeanPerturbedRSM[iDim][jDim] = 2.0 * turb_ke * (A_ij[iDim][jDim] + 1.0/3.0 * delta_ij); + MeanPerturbedRSM[iDim][jDim] = MeanReynoldsStress[iDim][jDim] + + uq_urlx*(MeanPerturbedRSM[iDim][jDim] - MeanReynoldsStress[iDim][jDim]); + } + } } /*! @@ -1473,73 +1548,37 @@ class CNumerics { */ inline void SetRadVarSource(const su2double *val_radvar_source) { RadVar_Source = val_radvar_source; } - /*! - * \brief Decomposes the symmetric matrix A_ij, into eigenvectors and eigenvalues - * \param[in] A_i: symmetric matrix to be decomposed - * \param[in] Eig_Vec: strores the eigenvectors - * \param[in] Eig_Val: stores the eigenvalues - * \param[in] n: order of matrix A_ij - */ - static void EigenDecomposition(su2double **A_ij, su2double **Eig_Vec, su2double *Eig_Val, unsigned short n); - - /*! - * \brief Recomposes the eigenvectors and eigenvalues into a matrix - * \param[in] A_ij: recomposed matrix - * \param[in] Eig_Vec: eigenvectors - * \param[in] Eig_Val: eigenvalues - * \param[in] n: order of matrix A_ij - */ - static void EigenRecomposition(su2double **A_ij, su2double **Eig_Vec, const su2double *Eig_Val, unsigned short n); - - /*! - * \brief tred2 - * \param[in] V: matrix that needs to be decomposed - * \param[in] d: holds eigenvalues - * \param[in] e: supplemental data structure - * \param[in] n: order of matrix V - */ - static void tred2(su2double **V, su2double *d, su2double *e, unsigned short n); - - /*! - * \brief tql2 - * \param[in] V: matrix that will hold the eigenvectors - * \param[in] d: array that will hold the ordered eigenvalues - * \param[in] e: supplemental data structure - * \param[in] n: order of matrix V - */ - static void tql2(su2double **V, su2double *d, su2double *e, unsigned short n); - /*! * \brief Set the pressure derivatives. - * \param[in] val_dPdU_i - pressure derivatives at i. + * \param[in] val_dPdU_i - pressure derivatives at i. * \param[in] val_dPdU_j - pressure derivatives at j. */ virtual inline void SetdPdU(su2double *val_dPdU_i, su2double *val_dPdU_j) { } /*! * \brief Set the temperature derivatives. - * \param[in] val_dTdU_i - temperature derivatives at i. + * \param[in] val_dTdU_i - temperature derivatives at i. * \param[in] val_dTdU_j - temperature derivatives at j. */ virtual inline void SetdTdU(su2double *val_dTdU_i, su2double *val_dTdU_j) { } /*! * \brief Set the vib-el temperture derivatives. - * \param[in] val_dTvedU_i - t_ve derivatives at i. + * \param[in] val_dTvedU_i - t_ve derivatives at i. * \param[in] val_dTvedU_j - t_ve derivatives at j. */ virtual inline void SetdTvedU(su2double *val_dTvedU_i, su2double *val_dTvedU_j) { } /*! * \brief Set the vib-elec energy. - * \param[in] val_Eve_i - vib-el energy at i. + * \param[in] val_Eve_i - vib-el energy at i. * \param[in] val_Eve_j - vib-el energy at j. */ virtual inline void SetEve(su2double *val_Eve_i, su2double *val_Eve_j) { } /*! * \brief Set the vib-elec specific heat. - * \param[in] val_Cvve_i - Cvve at i. + * \param[in] val_Cvve_i - Cvve at i. * \param[in] val_Cvve_j - Cvve at j. */ virtual inline void SetCvve(su2double *val_Cvve_i, su2double *val_Cvve_j) { } diff --git a/SU2_CFD/include/numerics/flow/convection/centered.hpp b/SU2_CFD/include/numerics/flow/convection/centered.hpp index af6024d3cb89..59ac6e83003b 100644 --- a/SU2_CFD/include/numerics/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics/flow/convection/centered.hpp @@ -30,196 +30,6 @@ #include "../../CNumerics.hpp" -/*! - * \class CCentBase_Flow - * \brief Intermediate class to define centered schemes. - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CCentBase_Flow : public CNumerics { - -protected: - unsigned short iDim, iVar, jVar; /*!< \brief Iteration on dimension and variables. */ - bool dynamic_grid; /*!< \brief Consider grid movement. */ - bool implicit; /*!< \brief Implicit calculation (compute Jacobians). */ - su2double fix_factor; /*!< \brief Fix factor for dissipation Jacobians (more diagonal dominance). */ - - su2double Velocity_i[MAXNDIM] = {0.0}; /*!< \brief Velocity at node i. */ - su2double Velocity_j[MAXNDIM] = {0.0}; /*!< \brief Velocity at node j. */ - su2double MeanVelocity[MAXNDIM] = {0.0}; /*!< \brief Mean velocity. */ - su2double ProjVelocity_i, ProjVelocity_j; /*!< \brief Velocities in the face normal direction. */ - su2double sq_vel_i, sq_vel_j; /*!< \brief Squared norm of the velocity vectors. */ - su2double Energy_i, Energy_j, MeanEnergy; /*!< \brief Energy at nodes i and j and mean. */ - su2double MeanDensity, MeanPressure, MeanEnthalpy; /*!< \brief Mean density, pressure, and enthalpy. */ - su2double *ProjFlux = nullptr; /*!< \brief "The" flux. */ - - su2double *Diff_U = nullptr, *Diff_Lapl = nullptr; /*!< \brief Differences of conservatives and undiv. Laplacians. */ - su2double Local_Lambda_i, Local_Lambda_j, MeanLambda; /*!< \brief Local eingenvalues. */ - su2double Param_p, Phi_i, Phi_j, StretchingFactor; /*!< \brief Streching parameters. */ - su2double cte_0, cte_1; /*!< \brief Constants for the scalar dissipation Jacobian. */ - - su2double ProjGridVel; /*!< \brief Projected grid velocity. */ - - su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */ - su2double** Jacobian_j = nullptr; /*!< \brief The Jacobian w.r.t. point j after computation. */ - - /*! - * \brief Hook method for derived classes to define preaccumulated variables, optional to implement. - * \return true if any variable was set as preacc. input, in which case the residual will be output. - */ - virtual bool SetPreaccInVars(void) {return false;} - - /*! - * \brief Derived classes must implement this method, called in ComputeResidual after inviscid part. - * \param[in,out] val_residual - Pointer to the convective flux contribution to the residual. - * \param[in,out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[in,out] val_Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - */ - virtual void DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) = 0; - - /*! - * \brief Add the contribution of a scalar dissipation term to the Jacobians. - * \param[in,out] val_Jacobian_i - Jacobian of the numerical method at node i. - * \param[in,out] val_Jacobian_j - Jacobian of the numerical method at node j. - */ - void ScalarDissipationJacobian(su2double **val_Jacobian_i, su2double **val_Jacobian_j); - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimension of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CCentBase_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - - /*! - * \brief Destructor of the class. - */ - ~CCentBase_Flow(void) override; - - /*! - * \brief Compute the flow residual using a centered method with artificial dissipation. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) final; - -}; - -/*! - * \class CCentLax_Flow - * \brief Class for computing the Lax-Friedrich centered scheme. - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CCentLax_Flow final : public CCentBase_Flow { -private: - su2double Param_Kappa_0; /*!< \brief Artificial dissipation parameter. */ - su2double sc0; /*!< \brief Streching parameter. */ - su2double Epsilon_0; /*!< \brief Artificial dissipation coefficient. */ - - /*! - * \brief Lax-Friedrich first order dissipation term. - * \param[in,out] val_residual - Pointer to the convective flux contribution to the residual. - * \param[in,out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[in,out] val_Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - */ - void DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) override; - - /*! - * \brief Set input variables for AD preaccumulation. - * \return true, as we will define inputs. - */ - bool SetPreaccInVars(void) override; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimension of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CCentLax_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - -}; - -/*! - * \class CCentJST_KE_Flow - * \brief Class for centered scheme - JST_KE (no 4th dissipation order term). - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CCentJST_KE_Flow final : public CCentBase_Flow { - -private: - su2double Param_Kappa_2; /*!< \brief Artificial dissipation parameter. */ - su2double sc2; /*!< \brief Streching parameter. */ - su2double Epsilon_2; /*!< \brief Artificial dissipation coefficient. */ - - /*! - * \brief JST_KE second order dissipation term. - * \param[in,out] val_residual - Pointer to the convective flux contribution to the residual. - * \param[in,out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[in,out] val_Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - */ - void DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) override; - - /*! - * \brief Set input variables for AD preaccumulation. - * \return true, as we will define inputs. - */ - bool SetPreaccInVars(void) override; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimension of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CCentJST_KE_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - -}; - -/*! - * \class CCentJST_Flow - * \brief Class for centered scheme - JST. - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CCentJST_Flow final : public CCentBase_Flow { - -private: - su2double Param_Kappa_2, Param_Kappa_4; /*!< \brief Artificial dissipation parameters. */ - su2double sc2, sc4; /*!< \brief Streching parameters. */ - su2double Epsilon_2, Epsilon_4; /*!< \brief Artificial dissipation coefficients. */ - - /*! - * \brief JST second and forth order dissipation terms. - * \param[in,out] val_residual - Pointer to the convective flux contribution to the residual. - * \param[in,out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[in,out] val_Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - */ - void DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) override; - - /*! - * \brief Set input variables for AD preaccumulation. - * \return true, as we will define inputs. - */ - bool SetPreaccInVars(void) override; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimension of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CCentJST_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); - -}; - /*! * \class CCentLaxInc_Flow * \brief Class for computing the Lax-Friedrich centered scheme (modified with incompressible preconditioning). diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index e36b8a2d01e2..75e722c4cc39 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -162,20 +162,6 @@ class CAvgGrad_Base : public CNumerics { su2double val_dist_ij_2, const unsigned short val_nPrimVar); - /*! - * \brief Initialize the Reynolds Stress Matrix - * \param[in] turb_ke turbulent kinetic energy of node - */ - void SetReynoldsStressMatrix(su2double turb_ke); - - /*! - * \brief Perturb the Reynolds stress tensor based on parameters - * \param[in] turb_ke: turbulent kinetic energy of the noce - * \param[in] Eig_Val_Comp: Defines type of eigenspace perturbation - * \param[in] beta_delta: Defines the amount of eigenvalue perturbation - */ - void SetPerturbedRSM(su2double turb_ke, const CConfig* config); - public: /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 70327066a197..44cb621479ec 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -322,21 +322,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics { bool sustaining_terms; /*! - * \brief Initialize the Reynolds Stress Matrix - * \param[in] turb_ke turbulent kinetic energy of node + * \brief A virtual member. Get strain magnitude based on perturbed reynolds stress matrix + * \param[in] turb_ke: turbulent kinetic energy of the node */ - void SetReynoldsStressMatrix(su2double turb_ke); - - /*! - * \brief Perturb the Reynolds stress tensor based on parameters - * \param[in] turb_ke: turbulent kinetic energy of the noce - * \param[in] config: config file - */ - void SetPerturbedRSM(su2double turb_ke, const CConfig* config); - /*! - * \brief A virtual member. Get strain magnitude based on perturbed reynolds stress matrix - * \param[in] turb_ke: turbulent kinetic energy of the node - */ void SetPerturbedStrainMag(su2double turb_ke); public: diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp index 83cd63c3b154..89242a55ef92 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp @@ -32,59 +32,109 @@ #include "flow/convection/centered.hpp" #include "flow/diffusion/viscous_fluxes.hpp" +namespace { + /*! - * \brief Generic factory implementation. + * \brief Upwind factory implementation for ideal gas. */ template -CNumericsSIMD* createNumerics(const CConfig& config, int iMesh) { +CNumericsSIMD* createUpwindIdealNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { + CNumericsSIMD* obj = nullptr; + switch (config.GetKind_Upwind_Flow()) { + case ROE: + obj = new CRoeScheme(config, iMesh, turbVars); + break; + } + return obj; +} + +/*! + * \brief Upwind factory implementation for real gas. + */ +template +CNumericsSIMD* createUpwindGeneralNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { + return nullptr; +} + +/*! + * \brief Centered factory implementation. + */ +template +CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { + CNumericsSIMD* obj = nullptr; + switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : LAX) { + case NO_CENTERED: + break; + case LAX: + obj = new CLaxScheme(config, iMesh, turbVars); + break; + case JST: + obj = new CJSTScheme(config, iMesh, turbVars); + break; + case JST_KE: + obj = new CJSTkeScheme(config, iMesh, turbVars); + break; + case JST_MAT: + obj = new CJSTmatScheme(config, iMesh, turbVars); + break; + } + return obj; +} + +/*! + * \brief Generic factory implementation. + */ +template +CNumericsSIMD* createNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) { CNumericsSIMD* obj = nullptr; + const bool ideal_gas = (config.GetKind_FluidModel() == STANDARD_AIR) || + (config.GetKind_FluidModel() == IDEAL_GAS); + switch (config.GetKind_ConvNumScheme_Flow()) { case SPACE_UPWIND: - switch (config.GetKind_Upwind_Flow()) { - case ROE: - obj = new CRoeScheme(config, iMesh); - break; + if (config.GetViscous()) { + if (ideal_gas) + obj = createUpwindIdealNumerics >(config, iMesh, turbVars); + else + obj = createUpwindGeneralNumerics >(config, iMesh, turbVars); + } + else { + if (ideal_gas) + obj = createUpwindIdealNumerics >(config, iMesh, turbVars); + else + obj = createUpwindGeneralNumerics >(config, iMesh, turbVars); } break; case SPACE_CENTERED: - switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : LAX) { - case NO_CENTERED: - break; - case LAX: - obj = new CLaxScheme(config, iMesh); - break; - case JST: - obj = new CJSTScheme(config, iMesh); - break; - case JST_KE: - obj = new CJSTkeScheme(config, iMesh); - break; - case JST_MAT: - obj = new CJSTmatScheme(config, iMesh); - break; + if (config.GetViscous()) { + if (ideal_gas) + obj = createCenteredNumerics >(config, iMesh, turbVars); + else + obj = createCenteredNumerics >(config, iMesh, turbVars); + } + else { + obj = createCenteredNumerics >(config, iMesh, turbVars); } break; } + return obj; } +} // namespace + /*! * \brief This function instantiates both 2D and 3D versions of the implementation in * createNumerics, which in turn instantiates the class templates of the different * numerical methods. */ -CNumericsSIMD* CNumericsSIMD::CreateNumerics(const CConfig& config, int nDim, int iMesh) { +CNumericsSIMD* CNumericsSIMD::CreateNumerics(const CConfig& config, int nDim, int iMesh, const CVariable* turbVars) { if ((Double::Size < 4) && (SU2_MPI::GetRank() == MASTER_NODE)) { cout << "WARNING: SU2 was not compiled for an AVX-capable architecture." << endl; } - CNumericsSIMD* obj = nullptr; - if (config.GetViscous()) { - if (nDim == 2) obj = createNumerics >(config, iMesh); - if (nDim == 3) obj = createNumerics >(config, iMesh); - } else { - if (nDim == 2) obj = createNumerics >(config, iMesh); - if (nDim == 3) obj = createNumerics >(config, iMesh); - } - return obj; + if (nDim == 2) return createNumerics<2>(config, iMesh, turbVars); + if (nDim == 3) return createNumerics<3>(config, iMesh, turbVars); + + return nullptr; } diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp index ad91b747d476..78c0f7355713 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp @@ -93,7 +93,8 @@ class CNumericsSIMD { * \param[in] config - Problem definitions. * \param[in] nDim - 2D or 3D. * \param[in] iMesh - Grid index. + * \param[in] turbVars - Turbulence variables. */ - static CNumericsSIMD* CreateNumerics(const CConfig& config, int nDim, int iMesh); + static CNumericsSIMD* CreateNumerics(const CConfig& config, int nDim, int iMesh, const CVariable* turbVars = nullptr); }; diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 03c770ece008..b074546722ad 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -30,6 +30,7 @@ #include "../../CNumericsSIMD.hpp" #include "../../util.hpp" #include "../variables.hpp" +#include "../../../numerics/CNumerics.hpp" /*! * \brief Average gradients at i/j points. @@ -65,14 +66,12 @@ FORCEINLINE void correctGradient(const PrimitiveType& V, } /*! - * \brief Compute the stress tensor (using the total viscosity). + * \brief Compute the stress tensor. * \note Second viscosity term ignored. */ -template -FORCEINLINE MatrixDbl stressTensor(const PrimitiveType& V, +template +FORCEINLINE MatrixDbl stressTensor(Double viscosity, const MatrixDbl grad) { - Double viscosity = V.laminarVisc() + V.eddyVisc(); - /*--- Hydrostatic term. ---*/ Double velDiv = 0.0; for (size_t iDim = 0; iDim < nDim; ++iDim) { @@ -91,6 +90,33 @@ FORCEINLINE MatrixDbl stressTensor(const PrimitiveType& V, return tau; } +/*! + * \brief Add perturbed stress tensor. + * \note Not inlined because it is not easy to vectorize properly, due to tred2 and tql2. + */ +template +NEVERINLINE void addPerturbedRSM(const PrimitiveType& V, + const MatrixType& grad, + const Double& turb_ke, + MatrixDbl& tau, + Ts... uq_args) { + /*--- Handle SIMD dimensions 1 by 1. ---*/ + for (size_t k = 0; k < Double::Size; ++k) { + su2double velgrad[nDim][nDim]; + for (size_t iVar = 0; iVar < nDim; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + velgrad[iVar][iDim] = grad(iVar+1,iDim)[k]; + + su2double rsm[3][3]; + CNumerics::ComputePerturbedRSM(nDim, uq_args..., velgrad, V.density()[k], + V.eddyVisc()[k], turb_ke[k], rsm); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + tau(iDim,jDim)[k] -= V.density()[k] * rsm[iDim][jDim]; + } +} + /*! * \brief SA-QCR2000 modification of the stress tensor. */ diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp index 3080990452b5..0dec388892f0 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -62,14 +62,13 @@ class CNoViscousFlux : public CNumericsSIMD { }; /*! - * \class CCompressibleViscousFlux - * \brief Decorator class to add viscous fluxes (compressible flow, ideal gas). + * \class CCompressibleViscousFluxBase + * \brief Decorator class to add viscous fluxes (compressible flow). */ -template -class CCompressibleViscousFlux : public CNumericsSIMD { +template +class CCompressibleViscousFluxBase : public CNumericsSIMD { protected: static constexpr size_t nDim = NDIM; - static constexpr size_t nPrimVar = nDim+7; static constexpr size_t nPrimVarGrad = nDim+1; const su2double gamma; @@ -79,19 +78,33 @@ class CCompressibleViscousFlux : public CNumericsSIMD { const su2double cp; const bool correct; const bool useSA_QCR; + const bool uq; + const bool uq_permute; + const size_t uq_eigval_comp; + const su2double uq_delta_b; + const su2double uq_urlx; + + const CVariable* turbVars; /*! * \brief Constructor, initialize constants and booleans. */ template - CCompressibleViscousFlux(const CConfig& config, int iMesh, Ts&...) : + CCompressibleViscousFluxBase(const CConfig& config, int iMesh, + const CVariable* turbVars_, Ts&...) : gamma(config.GetGamma()), gasConst(config.GetGas_ConstantND()), prandtlLam(config.GetPrandtl_Lam()), prandtlTurb(config.GetPrandtl_Turb()), cp(gamma * gasConst / (gamma - 1)), correct(iMesh == MESH_0), - useSA_QCR(config.GetQCR()) { + useSA_QCR(config.GetQCR()), + uq(config.GetUsing_UQ()), + uq_permute(config.GetUQ_Permute()), + uq_eigval_comp(config.GetEig_Val_Comp()), + uq_delta_b(config.GetUQ_Delta_B()), + uq_urlx(config.GetUQ_URLX()), + turbVars(turbVars_) { } /*! @@ -114,7 +127,11 @@ class CCompressibleViscousFlux : public CNumericsSIMD { MatrixDbl& jac_i, MatrixDbl& jac_j) const { - static_assert(PrimVarType::nVar <= nPrimVar,""); + static_assert(PrimVarType::nVar <= Derived::nPrimVar,""); + + /*--- Pointer on which to call the "compile-time virtual" methods. ---*/ + + const auto derived = static_cast(this); const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); @@ -130,14 +147,18 @@ class CCompressibleViscousFlux : public CNumericsSIMD { auto avgGrad = averageGradient(iPoint, jPoint, gradient); if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad); - /// TODO: Uncertainty quantification (needs a way to access tke, maybe in ctor). - /*--- Stress and heat flux tensors. ---*/ - auto tau = stressTensor(avgV, avgGrad); + auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad); if(useSA_QCR) addQCR(avgGrad, tau); + if(uq) { + Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + + gatherVariables(jPoint, turbVars->GetSolution())); + addPerturbedRSM(avgV, avgGrad, turb_ke, tau, + uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); + } - Double cond = cp * (avgV.laminarVisc()/prandtlLam + avgV.eddyVisc()/prandtlTurb); + Double cond = derived->thermalConductivity(avgV); VectorDbl heatFlux; for (size_t iDim = 0; iDim < nDim; ++iDim) { heatFlux(iDim) = cond * avgGrad(0,iDim); @@ -157,23 +178,9 @@ class CCompressibleViscousFlux : public CNumericsSIMD { Double dist_ij = sqrt(dist2_ij); auto dtau = stressTensorJacobian(avgV, unitNormal, dist_ij); - Double contraction = 0.0; - for (size_t iDim = 0; iDim < nDim; ++iDim) { - contraction += dtau(iDim,0) * avgV.velocity(iDim); - } /*--- Energy flux Jacobian. ---*/ - VectorDbl dEdU; - Double vel2 = 0.5 * squaredNorm(avgV.velocity()); - Double phi = (gamma-1) / avgV.density(); - Double RdTdrho = phi*vel2 - avgV.pressure() / pow(avgV.density(),2); - Double condOnRd = cond / (gasConst * dist_ij); - - dEdU(0) = area * (condOnRd * RdTdrho - contraction); - for (size_t iDim = 0; iDim < nDim; ++iDim) { - dEdU(iDim+1) = area * (condOnRd*phi*avgV.velocity(iDim) + dtau(iDim,0)); - } - dEdU(nDim+1) = area * condOnRd * phi; + auto dEdU = derived->energyJacobian(avgV, dtau, cond, area, dist_ij, iPoint, jPoint, solution); /*--- Update momentum and energy terms ("symmetric" part). ---*/ for (size_t iDim = 0; iDim < nDim; ++iDim) { @@ -234,3 +241,125 @@ class CCompressibleViscousFlux : public CNumericsSIMD { viscousTerms(iEdge, iPoint, jPoint, avgV, V, solution_, vector_ij, geometry, args...); } }; + +/*! + * \class CCompressibleViscousFlux + * \brief Decorator class to add viscous fluxes (compressible flow, ideal gas). + */ +template +class CCompressibleViscousFlux : public CCompressibleViscousFluxBase > { +public: + static constexpr size_t nPrimVar = NDIM+7; + using Base = CCompressibleViscousFluxBase >; + using Base::gamma; + using Base::gasConst; + using Base::prandtlLam; + using Base::prandtlTurb; + using Base::cp; + + /*! + * \brief Constructor, initialize constants and booleans. + */ + template + CCompressibleViscousFlux(Ts&... args) : Base(args...) {} + + /*! + * \brief Compute the thermal conductivity. + */ + template + FORCEINLINE Double thermalConductivity(const PrimitiveType& V) const { + return cp * (V.laminarVisc()/prandtlLam + V.eddyVisc()/prandtlTurb); + } + + /*! + * \brief Compute Jacobian of the energy flux, except the part due to the work of viscous forces. + */ + template + FORCEINLINE VectorDbl energyJacobian(const PrimitiveType& V, + const MatrixDbl& dtau, + Double thermalCond, + Double area, + Double dist_ij, + Ts&... args) const { + Double vel2 = 0.5 * squaredNorm(V.velocity()); + Double phi = (gamma-1) / V.density(); + Double RdTdrho = phi*vel2 - V.pressure() / pow(V.density(),2); + Double condOnRd = thermalCond / (gasConst * dist_ij); + Double contraction = 0.0; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + contraction += dtau(iDim,0) * V.velocity(iDim); + } + VectorDbl dEdU; + dEdU(0) = area * (condOnRd * RdTdrho - contraction); + for (size_t iDim = 0; iDim < nDim; ++iDim) { + dEdU(iDim+1) = area * (dtau(iDim,0) - condOnRd*phi*V.velocity(iDim)); + } + dEdU(nDim+1) = area * condOnRd * phi; + + return dEdU; + } +}; + +/*! + * \class CGeneralCompressibleViscousFlux + * \brief Decorator class to add viscous fluxes (compressible flow, real gas). + */ +template +class CGeneralCompressibleViscousFlux : public CCompressibleViscousFluxBase > { +public: + static constexpr size_t nPrimVar = NDIM+9; + static constexpr size_t nSecVar = 4; + using Base = CCompressibleViscousFluxBase >; + using Base::prandtlTurb; + + /*! + * \brief Constructor, initialize constants and booleans. + */ + template + CGeneralCompressibleViscousFlux(Ts&... args) : Base(args...) {} + + /*! + * \brief Compute the thermal conductivity. + */ + template + FORCEINLINE Double thermalConductivity(const PrimitiveType& V) const { + return V.thermalCond() + V.cp()*V.eddyVisc()/prandtlTurb; + } + + /*! + * \brief Compute Jacobian of the energy flux, except the part due to the work of viscous forces. + */ + template + FORCEINLINE VectorDbl energyJacobian(const PrimitiveType& V, + const MatrixDbl& dtau, + Double thermalCond, + Double area, + Double dist_ij, + Int iPoint, + Int jPoint, + const VariableType& solution) const { + Double vel2 = squaredNorm(V.velocity()); + Double contraction = 0.0; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + contraction += dtau(iDim,0) * V.velocity(iDim); + } + + auto secVar_i = gatherVariables(iPoint, solution.GetSecondary()); + auto secVar_j = gatherVariables(jPoint, solution.GetSecondary()); + + Double dTdrho_e = 0.5 * (secVar_i(2) + secVar_j(2)); + Double dTde_rho = 0.5 * (secVar_i(3) + secVar_j(3)) / V.density(); + + Double condOnDist = thermalCond / dist_ij; + Double dTdrho = dTdrho_e + dTde_rho*(vel2-V.enthalpy()+V.pressure()/V.density()); + + VectorDbl dEdU; + dEdU(0) = area * (condOnDist * dTdrho - contraction); + for (size_t iDim = 0; iDim < nDim; ++iDim) { + dEdU(iDim+1) = area * (dtau(iDim,0) - condOnDist*dTde_rho*V.velocity(iDim)); + } + dEdU(nDim+1) = area * condOnDist * dTde_rho; + + return dEdU; + } +}; diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index bb202fd99432..30d356bbc86e 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -50,13 +50,17 @@ struct CCompressiblePrimitives { FORCEINLINE const Double& velocity(size_t iDim) const { return all(iDim+1); } FORCEINLINE const Double* velocity() const { return &velocity(0); } - /*--- Un-reconstructed variables (not allocated by default). ---*/ + /*--- Un-reconstructed variables. ---*/ FORCEINLINE Double& speedSound() { return all(nDim+4); } FORCEINLINE Double& laminarVisc() { return all(nDim+5); } FORCEINLINE Double& eddyVisc() { return all(nDim+6); } + FORCEINLINE Double& thermalCond() { return all(nDim+7); } + FORCEINLINE Double& cp() { return all(nDim+8); } FORCEINLINE const Double& speedSound() const { return all(nDim+4); } FORCEINLINE const Double& laminarVisc() const { return all(nDim+5); } FORCEINLINE const Double& eddyVisc() const { return all(nDim+6); } + FORCEINLINE const Double& thermalCond() const { return all(nDim+7); } + FORCEINLINE const Double& cp() const { return all(nDim+8); } }; /*! diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 15526a90f807..f7f6af0ce5e8 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -278,6 +278,13 @@ class CEulerSolver : public CFVMFlowSolverBase { */ void SetCoefficient_Gradients(CConfig *config) const; + /*! + * \brief Instantiate a SIMD numerics object. + * \param[in] solvers - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) final; + public: /*! * \brief Constructor of the class. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 7ea8b5e651d6..068fc49c7e11 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -232,13 +232,18 @@ class CFVMFlowSolverBase : public CSolver { /*! * \brief Method to compute convective and viscous residual contribution using vectorized numerics. */ - void EdgeFluxResidual(const CGeometry *geometry, const CConfig *config); + void EdgeFluxResidual(const CGeometry *geometry, const CSolver* const* solvers, const CConfig *config); /*! * \brief Sum the edge fluxes for each cell to populate the residual vector, only used on coarse grids. */ void SumEdgeFluxes(const CGeometry* geometry); + /*! + * \brief Instantiate a SIMD numerics object. + */ + inline virtual void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) {} + /*! * \brief Destructor. */ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 812bc67a6f03..07e1577044ce 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1287,7 +1287,12 @@ void CFVMFlowSolverBase::BC_Custom(CGeometry* geometry, CSolver** solver_c } template -void CFVMFlowSolverBase::EdgeFluxResidual(const CGeometry *geometry, const CConfig *config) { +void CFVMFlowSolverBase::EdgeFluxResidual(const CGeometry *geometry, + const CSolver* const* solvers, + const CConfig *config) { + if (!edgeNumerics) { + InstantiateEdgeNumerics(solvers, config); + } /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 3bd79dae0b16..1277fd98989d 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -197,6 +197,10 @@ class CSolver { assert(base_nodes!=nullptr && "CSolver::base_nodes was not set properly, see brief for CSolver::SetBaseClassPointerToNodes()"); return base_nodes; } + inline const CVariable* GetNodes() const { + assert(base_nodes!=nullptr && "CSolver::base_nodes was not set properly, see brief for CSolver::SetBaseClassPointerToNodes()"); + return base_nodes; + } /*! * \brief Helper function to define the type and number of variables per point for each communication type. diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 3eee69cbd2cc..8f3de95ab4e0 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -291,24 +291,29 @@ class CEulerVariable : public CVariable { inline su2double *GetPrimitive(unsigned long iPoint) final {return Primitive[iPoint]; } /*! - * \brief Get the primitive variables. + * \brief Get all the secondary variables. + */ + inline const MatrixType& GetSecondary() const {return Secondary; } + + /*! + * \brief Get the secondary variables. * \param[in] iVar - Index of the variable. - * \return Value of the primitive variable for the index iVar. + * \return Value of the secondary variable for the index iVar. */ inline su2double GetSecondary(unsigned long iPoint, unsigned long iVar) const final {return Secondary(iPoint,iVar); } /*! - * \brief Set the value of the primitive variables. + * \brief Set the value of the secondary variables. * \param[in] iVar - Index of the variable. * \param[in] iVar - Index of the variable. - * \return Set the value of the primitive variable for the index iVar. + * \return Set the value of the secondary variable for the index iVar. */ inline void SetSecondary(unsigned long iPoint, unsigned long iVar, su2double val_secondary) final {Secondary(iPoint,iVar) = val_secondary; } /*! - * \brief Set the value of the primitive variables. + * \brief Set the value of the secondary variables. * \param[in] val_prim - Primitive variables. - * \return Set the value of the primitive variable for the index iVar. + * \return Set the value of the secondary variable for the index iVar. */ inline void SetSecondary(unsigned long iPoint, const su2double *val_secondary) final { for (unsigned long iVar = 0; iVar < nSecondaryVar; iVar++) @@ -316,8 +321,8 @@ class CEulerVariable : public CVariable { } /*! - * \brief Get the primitive variables of the problem. - * \return Pointer to the primitive variable vector. + * \brief Get the secondary variables of the problem. + * \return Pointer to the secondary variable vector. */ inline su2double *GetSecondary(unsigned long iPoint) final { return Secondary[iPoint]; } diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 159b7363c8b3..b103e17f4709 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -447,7 +447,7 @@ class CVariable { * \brief Get the entire solution of the problem. * \return Reference to the solution matrix. */ - inline const MatrixType& GetSolution(void) { return Solution; } + inline const MatrixType& GetSolution(void) const { return Solution; } /*! * \brief Get the solution of the problem. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index a71332346686..2a3448f1b256 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1629,24 +1629,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol case SPACE_CENTERED : if (compressible) { - /*--- Compressible flow ---*/ - switch (config->GetKind_Centered_Flow()) { - case LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLax_Flow(nDim, nVar_Flow, config); break; - case JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJST_Flow(nDim, nVar_Flow, config); break; - case JST_KE : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJST_KE_Flow(nDim, nVar_Flow, config); break; - case JST_MAT : - if (!config->GetUseVectorization()) { - SU2_MPI::Error("JST with matrix dissipation requires USE_VECTORIZATION=YES.", CURRENT_FUNCTION); - } - break; - default: - SU2_OMP_MASTER - SU2_MPI::Error("Invalid centered scheme or not implemented.", CURRENT_FUNCTION); - break; - } - - for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][conv_term] = new CCentLax_Flow(nDim, nVar_Flow, config); + /*--- "conv_term" is not instantiated as all compressible centered schemes are vectorized. ---*/ /*--- Definition of the boundary condition method ---*/ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 4e37f06ce75b..ccfd713b03d3 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -40,9 +40,7 @@ CNumerics::CNumerics(void) { Proj_Flux_Tensor = nullptr; Flux_Tensor = nullptr; - tau = nullptr; - delta = nullptr; - delta3 = nullptr; + tau = nullptr; Vector = nullptr; @@ -82,18 +80,6 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, for (iDim = 0; iDim < nDim; iDim++) tau[iDim] = new su2double [nDim] (); - delta = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) { - delta[iDim] = new su2double [nDim] (); - delta[iDim][iDim] = 1.0; - } - - delta3 = new su2double* [3]; - for (iDim = 0; iDim < 3; iDim++) { - delta3[iDim] = new su2double [3] (); - delta3[iDim][iDim] = 1.0; - } - Proj_Flux_Tensor = new su2double [nVar] (); turb_ke_i = 0.0; @@ -108,40 +94,10 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, /* --- Initializing variables for the UQ methodology --- */ using_uq = config->GetUsing_UQ(); - if (using_uq){ - MeanReynoldsStress = new su2double* [3]; - MeanPerturbedRSM = new su2double* [3]; - A_ij = new su2double* [3]; - newA_ij = new su2double* [3]; - Eig_Vec = new su2double* [3]; - New_Eig_Vec = new su2double* [3]; - Corners = new su2double* [3]; - Eig_Val = new su2double [3]; - Barycentric_Coord = new su2double [2]; - New_Coord = new su2double [2]; - for (iDim = 0; iDim < 3; iDim++){ - MeanReynoldsStress[iDim] = new su2double [3]; - MeanPerturbedRSM[iDim] = new su2double [3]; - A_ij[iDim] = new su2double [3]; - newA_ij[iDim] = new su2double [3]; - Eig_Vec[iDim] = new su2double [3]; - New_Eig_Vec[iDim] = new su2double [3]; - Corners[iDim] = new su2double [2]; - Eig_Val[iDim] = 0; - } - Eig_Val_Comp = config->GetEig_Val_Comp(); - uq_delta_b = config->GetUQ_Delta_B(); - uq_urlx = config->GetUQ_URLX(); - uq_permute = config->GetUQ_Permute(); - - /* define barycentric triangle corner points */ - Corners[0][0] = 1.0; - Corners[0][1] = 0.0; - Corners[1][0] = 0.0; - Corners[1][1] = 0.0; - Corners[2][0] = 0.5; - Corners[2][1] = 0.866025; - } + Eig_Val_Comp = config->GetEig_Val_Comp(); + uq_delta_b = config->GetUQ_Delta_B(); + uq_urlx = config->GetUQ_URLX(); + uq_permute = config->GetUQ_Permute(); } @@ -165,44 +121,10 @@ CNumerics::~CNumerics(void) { delete [] tau; } - if (delta) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - delete [] delta[iDim]; - delete [] delta; - } - - if (delta3) { - for (unsigned short iDim = 0; iDim < 3; iDim++) - delete [] delta3[iDim]; - delete [] delta3; - } - delete [] Vector; delete [] l; delete [] m; - - if (using_uq) { - for (unsigned short iDim = 0; iDim < 3; iDim++){ - delete [] MeanReynoldsStress[iDim]; - delete [] MeanPerturbedRSM[iDim]; - delete [] A_ij[iDim]; - delete [] newA_ij[iDim]; - delete [] Eig_Vec[iDim]; - delete [] New_Eig_Vec[iDim]; - delete [] Corners[iDim]; - } - delete [] MeanReynoldsStress; - delete [] MeanPerturbedRSM; - delete [] A_ij; - delete [] newA_ij; - delete [] Eig_Vec; - delete [] New_Eig_Vec; - delete [] Corners; - delete [] Eig_Val; - delete [] Barycentric_Coord; - delete [] New_Coord; - } } void CNumerics::GetInviscidFlux(su2double val_density, const su2double *val_velocity, @@ -1757,7 +1679,7 @@ void CNumerics::GetPrimitive2Conservative (const su2double *val_Mean_PrimVar, } void CNumerics::CreateBasis(const su2double *val_Normal) { - + unsigned short iDim; su2double modm, modl; @@ -1880,348 +1802,3 @@ su2double CNumerics::GetRoe_Dissipation(const su2double Dissipation_i, } return Dissipation_ij; } - -void CNumerics::EigenDecomposition(su2double **A_ij, su2double **Eig_Vec, su2double *Eig_Val, unsigned short n){ - int iDim,jDim; - su2double *e = new su2double [n]; - for (iDim= 0; iDim< n; iDim++){ - e[iDim] = 0; - for (jDim = 0; jDim < n; jDim++){ - Eig_Vec[iDim][jDim] = A_ij[iDim][jDim]; - } - } - tred2(Eig_Vec, Eig_Val, e, n); - tql2(Eig_Vec, Eig_Val, e, n); - - delete [] e; -} - -void CNumerics::EigenRecomposition(su2double **A_ij, su2double **Eig_Vec, const su2double *Eig_Val, unsigned short n){ - unsigned short i,j,k; - su2double **tmp = new su2double* [n]; - su2double **deltaN = new su2double* [n]; - - for (i= 0; i< n; i++){ - tmp[i] = new su2double [n]; - deltaN[i] = new su2double [n]; - } - - for (i = 0; i < n; i++) { - for (j = 0; j < n; j++) { - if (i == j) deltaN[i][j] = 1.0; - else deltaN[i][j]=0.0; - } - } - - for (i= 0; i< n; i++){ - for (j = 0; j < n; j++){ - tmp[i][j] = 0.0; - for (k = 0; k < n; k++){ - tmp[i][j] += Eig_Vec[i][k] * Eig_Val[k] * deltaN[k][j]; - } - } - } - - for (i= 0; i< n; i++){ - for (j = 0; j < n; j++){ - A_ij[i][j] = 0.0; - for (k = 0; k < n; k++){ - A_ij[i][j] += tmp[i][k] * Eig_Vec[j][k]; - } - } - } - - for (i = 0; i < n; i++){ - delete [] tmp[i]; - delete [] deltaN[i]; - } - delete [] tmp; - delete [] deltaN; -} - -void CNumerics::tred2(su2double **V, su2double *d, su2double *e, unsigned short n) { -/* Author: - - * Original FORTRAN77 version by Smith, Boyle, Dongarra, Garbow, Ikebe, - * Klema, Moler. - * C++ version by Aashwin Mishra and Jayant Mukhopadhaya. - - * Reference: - - * Martin, Reinsch, Wilkinson, - * TRED2, - * Numerische Mathematik, - * Volume 11, pages 181-195, 1968. - - * James Wilkinson, Christian Reinsch, - * Handbook for Automatic Computation, - * Volume II, Linear Algebra, Part 2, - * Springer, 1971, - * ISBN: 0387054146, - * LC: QA251.W67. - - * Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, - * Yasuhiko Ikebe, Virginia Klema, Cleve Moler, - * Matrix Eigensystem Routines, EISPACK Guide, - * Lecture Notes in Computer Science, Volume 6, - * Springer Verlag, 1976, - * ISBN13: 978-3540075462, - * LC: QA193.M37 - -*/ - - unsigned short i,j,k; - - for (j = 0; j < n; j++) { - d[j] = V[n-1][j]; - } - - /* Householder reduction to tridiagonal form. */ - - for (i = n-1; i > 0; i--) { - - /* Scale to avoid under/overflow. */ - - su2double scale = 0.0; - su2double h = 0.0; - for (k = 0; k < i; k++) { - scale = scale + fabs(d[k]); - } - if (scale == 0.0) { - e[i] = d[i-1]; - for (j = 0; j < i; j++) { - d[j] = V[i-1][j]; - V[i][j] = 0.0; - V[j][i] = 0.0; - } - } - else { - - /* Generate Householder vector. */ - - for (k = 0; k < i; k++) { - d[k] /= scale; - h += d[k] * d[k]; - } - su2double f = d[i-1]; - su2double g = sqrt(h); - if (f > 0) { - g = -g; - } - e[i] = scale * g; - h = h - f * g; - d[i-1] = f - g; - for (j = 0; j < i; j++) { - e[j] = 0.0; - } - - /* Apply similarity transformation to remaining columns. */ - - for (j = 0; j < i; j++) { - f = d[j]; - V[j][i] = f; - g = e[j] + V[j][j] * f; - for (k = j+1; k <= i-1; k++) { - g += V[k][j] * d[k]; - e[k] += V[k][j] * f; - } - e[j] = g; - } - f = 0.0; - for (j = 0; j < i; j++) { - e[j] /= h; - f += e[j] * d[j]; - } - su2double hh = f / (h + h); - for (j = 0; j < i; j++) { - e[j] -= hh * d[j]; - } - for (j = 0; j < i; j++) { - f = d[j]; - g = e[j]; - for (k = j; k <= i-1; k++) { - V[k][j] -= (f * e[k] + g * d[k]); - } - d[j] = V[i-1][j]; - V[i][j] = 0.0; - } - } - d[i] = h; - } - - /* Accumulate transformations. */ - - for (i = 0; i < n-1; i++) { - V[n-1][i] = V[i][i]; - V[i][i] = 1.0; - su2double h = d[i+1]; - if (h != 0.0) { - for (k = 0; k <= i; k++) { - d[k] = V[k][i+1] / h; - } - for (j = 0; j <= i; j++) { - su2double g = 0.0; - for (k = 0; k <= i; k++) { - g += V[k][i+1] * V[k][j]; - } - for (k = 0; k <= i; k++) { - V[k][j] -= g * d[k]; - } - } - } - for (k = 0; k <= i; k++) { - V[k][i+1] = 0.0; - } - } - for (j = 0; j < n; j++) { - d[j] = V[n-1][j]; - V[n-1][j] = 0.0; - } - V[n-1][n-1] = 1.0; - e[0] = 0.0; -} - -void CNumerics::tql2(su2double **V, su2double *d, su2double *e, unsigned short n) { - -/* Author: - - * Original FORTRAN77 version by Smith, Boyle, Dongarra, Garbow, Ikebe, - * Klema, Moler. - * C++ version by Aashwin Mishra and Jayant Mukhopadhaya. - - * Reference: - - * Bowdler, Martin, Reinsch, Wilkinson, - * TQL2, - * Numerische Mathematik, - * Volume 11, pages 293-306, 1968. - - * James Wilkinson, Christian Reinsch, - * Handbook for Automatic Computation, - * Volume II, Linear Algebra, Part 2, - * Springer, 1971, - * ISBN: 0387054146, - * LC: QA251.W67. - - * Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, - * Yasuhiko Ikebe, Virginia Klema, Cleve Moler, - * Matrix Eigensystem Routines, EISPACK Guide, - * Lecture Notes in Computer Science, Volume 6, - * Springer Verlag, 1976, - * ISBN13: 978-3540075462, - * LC: QA193.M37 - -*/ - - int i,j,k,l; - for (i = 1; i < n; i++) { - e[i-1] = e[i]; - } - e[n-1] = 0.0; - - su2double f = 0.0; - su2double tst1 = 0.0; - su2double eps = pow(2.0,-52.0); - for (l = 0; l < n; l++) { - - /* Find small subdiagonal element */ - - tst1 = max(tst1,(fabs(d[l]) + fabs(e[l]))); - int m = l; - while (m < n) { - if (fabs(e[m]) <= eps*tst1) { - break; - } - m++; - } - - /* If m == l, d[l] is an eigenvalue, */ - /* otherwise, iterate. */ - - if (m > l) { - int iter = 0; - do { - iter = iter + 1; /* (Could check iteration count here.) */ - - /* Compute implicit shift */ - - su2double g = d[l]; - su2double p = (d[l+1] - g) / (2.0 * e[l]); - su2double r = sqrt(p*p+1.0); - if (p < 0) { - r = -r; - } - d[l] = e[l] / (p + r); - d[l+1] = e[l] * (p + r); - su2double dl1 = d[l+1]; - su2double h = g - d[l]; - for (i = l+2; i < n; i++) { - d[i] -= h; - } - f = f + h; - - /* Implicit QL transformation. */ - - p = d[m]; - su2double c = 1.0; - su2double c2 = c; - su2double c3 = c; - su2double el1 = e[l+1]; - su2double s = 0.0; - su2double s2 = 0.0; - for (i = m-1; i >= l; i--) { - c3 = c2; - c2 = c; - s2 = s; - g = c * e[i]; - h = c * p; - r = sqrt(p*p+e[i]*e[i]); - e[i+1] = s * r; - s = e[i] / r; - c = p / r; - p = c * d[i] - s * g; - d[i+1] = h + s * (c * g + s * d[i]); - - /* Accumulate transformation. */ - - for (k = 0; k < n; k++) { - h = V[k][i+1]; - V[k][i+1] = s * V[k][i] + c * h; - V[k][i] = c * V[k][i] - s * h; - } - } - p = -s * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; - - /* Check for convergence. */ - - } while (fabs(e[l]) > eps*tst1); - } - d[l] = d[l] + f; - e[l] = 0.0; - } - - /* Sort eigenvalues and corresponding vectors. */ - - for (i = 0; i < n-1; i++) { - k = i; - su2double p = d[i]; - for (j = i+1; j < n; j++) { - if (d[j] < p) { - k = j; - p = d[j]; - } - } - if (k != i) { - d[k] = d[i]; - d[i] = p; - for (j = 0; j < n; j++) { - p = V[j][i]; - V[j][i] = V[j][k]; - V[j][k] = p; - } - } - } -} - diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index 0e7fd1532da7..d95c9b88b8b5 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -28,316 +28,6 @@ #include "../../../../include/numerics/flow/convection/centered.hpp" #include "../../../../../Common/include/toolboxes/geometry_toolbox.hpp" -CCentBase_Flow::CCentBase_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : - CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Flow() == 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(); - fix_factor = config->GetCent_Jac_Fix_Factor(); - - Gamma = config->GetGamma(); - Gamma_Minus_One = Gamma - 1.0; - - /*--- Allocate required structures ---*/ - Diff_U = new su2double [nVar]; - Diff_Lapl = new su2double [nVar]; - ProjFlux = new su2double [nVar]; - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - } -} - -CCentBase_Flow::~CCentBase_Flow(void) { - delete [] Diff_U; - delete [] Diff_Lapl; - delete [] ProjFlux; - - if (Jacobian_i != nullptr) { - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - } - delete [] Jacobian_i; - delete [] Jacobian_j; - } -} - -CNumerics::ResidualType<> CCentBase_Flow::ComputeResidual(const CConfig* config) { - - implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - su2double U_i[5] = {0.0}, U_j[5] = {0.0}; - - bool preacc = SetPreaccInVars(); - - if (preacc) { - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(V_i, nDim+5); AD::SetPreaccIn(V_j, nDim+5); - AD::SetPreaccIn(Lambda_i, Lambda_j); - if (dynamic_grid) { - AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); - } - } - - /*--- Pressure, density, enthalpy, energy, and velocity at points i and j ---*/ - - Pressure_i = V_i[nDim+1]; Pressure_j = V_j[nDim+1]; - Density_i = V_i[nDim+2]; Density_j = V_j[nDim+2]; - Enthalpy_i = V_i[nDim+3]; Enthalpy_j = V_j[nDim+3]; - SoundSpeed_i = V_i[nDim+4]; SoundSpeed_j = V_j[nDim+4]; - Energy_i = Enthalpy_i - Pressure_i/Density_i; Energy_j = Enthalpy_j - Pressure_j/Density_j; - - sq_vel_i = 0.0; sq_vel_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - Velocity_j[iDim] = V_j[iDim+1]; - sq_vel_i += 0.5*Velocity_i[iDim]*Velocity_i[iDim]; - sq_vel_j += 0.5*Velocity_j[iDim]*Velocity_j[iDim]; - } - - /*--- Recompute conservative variables ---*/ - - U_i[0] = Density_i; U_j[0] = Density_j; - for (iDim = 0; iDim < nDim; iDim++) { - U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; - } - U_i[nDim+1] = Density_i*Energy_i; U_j[nDim+1] = Density_j*Energy_j; - - /*--- Compute mean values of the variables ---*/ - - MeanDensity = 0.5*(Density_i+Density_j); - MeanPressure = 0.5*(Pressure_i+Pressure_j); - MeanEnthalpy = 0.5*(Enthalpy_i+Enthalpy_j); - for (iDim = 0; iDim < nDim; iDim++) - MeanVelocity[iDim] = 0.5*(Velocity_i[iDim]+Velocity_j[iDim]); - MeanEnergy = 0.5*(Energy_i+Energy_j); - - /*--- Residual of the inviscid flux ---*/ - - GetInviscidProjFlux(&MeanDensity, MeanVelocity, &MeanPressure, &MeanEnthalpy, Normal, ProjFlux); - - /*--- Jacobians of the inviscid flux, scale = 0.5 because ProjFlux ~ 0.5*(fc_i+fc_j)*Normal ---*/ - - if (implicit) { - GetInviscidProjJac(MeanVelocity, &MeanEnergy, Normal, 0.5, Jacobian_i); - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_j[iVar][jVar] = Jacobian_i[iVar][jVar]; - } - - /*--- Adjustment due to grid motion ---*/ - - if (dynamic_grid) { - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - - for (iVar = 0; iVar < nVar; iVar++) { - ProjFlux[iVar] -= ProjGridVel * 0.5*(U_i[iVar] + U_j[iVar]); - if (implicit) { - Jacobian_i[iVar][iVar] -= 0.5*ProjGridVel; - Jacobian_j[iVar][iVar] -= 0.5*ProjGridVel; - } - } - } - - /*--- Compute the local spectral radius and the stretching factor ---*/ - - ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjVelocity_i += Velocity_i[iDim]*Normal[iDim]; - ProjVelocity_j += Velocity_j[iDim]*Normal[iDim]; - } - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Adjustment due to mesh motion ---*/ - - if (dynamic_grid) { - ProjVelocity_i -= ProjGridVel; - ProjVelocity_j -= ProjGridVel; - } - - /*--- Dissipation term ---*/ - - Local_Lambda_i = (fabs(ProjVelocity_i)+SoundSpeed_i*Area); - Local_Lambda_j = (fabs(ProjVelocity_j)+SoundSpeed_j*Area); - MeanLambda = 0.5*(Local_Lambda_i+Local_Lambda_j); - - Phi_i = pow(Lambda_i/(4.0*MeanLambda), Param_p); - Phi_j = pow(Lambda_j/(4.0*MeanLambda), Param_p); - StretchingFactor = 4.0*Phi_i*Phi_j/(Phi_i+Phi_j); - - /*--- Compute differences btw. conservative variables, with a correction for enthalpy ---*/ - - for (iVar = 0; iVar < nVar-1; iVar++) { - Diff_U[iVar] = U_i[iVar]-U_j[iVar]; - } - Diff_U[nVar-1] = Density_i*Enthalpy_i-Density_j*Enthalpy_j; - - DissipationTerm(ProjFlux, Jacobian_i, Jacobian_j); - - if (preacc) { - AD::SetPreaccOut(ProjFlux, nVar); - AD::EndPreacc(); - } - - return ResidualType<>(ProjFlux, Jacobian_i, Jacobian_j); - -} - -void CCentBase_Flow::ScalarDissipationJacobian(su2double **val_Jacobian_i, su2double **val_Jacobian_j) { - - /*--- n-1 diagonal entries ---*/ - - for (iVar = 0; iVar < (nVar-1); iVar++) { - val_Jacobian_i[iVar][iVar] += fix_factor*cte_0; - val_Jacobian_j[iVar][iVar] -= fix_factor*cte_1; - } - - /*--- Last row of Jacobian_i ---*/ - - val_Jacobian_i[nVar-1][0] += fix_factor*cte_0*Gamma_Minus_One*sq_vel_i; - for (iDim = 0; iDim < nDim; iDim++) - val_Jacobian_i[nVar-1][iDim+1] -= fix_factor*cte_0*Gamma_Minus_One*Velocity_i[iDim]; - val_Jacobian_i[nVar-1][nVar-1] += fix_factor*cte_0*Gamma; - - /*--- Last row of Jacobian_j ---*/ - - val_Jacobian_j[nVar-1][0] -= fix_factor*cte_1*Gamma_Minus_One*sq_vel_j; - for (iDim = 0; iDim < nDim; iDim++) - val_Jacobian_j[nVar-1][iDim+1] += fix_factor*cte_1*Gamma_Minus_One*Velocity_j[iDim]; - val_Jacobian_j[nVar-1][nVar-1] -= fix_factor*cte_1*Gamma; - -} - -CCentLax_Flow::CCentLax_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : - CCentBase_Flow(val_nDim, val_nVar, config) { - - /*--- Artifical dissipation parameters ---*/ - Param_p = 0.3; - Param_Kappa_0 = config->GetKappa_1st_Flow(); - -} - -void CCentLax_Flow::DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) { - - /*--- Compute dissipation coefficient ---*/ - - sc0 = 3.0*(su2double(Neighbor_i)+su2double(Neighbor_j))/(su2double(Neighbor_i)*su2double(Neighbor_j)); - Epsilon_0 = Param_Kappa_0*sc0*su2double(nDim)/3.0; - - /*--- Compute viscous part of the residual ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - val_residual[iVar] += Epsilon_0*Diff_U[iVar]*StretchingFactor*MeanLambda; - - /*--- Jacobian computation ---*/ - - if (implicit) { - - cte_0 = Epsilon_0*StretchingFactor*MeanLambda; - cte_1 = cte_0; - - ScalarDissipationJacobian(val_Jacobian_i, val_Jacobian_j); - } -} - -bool CCentLax_Flow::SetPreaccInVars(void) { - AD::StartPreacc(); - return true; -} - -CCentJST_KE_Flow::CCentJST_KE_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : - CCentBase_Flow(val_nDim, val_nVar, config) { - - /*--- Artifical dissipation parameters ---*/ - Param_p = 0.3; - Param_Kappa_2 = config->GetKappa_2nd_Flow(); - -} - -void CCentJST_KE_Flow::DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) { - - /*--- Compute dissipation coefficient ---*/ - - sc2 = 3.0*(su2double(Neighbor_i)+su2double(Neighbor_j))/(su2double(Neighbor_i)*su2double(Neighbor_j)); - Epsilon_2 = Param_Kappa_2*0.5*(Sensor_i+Sensor_j)*sc2; - - /*--- Compute viscous part of the residual ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - val_residual[iVar] += Epsilon_2*(Diff_U[iVar])*StretchingFactor*MeanLambda; - - /*--- Jacobian computation ---*/ - - if (implicit) { - - cte_0 = Epsilon_2*StretchingFactor*MeanLambda; - cte_1 = cte_0; - - ScalarDissipationJacobian(val_Jacobian_i, val_Jacobian_j); - } -} - -bool CCentJST_KE_Flow::SetPreaccInVars(void) { - AD::StartPreacc(); - AD::SetPreaccIn(Sensor_i); AD::SetPreaccIn(Sensor_j); - return true; -} - -CCentJST_Flow::CCentJST_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : - CCentBase_Flow(val_nDim, val_nVar, config) { - - /*--- Artifical dissipation parameters ---*/ - Param_p = 0.3; - Param_Kappa_2 = config->GetKappa_2nd_Flow(); - Param_Kappa_4 = config->GetKappa_4th_Flow(); - -} - -void CCentJST_Flow::DissipationTerm(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j) { - - /*--- Compute differences btw. Laplacians ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Diff_Lapl[iVar] = Und_Lapl_i[iVar]-Und_Lapl_j[iVar]; - } - - /*--- Compute dissipation coefficients ---*/ - - sc2 = 3.0*(su2double(Neighbor_i)+su2double(Neighbor_j))/(su2double(Neighbor_i)*su2double(Neighbor_j)); - sc4 = sc2*sc2/4.0; - - Epsilon_2 = Param_Kappa_2*0.5*(Sensor_i+Sensor_j)*sc2; - Epsilon_4 = max(0.0, Param_Kappa_4-Epsilon_2)*sc4; - - /*--- Compute viscous part of the residual ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - val_residual[iVar] += (Epsilon_2*Diff_U[iVar] - Epsilon_4*Diff_Lapl[iVar])*StretchingFactor*MeanLambda; - - /*--- Jacobian computation ---*/ - - if (implicit) { - - cte_0 = (Epsilon_2 + Epsilon_4*su2double(Neighbor_i+1))*StretchingFactor*MeanLambda; - cte_1 = (Epsilon_2 + Epsilon_4*su2double(Neighbor_j+1))*StretchingFactor*MeanLambda; - - ScalarDissipationJacobian(val_Jacobian_i, val_Jacobian_j); - } -} - -bool CCentJST_Flow::SetPreaccInVars(void) { - AD::StartPreacc(); - AD::SetPreaccIn(Sensor_i); AD::SetPreaccIn(Und_Lapl_i, nVar); - AD::SetPreaccIn(Sensor_j); AD::SetPreaccIn(Und_Lapl_j, nVar); - return true; -} - CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index d7ff4ab9e51b..570a0f595a6e 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -139,7 +139,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } else { // compute both parts in one step const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; - ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, 0.0); // TODO why ignore turb_ke? + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0)); // TODO why ignore turb_ke? } } @@ -183,116 +183,6 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal, tau[iDim][jDim] = tau[iDim][jDim]*(val_tau_wall/WallShearStress); } -void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){ - su2double meandensity = Mean_PrimVar[nDim+2]; - ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar+1, Mean_Eddy_Viscosity, meandensity, turb_ke, true); - for(unsigned short iDim=0; iDim<3; iDim++){ - for(unsigned short jDim=0; jDim<3; jDim++){ - MeanReynoldsStress[iDim][jDim] /= (-meandensity); - } - } -} - -void CAvgGrad_Base::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ - - unsigned short iDim,jDim; - - /* --- Calculate anisotropic part of Reynolds Stress tensor --- */ - - for (iDim = 0; iDim< 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - A_ij[iDim][jDim] = .5 * MeanReynoldsStress[iDim][jDim] / turb_ke - delta3[iDim][jDim] / 3.0; - Eig_Vec[iDim][jDim] = A_ij[iDim][jDim]; - } - } - - /* --- Get ordered eigenvectors and eigenvalues of A_ij --- */ - - EigenDecomposition(A_ij, Eig_Vec, Eig_Val, 3); - - /* compute convex combination coefficients */ - su2double c1c = Eig_Val[2] - Eig_Val[1]; - su2double c2c = 2.0 * (Eig_Val[1] - Eig_Val[0]); - su2double c3c = 3.0 * Eig_Val[0] + 1.0; - - /* define barycentric traingle corner points */ - Corners[0][0] = 1.0; - Corners[0][1] = 0.0; - Corners[1][0] = 0.0; - Corners[1][1] = 0.0; - Corners[2][0] = 0.5; - Corners[2][1] = 0.866025; - - /* define barycentric coordinates */ - Barycentric_Coord[0] = Corners[0][0] * c1c + Corners[1][0] * c2c + Corners[2][0] * c3c; - Barycentric_Coord[1] = Corners[0][1] * c1c + Corners[1][1] * c2c + Corners[2][1] * c3c; - - if (Eig_Val_Comp == 1) { - /* 1C turbulence */ - New_Coord[0] = Corners[0][0]; - New_Coord[1] = Corners[0][1]; - } - else if (Eig_Val_Comp== 2) { - /* 2C turbulence */ - New_Coord[0] = Corners[1][0]; - New_Coord[1] = Corners[1][1]; - } - else if (Eig_Val_Comp == 3) { - /* 3C turbulence */ - New_Coord[0] = Corners[2][0]; - New_Coord[1] = Corners[2][1]; - } - else { - /* 2C turbulence */ - New_Coord[0] = Corners[1][0]; - New_Coord[1] = Corners[1][1]; - } - - /* calculate perturbed barycentric coordinates */ - Barycentric_Coord[0] = Barycentric_Coord[0] + (uq_delta_b) * (New_Coord[0] - Barycentric_Coord[0]); - Barycentric_Coord[1] = Barycentric_Coord[1] + (uq_delta_b) * (New_Coord[1] - Barycentric_Coord[1]); - - /* rebuild c1c,c2c,c3c based on perturbed barycentric coordinates */ - c3c = Barycentric_Coord[1] / Corners[2][1]; - c1c = Barycentric_Coord[0] - Corners[2][0] * c3c; - c2c = 1 - c1c - c3c; - - /* build new anisotropy eigenvalues */ - Eig_Val[0] = (c3c - 1) / 3.0; - Eig_Val[1] = 0.5 *c2c + Eig_Val[0]; - Eig_Val[2] = c1c + Eig_Val[1]; - - /* permute eigenvectors if required */ - if (uq_permute) { - for (iDim=0; iDim<3; iDim++) { - for (jDim=0; jDim<3; jDim++) { - New_Eig_Vec[iDim][jDim] = Eig_Vec[2-iDim][jDim]; - } - } - } - - else { - for (iDim=0; iDim<3; iDim++) { - for (jDim=0; jDim<3; jDim++) { - New_Eig_Vec[iDim][jDim] = Eig_Vec[iDim][jDim]; - } - } - } - - EigenRecomposition(newA_ij, New_Eig_Vec, Eig_Val, 3); - - /* compute perturbed Reynolds stress matrix; use under-relaxation factor (uq_urlx)*/ - for (iDim = 0; iDim< 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - MeanPerturbedRSM[iDim][jDim] = 2.0 * turb_ke * (newA_ij[iDim][jDim] + 1.0/3.0 * delta3[iDim][jDim]); - MeanPerturbedRSM[iDim][jDim] = MeanReynoldsStress[iDim][jDim] + - uq_urlx*(MeanPerturbedRSM[iDim][jDim] - MeanReynoldsStress[iDim][jDim]); - } - } - -} - - void CAvgGrad_Base::SetTauJacobian(const su2double *val_Mean_PrimVar, const su2double val_laminar_viscosity, const su2double val_eddy_viscosity, @@ -556,8 +446,9 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) /* --- If using UQ methodology, set Reynolds Stress tensor and perform perturbation--- */ if (using_uq){ - SetReynoldsStressMatrix(Mean_turb_ke); - SetPerturbedRSM(Mean_turb_ke, config); + ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, + Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, + Mean_turb_ke, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ @@ -1042,17 +933,27 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c dist_ij_2, nDim+1); } + /*--- Wall shear stress values (wall functions) ---*/ + + if (TauWall_i > 0.0 && TauWall_j > 0.0) Mean_TauWall = 0.5*(TauWall_i + TauWall_j); + else if (TauWall_i > 0.0) Mean_TauWall = TauWall_i; + else if (TauWall_j > 0.0) Mean_TauWall = TauWall_j; + else Mean_TauWall = -1.0; + /* --- If using UQ methodology, set Reynolds Stress tensor and perform perturbation--- */ if (using_uq){ - SetReynoldsStressMatrix(Mean_turb_ke); - SetPerturbedRSM(Mean_turb_ke, config); + ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, + Mean_GradPrimVar+1, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, + Mean_turb_ke, MeanPerturbedRSM); } /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + if (config->GetQCR()) AddQCR(nDim, &Mean_GradPrimVar[1], tau); + if (Mean_TauWall > 0) AddTauWall(Normal, Mean_TauWall); SetHeatFluxVector(Mean_GradPrimVar, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, Mean_Thermal_Conductivity, Mean_Cp); diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp index 73867bef1406..8940306edf33 100644 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp @@ -835,8 +835,9 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi /* if using UQ methodolgy, calculate production using perturbed Reynolds stress matrix */ if (using_uq){ - SetReynoldsStressMatrix(TurbVar_i[0]); - SetPerturbedRSM(TurbVar_i[0], config); + ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, + PrimVar_Grad_i+1, Density_i, Eddy_Viscosity_i, + TurbVar_i[0], MeanPerturbedRSM); SetPerturbedStrainMag(TurbVar_i[0]); pk = Eddy_Viscosity_i*PerturbedStrainMag*PerturbedStrainMag - 2.0/3.0*Density_i*TurbVar_i[0]*diverg; @@ -905,152 +906,19 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi } -void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){ - ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i+1, Eddy_Viscosity_i, Density_i, turb_ke, true); - for(unsigned short iDim=0; iDim<3; iDim++){ - for(unsigned short jDim=0; jDim<3; jDim++){ - MeanReynoldsStress[iDim][jDim] /= (-Density_i); - } - } -} - -void CSourcePieceWise_TurbSST::SetPerturbedRSM(su2double turb_ke, const CConfig* config){ - - unsigned short iDim,jDim; - - /* --- Calculate anisotropic part of Reynolds Stress tensor --- */ - - for (iDim = 0; iDim< 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - A_ij[iDim][jDim] = .5 * MeanReynoldsStress[iDim][jDim] / turb_ke - delta3[iDim][jDim] / 3.0; - Eig_Vec[iDim][jDim] = A_ij[iDim][jDim]; - } - } - - /* --- Get ordered eigenvectors and eigenvalues of A_ij --- */ - - EigenDecomposition(A_ij, Eig_Vec, Eig_Val, 3); - - /* compute convex combination coefficients */ - su2double c1c = Eig_Val[2] - Eig_Val[1]; - su2double c2c = 2.0 * (Eig_Val[1] - Eig_Val[0]); - su2double c3c = 3.0 * Eig_Val[0] + 1.0; - - /* define barycentric traingle corner points */ - Corners[0][0] = 1.0; - Corners[0][1] = 0.0; - Corners[1][0] = 0.0; - Corners[1][1] = 0.0; - Corners[2][0] = 0.5; - Corners[2][1] = 0.866025; - - /* define barycentric coordinates */ - Barycentric_Coord[0] = Corners[0][0] * c1c + Corners[1][0] * c2c + Corners[2][0] * c3c; - Barycentric_Coord[1] = Corners[0][1] * c1c + Corners[1][1] * c2c + Corners[2][1] * c3c; - - if (Eig_Val_Comp == 1) { - /* 1C turbulence */ - New_Coord[0] = Corners[0][0]; - New_Coord[1] = Corners[0][1]; - } - else if (Eig_Val_Comp == 2) { - /* 2C turbulence */ - New_Coord[0] = Corners[1][0]; - New_Coord[1] = Corners[1][1]; - } - else if (Eig_Val_Comp == 3) { - /* 3C turbulence */ - New_Coord[0] = Corners[2][0]; - New_Coord[1] = Corners[2][1]; - } - else { - /* 2C turbulence */ - New_Coord[0] = Corners[1][0]; - New_Coord[1] = Corners[1][1]; - } - /* calculate perturbed barycentric coordinates */ - - Barycentric_Coord[0] = Barycentric_Coord[0] + (uq_delta_b) * (New_Coord[0] - Barycentric_Coord[0]); - Barycentric_Coord[1] = Barycentric_Coord[1] + (uq_delta_b) * (New_Coord[1] - Barycentric_Coord[1]); - - /* rebuild c1c,c2c,c3c based on new barycentric coordinates */ - c3c = Barycentric_Coord[1] / Corners[2][1]; - c1c = Barycentric_Coord[0] - Corners[2][0] * c3c; - c2c = 1 - c1c - c3c; - - /* build new anisotropy eigenvalues */ - Eig_Val[0] = (c3c - 1) / 3.0; - Eig_Val[1] = 0.5 *c2c + Eig_Val[0]; - Eig_Val[2] = c1c + Eig_Val[1]; - - /* permute eigenvectors if required */ - if (uq_permute) { - for (iDim=0; iDim<3; iDim++) { - for (jDim=0; jDim<3; jDim++) { - New_Eig_Vec[iDim][jDim] = Eig_Vec[2-iDim][jDim]; - } - } - } - - else { - for (iDim=0; iDim<3; iDim++) { - for (jDim=0; jDim<3; jDim++) { - New_Eig_Vec[iDim][jDim] = Eig_Vec[iDim][jDim]; - } - } - } - - EigenRecomposition(newA_ij, New_Eig_Vec, Eig_Val, 3); - - /* compute perturbed Reynolds stress matrix; use under-relaxation factor (urlx)*/ - for (iDim = 0; iDim< 3; iDim++){ - for (jDim = 0; jDim < 3; jDim++){ - MeanPerturbedRSM[iDim][jDim] = 2.0 * turb_ke * (newA_ij[iDim][jDim] + 1.0/3.0 * delta3[iDim][jDim]); - MeanPerturbedRSM[iDim][jDim] = MeanReynoldsStress[iDim][jDim] + - uq_urlx*(MeanPerturbedRSM[iDim][jDim] - MeanReynoldsStress[iDim][jDim]); - } - } +void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke){ -} + /*--- Compute norm of perturbed strain rate tensor. ---*/ -void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke){ - unsigned short iDim, jDim; PerturbedStrainMag = 0; - su2double **StrainRate = new su2double* [nDim]; - for (iDim= 0; iDim< nDim; iDim++){ - StrainRate[iDim] = new su2double [nDim]; - } + for (unsigned short iDim = 0; iDim < nDim; iDim++){ + for (unsigned short jDim = 0; jDim < nDim; jDim++){ + su2double StrainRate_ij = MeanPerturbedRSM[iDim][jDim] - TWO3 * turb_ke * delta[iDim][jDim]; + StrainRate_ij = - StrainRate_ij * Density_i / (2 * Eddy_Viscosity_i); - /* compute perturbed strain rate tensor */ - - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim =0; jDim < nDim; jDim++){ - StrainRate[iDim][jDim] = MeanPerturbedRSM[iDim][jDim] - - TWO3 * turb_ke * delta[iDim][jDim]; - StrainRate[iDim][jDim] = - StrainRate[iDim][jDim] * Density_i / (2 * Eddy_Viscosity_i); + PerturbedStrainMag += pow(StrainRate_ij, 2.0); } } - - /*--- Add diagonal part ---*/ - - for (iDim = 0; iDim < nDim; iDim++) { - PerturbedStrainMag += pow(StrainRate[iDim][iDim], 2.0); - } - - /*--- Add off diagonals ---*/ - - PerturbedStrainMag += 2.0*pow(StrainRate[1][0], 2.0); - - if (nDim == 3) { - PerturbedStrainMag += 2.0*pow(StrainRate[0][2], 2.0); - PerturbedStrainMag += 2.0*pow(StrainRate[1][2], 2.0); - } - PerturbedStrainMag = sqrt(2.0*PerturbedStrainMag); - for (iDim= 0; iDim< nDim; iDim++){ - delete [] StrainRate[iDim]; - } - - delete [] StrainRate; } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index c8b2d5f1e0c3..467f4e9b9fa3 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -349,26 +349,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, /*--- Add the solver name (max 8 characters). ---*/ SolverName = "C.FLOW"; - /*--- Vectorized numerics. ---*/ - if (config->GetUseVectorization()) { - const bool uncertain = config->GetUsing_UQ(); - const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || - (config->GetKind_FluidModel() == IDEAL_GAS); - const bool low_mach_corr = config->Low_Mach_Correction(); - - if (uncertain || !ideal_gas || low_mach_corr) { - SU2_MPI::Error("Some of the requested features are not yet " - "supported with vectorization.", CURRENT_FUNCTION); - } - - edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, iMesh); - - if (!edgeNumerics) { - SU2_MPI::Error("The numerical scheme in use does not " - "support vectorization.", CURRENT_FUNCTION); - } - } - /*--- Finally, check that the static arrays will be large enough (keep this * check at the bottom to make sure we consider the "final" values). ---*/ if((nDim > MAXNDIM) || (nPrimVar > MAXNVAR) || (nSecondaryVar > MAXNVAR)) @@ -688,6 +668,27 @@ CEulerSolver::~CEulerSolver(void) { } +void CEulerSolver::InstantiateEdgeNumerics(const CSolver* const* solver_container, const CConfig* config) { + + SU2_OMP_BARRIER + SU2_OMP_MASTER { + + if (config->Low_Mach_Correction()) + SU2_MPI::Error("Low-Mach correction is not supported with vectorization.", CURRENT_FUNCTION); + + if (solver_container[TURB_SOL]) + edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, MGLevel, solver_container[TURB_SOL]->GetNodes()); + else + edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, MGLevel); + + if (!edgeNumerics) + SU2_MPI::Error("The numerical scheme + gas model in use do not " + "support vectorization.", CURRENT_FUNCTION); + + } + SU2_OMP_BARRIER +} + void CEulerSolver::InitTurboContainers(CGeometry *geometry, CConfig *config){ unsigned short iMarker, iSpan, iVar; nSpanMax = config->GetnSpanMaxAllZones(); @@ -2645,95 +2646,16 @@ void CEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep) { - /*--- If possible use the vectorized numerics instead. ---*/ - if (edgeNumerics) { EdgeFluxResidual(geometry, config); return; } - - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool jst_scheme = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); - const bool jst_ke_scheme = (config->GetKind_Centered_Flow() == JST_KE) && (iMesh == MESH_0); - - /*--- Pick one numerics object per thread. ---*/ - CNumerics* numerics = numerics_container[CONV_TERM + omp_get_thread_num()*MAX_TERMS]; - - /*--- Loop over edge colors. ---*/ - for (auto color : EdgeColoring) - { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) - for(auto k = 0ul; k < color.size; ++k) { - - auto iEdge = color.indices[k]; - - /*--- Points in edge, set normal vectors, and number of neighbors ---*/ - - auto iPoint = geometry->edges->GetNode(iEdge,0); - auto jPoint = geometry->edges->GetNode(iEdge,1); - - numerics->SetNormal(geometry->edges->GetNormal(iEdge)); - numerics->SetNeighbor(geometry->nodes->GetnNeighbor(iPoint), geometry->nodes->GetnNeighbor(jPoint)); - - /*--- Set primitive variables w/o reconstruction ---*/ - - numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint)); - - /*--- Set the largest convective eigenvalue ---*/ - - numerics->SetLambda(nodes->GetLambda(iPoint), nodes->GetLambda(jPoint)); - - /*--- Set undivided laplacian an pressure based sensor ---*/ - - if (jst_scheme) { - numerics->SetUndivided_Laplacian(nodes->GetUndivided_Laplacian(iPoint), - nodes->GetUndivided_Laplacian(jPoint)); - } - if (jst_scheme || jst_ke_scheme) { - numerics->SetSensor(nodes->GetSensor(iPoint), nodes->GetSensor(jPoint)); - } - - /*--- Grid movement ---*/ - - if (dynamic_grid) { - numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(jPoint)); - } - - /*--- Compute residuals, and Jacobians ---*/ - - auto residual = numerics->ComputeResidual(config); - - /*--- Update convective and artificial dissipation residuals. ---*/ - - if (ReducerStrategy) { - EdgeFluxes.SetBlock(iEdge, residual); - if (implicit) - Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); - } - else { - LinSysRes.AddBlock(iPoint, residual); - LinSysRes.SubtractBlock(jPoint, residual); - if (implicit) - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); - } - - /*--- Viscous contribution. ---*/ - - Viscous_Residual(iEdge, geometry, solver_container, - numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); - } - } // end color loop - - if (ReducerStrategy) { - SumEdgeFluxes(geometry); - if (implicit) - Jacobian.SetDiagonalAsColumnSum(); - } - + EdgeFluxResidual(geometry, solver_container, config); } void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { - /*--- If possible use the vectorized numerics instead. ---*/ - if (edgeNumerics) { EdgeFluxResidual(geometry, config); return; } + if (config->GetUseVectorization()) { + EdgeFluxResidual(geometry, solver_container, config); + return; + } const auto InnerIter = config->GetInnerIter(); const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); diff --git a/TestCases/cont_adj_euler/naca0012/of_grad_cd.dat.ref b/TestCases/cont_adj_euler/naca0012/of_grad_cd.dat.ref index 37151f9b51b8..2f88439a0b9a 100644 --- a/TestCases/cont_adj_euler/naca0012/of_grad_cd.dat.ref +++ b/TestCases/cont_adj_euler/naca0012/of_grad_cd.dat.ref @@ -1,39 +1,39 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , 0.00535089 , 0.0001 - 1 , 0.0242911 , 0.0001 - 2 , 0.0299764 , 0.0001 - 3 , 0.0234387 , 0.0001 - 4 , 0.00340296 , 0.0001 - 5 , -0.0308074 , 0.0001 - 6 , -0.0775488 , 0.0001 - 7 , -0.132036 , 0.0001 - 8 , -0.186829 , 0.0001 - 9 , -0.23393 , 0.0001 - 10 , -0.268075 , 0.0001 - 11 , -0.28975 , 0.0001 - 12 , -0.305723 , 0.0001 - 13 , -0.325869 , 0.0001 - 14 , -0.358213 , 0.0001 - 15 , -0.407235 , 0.0001 - 16 , -0.479256 , 0.0001 - 17 , -0.595829 , 0.0001 - 18 , -0.849999 , 0.0001 - 19 , 0.56594 , 0.0001 - 20 , 1.14452 , 0.0001 - 21 , 1.50242 , 0.0001 - 22 , 1.6287 , 0.0001 - 23 , 1.54945 , 0.0001 - 24 , 1.3156 , 0.0001 - 25 , 0.989938 , 0.0001 - 26 , 0.634635 , 0.0001 - 27 , 0.301857 , 0.0001 - 28 , 0.0291461 , 0.0001 - 29 , -0.160219 , 0.0001 - 30 , -0.253768 , 0.0001 - 31 , -0.247224 , 0.0001 - 32 , -0.147422 , 0.0001 - 33 , 0.0181747 , 0.0001 - 34 , 0.197239 , 0.0001 - 35 , 0.341459 , 0.0001 + 0 , 0.00534068 , 0.0001 + 1 , 0.024277 , 0.0001 + 2 , 0.0299625 , 0.0001 + 3 , 0.0234274 , 0.0001 + 4 , 0.00339532 , 0.0001 + 5 , -0.0308111 , 0.0001 + 6 , -0.0775489 , 0.0001 + 7 , -0.132033 , 0.0001 + 8 , -0.186824 , 0.0001 + 9 , -0.233924 , 0.0001 + 10 , -0.268068 , 0.0001 + 11 , -0.289743 , 0.0001 + 12 , -0.305716 , 0.0001 + 13 , -0.325863 , 0.0001 + 14 , -0.358209 , 0.0001 + 15 , -0.407234 , 0.0001 + 16 , -0.47926 , 0.0001 + 17 , -0.595838 , 0.0001 + 18 , -0.850012 , 0.0001 + 19 , 0.565996 , 0.0001 + 20 , 1.14457 , 0.0001 + 21 , 1.50247 , 0.0001 + 22 , 1.62875 , 0.0001 + 23 , 1.54948 , 0.0001 + 24 , 1.31563 , 0.0001 + 25 , 0.989959 , 0.0001 + 26 , 0.634653 , 0.0001 + 27 , 0.301875 , 0.0001 + 28 , 0.0291654 , 0.0001 + 29 , -0.160196 , 0.0001 + 30 , -0.253741 , 0.0001 + 31 , -0.247195 , 0.0001 + 32 , -0.147396 , 0.0001 + 33 , 0.0181901 , 0.0001 + 34 , 0.197238 , 0.0001 + 35 , 0.341451 , 0.0001 36 , 0.467151 , 0.0001 - 37 , 0.682813 , 0.0001 + 37 , 0.682825 , 0.0001 diff --git a/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref b/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref index 9900d6061898..b471e28d3149 100644 --- a/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref +++ b/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref @@ -1,39 +1,39 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , -2524.84 , 0.001 - 1 , -12709.5 , 0.001 - 2 , -21482.2 , 0.001 - 3 , -27443.5 , 0.001 - 4 , -30308.7 , 0.001 - 5 , -30426.4 , 0.001 - 6 , -28432.2 , 0.001 - 7 , -24977.3 , 0.001 - 8 , -20586.8 , 0.001 - 9 , -15680.0 , 0.001 - 10 , -10727.7 , 0.001 - 11 , -6467.93 , 0.001 - 12 , -4048.9 , 0.001 - 13 , -4948.64 , 0.001 - 14 , -10478.5 , 0.001 - 15 , -20444.6 , 0.001 - 16 , -29988.8 , 0.001 - 17 , -26118.8 , 0.001 - 18 , -61059.9 , 0.001 - 19 , -3516.44 , 0.001 - 20 , -1725.14 , 0.001 - 21 , -1598.14 , 0.001 - 22 , -2571.41 , 0.001 - 23 , -6220.88 , 0.001 - 24 , -13796.7 , 0.001 - 25 , -25288.2 , 0.001 - 26 , -39170.3 , 0.001 - 27 , -52448.2 , 0.001 - 28 , -60910.8 , 0.001 - 29 , -59787.0 , 0.001 - 30 , -45235.2 , 0.001 - 31 , -17007.3 , 0.001 - 32 , 18476.1 , 0.001 - 33 , 48243.5 , 0.001 - 34 , 61878.7 , 0.001 - 35 , 63931.9 , 0.001 - 36 , 54760.0 , 0.001 - 37 , 62319.5 , 0.001 + 0 , -2654.17 , 0.001 + 1 , -12995.5 , 0.001 + 2 , -21781.7 , 0.001 + 3 , -27675.9 , 0.001 + 4 , -30438.6 , 0.001 + 5 , -30444.6 , 0.001 + 6 , -28343.8 , 0.001 + 7 , -24794.9 , 0.001 + 8 , -20330.1 , 0.001 + 9 , -15377.4 , 0.001 + 10 , -10418.8 , 0.001 + 11 , -6202.42 , 0.001 + 12 , -3878.85 , 0.001 + 13 , -4911.73 , 0.001 + 14 , -10574.3 , 0.001 + 15 , -20611.7 , 0.001 + 16 , -30111.5 , 0.001 + 17 , -26295.8 , 0.001 + 18 , -62395.3 , 0.001 + 19 , -2807.68 , 0.001 + 20 , -838.255 , 0.001 + 21 , -672.494 , 0.001 + 22 , -1776.6 , 0.001 + 23 , -5721.22 , 0.001 + 24 , -13719.0 , 0.001 + 25 , -25698.8 , 0.001 + 26 , -40060.1 , 0.001 + 27 , -53726.4 , 0.001 + 28 , -62410.5 , 0.001 + 29 , -61286.9 , 0.001 + 30 , -46498.5 , 0.001 + 31 , -17834.2 , 0.001 + 32 , 18190.2 , 0.001 + 33 , 48451.4 , 0.001 + 34 , 62322.0 , 0.001 + 35 , 64135.8 , 0.001 + 36 , 54563.9 , 0.001 + 37 , 64854.5 , 0.001 diff --git a/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref b/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref index 1c0e0f22f813..4992e46a279c 100644 --- a/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref +++ b/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref @@ -1,4 +1,4 @@ VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE" - 0 , 0.2285012284 , -106.6962871 , 0.2620099206 , -1.533188961 , 0.0 , -1.53853982 , 0.0 , 0.0 , 1.203388761 , 0.0 - 1 , 0.3866276739 , -174.5390955 , 0.439411593 , -2.414835143 , 0.0 , -2.423846191 , 0.0 , 0.0 , 1.053071502 , 0.0 - 2 , 0.5178161379 , -230.224511 , 0.5861998797 , -3.128333793 , 0.0 , -3.140377217 , 0.0 , 0.0 , 0.6530524517 , 0.0 + 0 , 0.2253591473 , -105.6097088 , 0.2588459007 , -1.5322178 , 0.0 , -1.537499867 , 0.0 , 0.0 , 1.202899757 , 0.0 + 1 , 0.3835809166 , -173.3502205 , 0.4363886002 , -2.415957492 , 0.0 , -2.424902327 , 0.0 , 0.0 , 1.053347497 , 0.0 + 2 , 0.5151776249 , -228.9760041 , 0.5835870252 , -3.129538494 , 0.0 , -3.141524632 , 0.0 , 0.0 , 0.6540715539 , 0.0 diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 054a7ba1eb48..4dc0b6d505c1 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -47,7 +47,7 @@ def main(): channel.cfg_dir = "euler/channel" channel.cfg_file = "inv_channel_RK.cfg" channel.test_iter = 20 - channel.test_vals = [-2.667326, 2.797439, 0.018717, 0.006906] + channel.test_vals = [-2.667328, 2.797437, 0.018714, 0.006906] test_list.append(channel) # NACA0012 @@ -71,7 +71,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [0.281704, 0.011821] + oneram6.test_vals = [0.281703, 0.011821] test_list.append(oneram6) # Fixed CL NACA0012 @@ -79,7 +79,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-12.130189, -6.702728, 0.300000, 0.019470] + fixedCL_naca0012.test_vals = [-7.374790, -1.872333, 0.300000, 0.019471] test_list.append(fixedCL_naca0012) # HYPERSONIC FLOW PAST BLUNT BODY @@ -143,7 +143,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.021224, -5.268445, 0.807582, 0.060731] + rae2822_sa.test_vals = [-2.020123, -5.269330, 0.807147, 0.060499] test_list.append(rae2822_sa) # RAE2822 SST @@ -151,7 +151,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-0.510637, 4.876603, 0.812485, 0.061969] + rae2822_sst.test_vals = [-0.510633, 4.871233, 0.811923, 0.061627] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -159,7 +159,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.429813, 4.876602, 0.812485, 0.061969] + rae2822_sst_sust.test_vals = [-2.430689, 4.871233, 0.811923, 0.061627] test_list.append(rae2822_sst_sust) # Flat plate @@ -183,7 +183,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.537781, -14.899750, 1.064330, 0.019756] + turb_naca0012_sa.test_vals = [-11.531286, -14.899968, 1.064330, 0.019756] test_list.append(turb_naca0012_sa) # NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253) @@ -191,7 +191,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 = [-12.797090, -5.872763, 1.049989, 0.019163] + turb_naca0012_sst.test_vals = [-12.797872, -5.863656, 1.049989, 0.019163] test_list.append(turb_naca0012_sst) # NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253) @@ -199,7 +199,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.640091, -5.751854, 1.005233, 0.019017] + turb_naca0012_sst_sust.test_vals = [-12.640670, -5.746919, 1.005233, 0.019017] test_list.append(turb_naca0012_sst_sust) # PROPELLER @@ -220,7 +220,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-7.653296, -7.729472, -1.981061, -0.000016, 0.079062] + turb_naca0012_sst_restart_mg.test_vals = [-7.652983, -7.729472, -1.981061, -0.000015, 0.079061] test_list.append(turb_naca0012_sst_restart_mg) ############################# @@ -232,7 +232,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.978913, 1.140283, 1.211887, 0.194208] + turb_naca0012_1c.test_vals = [-4.979967, 1.140048, 1.215111, 0.195366] test_list.append(turb_naca0012_1c) # NACA0012 2c @@ -276,7 +276,7 @@ def main(): harmonic_balance.cfg_dir = "harmonic_balance" harmonic_balance.cfg_file = "HB.cfg" harmonic_balance.test_iter = 25 - harmonic_balance.test_vals = [-1.589755, 3.922207, 0.006725, 0.099455] + harmonic_balance.test_vals = [-1.589740, 3.922578, 0.006703, 0.099632] harmonic_balance.new_output = False test_list.append(harmonic_balance) @@ -285,7 +285,7 @@ def main(): hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 - hb_rans_preconditioning.test_vals = [-1.902391, -5.950120, 0.007786, 0.128110] + hb_rans_preconditioning.test_vals = [-1.902111, -5.949291, 0.007768, 0.128060] hb_rans_preconditioning.new_output = False test_list.append(hb_rans_preconditioning) @@ -327,7 +327,7 @@ def main(): sine_gust.cfg_dir = "gust" sine_gust.cfg_file = "inv_gust_NACA0012.cfg" sine_gust.test_iter = 5 - sine_gust.test_vals = [-1.977520, 3.481804, -0.012404, -0.007452] + sine_gust.test_vals = [-1.977520, 3.481804, -0.012403, -0.007453] sine_gust.unsteady = True test_list.append(sine_gust) @@ -336,7 +336,7 @@ def main(): aeroelastic.cfg_dir = "aeroelastic" aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" aeroelastic.test_iter = 2 - aeroelastic.test_vals = [0.074525, 0.033127, -0.001650, -0.000127] + aeroelastic.test_vals = [0.074447, 0.033116, -0.001649, -0.000127] aeroelastic.unsteady = True test_list.append(aeroelastic) @@ -414,7 +414,7 @@ def main(): transonic_stator_rst.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_rst.cfg_file = "transonic_stator_rst.cfg" transonic_stator_rst.test_iter = 20 - transonic_stator_rst.test_vals = [-6.621626, -0.614379, 5.002986, 0.002951] + transonic_stator_rst.test_vals = [-6.621626, -0.614366, 5.002986, 0.002951] transonic_stator_rst.new_output = False test_list.append(transonic_stator_rst) @@ -447,7 +447,7 @@ def main(): channel_3D.cfg_dir = "sliding_interface/channel_3D" channel_3D.cfg_file = "channel_3D_WA.cfg" channel_3D.test_iter = 2 - channel_3D.test_vals = [2.000000, 0.000000, 0.620151, 0.505157, 0.415249] + channel_3D.test_vals = [2.000000, 0.000000, 0.620151, 0.505156, 0.415292] channel_3D.unsteady = True channel_3D.multizone = True test_list.append(channel_3D) @@ -487,7 +487,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.619179, -1.564701] + bars_SST_2D.test_vals = [13.000000, -0.619686, -1.564594] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -500,7 +500,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-2.378370, -1.585252, -2.028505, 6.4359e+04] + statbeam3d.test_vals = [-2.378370, -1.585252, -2.028505, 64359.000000] test_list.append(statbeam3d) # Dynamic beam, 2d diff --git a/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref b/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref index 37e5fff513f6..42013415280d 100644 --- a/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref +++ b/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref @@ -1,3 +1,3 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , 0.0770973 , 0.001 - 1 , -0.113478 , 0.001 + 0 , 0.076707 , 0.001 + 1 , -0.113016 , 0.001 diff --git a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref index e279f1f88b94..5e90b70e61f9 100644 --- a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref +++ b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref @@ -1,3 +1,3 @@ VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE" - 0 , 0.05698141666 , -0.4016762951 , 0.05429379682 , 0.1237934687 , 0.0 , 0.1225795963 , 0.0 , 0.0 , 0.009159910311 , 0.0 - 1 , -0.08029896733 , 8.980081481 , -0.08626576809 , 0.2725786469 , 0.0 , 0.2743956584 , 0.0 , 0.0 , 0.008211157534 , 0.0 + 0 , 0.05682885157 , -0.4029984608 , 0.05415019059 , 0.123381226 , 0.0 , 0.1221705845 , 0.0 , 0.0 , 0.009249712749 , 0.0 + 1 , -0.08020054725 , 8.95834958 , -0.08614128407 , 0.2713852273 , 0.0 , 0.2731998072 , 0.0 , 0.0 , 0.00847998417 , 0.0 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 7db7680864fe..1375f196c901 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -99,7 +99,7 @@ def main(): channel.cfg_dir = "euler/channel" channel.cfg_file = "inv_channel_RK.cfg" channel.test_iter = 20 - channel.test_vals = [-2.647975, 2.818091, 0.022281, 0.004645] + channel.test_vals = [-2.647975, 2.818090, 0.022280, 0.004644] channel.su2_exec = "parallel_computation.py -f" channel.timeout = 1600 channel.tol = 0.00001 @@ -132,7 +132,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-10.198925, -9.616270, 0.281704, 0.011821] + oneram6.test_vals = [-9.277150, -8.694005, 0.281703, 0.011821] oneram6.su2_exec = "parallel_computation.py -f" oneram6.timeout = 3200 oneram6.tol = 0.00001 @@ -143,7 +143,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-12.134366, -6.700512, 0.300000, 0.019470] + fixedCL_naca0012.test_vals = [-7.379831, -1.886302, 0.300000, 0.019471] fixedCL_naca0012.su2_exec = "parallel_computation.py -f" fixedCL_naca0012.timeout = 1600 fixedCL_naca0012.tol = 0.00001 @@ -155,7 +155,7 @@ def main(): polar_naca0012.cfg_file = "inv_NACA0012.cfg" polar_naca0012.polar = True polar_naca0012.test_iter = 10 - polar_naca0012.test_vals = [-1.220449, 4.250991, 0.008963, 0.016640] + polar_naca0012.test_vals = [-1.217981, 4.256386, 0.009084, 0.016823] polar_naca0012.su2_exec = "compute_polar.py -i 11" polar_naca0012.timeout = 1600 polar_naca0012.tol = 0.00001 @@ -240,7 +240,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.005647, -5.264972, 0.809933, 0.062257] + rae2822_sa.test_vals = [-2.004689, -5.265793, 0.809463, 0.062016] rae2822_sa.su2_exec = "parallel_computation.py -f" rae2822_sa.timeout = 1600 rae2822_sa.tol = 0.00001 @@ -251,7 +251,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-0.510644, 4.875481, 0.814472, 0.062806] + rae2822_sst.test_vals = [-0.510641, 4.870022, 0.813722, 0.062439] rae2822_sst.su2_exec = "parallel_computation.py -f" rae2822_sst.timeout = 1600 rae2822_sst.tol = 0.00001 @@ -262,7 +262,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.435061, 4.875481, 0.814472, 0.062806] + rae2822_sst_sust.test_vals = [-2.435890, 4.870022, 0.813722, 0.062439] rae2822_sst_sust.su2_exec = "parallel_computation.py -f" rae2822_sst_sust.timeout = 1600 rae2822_sst_sust.tol = 0.00001 @@ -295,7 +295,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.155953, -14.468619, 1.064330, 0.019756] #last 4 columns + turb_naca0012_sa.test_vals = [-11.147929, -14.466781, 1.064330, 0.019756] turb_naca0012_sa.su2_exec = "parallel_computation.py -f" turb_naca0012_sa.timeout = 3200 turb_naca0012_sa.tol = 0.00001 @@ -306,7 +306,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 = [-12.799245, -5.875128, 1.049989, 0.019163] #last 4 columns + turb_naca0012_sst.test_vals = [-12.800055, -5.865784, 1.049989, 0.019163] turb_naca0012_sst.su2_exec = "parallel_computation.py -f" turb_naca0012_sst.timeout = 3200 turb_naca0012_sst.tol = 0.00001 @@ -317,7 +317,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.641087, -5.753486, 1.005233, 0.019017] #last 4 columns + turb_naca0012_sst_sust.test_vals = [-12.641676, -5.748419, 1.005233, 0.019017] turb_naca0012_sst_sust.su2_exec = "parallel_computation.py -f" turb_naca0012_sst_sust.timeout = 3200 turb_naca0012_sst_sust.tol = 0.00001 @@ -344,7 +344,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-7.620278, -7.729499, -1.981040, -0.000016, 0.079062] + turb_naca0012_sst_restart_mg.test_vals = [-7.619889, -7.729499, -1.981039, -0.000016, 0.079062] turb_naca0012_sst_restart_mg.su2_exec = "parallel_computation.py -f" turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 @@ -594,7 +594,7 @@ def main(): contadj_fixed_CL_naca0012.cfg_dir = "fixed_cl/naca0012" contadj_fixed_CL_naca0012.cfg_file = "inv_NACA0012_ContAdj.cfg" contadj_fixed_CL_naca0012.test_iter = 100 - contadj_fixed_CL_naca0012.test_vals = [0.275868, -5.200487, 0.342720, 0.000105] + contadj_fixed_CL_naca0012.test_vals = [0.275856, -5.200511, 0.342710, 0.000105] contadj_fixed_CL_naca0012.su2_exec = "parallel_computation.py -f" contadj_fixed_CL_naca0012.timeout = 1600 contadj_fixed_CL_naca0012.tol = 0.00001 @@ -683,7 +683,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.973120, 1.141760, 0.861180, 0.014211] + turb_naca0012_1c.test_vals = [-4.981674, 1.138300, 0.867182, 0.015761] turb_naca0012_1c.su2_exec = "parallel_computation.py -f" turb_naca0012_1c.timeout = 1600 turb_naca0012_1c.tol = 0.00001 @@ -742,7 +742,7 @@ def main(): harmonic_balance.cfg_dir = "harmonic_balance" harmonic_balance.cfg_file = "HB.cfg" harmonic_balance.test_iter = 25 - harmonic_balance.test_vals = [-1.589755, 3.922208, 0.006724, 0.099454] #last 4 columns + harmonic_balance.test_vals = [-1.589739, 3.922579, 0.006702, 0.099632] harmonic_balance.su2_exec = "parallel_computation.py -f" harmonic_balance.timeout = 1600 harmonic_balance.tol = 0.00001 @@ -754,7 +754,7 @@ def main(): hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 - hb_rans_preconditioning.test_vals = [-1.902361, -5.950093, 0.007786, 0.128110] #last 4 columns + hb_rans_preconditioning.test_vals = [-1.902099, -5.949279, 0.007768, 0.128062] hb_rans_preconditioning.su2_exec = "parallel_computation.py -f" hb_rans_preconditioning.timeout = 1600 hb_rans_preconditioning.tol = 0.00001 @@ -770,7 +770,7 @@ def main(): rot_naca0012.cfg_dir = "rotating/naca0012" rot_naca0012.cfg_file = "rot_NACA0012.cfg" rot_naca0012.test_iter = 25 - rot_naca0012.test_vals = [-2.668380, 2.875334, -0.079260, 0.002252] + rot_naca0012.test_vals = [-2.698005, 2.845328, -0.079439, 0.002128] rot_naca0012.su2_exec = "parallel_computation.py -f" rot_naca0012.timeout = 1600 rot_naca0012.tol = 0.00001 @@ -831,7 +831,7 @@ def main(): aeroelastic.cfg_dir = "aeroelastic" aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" aeroelastic.test_iter = 2 - aeroelastic.test_vals = [0.076565, 0.033057, -0.001650, -0.000127] + aeroelastic.test_vals = [0.076550, 0.033042, -0.001650, -0.000127] aeroelastic.su2_exec = "parallel_computation.py -f" aeroelastic.timeout = 1600 aeroelastic.tol = 0.00001 @@ -1039,7 +1039,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.619179, -1.564701] #last 4 columns + bars_SST_2D.test_vals = [13.000000, -0.619686, -1.564594] bars_SST_2D.su2_exec = "SU2_CFD" bars_SST_2D.timeout = 1600 bars_SST_2D.tol = 0.00001 @@ -1127,7 +1127,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.379832, -4.005999, 0.000000, 0.000000] #last 4 columns + dyn_fsi.test_vals = [-4.379832, -4.006001, 0.000000, 0.000000] dyn_fsi.multizone = True dyn_fsi.unsteady = True dyn_fsi.su2_exec = "mpirun -n 2 SU2_CFD" @@ -1197,7 +1197,7 @@ def main(): cht_compressible.cfg_dir = "coupled_cht/comp_2d" cht_compressible.cfg_file = "cht_2d_3cylinders.cfg" cht_compressible.test_iter = 10 - cht_compressible.test_vals = [-4.256303, -0.532538, -0.532538, -0.532537] #last 4 columns + cht_compressible.test_vals = [-4.256032, -0.532728, -0.532729, -0.532728] cht_compressible.su2_exec = "SU2_CFD" cht_compressible.timeout = 1600 cht_compressible.multizone = True @@ -1224,7 +1224,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 = [-12.799245, -5.875128, 1.049989, 0.019163] #last 4 columns + pywrapper_turb_naca0012_sst.test_vals = [-12.800055, -5.865784, 1.049989, 0.019163] pywrapper_turb_naca0012_sst.su2_exec = "mpirun -np 2 SU2_CFD.py --parallel -f" pywrapper_turb_naca0012_sst.timeout = 3200 pywrapper_turb_naca0012_sst.tol = 0.00001 @@ -1247,7 +1247,7 @@ def main(): pywrapper_aeroelastic.cfg_dir = "aeroelastic" pywrapper_aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" pywrapper_aeroelastic.test_iter = 2 - pywrapper_aeroelastic.test_vals = [0.076565, 0.033057, -0.001650, -0.000127] #last 4 columns + pywrapper_aeroelastic.test_vals = [0.076550, 0.033042, -0.001650, -0.000127] pywrapper_aeroelastic.su2_exec = "mpirun -np 2 SU2_CFD.py --parallel -f" pywrapper_aeroelastic.timeout = 1600 pywrapper_aeroelastic.tol = 0.00001 @@ -1272,7 +1272,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.245726, -0.001238, 0.175705] #last 4 columns + pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.245726, -0.001240, 0.175715] pywrapper_unsteadyCHT.su2_exec = "mpirun -np 2 python launch_unsteady_CHT_FlatPlate.py --parallel -f" pywrapper_unsteadyCHT.timeout = 1600 pywrapper_unsteadyCHT.tol = 0.00001 @@ -1285,7 +1285,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614165, 2.242643, -0.038313, 0.173856] #last 4 columns + pywrapper_rigidMotion.test_vals = [-1.614165, 2.242641, -0.038307, 0.173866] pywrapper_rigidMotion.su2_exec = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f" pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index d5485d16215a..f6b6f4f7cacd 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -47,7 +47,7 @@ def main(): discadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" discadj_naca0012.cfg_file = "inv_NACA0012_discadj.cfg" discadj_naca0012.test_iter = 100 - discadj_naca0012.test_vals = [-3.559002, -8.926022, -0.000000, 0.005588] #last 4 columns + discadj_naca0012.test_vals = [-3.561506, -8.926634, -0.000000, 0.005587] discadj_naca0012.su2_exec = "parallel_computation.py -f" discadj_naca0012.timeout = 1600 discadj_naca0012.tol = 0.00001 @@ -58,7 +58,7 @@ def main(): discadj_cylinder3D.cfg_dir = "disc_adj_euler/cylinder3D" discadj_cylinder3D.cfg_file = "inv_cylinder3D.cfg" discadj_cylinder3D.test_iter = 5 - discadj_cylinder3D.test_vals = [-3.756465, -3.861233, 0.000000, 0.000000] + discadj_cylinder3D.test_vals = [-3.756779, -3.861913, 0.000000, 0.000000] discadj_cylinder3D.su2_exec = "parallel_computation.py -f" discadj_cylinder3D.timeout = 1600 discadj_cylinder3D.tol = 0.00001 @@ -69,7 +69,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [2.190287, 1.635972, 47258.000000, 0.000000] + discadj_arina2k.test_vals = [2.189902, 1.635938, 47258.000000, 0.000000] discadj_arina2k.su2_exec = "parallel_computation.py -f" discadj_arina2k.timeout = 8400 discadj_arina2k.tol = 0.00001 @@ -95,7 +95,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.221231, -0.491824, 0.557480, 0.000027] #last 4 columns + discadj_rans_naca0012_sst.test_vals = [-2.221040, -0.491810, 0.557480, 0.000027] discadj_rans_naca0012_sst.su2_exec = "parallel_computation.py -f" discadj_rans_naca0012_sst.timeout = 1600 discadj_rans_naca0012_sst.tol = 0.00001 @@ -230,7 +230,7 @@ def main(): discadj_pitchingNACA0012.cfg_dir = "disc_adj_euler/naca0012_pitching" discadj_pitchingNACA0012.cfg_file = "inv_NACA0012_pitching.cfg" discadj_pitchingNACA0012.test_iter = 4 - discadj_pitchingNACA0012.test_vals = [-1.228572, -1.640409, -0.007617, 0.000013] + discadj_pitchingNACA0012.test_vals = [-1.223480, -1.639387, -0.007591, 0.000013] discadj_pitchingNACA0012.su2_exec = "parallel_computation.py -f" discadj_pitchingNACA0012.timeout = 1600 discadj_pitchingNACA0012.tol = 0.00001 @@ -302,7 +302,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-5.323990, -2.4380e-13] + discadj_fsi2.test_vals = [-5.318452, -2.4380e-13] discadj_fsi2.su2_exec = "mpirun -n 2 SU2_CFD_AD" discadj_fsi2.timeout = 1600 discadj_fsi2.tol = 1e-16 @@ -389,7 +389,7 @@ def main(): naca_restart_shape_opt.cfg_dir = "optimization_rans/naca0012" naca_restart_shape_opt.cfg_file = "naca0012.cfg" naca_restart_shape_opt.test_iter = 1 - naca_restart_shape_opt.test_vals = [1.000000, 1.000000, 0.007046, 0.196883] #last 4 columns + naca_restart_shape_opt.test_vals = [1.000000, 1.000000, 0.007046, 0.196671] naca_restart_shape_opt.su2_exec = "shape_optimization.py -f" naca_restart_shape_opt.timeout = 1600 naca_restart_shape_opt.tol = 0.00001 diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg index 9ed891401238..71e181c873b6 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg @@ -157,6 +157,7 @@ MG_DAMP_PROLONGATION= 0.75 % Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, % TURKEL_PREC, MSW) CONV_NUM_METHOD_FLOW= ROE +USE_VECTORIZATION= YES % % Spatial numerical order integration (1ST_ORDER, 2ND_ORDER, 2ND_ORDER_LIMITER) MUSCL_FLOW= YES diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 199bf0c425a8..3e7399b06b59 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -105,7 +105,7 @@ def main(): channel.cfg_dir = "euler/channel" channel.cfg_file = "inv_channel_RK.cfg" channel.test_iter = 10 - channel.test_vals = [-2.475874, 3.046368, -0.203968, 0.036020] #last 4 columns + channel.test_vals = [-2.475872, 3.046370, -0.203974, 0.036018] channel.su2_exec = "SU2_CFD" channel.timeout = 1600 channel.new_output = True @@ -141,7 +141,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-10.200724, -9.619291, 0.281704, 0.011821] #last 4 columns + oneram6.test_vals = [-9.279396, -8.697739, 0.281703, 0.011821] oneram6.su2_exec = "SU2_CFD" oneram6.timeout = 9600 oneram6.new_output = True @@ -153,7 +153,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-12.128275, -6.700329, 0.300000, 0.019470] #last 4 columns + fixedCL_naca0012.test_vals = [-7.382410, -1.879887, 0.300000, 0.019471] fixedCL_naca0012.su2_exec = "SU2_CFD" fixedCL_naca0012.new_output = True fixedCL_naca0012.timeout = 1600 @@ -167,7 +167,7 @@ def main(): polar_naca0012.polar = True polar_naca0012.new_output= True polar_naca0012.test_iter = 10 - polar_naca0012.test_vals = [-1.244014, 4.220468, 0.016275, 0.015988] #last 4 columns + polar_naca0012.test_vals = [-1.243326, 4.224483, 0.016432, 0.016145] polar_naca0012.su2_exec = "compute_polar.py -n 1 -i 11" polar_naca0012.timeout = 1600 polar_naca0012.tol = 0.00001 @@ -226,7 +226,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.850123, -1.388088, -0.056090, 108.140177] #last 4 columns + cylinder_lowmach.test_vals = [-6.850123, -1.388088, -0.056090, 108.140176] cylinder_lowmach.su2_exec = "SU2_CFD" cylinder_lowmach.new_output = True cylinder_lowmach.timeout = 1600 @@ -274,7 +274,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.021224, -5.268445, 0.807582, 0.060731] #last 4 columns + rae2822_sa.test_vals = [-2.020123, -5.269330, 0.807147, 0.060499] rae2822_sa.su2_exec = "SU2_CFD" rae2822_sa.timeout = 1600 rae2822_sa.new_output = True @@ -286,7 +286,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-0.510641, 4.877423, 0.813358, 0.061455] #last 4 columns + rae2822_sst.test_vals = [-0.510639, 4.872266, 0.812659, 0.061095] rae2822_sst.su2_exec = "SU2_CFD" rae2822_sst.new_output = True rae2822_sst.timeout = 1600 @@ -298,7 +298,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.430811, 4.877422, 0.813358, 0.061455] #last 4 columns + rae2822_sst_sust.test_vals = [-2.431741, 4.872266, 0.812658, 0.061095] rae2822_sst_sust.su2_exec = "SU2_CFD" rae2822_sst_sust.timeout = 1600 rae2822_sst_sust.tol = 0.00001 @@ -333,7 +333,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-11.141831, -14.498856, 1.064330, 0.019756] #last 4 columns + turb_naca0012_sa.test_vals = [-11.133933, -14.498178, 1.064330, 0.019756] turb_naca0012_sa.su2_exec = "SU2_CFD" turb_naca0012_sa.new_output = True turb_naca0012_sa.timeout = 3200 @@ -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 = [-12.797476, -5.873045, 1.049989, 0.019163] #last 4 columns + turb_naca0012_sst.test_vals = [-12.798258, -5.863895, 1.049989, 0.019163] turb_naca0012_sst.su2_exec = "SU2_CFD" turb_naca0012_sst.new_output = True turb_naca0012_sst.timeout = 3200 @@ -357,7 +357,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.640277, -5.752224, 1.005233, 0.019017] #last 4 columns + turb_naca0012_sst_sust.test_vals = [-12.640857, -5.747260, 1.005233, 0.019017] turb_naca0012_sst_sust.su2_exec = "SU2_CFD" turb_naca0012_sst_sust.timeout = 3200 turb_naca0012_sst_sust.tol = 0.00001 @@ -385,7 +385,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 50 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-7.653492, -7.729550, -1.981855, -0.000015, 0.079062] #last 5 columns + turb_naca0012_sst_restart_mg.test_vals = [-7.653235, -7.729550, -1.981855, -0.000015, 0.079061] turb_naca0012_sst_restart_mg.su2_exec = "SU2_CFD" turb_naca0012_sst_restart_mg.new_output = True turb_naca0012_sst_restart_mg.timeout = 3200 @@ -702,7 +702,7 @@ def main(): contadj_fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" contadj_fixedCL_naca0012.cfg_file = "inv_NACA0012_ContAdj.cfg" contadj_fixedCL_naca0012.test_iter = 100 - contadj_fixedCL_naca0012.test_vals = [0.293257, -5.201691, 0.360610, -0.000021] #last 4 columns + contadj_fixedCL_naca0012.test_vals = [0.293213, -5.201710, 0.360590, -0.000022] contadj_fixedCL_naca0012.su2_exec = "SU2_CFD" contadj_fixedCL_naca0012.new_output= True contadj_fixedCL_naca0012.timeout = 1600 @@ -790,7 +790,7 @@ def main(): contadj_rans_rae2822.cfg_dir = "cont_adj_rans/rae2822" contadj_rans_rae2822.cfg_file = "turb_SA_RAE2822.cfg" contadj_rans_rae2822.test_iter = 20 - contadj_rans_rae2822.test_vals = [-5.369688, -10.872209, -0.212470, 0.005448] #last 4 columns + contadj_rans_rae2822.test_vals = [-5.369688, -10.872211, -0.212470, 0.005448] contadj_rans_rae2822.su2_exec = "SU2_CFD" contadj_rans_rae2822.new_output = True contadj_rans_rae2822.timeout = 1600 @@ -806,7 +806,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.978068, 1.139123, 0.806952, 0.064685] #last 4 columns + turb_naca0012_1c.test_vals = [-4.977752, 1.139722, 0.804826, 0.063918] turb_naca0012_1c.su2_exec = "SU2_CFD" turb_naca0012_1c.new_output = True turb_naca0012_1c.timeout = 1600 @@ -870,7 +870,7 @@ def main(): harmonic_balance.cfg_dir = "harmonic_balance" harmonic_balance.cfg_file = "HB.cfg" harmonic_balance.test_iter = 25 - harmonic_balance.test_vals = [-1.589755, 3.922208, 0.006724, 0.099454] #last 4 columns + harmonic_balance.test_vals = [-1.589739, 3.922579, 0.006702, 0.099632] harmonic_balance.su2_exec = "SU2_CFD" harmonic_balance.new_output = False harmonic_balance.timeout = 1600 @@ -882,7 +882,7 @@ def main(): hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 - hb_rans_preconditioning.test_vals = [-1.902362, -5.950093, 0.007786, 0.128110] #last 4 columns + hb_rans_preconditioning.test_vals = [-1.902098, -5.949277, 0.007768, 0.128061] hb_rans_preconditioning.su2_exec = "SU2_CFD" hb_rans_preconditioning.new_output= False hb_rans_preconditioning.timeout = 1600 @@ -898,7 +898,7 @@ def main(): rot_naca0012.cfg_dir = "rotating/naca0012" rot_naca0012.cfg_file = "rot_NACA0012.cfg" rot_naca0012.test_iter = 25 - rot_naca0012.test_vals = [-2.657250, 2.889308, -0.079043, 0.002261] #last 4 columns + rot_naca0012.test_vals = [-2.688979, 2.857521, -0.079219, 0.002135] rot_naca0012.su2_exec = "SU2_CFD" rot_naca0012.timeout = 1600 rot_naca0012.tol = 0.00001 @@ -909,7 +909,7 @@ def main(): cavity.cfg_dir = "moving_wall/cavity" cavity.cfg_file = "lam_cavity.cfg" cavity.test_iter = 25 - cavity.test_vals = [-5.627934, -0.164470, 0.051972, 2.547034] #last 4 columns + cavity.test_vals = [-5.627934, -0.164470, 0.051972, 2.547039] cavity.su2_exec = "SU2_CFD" cavity.new_output = True cavity.timeout = 1600 @@ -950,7 +950,7 @@ def main(): sine_gust.cfg_dir = "gust" sine_gust.cfg_file = "inv_gust_NACA0012.cfg" sine_gust.test_iter = 5 - sine_gust.test_vals = [-1.977520, 3.481804, -0.012277, -0.007308] #last 4 columns + sine_gust.test_vals = [-1.977520, 3.481804, -0.012277, -0.007309] sine_gust.su2_exec = "SU2_CFD" sine_gust.timeout = 1600 sine_gust.tol = 0.00001 @@ -963,7 +963,7 @@ def main(): aeroelastic.cfg_dir = "aeroelastic" aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" aeroelastic.test_iter = 2 - aeroelastic.test_vals = [0.074885, 0.033116, -0.001650, -0.000127] #last 4 columns + aeroelastic.test_vals = [0.074836, 0.033102, -0.001650, -0.000127] aeroelastic.su2_exec = "SU2_CFD" aeroelastic.timeout = 1600 aeroelastic.tol = 0.00001 @@ -1209,7 +1209,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.619179, -1.564701] #last 3 columns + bars_SST_2D.test_vals = [13.000000, -0.619686, -1.564594] bars_SST_2D.su2_exec = "SU2_CFD" bars_SST_2D.timeout = 1600 bars_SST_2D.tol = 0.00001 @@ -1343,7 +1343,7 @@ def main(): airfoilRBF.cfg_dir = "fea_fsi/Airfoil_RBF" airfoilRBF.cfg_file = "config.cfg" airfoilRBF.test_iter = 1 - airfoilRBF.test_vals = [1.000000, -2.791154, -4.961536] + airfoilRBF.test_vals = [1.000000, -2.786185, -4.977948] airfoilRBF.su2_exec = "SU2_CFD" airfoilRBF.timeout = 1600 airfoilRBF.multizone = True @@ -1408,7 +1408,7 @@ def main(): cht_incompressible.cfg_dir = "coupled_cht/comp_2d" cht_incompressible.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible.test_iter = 10 - cht_incompressible.test_vals = [-4.256303, -0.532538, -0.532538, -0.532537] #last 4 columns + cht_incompressible.test_vals = [-4.256032, -0.532728, -0.532729, -0.532728] cht_incompressible.su2_exec = "SU2_CFD" cht_incompressible.timeout = 1600 cht_incompressible.multizone = True @@ -1795,7 +1795,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 = [-12.797476, -5.873045, 1.049989, 0.019163] #last 4 columns + pywrapper_turb_naca0012_sst.test_vals = [-12.798258, -5.863895, 1.049989, 0.019163] pywrapper_turb_naca0012_sst.su2_exec = "SU2_CFD.py -f" pywrapper_turb_naca0012_sst.new_output = True pywrapper_turb_naca0012_sst.timeout = 3200 @@ -1808,7 +1808,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 3 - pywrapper_square_cylinder.test_vals = [-1.162561, 0.066414, 1.399788, 2.220411] #last 4 columns + pywrapper_square_cylinder.test_vals = [-1.162560, 0.066414, 1.399788, 2.220411] pywrapper_square_cylinder.su2_exec = "SU2_CFD.py -f" pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 @@ -1821,7 +1821,7 @@ def main(): pywrapper_aeroelastic.cfg_dir = "aeroelastic" pywrapper_aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" pywrapper_aeroelastic.test_iter = 2 - pywrapper_aeroelastic.test_vals = [0.074885, 0.033116, -0.001650, -0.000127] #last 4 columns + pywrapper_aeroelastic.test_vals = [0.074836, 0.033102, -0.001650, -0.000127] pywrapper_aeroelastic.su2_exec = "SU2_CFD.py -f" pywrapper_aeroelastic.new_output = True pywrapper_aeroelastic.timeout = 1600 @@ -1850,7 +1850,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.245717, 0.000766, 0.175707] #last 4 columns + pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.245716, 0.000766, 0.175719] pywrapper_unsteadyCHT.su2_exec = "python launch_unsteady_CHT_FlatPlate.py -f" pywrapper_unsteadyCHT.timeout = 1600 pywrapper_unsteadyCHT.new_output = True @@ -1864,7 +1864,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614167, 2.242633, -0.037878, 0.173901] #last 4 columns + pywrapper_rigidMotion.test_vals = [-1.614167, 2.242632, -0.037871, 0.173912] pywrapper_rigidMotion.su2_exec = "python launch_flatPlate_rigidMotion.py -f" pywrapper_rigidMotion.new_output = True pywrapper_rigidMotion.timeout = 1600 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 183ad99c2196..f6aa6be6d034 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -47,7 +47,7 @@ def main(): discadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" discadj_naca0012.cfg_file = "inv_NACA0012_discadj.cfg" discadj_naca0012.test_iter = 100 - discadj_naca0012.test_vals = [-3.559002, -8.926022, -0.000000, 0.005588] #last 4 columns + discadj_naca0012.test_vals = [-3.561506, -8.926634, -0.000000, 0.005587] discadj_naca0012.su2_exec = "SU2_CFD_AD" discadj_naca0012.timeout = 1600 discadj_naca0012.tol = 0.00001 @@ -58,7 +58,7 @@ def main(): discadj_cylinder3D.cfg_dir = "disc_adj_euler/cylinder3D" discadj_cylinder3D.cfg_file = "inv_cylinder3D.cfg" discadj_cylinder3D.test_iter = 5 - discadj_cylinder3D.test_vals = [-3.759637, -3.864023, -0.000000, 0.000000] #last 4 columns + discadj_cylinder3D.test_vals = [-3.759951, -3.864587, -0.000000, 0.000000] discadj_cylinder3D.su2_exec = "SU2_CFD_AD" discadj_cylinder3D.timeout = 1600 discadj_cylinder3D.tol = 0.00001 @@ -69,7 +69,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [2.108283, 1.574295, 47250.0, 0.0]#last 4 columns + discadj_arina2k.test_vals = [2.107800, 1.574246, 47250.000000, 0.000000] discadj_arina2k.su2_exec = "SU2_CFD_AD" discadj_arina2k.timeout = 8400 discadj_arina2k.tol = 0.00001 @@ -95,7 +95,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.221232, -0.492215, 0.557480, 0.000027] #last 4 columns + discadj_rans_naca0012_sst.test_vals = [-2.221040, -0.492202, 0.557470, 0.000027] discadj_rans_naca0012_sst.su2_exec = "SU2_CFD_AD" discadj_rans_naca0012_sst.timeout = 1600 discadj_rans_naca0012_sst.tol = 0.00001 @@ -110,7 +110,7 @@ def main(): discadj_incomp_NACA0012.cfg_dir = "disc_adj_incomp_euler/naca0012" discadj_incomp_NACA0012.cfg_file = "incomp_NACA0012_disc.cfg" discadj_incomp_NACA0012.test_iter = 20 - discadj_incomp_NACA0012.test_vals = [20.0, -4.092007, -2.652750, 0.0] #last 4 columns + discadj_incomp_NACA0012.test_vals = [20.000000, -4.092007, -2.652750, 0.000000] discadj_incomp_NACA0012.su2_exec = "SU2_CFD_AD" discadj_incomp_NACA0012.timeout = 1600 discadj_incomp_NACA0012.tol = 0.00001 @@ -198,7 +198,7 @@ def main(): discadj_pitchingNACA0012.cfg_dir = "disc_adj_euler/naca0012_pitching" discadj_pitchingNACA0012.cfg_file = "inv_NACA0012_pitching.cfg" discadj_pitchingNACA0012.test_iter = 4 - discadj_pitchingNACA0012.test_vals = [-1.223509, -1.646090, -0.007671, 0.000013] #last 4 columns + discadj_pitchingNACA0012.test_vals = [-1.218846, -1.645199, -0.007645, 0.000013] discadj_pitchingNACA0012.su2_exec = "SU2_CFD_AD" discadj_pitchingNACA0012.timeout = 1600 discadj_pitchingNACA0012.tol = 0.00001 @@ -210,7 +210,7 @@ def main(): unst_deforming_naca0012.cfg_dir = "disc_adj_euler/naca0012_pitching_def" unst_deforming_naca0012.cfg_file = "inv_NACA0012_pitching_deform_ad.cfg" unst_deforming_naca0012.test_iter = 4 - unst_deforming_naca0012.test_vals = [-1.867980, -1.741743, 1090.200000, 0.000006] #last 4 columns + unst_deforming_naca0012.test_vals = [-1.864178, -1.742176, 1090.300000, 0.000006] unst_deforming_naca0012.su2_exec = "SU2_CFD_AD" unst_deforming_naca0012.timeout = 1600 unst_deforming_naca0012.tol = 0.00001 diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index e59f7a9b4763..8a452492b98b 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -49,7 +49,7 @@ def main(): tutorial_inv_bump.cfg_dir = "../Tutorials/compressible_flow/Inviscid_Bump" tutorial_inv_bump.cfg_file = "inv_channel.cfg" tutorial_inv_bump.test_iter = 0 - tutorial_inv_bump.test_vals = [-1.437425, 4.075857, 0.005477, 0.012933] + tutorial_inv_bump.test_vals = [-1.437425, 4.075857, 0.005439, 0.012998] tutorial_inv_bump.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_inv_bump.timeout = 1600 tutorial_inv_bump.tol = 0.00001 @@ -73,7 +73,7 @@ def main(): tutorial_inv_onera.cfg_dir = "../Tutorials/compressible_flow/Inviscid_ONERAM6" tutorial_inv_onera.cfg_file = "inv_ONERAM6.cfg" tutorial_inv_onera.test_iter = 0 - tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.247505, 0.085573] + tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.247451, 0.085770] tutorial_inv_onera.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_inv_onera.timeout = 1600 tutorial_inv_onera.tol = 0.00001 @@ -133,7 +133,7 @@ def main(): tutorial_turb_oneram6.cfg_dir = "../Tutorials/compressible_flow/Turbulent_ONERAM6" tutorial_turb_oneram6.cfg_file = "turb_ONERAM6.cfg" tutorial_turb_oneram6.test_iter = 0 - tutorial_turb_oneram6.test_vals = [-4.564441, -11.523120, 0.333773, 0.098351] + tutorial_turb_oneram6.test_vals = [-4.564441, -11.524277, 0.327954, 0.097349] tutorial_turb_oneram6.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_turb_oneram6.timeout = 1600 tutorial_turb_oneram6.tol = 0.00001 @@ -169,7 +169,7 @@ def main(): propeller_var_load.cfg_dir = "../Tutorials/compressible_flow/ActuatorDisk_VariableLoad" propeller_var_load.cfg_file = "propeller_variable_load.cfg" propeller_var_load.test_iter = 20 - propeller_var_load.test_vals = [-1.837780, -4.535082, -0.000315, 0.171887] + propeller_var_load.test_vals = [-1.830252, -4.535038, -0.000323, 0.171648] propeller_var_load.su2_exec = "mpirun -np 2 SU2_CFD" propeller_var_load.timeout = 3200 propeller_var_load.tol = 0.00001 @@ -182,7 +182,7 @@ def main(): tutorial_design_inv_naca0012.cfg_dir = "../Tutorials/design/Inviscid_2D_Unconstrained_NACA0012" tutorial_design_inv_naca0012.cfg_file = "inv_NACA0012_basic.cfg" tutorial_design_inv_naca0012.test_iter = 0 - tutorial_design_inv_naca0012.test_vals = [-3.585391, -2.989014, 0.135675, 0.209070] + tutorial_design_inv_naca0012.test_vals = [-3.585391, -2.989014, 0.135070, 0.208565] tutorial_design_inv_naca0012.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_design_inv_naca0012.timeout = 1600 tutorial_design_inv_naca0012.tol = 0.00001 @@ -194,7 +194,7 @@ def main(): tutorial_design_turb_rae2822.cfg_dir = "../Tutorials/design/Turbulent_2D_Constrained_RAE2822" tutorial_design_turb_rae2822.cfg_file = "turb_SA_RAE2822.cfg" tutorial_design_turb_rae2822.test_iter = 0 - tutorial_design_turb_rae2822.test_vals = [-1.700114, -4.941261, 0.218432, 0.190639] #last 4 columns + tutorial_design_turb_rae2822.test_vals = [-1.700114, -4.941305, 0.218348, 0.190357] tutorial_design_turb_rae2822.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_design_turb_rae2822.timeout = 1600 tutorial_design_turb_rae2822.tol = 0.00001