diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 949a094969ed..9f170cac9b83 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -182,7 +182,6 @@ class CNumerics { 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 */ @@ -689,17 +688,15 @@ class CNumerics { /*! * \brief Set the value of the second blending function. - * \param[in] val_F1_i - Value of the second Menter blending function at point i. - * \param[in] val_F1_j - Value of the second Menter blending function at point j. + * \param[in] val_F2_i - Value of the second Menter blending function at point i. */ - virtual void SetF2blending(su2double val_F1_i, su2double val_F1_j) {/* empty */}; + virtual void SetF2blending(su2double val_F2_i) {/* empty */}; /*! * \brief Set the value of the cross diffusion for the SST model. * \param[in] val_CDkw_i - Value of the cross diffusion at point i. - * \param[in] val_CDkw_j - Value of the cross diffusion at point j. */ - virtual void SetCrossDiff(su2double val_CDkw_i, su2double val_CDkw_j) {/* empty */}; + virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */}; /*! * \brief Set the gradient of the auxiliary variables. @@ -1425,47 +1422,6 @@ class CNumerics { return ERR; } - /*! - * \brief Set intermittency for numerics (used in SA with LM transition model) - */ - inline virtual void SetIntermittency(su2double intermittency_in) { } - - /*! - * \brief Residual for source term integration. - * \param[in] val_production - Value of the Production. - */ - inline virtual void SetProduction(su2double val_production) { } - - /*! - * \brief Residual for source term integration. - * \param[in] val_destruction - Value of the Destruction. - */ - inline virtual void SetDestruction(su2double val_destruction) { } - - /*! - * \brief Residual for source term integration. - * \param[in] val_crossproduction - Value of the CrossProduction. - */ - inline virtual void SetCrossProduction(su2double val_crossproduction) { } - - /*! - * \brief Residual for source term integration. - * \param[in] val_production - Value of the Production. - */ - inline virtual su2double GetProduction(void) const { return 0; } - - /*! - * \brief Residual for source term integration. - * \param[in] val_destruction - Value of the Destruction. - */ - inline virtual su2double GetDestruction(void) const { return 0; } - - /*! - * \brief Residual for source term integration. - * \param[in] val_crossproduction - Value of the CrossProduction. - */ - inline virtual su2double GetCrossProduction(void) const { return 0; } - /*! * \brief A virtual member. */ diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 522f6055e36d..9e778ea70b60 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -1,8 +1,6 @@ /*! * \file turb_sources.hpp - * \brief Delarations of numerics classes for integration of source - * terms in turbulence problems. - * \author F. Palacios, T. Economon + * \brief Numerics classes for integration of source terms in turbulence problems. * \version 7.3.0 "Blackbird" * * SU2 Project Website: https://su2code.github.io @@ -28,400 +26,658 @@ #pragma once +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../scalar/scalar_sources.hpp" /*! - * \class CSourcePieceWise_TurbSA + * \class CSAVariables + * \brief Structure with SA common auxiliary functions and constants. + */ +struct CSAVariables { + /*--- List of constants ---*/ + const su2double cv1_3 = pow(7.1, 3); + const su2double k2 = pow(0.41, 2); + const su2double cb1 = 0.1355; + const su2double cw2 = 0.3; + const su2double ct3 = 1.2; + const su2double ct4 = 0.5; + const su2double cw3_6 = pow(2, 6); + const su2double sigma = 2.0 / 3.0; + const su2double cb2 = 0.622; + const su2double cb2_sigma = cb2 / sigma; + const su2double cw1 = cb1 / k2 + (1 + cb2) / sigma; + const su2double cr1 = 0.5; + + /*--- List of auxiliary functions ---*/ + su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2; + + /*--- List of helpers ---*/ + su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad, gamma_bc; +}; + +/*! + * \class CSourceBase_TurbSA * \brief Class for integrating the source terms of the Spalart-Allmaras turbulence model equation. - * \ingroup SourceDiscr - * \author A. Bueno. + * The variables that are subject to change in each variation/correction have their own class. + * \note Additional source terms (e.g. compressibility) are implemented as decorators. */ +template class CSourceBase_TurbSA : public CNumerics { -protected: - su2double cv1_3; - su2double k2; - su2double cb1; - su2double cw2; - su2double ct3; - su2double ct4; - su2double cw3_6; - su2double cb2_sigma; - su2double sigma; - su2double cb2; - su2double cw1; - su2double cr1; - + protected: su2double Gamma_BC = 0.0; - su2double intermittency; - su2double Production, Destruction, CrossProduction; + /*--- Residual and Jacobian ---*/ su2double Residual, *Jacobian_i; -private: - su2double Jacobian_Buffer; /// Static storage for the Jacobian (which needs to be pointer for return type). + su2double Jacobian_Buffer; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */ -protected: - const bool incompressible = false, rotating_frame = false; - bool roughwall = false; + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + const bool rotating_frame = false; + const bool transition = false; -public: + public: /*! * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. + * \param[in] nDim - Number of dimensions of the problem. * \param[in] config - Definition of the particular problem. */ - CSourceBase_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); + CSourceBase_TurbSA(unsigned short nDim, unsigned short, const CConfig* config) + : CNumerics(nDim, 1, config), + idx(nDim, config->GetnSpecies()), + rotating_frame(config->GetRotating_Frame()), + transition(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) { + /*--- Setup the Jacobian pointer, we need to return su2double** but we know + * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/ + Jacobian_i = &Jacobian_Buffer; + } /*! - * \brief Residual for source term integration. - * \param[in] intermittency_in - Value of the intermittency. + * \brief Get the intermittency for the BC transition model. + * \return Value of the intermittency. */ - inline void SetIntermittency(su2double intermittency_in) final { intermittency = intermittency_in; } + su2double GetGammaBC() const final { return Gamma_BC; } /*! * \brief Residual for source term integration. - * \param[in] val_production - Value of the Production. + * \param[in] config - Definition of the particular problem. + * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ - inline void SetProduction(su2double val_production) final { Production = val_production; } + ResidualType<> ComputeResidual(const CConfig* config) override { + const auto& density = V_i[idx.Density()]; + const auto& laminar_viscosity = V_i[idx.LaminarViscosity()]; + + AD::StartPreacc(); + AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, ScalarVar_i[0], Volume, dist_i, roughness_i); + AD::SetPreaccIn(Vorticity_i, 3); + AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim); + AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); - /*! - * \brief Residual for source term integration. - * \param[in] val_destruction - Value of the Destruction. - */ - inline void SetDestruction(su2double val_destruction) final { Destruction = val_destruction; } + /*--- Common auxiliary variables and constants of the model. ---*/ + CSAVariables var; - /*! - * \brief Residual for source term integration. - * \param[in] val_crossproduction - Value of the CrossProduction. - */ - inline void SetCrossProduction(su2double val_crossproduction) final { CrossProduction = val_crossproduction; } + Residual = 0.0; + Jacobian_i[0] = 0.0; - /*! - * \brief ______________. - */ - inline su2double GetProduction(void) const final { return Production; } + /*--- Evaluate Omega with a rotational correction term. ---*/ - /*! - * \brief Get the intermittency for the BC trans. model. - * \return Value of the intermittency. - */ - inline su2double GetGammaBC(void) const final { return Gamma_BC; } + Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var); - /*! - * \brief ______________. - */ - inline su2double GetDestruction(void) const final { return Destruction; } + /*--- Dacles-Mariani et. al. rotation correction ("-R"), this is applied by + * default for rotating frame, but should be controled in the config. ---*/ + if (rotating_frame) { + var.Omega += 2.0 * min(0.0, StrainMag_i - var.Omega); + } - /*! - * \brief ______________. - */ - inline su2double GetCrossProduction(void) const final { return CrossProduction; } + if (dist_i > 1e-10) { + /*--- Vorticity ---*/ + var.S = var.Omega; + + var.dist_i_2 = pow(dist_i, 2); + const su2double nu = laminar_viscosity / density; + var.inv_k2_d2 = 1.0 / (var.k2 * var.dist_i_2); + + /*--- Modified values for roughness, roughness_i = 0 for smooth walls and Ji remains the same. + * Ref: Aupoix, B. and Spalart, P. R., "Extensions of the Spalart-Allmaras Turbulence Model to Account for Wall + * Roughness," International Journal of Heat and Fluid Flow, Vol. 24, 2003, pp. 454-462. + * See https://turbmodels.larc.nasa.gov/spalart.html#sarough for detailed explanation. ---*/ + var.Ji = ScalarVar_i[0] / nu + var.cr1 * (roughness_i / (dist_i + EPS)); + var.d_Ji = 1.0 / nu; + + const su2double Ji_2 = pow(var.Ji, 2); + const su2double Ji_3 = Ji_2 * var.Ji; + + var.fv1 = Ji_3 / (Ji_3 + var.cv1_3); + var.d_fv1 = 3 * Ji_2 * var.cv1_3 / (nu * pow(Ji_3 + var.cv1_3, 2)); + + /*--- Using a modified relation so as to not change the Shat that depends on fv2. + * From NASA turb modeling resource and 2003 paper. ---*/ + var.fv2 = 1 - ScalarVar_i[0] / (nu + ScalarVar_i[0] * var.fv1); + var.d_fv2 = -(1 / nu - Ji_2 * var.d_fv1) / pow(1 + var.Ji * var.fv1, 2); + + /*--- Compute ft2 term ---*/ + ft2::get(var); + + /*--- Compute modified vorticity ---*/ + ModVort::get(ScalarVar_i[0], nu, var); + var.inv_Shat = 1.0 / var.Shat; + + /*--- Compute auxiliary function r ---*/ + rFunc::get(ScalarVar_i[0], var); + + var.g = var.r + var.cw2 * (pow(var.r, 6) - var.r); + var.g_6 = pow(var.g, 6); + var.glim = pow((1 + var.cw3_6) / (var.g_6 + var.cw3_6), 1.0 / 6.0); + var.fw = var.g * var.glim; + + var.d_g = var.d_r * (1 + var.cw2 * (6 * pow(var.r, 5) - 1)); + var.d_fw = var.d_g * var.glim * (1 - var.g_6 / (var.g_6 + var.cw3_6)); + + var.norm2_Grad = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0]); + + if (transition) { + /*--- BC transition model (2020 revision). This should only be used with SA-noft2. + * TODO: Consider making this part of the "SourceTerms" template. ---*/ + const su2double chi_1 = 0.002; + const su2double chi_2 = 50.0; + + /*--- Turbulence intensity is u'/U so we multiply by 100 to get percentage. ---*/ + const su2double tu = 100.0 * config->GetTurbulenceIntensity_FreeStream(); + const su2double nu_t = ScalarVar_i[0] * var.fv1; + + const su2double re_v = density * var.dist_i_2 / laminar_viscosity * var.Omega; + const su2double re_theta = re_v / 2.193; + /*--- Menter correlation. ---*/ + const su2double re_theta_t = 803.73 * pow(tu + 0.6067, -1.027); + + const su2double term1 = sqrt(max(re_theta - re_theta_t, 0.0) / (chi_1 * re_theta_t)); + const su2double term2 = sqrt(max((nu_t * chi_2) / nu, 0.0)); + const su2double term_exponential = (term1 + term2); + + Gamma_BC = 1.0 - exp(-term_exponential); + var.gamma_bc = Gamma_BC; + } else { + /*--- Do not modify the production. ---*/ + var.gamma_bc = 1.0; + } + + /*--- Compute production, destruction and cross production and jacobian ---*/ + su2double Production = 0.0, Destruction = 0.0, CrossProduction = 0.0; + SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, CrossProduction, Jacobian_i[0]); + + Residual = (Production - Destruction + CrossProduction) * Volume; + Jacobian_i[0] *= Volume; + } + + AD::SetPreaccOut(Residual); + AD::EndPreacc(); + + return ResidualType<>(&Residual, &Jacobian_i, nullptr); + } }; +namespace detail { + +/* ============================================================================= + * SPALART-ALLMARAS VARIATIONS AND CORRECTIONS + * ============================================================================*/ /*! - * \class CSourcePieceWise_TurbSA - * \brief Class for integrating the source terms of the Spalart-Allmaras turbulence model equation. - * \ingroup SourceDiscr - * \author A. Bueno. + * \brief Strain rate classes. + * \param[in] vorticity: Vorticity array. + * \param[in] nDim: Problem dimension. + * \param[in] velocity_grad: Velocity gradients. + * \param[out] var: Common SA variables struct (to set Omega). */ -template -class CSourcePieceWise_TurbSA final : public CSourceBase_TurbSA { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - - su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; - su2double r, g, g_6, glim, fw; - su2double norm2_Grad; - su2double dfv1, dfv2, dShat; - su2double dr, dg, dfw; - unsigned short iDim; - bool transition; - bool axisymmetric; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CSourcePieceWise_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); +namespace Omega { - /*! - * \brief Residual for source term integration. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) override; +/*! \brief Baseline. */ +struct Bsl { + template + static void get(const su2double* vorticity, unsigned short, const MatrixType&, CSAVariables& var) { + var.Omega = GeometryToolbox::Norm(3, vorticity); + } +}; +/*! \brief Edward. Here Omega is the Strain Rate. */ +struct Edw { + template + static void get(const su2double*, unsigned short nDim, const MatrixType& velocity_grad, CSAVariables& var) { + su2double Sbar = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + for (unsigned short jDim = 0; jDim < nDim; ++jDim) { + Sbar += (velocity_grad[iDim][jDim] + velocity_grad[jDim][iDim]) * velocity_grad[iDim][jDim]; + } + } + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + Sbar -= (2.0 / 3.0) * pow(velocity_grad[iDim][iDim], 2); + } + var.Omega = sqrt(max(Sbar, 0.0)); + } }; +} // namespace Omega /*! - * \class CSourcePieceWise_TurbSA_COMP - * \brief Class for integrating the source terms of the Spalart-Allmaras CC modification turbulence model equation. - * \ingroup SourceDiscr - * \author E.Molina, A. Bueno. - * \version 7.3.0 "Blackbird" + * \brief Classes to set the ft2 term and its derivative. + * \param[in,out] var: Common SA variables struct. */ -template -class CSourcePieceWise_TurbSA_COMP final : public CSourceBase_TurbSA { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - - su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; - su2double r, g, g_6, glim, fw; - su2double norm2_Grad; - su2double dfv1, dfv2, dShat; - su2double dr, dg, dfw; - su2double aux_cc, CompCorrection, c5; - unsigned short iDim, jDim; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CSourcePieceWise_TurbSA_COMP(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); +namespace ft2 { - /*! - * \brief Residual for source term integration. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) override; +/*! \brief No-ft2 term. */ +struct Zero { + static void get(CSAVariables& var) { + var.ft2 = 0.0; + var.d_ft2 = 0.0; + } +}; +/*! \brief Non-zero ft2 term according to the literature. */ +struct Nonzero { + static void get(CSAVariables& var) { + const su2double xsi2 = pow(var.Ji, 2); + var.ft2 = var.ct3 * exp(-var.ct4 * xsi2); + var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; + } }; +} // namespace ft2 /*! - * \class CSourcePieceWise_TurbSA_E - * \brief Class for integrating the source terms of the Spalart-Allmaras Edwards modification turbulence model equation. - * \ingroup SourceDiscr - * \author E.Molina, A. Bueno. - * \version 7.3.0 "Blackbird" + * \brief Classes to compute the modified vorticity (\tilde{S}) and its derivative. + * \param[in] nue: SA variable. + * \param[in] nu: Laminar viscosity. + * \param[in,out] var: Common SA variables struct. */ -template -class CSourcePieceWise_TurbSA_E final : public CSourceBase_TurbSA { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - - su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; - su2double r, g, g_6, glim, fw; - su2double norm2_Grad; - su2double dfv1, dfv2, dShat; - su2double dr, dg, dfw; - su2double Sbar; - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CSourcePieceWise_TurbSA_E(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); +namespace ModVort { + +/*! \brief Baseline. */ +struct Bsl { + static void get(const su2double& nue, const su2double& nu, CSAVariables& var) { + const su2double Sbar = nue * var.fv2 * var.inv_k2_d2; + var.Shat = var.S + Sbar; + var.Shat = max(var.Shat, 1.0e-10); + if (var.Shat <= 1.0e-10) { + var.d_Shat = 0.0; + } else { + var.d_Shat = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2; + } + } +}; - /*! - * \brief Residual for source term integration. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) override; +/*! \brief Edward. */ +struct Edw { + static void get(const su2double& nue, const su2double& nu, CSAVariables& var) { + var.Shat = max(var.S * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16); + var.Shat = max(var.Shat, 1.0e-10); + if (var.Shat <= 1.0e-10) { + var.d_Shat = 0.0; + } else { + var.d_Shat = -var.S * pow(var.Ji, -2) / nu + var.S * var.d_fv1; + } + } +}; + +/*! \brief Negative. */ +struct Neg { + static void get(const su2double& nue, const su2double& nu, CSAVariables& var) { + if (nue > 0.0) { + // Baseline solution + Bsl::get(nue, nu, var); + } else { + var.Shat = 1.0e-10; + var.d_Shat = 0.0; + } + /*--- Don't check whether Sbar <>= -cv2*S. + * Steven R. Allmaras, Forrester T. Johnson and Philippe R. Spalart - + * "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model" eq. 12 + * No need for Sbar ---*/ + } }; +} // namespace ModVort /*! - * \class CSourcePieceWise_TurbSA_E_COMP - * \brief Class for integrating the source terms of the Spalart-Allmaras Edwards modification with CC turbulence model equation. - * \ingroup SourceDiscr - * \author E.Molina, A. Bueno. - * \version 7.3.0 "Blackbird" + * \brief Auxiliary function r and its derivative. + * \param[in] nue: SA variable. + * \param[in,out] var: Common SA variables struct. */ -template -class CSourcePieceWise_TurbSA_E_COMP : public CSourceBase_TurbSA { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - - su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; - su2double r, g, g_6, glim, fw; - su2double norm2_Grad; - su2double dfv1, dfv2, dShat; - su2double dr, dg, dfw; - su2double Sbar; - su2double aux_cc, CompCorrection, c5; - -public: +namespace r { + +/*! \brief Baseline. */ +struct Bsl { + static void get(const su2double& nue, CSAVariables& var) { + var.r = min(nue * var.inv_Shat * var.inv_k2_d2, 10.0); + var.d_r = (var.Shat - nue * var.d_Shat) * pow(var.inv_Shat, 2) * var.inv_k2_d2; + if (var.r >= 10.0) var.d_r = 0.0; + } +}; + +/*! \brief Edward. */ +struct Edw { + static void get(const su2double& nue, CSAVariables& var) { + var.r = min(nue * var.inv_Shat * var.inv_k2_d2, 10.0); + var.r = tanh(var.r) / tanh(1.0); + + var.d_r = (var.Shat - nue * var.d_Shat) * pow(var.inv_Shat, 2) * var.inv_k2_d2; + var.d_r = (1 - pow(tanh(var.r), 2.0)) * (var.d_r) / tanh(1.0); + } +}; +} // namespace r + +/*! + * \brief Source terms classes: production, destruction and cross-productions term and their derivative. + * \param[in] nue: SA variable. + * \param[in] var: Common SA variables struct. + * \param[out] production: Production term. + * \param[out] destruction: Destruction term. + * \param[out] cross_production: CrossProduction term. + * \param[out] jacobian: Derivative of the combined source term wrt nue. + */ +namespace SourceTerms { + +/*! \brief Baseline (Original SA model). */ +struct Bsl { + static void get(const su2double& nue, const CSAVariables& var, su2double& production, su2double& destruction, + su2double& cross_production, su2double& jacobian) { + ComputeProduction(nue, var, production, jacobian); + ComputeDestruction(nue, var, destruction, jacobian); + ComputeCrossProduction(nue, var, cross_production, jacobian); + } + + static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production, + su2double& jacobian) { + const su2double factor = var.gamma_bc * var.cb1; + production = factor * (1.0 - var.ft2) * var.Shat * nue; + jacobian += factor * (-var.Shat * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Shat)); + } + + static void ComputeDestruction(const su2double& nue, const CSAVariables& var, su2double& destruction, + su2double& jacobian) { + const su2double cb1_k2 = var.cb1 / var.k2; + const su2double factor = var.cw1 * var.fw - cb1_k2 * var.ft2; + destruction = factor * pow(nue, 2) / var.dist_i_2; + jacobian -= ((var.cw1 * var.d_fw - cb1_k2 * var.d_ft2) * pow(nue, 2) + factor * 2 * nue) / var.dist_i_2; + } + + static void ComputeCrossProduction(const su2double& nue, const CSAVariables& var, su2double& cross_production, + su2double&) { + cross_production = var.cb2_sigma * var.norm2_Grad; + /*--- No contribution to the jacobian. ---*/ + } +}; + +/*! \brief Negative. */ +struct Neg { + static void get(const su2double& nue, const CSAVariables& var, su2double& production, su2double& destruction, + su2double& cross_production, su2double& jacobian) { + if (nue > 0.0) { + Bsl::get(nue, var, production, destruction, cross_production, jacobian); + } else { + ComputeProduction(nue, var, production, jacobian); + ComputeDestruction(nue, var, destruction, jacobian); + ComputeCrossProduction(nue, var, cross_production, jacobian); + } + } + + static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production, + su2double& jacobian) { + const su2double dP_dnu = var.gamma_bc * var.cb1 * (1.0 - var.ct3) * var.S; + production = dP_dnu * nue; + jacobian += dP_dnu; + } + + static void ComputeDestruction(const su2double& nue, const CSAVariables& var, su2double& destruction, + su2double& jacobian) { + /*--- The destruction when nue < 0 is added instead of the usual subtraction, hence the negative sign. ---*/ + const su2double dD_dnu = -var.cw1 * nue / var.dist_i_2; + destruction = dD_dnu * nue; + jacobian -= 2 * dD_dnu; + } + + static void ComputeCrossProduction(const su2double& nue, const CSAVariables& var, su2double& cross_production, + su2double& jacobian) { + Bsl::ComputeCrossProduction(nue, var, cross_production, jacobian); + } +}; +} // namespace SourceTerms + +/* ============================================================================= + * SPALART-ALLMARAS ADDITIONAL SOURCE TERMS DECORATORS + * ============================================================================*/ + +/*! + * \class CCompressibilityCorrection + * \brief Mixing Layer Compressibility Correction (SA-comp). + */ +template +class CCompressibilityCorrection final : public ParentClass { + private: + using ParentClass::Gamma; + using ParentClass::idx; + using ParentClass::nDim; + using ParentClass::PrimVar_Grad_i; + using ParentClass::ScalarVar_i; + using ParentClass::V_i; + using ParentClass::Volume; + + using ResidualType = typename ParentClass::template ResidualType<>; + + const su2double c5 = 3.5; + + public: /*! * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. + * \param[in] nDim - Number of dimensions of the problem. * \param[in] config - Definition of the particular problem. */ - CSourcePieceWise_TurbSA_E_COMP(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); + CCompressibilityCorrection(unsigned short nDim, unsigned short, const CConfig* config) + : ParentClass(nDim, 0, config) {} /*! * \brief Residual for source term integration. * \param[in] config - Definition of the particular problem. * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ - ResidualType<> ComputeResidual(const CConfig* config) override; + ResidualType ComputeResidual(const CConfig* config) override { + /*--- Residual from standard SA ---*/ + ParentClass::ComputeResidual(config); + + /*--- Compressibility Correction term ---*/ + const auto& pressure = V_i[idx.Pressure()]; + const auto& density = V_i[idx.Density()]; + const su2double sound_speed = sqrt(pressure * Gamma / density); + su2double aux_cc = 0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + for (unsigned short jDim = 0; jDim < nDim; ++jDim) { + aux_cc += pow(PrimVar_Grad_i[idx.Velocity() + iDim][jDim], 2); + } + } + const su2double d_CompCorrection = 2.0 * c5 * ScalarVar_i[0] / pow(sound_speed, 2) * aux_cc * Volume; + const su2double CompCorrection = 0.5 * ScalarVar_i[0] * d_CompCorrection; + + this->Residual -= CompCorrection; + this->Jacobian_i[0] -= d_CompCorrection; + + return ResidualType(&this->Residual, &this->Jacobian_i, nullptr); + } }; +} // namespace detail + +/* ============================================================================= + * SPALART-ALLMARAS CLASSES + * ============================================================================*/ + +/// TODO: Factory method to create combinations of the different variations based on the config. +/// See PR #1413, the combinations exposed should follow https://turbmodels.larc.nasa.gov/spalart.html + /*! - * \class CSourcePieceWise_TurbSA_Neg + * \class CSourcePieceWise_TurbSA * \brief Class for integrating the source terms of the Spalart-Allmaras turbulence model equation. - * \ingroup SourceDiscr - * \author F. Palacios */ template -class CSourcePieceWise_TurbSA_Neg : public CSourceBase_TurbSA { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ +using CSourcePieceWise_TurbSA = CSourceBase_TurbSA; - su2double nu, Ji, fv1, fv2, ft2, Omega, S, Shat, inv_Shat, dist_i_2, Ji_2, Ji_3, inv_k2_d2; - su2double r, g, g_6, glim, fw; - su2double norm2_Grad; - su2double dfv1, dfv2, dShat; - su2double dr, dg, dfw; +/*! + * \class CSourcePieceWise_TurbSA_COMP + * \brief Class for integrating the source terms of the Spalart-Allmaras model with compressibility correction. + */ +template +using CSourcePieceWise_TurbSA_COMP = detail::CCompressibilityCorrection>; -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CSourcePieceWise_TurbSA_Neg(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); +/*! + * \class CSourcePieceWise_TurbSA_E + * \brief Class for integrating the source terms of the Spalart-Allmaras model with Edwards modification. + * \note From NASA Turbulence model site. http://turbmodels.larc.nasa.gov/spalart.html + * This form was developed primarily to improve the near-wall numerical behavior of the model (i.e. the goal was to + * improve the convergence behavior). The reference is: Edwards, J. R. and Chandra, S. "Comparison of Eddy + * Viscosity-Transport Turbulence Models for Three-Dimensional, Shock-Separated Flowfields," AIAA Journal, Vol. 34, No. + * 4, 1996, pp. 756-763. + */ +template +using CSourcePieceWise_TurbSA_E = CSourceBase_TurbSA; - /*! - * \brief Residual for source term integration. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) override; +/*! + * \class CSourcePieceWise_TurbSA_E_COMP + * \brief Class for integrating the source terms of the Spalart-Allmaras model with Edwards modification + * and compressibility correction. + */ +template +using CSourcePieceWise_TurbSA_E_COMP = detail::CCompressibilityCorrection>; -}; +/*! + * \class CSourcePieceWise_TurbSA_Neg + * \brief Class for integrating the source terms of the negative Spalart-Allmaras model. + */ +template +using CSourcePieceWise_TurbSA_Neg = CSourceBase_TurbSA; /*! * \class CSourcePieceWise_TurbSST * \brief Class for integrating the source terms of the Menter SST turbulence model equations. - * \ingroup SourceDiscr - * \author A. Campos. */ template class CSourcePieceWise_TurbSST final : public CNumerics { -private: - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - - su2double F1_i, - F1_j, - F2_i, - F2_j; - - su2double alfa_1, - alfa_2, - beta_1, - beta_2, - sigma_k_1, - sigma_k_2, - sigma_w_1, - sigma_w_2, - beta_star, - a1; + private: + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + const bool sustaining_terms = false; + const bool axisymmetric = false; - su2double CDkw_i, CDkw_j; + /*--- Closure constants ---*/ + const su2double sigma_k_1, sigma_k_2, sigma_w_1, sigma_w_2, beta_1, beta_2, beta_star, a1, alfa_1, alfa_2; - su2double kAmb, omegaAmb; + /*--- Ambient values for SST-SUST. ---*/ + const su2double kAmb, omegaAmb; + su2double F1_i, F2_i, CDkw_i; su2double Residual[2]; su2double* Jacobian_i[2]; - su2double Jacobian_Buffer[4]; /// Static storage for the Jacobian (which needs to be pointer for return type). - - bool incompressible; - bool sustaining_terms; - bool axisymmetric; + su2double Jacobian_Buffer[4]; /// Static storage for the Jacobian (which needs to be pointer for return type). /*! - * \brief A virtual member. Get strain magnitude based on perturbed reynolds stress matrix - * \param[in] turb_ke: turbulent kinetic energy of the node + * \brief Get strain magnitude based on perturbed reynolds stress matrix. + * \param[in] turb_ke: turbulent kinetic energy of the node. */ - void SetPerturbedStrainMag(su2double turb_ke); + inline su2double PerturbedStrainMag(su2double turb_ke) const { + /*--- Compute norm of perturbed strain rate tensor. ---*/ + + su2double perturbedStrainMag = 0; + 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); + + perturbedStrainMag += pow(StrainRate_ij, 2.0); + } + } + return sqrt(2.0 * perturbedStrainMag); + } /*! * \brief Add contribution due to axisymmetric formulation to 2D residual */ inline void ResidualAxisymmetric(su2double alfa_blended, su2double zeta) { - if (Coord_i[1] < EPS) return; - su2double yinv, rhov, k, w; - su2double sigma_k_i, sigma_w_i; - su2double pk_axi, pw_axi, cdk_axi, cdw_axi; - AD::SetPreaccIn(Coord_i[1]); - yinv = 1.0/Coord_i[1]; - rhov = Density_i*V_i[2]; - k = ScalarVar_i[0]; - w = ScalarVar_i[1]; + const su2double yinv = 1.0 / Coord_i[1]; + const su2double rhov = Density_i * V_i[2]; + const su2double& k = ScalarVar_i[0]; + const su2double& w = ScalarVar_i[1]; /*--- Compute blended constants ---*/ - sigma_k_i = F1_i*sigma_k_1+(1.0-F1_i)*sigma_k_2; - sigma_w_i = F1_i*sigma_w_1+(1.0-F1_i)*sigma_w_2; + const su2double sigma_k_i = F1_i * sigma_k_1 + (1.0 - F1_i) * sigma_k_2; + const su2double sigma_w_i = F1_i * sigma_w_1 + (1.0 - F1_i) * sigma_w_2; /*--- Production ---*/ - pk_axi = max(0.0,2.0/3.0*rhov*k*((2.0*yinv*V_i[2]-PrimVar_Grad_i[2][1]-PrimVar_Grad_i[1][0])/zeta-1.0)); - pw_axi = alfa_blended*zeta/k*pk_axi; + const su2double pk_axi = max( + 0.0, 2.0 / 3.0 * rhov * k * ((2.0 * yinv * V_i[2] - PrimVar_Grad_i[2][1] - PrimVar_Grad_i[1][0]) / zeta - 1.0)); + const su2double pw_axi = alfa_blended * zeta / k * pk_axi; /*--- Convection-Diffusion ---*/ - cdk_axi = rhov*k-(Laminar_Viscosity_i+sigma_k_i*Eddy_Viscosity_i)*ScalarVar_Grad_i[0][1]; - cdw_axi = rhov*w-(Laminar_Viscosity_i+sigma_w_i*Eddy_Viscosity_i)*ScalarVar_Grad_i[1][1]; + const su2double cdk_axi = rhov * k - (Laminar_Viscosity_i + sigma_k_i * Eddy_Viscosity_i) * ScalarVar_Grad_i[0][1]; + const su2double cdw_axi = rhov * w - (Laminar_Viscosity_i + sigma_w_i * Eddy_Viscosity_i) * ScalarVar_Grad_i[1][1]; /*--- Add terms to the residuals ---*/ - Residual[0] += yinv*Volume*(pk_axi-cdk_axi); - Residual[1] += yinv*Volume*(pw_axi-cdw_axi); - + Residual[0] += yinv * Volume * (pk_axi - cdk_axi); + Residual[1] += yinv * Volume * (pw_axi - cdw_axi); } -public: + public: /*! * \brief Constructor of the class. * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. + * \param[in] constants - SST model constants. + * \param[in] val_kine_Inf - Freestream k, for SST with sustaining terms. + * \param[in] val_omega_Inf - Freestream w, for SST with sustaining terms. * \param[in] config - Definition of the particular problem. */ - CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const su2double* constants, - su2double val_kine_Inf, su2double val_omega_Inf, const CConfig* config); + CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short, const su2double* constants, su2double val_kine_Inf, + su2double val_omega_Inf, const CConfig* config) + : CNumerics(val_nDim, 2, config), + idx(val_nDim, config->GetnSpecies()), + sustaining_terms(config->GetKind_Turb_Model() == TURB_MODEL::SST_SUST), + axisymmetric(config->GetAxisymmetric()), + sigma_k_1(constants[0]), + sigma_k_2(constants[1]), + sigma_w_1(constants[2]), + sigma_w_2(constants[3]), + beta_1(constants[4]), + beta_2(constants[5]), + beta_star(constants[6]), + a1(constants[7]), + alfa_1(constants[8]), + alfa_2(constants[9]), + kAmb(val_kine_Inf), + omegaAmb(val_omega_Inf) { + /*--- "Allocate" the Jacobian using the static buffer. ---*/ + Jacobian_i[0] = Jacobian_Buffer; + Jacobian_i[1] = Jacobian_Buffer + 2; + } /*! * \brief Set the value of the first blending function. * \param[in] val_F1_i - Value of the first blending function at point i. - * \param[in] val_F1_j - Value of the first blending function at point j. + * \param[in] Not used. */ - inline void SetF1blending(su2double val_F1_i, su2double val_F1_j) override { + inline void SetF1blending(su2double val_F1_i, su2double) override { F1_i = val_F1_i; - F1_j = val_F1_j; } /*! * \brief Set the value of the second blending function. * \param[in] val_F2_i - Value of the second blending function at point i. - * \param[in] val_F2_j - Value of the second blending function at point j. */ - inline void SetF2blending(su2double val_F2_i, su2double val_F2_j) override { + inline void SetF2blending(su2double val_F2_i) override { F2_i = val_F2_i; - F2_j = val_F2_j; } /*! * \brief Set the value of the cross diffusion for the SST model. * \param[in] val_CDkw_i - Value of the cross diffusion at point i. - * \param[in] val_CDkw_j - Value of the cross diffusion at point j. */ - inline void SetCrossDiff(su2double val_CDkw_i, su2double val_CDkw_j) override { + inline void SetCrossDiff(su2double val_CDkw_i) override { CDkw_i = val_CDkw_i; - CDkw_j = val_CDkw_j; } /*! @@ -429,6 +685,103 @@ class CSourcePieceWise_TurbSST final : public CNumerics { * \param[in] config - Definition of the particular problem. * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ - ResidualType<> ComputeResidual(const CConfig* config) override; + ResidualType<> ComputeResidual(const CConfig* config) override { + AD::StartPreacc(); + AD::SetPreaccIn(StrainMag_i); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(Volume); + AD::SetPreaccIn(dist_i); + AD::SetPreaccIn(F1_i); + AD::SetPreaccIn(F2_i); + AD::SetPreaccIn(CDkw_i); + AD::SetPreaccIn(PrimVar_Grad_i, nDim + idx.Velocity(), nDim); + AD::SetPreaccIn(Vorticity_i, 3); + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; + Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; + + Residual[0] = 0.0; + Residual[1] = 0.0; + Jacobian_i[0][0] = 0.0; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = 0.0; + + /*--- Computation of blended constants for the source terms ---*/ + + const su2double alfa_blended = F1_i * alfa_1 + (1.0 - F1_i) * alfa_2; + const su2double beta_blended = F1_i * beta_1 + (1.0 - F1_i) * beta_2; + + if (dist_i > 1e-10) { + /*--- Production ---*/ + + su2double diverg = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + diverg += PrimVar_Grad_i[iDim + idx.Velocity()][iDim]; + + /*--- If using UQ methodolgy, calculate production using perturbed Reynolds stress matrix ---*/ + + su2double StrainMag = StrainMag_i; + + if (using_uq) { + ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, PrimVar_Grad_i + idx.Velocity(), + Density_i, Eddy_Viscosity_i, ScalarVar_i[0], MeanPerturbedRSM); + StrainMag = PerturbedStrainMag(ScalarVar_i[0]); + } + + su2double pk = Eddy_Viscosity_i * pow(StrainMag, 2) - 2.0 / 3.0 * Density_i * ScalarVar_i[0] * diverg; + pk = max(0.0, min(pk, 20.0 * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0])); + + const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity_i); + const su2double zeta = max(ScalarVar_i[1], VorticityMag * F2_i / a1); + su2double pw = alfa_blended * Density_i * max(pow(StrainMag, 2) - 2.0 / 3.0 * zeta * diverg, 0.0); + + /*--- Sustaining terms, if desired. Note that if the production terms are + larger equal than the sustaining terms, the original formulation is + obtained again. This is in contrast to the version in literature + where the sustaining terms are simply added. This latter approach could + lead to problems for very big values of the free-stream turbulence + intensity. ---*/ + + if (sustaining_terms) { + const su2double sust_k = beta_star * Density_i * kAmb * omegaAmb; + const su2double sust_w = beta_blended * Density_i * omegaAmb * omegaAmb; + pk = max(pk, sust_k); + pw = max(pw, sust_w); + } + + /*--- Add the production terms to the residuals. ---*/ + + Residual[0] += pk * Volume; + Residual[1] += pw * Volume; + + /*--- Dissipation ---*/ + + Residual[0] -= beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Volume; + Residual[1] -= beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * Volume; + + /*--- Cross diffusion ---*/ + + Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; + + /*--- Contribution due to 2D axisymmetric formulation ---*/ + + if (axisymmetric) ResidualAxisymmetric(alfa_blended, zeta); + + /*--- Implicit part ---*/ + + Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume; + Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume; + } + + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(Residual, Jacobian_i, nullptr); + } }; diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 7af3c2d9e585..2b031040559c 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -108,7 +108,6 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/numerics/continuous_adjoint/adj_sources.cpp \ ../src/numerics/scalar/scalar_sources.cpp \ ../src/numerics/species/species_sources.cpp \ - ../src/numerics/turbulent/turb_sources.cpp \ ../src/numerics/elasticity/CFEAElasticity.cpp \ ../src/numerics/elasticity/CFEALinearElasticity.cpp \ ../src/numerics/elasticity/CFEANonlinearElasticity.cpp \ diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index a99192bc4ecb..39934aa7926e 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -141,7 +141,6 @@ su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/continuous_adjoint/adj_diffusion.cpp', 'numerics/continuous_adjoint/adj_sources.cpp', 'numerics/scalar/scalar_sources.cpp', - 'numerics/turbulent/turb_sources.cpp', 'numerics/species/species_sources.cpp', 'numerics/elasticity/CFEAElasticity.cpp', 'numerics/elasticity/CFEALinearElasticity.cpp', diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index dbad754d5d0e..dc9a68d9082e 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -131,15 +131,16 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, * parts of tau can be computed with the total viscosity. --- */ if (using_uq){ - ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity); // laminar part + // laminar part + ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity); // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) tau[iDim][jDim] += (-Density) * MeanPerturbedRSM[iDim][jDim]; } else { - // 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, su2double(0.0)); // TODO why ignore turb_ke? + // turb_ke is not considered in the stress tensor, see #797 + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0)); } } diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp deleted file mode 100644 index 31e397f20a40..000000000000 --- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp +++ /dev/null @@ -1,952 +0,0 @@ -/*! - * \file turb_sources.cpp - * \brief Implementation of numerics classes for integration of - * turbulence source-terms. - * \author F. Palacios, T. Economon - * \version 7.3.0 "Blackbird" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) - * - * SU2 is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * SU2 is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with SU2. If not, see . - */ - -#include "../../../include/numerics/turbulent/turb_sources.hpp" - -#include "../../../include/variables/CEulerVariable.hpp" -#include "../../../include/variables/CIncEulerVariable.hpp" -#include "../../../include/variables/CNEMOEulerVariable.hpp" - -CSourceBase_TurbSA::CSourceBase_TurbSA(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CNumerics(val_nDim, val_nVar, config), - incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), - rotating_frame(config->GetRotating_Frame()) -{ - /*--- Spalart-Allmaras closure constants ---*/ - - cv1_3 = pow(7.1, 3.0); - k2 = pow(0.41, 2.0); - cb1 = 0.1355; - cw2 = 0.3; - ct3 = 1.2; - ct4 = 0.5; - cw3_6 = pow(2.0, 6.0); - sigma = 2./3.; - cb2 = 0.622; - cb2_sigma = cb2/sigma; - cw1 = cb1/k2+(1.0+cb2)/sigma; - cr1 = 0.5; - - /*--- Setup the Jacobian pointer, we need to return su2double** but - * we know the Jacobian is 1x1 so we use this "trick" to avoid - * having to dynamically allocate. ---*/ - - Jacobian_i = &Jacobian_Buffer; - -} - -template -CSourcePieceWise_TurbSA::CSourcePieceWise_TurbSA(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()) { - - transition = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC); -} - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig* config) { - -// AD::StartPreacc(); -// AD::SetPreaccIn(V_i, nDim+6); -// AD::SetPreaccIn(Vorticity_i, nDim); -// AD::SetPreaccIn(StrainMag_i); -// AD::SetPreaccIn(ScalarVar_i[0]); -// AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); -// AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - - // Set the boolean here depending on whether the point is closest to a rough wall or not. - roughwall = (roughness_i > 0.0); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - - Residual = 0.0; - Production = 0.0; - Destruction = 0.0; - CrossProduction = 0.0; - Jacobian_i[0] = 0.0; - - /*--- Evaluate Omega ---*/ - - Omega = sqrt(Vorticity_i[0]*Vorticity_i[0] + Vorticity_i[1]*Vorticity_i[1] + Vorticity_i[2]*Vorticity_i[2]); - - /*--- Rotational correction term ---*/ - - if (rotating_frame) { Omega += 2.0*min(0.0, StrainMag_i-Omega); } - - if (dist_i > 1e-10) { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - nu = Laminar_Viscosity_i/Density_i; - - /*--- Modified values for roughness ---*/ - /*--- Ref: Aupoix, B. and Spalart, P. R., "Extensions of the Spalart-Allmaras Turbulence Model to Account for Wall Roughness," - * International Journal of Heat and Fluid Flow, Vol. 24, 2003, pp. 454-462. ---*/ - /* --- See https://turbmodels.larc.nasa.gov/spalart.html#sarough for detailed explanation. ---*/ - - Ji = ScalarVar_i[0]/nu + cr1*(roughness_i/(dist_i+EPS)); //roughness_i = 0 for smooth walls and Ji remains the same, changes only if roughness is specified. - Ji_2 = Ji*Ji; - Ji_3 = Ji_2*Ji; - fv1 = Ji_3/(Ji_3+cv1_3); - - /*--- Using a modified relation so as to not change the Shat that depends on fv2. ---*/ - fv2 = 1.0 - ScalarVar_i[0]/(nu+ScalarVar_i[0]*fv1); // From NASA turb modeling resource and 2003 paper - - ft2 = ct3*exp(-ct4*Ji_2); - S = Omega; - inv_k2_d2 = 1.0/(k2*dist_i_2); - - Shat = S + ScalarVar_i[0]*fv2*inv_k2_d2; - Shat = max(Shat, 1.0e-10); - inv_Shat = 1.0/Shat; - -// Original SA model -// Production = cb1*(1.0-ft2)*Shat*ScalarVar_i[0]*Volume; - - if (transition) { - - /*--- BC model constants (2020 revision). ---*/ - const su2double chi_1 = 0.002; - const su2double chi_2 = 50.0; - - /*--- turbulence intensity is u'/U so we multiply by 100 to get percentage ---*/ - su2double tu = 100.0 * config->GetTurbulenceIntensity_FreeStream(); - - su2double nu_t = (ScalarVar_i[0]*fv1); //S-A variable - - su2double re_v = ((Density_i*pow(dist_i,2.))/(Laminar_Viscosity_i))*Omega; - su2double re_theta = re_v/2.193; - su2double re_theta_t = (803.73 * pow((tu + 0.6067),-1.027)); //MENTER correlation - //re_theta_t = 163.0 + exp(6.91-tu); //ABU-GHANNAM & SHAW correlation - - su2double term1 = sqrt(max(re_theta-re_theta_t,0.)/(chi_1*re_theta_t)); - su2double term2 = sqrt(max((nu_t*chi_2)/nu,0.)); - su2double term_exponential = (term1 + term2); - - Gamma_BC = 1.0 - exp(-term_exponential); - - Production = Gamma_BC*cb1*Shat*ScalarVar_i[0]*Volume; - } - else { - Production = cb1*Shat*ScalarVar_i[0]*Volume; - } - - /*--- Destruction term ---*/ - - r = min(ScalarVar_i[0]*inv_Shat*inv_k2_d2,10.0); - g = r + cw2*(pow(r,6.0)-r); - g_6 = pow(g,6.0); - glim = pow((1.0+cw3_6)/(g_6+cw3_6),1.0/6.0); - fw = g*glim; - -// Original SA model -// Destruction = (cw1*fw-cb1*ft2/k2)*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - Destruction = cw1*fw*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production - Destruction + CrossProduction; - - /*--- Implicit part, production term ---*/ - - dfv1 = 3.0*Ji_2*cv1_3/(nu*pow(Ji_3+cv1_3,2.)); - dfv2 = -(1/nu-Ji_2*dfv1)/pow(1.+Ji*fv1,2.); - if ( Shat <= 1.0e-10 ) dShat = 0.0; - else dShat = (fv2+ScalarVar_i[0]*dfv2)*inv_k2_d2; - - if (transition) { - Jacobian_i[0] += Gamma_BC*cb1*(ScalarVar_i[0]*dShat+Shat)*Volume; - } - else { - Jacobian_i[0] += cb1*(ScalarVar_i[0]*dShat+Shat)*Volume; - } - - /*--- Implicit part, destruction term ---*/ - - dr = (Shat-ScalarVar_i[0]*dShat)*inv_Shat*inv_Shat*inv_k2_d2; - if (r == 10.0) dr = 0.0; - dg = dr*(1.+cw2*(6.0*pow(r,5.0)-1.0)); - dfw = dg*glim*(1.-g_6/(g_6+cw3_6)); - Jacobian_i[0] -= cw1*(dfw*ScalarVar_i[0] + 2.0*fw)*ScalarVar_i[0]/dist_i_2*Volume; - - } - -// AD::SetPreaccOut(Residual); -// AD::EndPreacc(); - - return ResidualType<>(&Residual, &Jacobian_i, nullptr); - -} - -template -CSourcePieceWise_TurbSA_COMP::CSourcePieceWise_TurbSA_COMP(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()), - c5(3.5) { } - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_COMP::ComputeResidual(const CConfig* config) { - - // AD::StartPreacc(); - // AD::SetPreaccIn(V_i, nDim+6); - // AD::SetPreaccIn(Vorticity_i, nDim); - // AD::SetPreaccIn(StrainMag_i); - // AD::SetPreaccIn(ScalarVar_i[0]); - // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); - // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - - Residual = 0.0; - Production = 0.0; - Destruction = 0.0; - CrossProduction = 0.0; - Jacobian_i[0] = 0.0; - - /*--- Evaluate Omega ---*/ - - Omega = sqrt(Vorticity_i[0]*Vorticity_i[0] + Vorticity_i[1]*Vorticity_i[1] + Vorticity_i[2]*Vorticity_i[2]); - - /*--- Rotational correction term ---*/ - - if (rotating_frame) { Omega += 2.0*min(0.0, StrainMag_i-Omega); } - - if (dist_i > 1e-10) { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - nu = Laminar_Viscosity_i/Density_i; - Ji = ScalarVar_i[0]/nu; - Ji_2 = Ji*Ji; - Ji_3 = Ji_2*Ji; - fv1 = Ji_3/(Ji_3+cv1_3); - fv2 = 1.0 - Ji/(1.0+Ji*fv1); - ft2 = ct3*exp(-ct4*Ji_2); - S = Omega; - inv_k2_d2 = 1.0/(k2*dist_i_2); - - Shat = S + ScalarVar_i[0]*fv2*inv_k2_d2; - Shat = max(Shat, 1.0e-10); - inv_Shat = 1.0/Shat; - - /*--- Production term ---*/; - - Production = cb1*Shat*ScalarVar_i[0]*Volume; - - /*--- Destruction term ---*/ - - r = min(ScalarVar_i[0]*inv_Shat*inv_k2_d2,10.0); - g = r + cw2*(pow(r,6.0)-r); - g_6 = pow(g,6.0); - glim = pow((1.0+cw3_6)/(g_6+cw3_6),1.0/6.0); - fw = g*glim; - - Destruction = cw1*fw*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production - Destruction + CrossProduction; - - /*--- Compressibility Correction term ---*/ - Pressure_i = V_i[idx.Pressure()]; - SoundSpeed_i = sqrt(Pressure_i*Gamma/Density_i); - aux_cc=0; - for(iDim=0;iDim(&Residual, &Jacobian_i, nullptr); - -} - -template -CSourcePieceWise_TurbSA_E::CSourcePieceWise_TurbSA_E(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()) { } - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E::ComputeResidual(const CConfig* config) { - - unsigned short iDim, jDim; - - // AD::StartPreacc(); - // AD::SetPreaccIn(V_i, nDim+6); - // AD::SetPreaccIn(Vorticity_i, nDim); - // AD::SetPreaccIn(StrainMag_i); - // AD::SetPreaccIn(ScalarVar_i[0]); - // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); - // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - - Residual = 0.0; - Production = 0.0; - Destruction = 0.0; - CrossProduction = 0.0; - Jacobian_i[0] = 0.0; - - /* - From NASA Turbulence model site. http://turbmodels.larc.nasa.gov/spalart.html - This form was developed primarily to improve the near-wall numerical behavior of the model (i.e., the goal was to improve the convergence behavior). The reference is: - Edwards, J. R. and Chandra, S. "Comparison of Eddy Viscosity-Transport Turbulence Models for Three-Dimensional, Shock-Separated Flowfields," AIAA Journal, Vol. 34, No. 4, 1996, pp. 756-763. - In this modificaton Omega is replaced by Strain Rate - */ - - /*--- Evaluate Omega, here Omega is the Strain Rate ---*/ - - Sbar = 0.0; - for(iDim=0;iDim 1e-10) { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - nu = Laminar_Viscosity_i/Density_i; - Ji = ScalarVar_i[0]/nu; - Ji_2 = Ji*Ji; - Ji_3 = Ji_2*Ji; - fv1 = Ji_3/(Ji_3+cv1_3); - fv2 = 1.0 - Ji/(1.0+Ji*fv1); - ft2 = ct3*exp(-ct4*Ji_2); - S = Omega; - inv_k2_d2 = 1.0/(k2*dist_i_2); - - //Shat = S + ScalarVar_i[0]*fv2*inv_k2_d2; - Shat = max(S*((1.0/max(Ji,1.0e-16))+fv1),1.0e-16); - - Shat = max(Shat, 1.0e-10); - inv_Shat = 1.0/Shat; - - /*--- Production term ---*/; - - Production = cb1*Shat*ScalarVar_i[0]*Volume; - - /*--- Destruction term ---*/ - - r = min(ScalarVar_i[0]*inv_Shat*inv_k2_d2,10.0); - r=tanh(r)/tanh(1.0); - - g = r + cw2*(pow(r,6.0)-r); - g_6 = pow(g,6.0); - glim = pow((1.0+cw3_6)/(g_6+cw3_6),1.0/6.0); - fw = g*glim; - - Destruction = cw1*fw*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production - Destruction + CrossProduction; - - /*--- Implicit part, production term ---*/ - - dfv1 = 3.0*Ji_2*cv1_3/(nu*pow(Ji_3+cv1_3,2.)); - dfv2 = -(1/nu-Ji_2*dfv1)/pow(1.+Ji*fv1,2.); - - if ( Shat <= 1.0e-10 ) dShat = 0.0; - else dShat = -S*pow(Ji,-2.0)/nu + S*dfv1; - Jacobian_i[0] += cb1*(ScalarVar_i[0]*dShat+Shat)*Volume; - - /*--- Implicit part, destruction term ---*/ - - dr = (Shat-ScalarVar_i[0]*dShat)*inv_Shat*inv_Shat*inv_k2_d2; - dr=(1-pow(tanh(r),2.0))*(dr)/tanh(1.0); - dg = dr*(1.+cw2*(6.0*pow(r,5.0)-1.0)); - dfw = dg*glim*(1.-g_6/(g_6+cw3_6)); - Jacobian_i[0] -= cw1*(dfw*ScalarVar_i[0] + 2.0*fw)*ScalarVar_i[0]/dist_i_2*Volume; - - } - - // AD::SetPreaccOut(Residual); - // AD::EndPreacc(); - - return ResidualType<>(&Residual, &Jacobian_i, nullptr); - -} - -template -CSourcePieceWise_TurbSA_E_COMP::CSourcePieceWise_TurbSA_E_COMP(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()) { } - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_E_COMP::ComputeResidual(const CConfig* config) { - - unsigned short iDim, jDim; - - // AD::StartPreacc(); - // AD::SetPreaccIn(V_i, nDim+6); - // AD::SetPreaccIn(Vorticity_i, nDim); - // AD::SetPreaccIn(StrainMag_i); - // AD::SetPreaccIn(ScalarVar_i[0]); - // AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); - // AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - - Residual = 0.0; - Production = 0.0; - Destruction = 0.0; - CrossProduction = 0.0; - Jacobian_i[0] = 0.0; - - /* - From NASA Turbulence model site. http://turbmodels.larc.nasa.gov/spalart.html - This form was developed primarily to improve the near-wall numerical behavior of the model (i.e., the goal was to improve the convergence behavior). The reference is: - Edwards, J. R. and Chandra, S. "Comparison of Eddy Viscosity-Transport Turbulence Models for Three-Dimensional, Shock-Separated Flowfields," AIAA Journal, Vol. 34, No. 4, 1996, pp. 756-763. - In this modificaton Omega is replaced by Strain Rate - */ - - /*--- Evaluate Omega, here Omega is the Strain Rate ---*/ - - Sbar = 0.0; - for(iDim=0;iDim 1e-10) { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - nu = Laminar_Viscosity_i/Density_i; - Ji = ScalarVar_i[0]/nu; - Ji_2 = Ji*Ji; - Ji_3 = Ji_2*Ji; - fv1 = Ji_3/(Ji_3+cv1_3); - fv2 = 1.0 - Ji/(1.0+Ji*fv1); - ft2 = ct3*exp(-ct4*Ji_2); - S = Omega; - inv_k2_d2 = 1.0/(k2*dist_i_2); - - Shat = max(S*((1.0/max(Ji,1.0e-16))+fv1),1.0e-16); - - Shat = max(Shat, 1.0e-10); - inv_Shat = 1.0/Shat; - - /*--- Production term ---*/; - - Production = cb1*Shat*ScalarVar_i[0]*Volume; - - /*--- Destruction term ---*/ - - r = min(ScalarVar_i[0]*inv_Shat*inv_k2_d2,10.0); - r=tanh(r)/tanh(1.0); - - g = r + cw2*(pow(r,6.0)-r); - g_6 = pow(g,6.0); - glim = pow((1.0+cw3_6)/(g_6+cw3_6),1.0/6.0); - fw = g*glim; - - Destruction = cw1*fw*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production - Destruction + CrossProduction; - - /*--- Compressibility Correction term ---*/ - Pressure_i = V_i[idx.Pressure()]; - SoundSpeed_i = sqrt(Pressure_i*Gamma/Density_i); - aux_cc=0; - for(iDim=0;iDim(&Residual, &Jacobian_i, nullptr); - -} - -template -CSourcePieceWise_TurbSA_Neg::CSourcePieceWise_TurbSA_Neg(unsigned short val_nDim, - unsigned short val_nVar, - const CConfig* config) : - CSourceBase_TurbSA(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()) { } - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSA_Neg::ComputeResidual(const CConfig* config) { - - unsigned short iDim; - -// AD::StartPreacc(); -// AD::SetPreaccIn(V_i, nDim+6); -// AD::SetPreaccIn(Vorticity_i, nDim); -// AD::SetPreaccIn(StrainMag_i); -// AD::SetPreaccIn(ScalarVar_i[0]); -// AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); -// AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - - Residual = 0.0; - Production = 0.0; - Destruction = 0.0; - CrossProduction = 0.0; - Jacobian_i[0] = 0.0; - - /*--- Evaluate Omega ---*/ - - Omega = sqrt(Vorticity_i[0]*Vorticity_i[0] + Vorticity_i[1]*Vorticity_i[1] + Vorticity_i[2]*Vorticity_i[2]); - - /*--- Rotational correction term ---*/ - - if (rotating_frame) { Omega += 2.0*min(0.0, StrainMag_i-Omega); } - - if (dist_i > 1e-10) { - - if (ScalarVar_i[0] > 0.0) { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - nu = Laminar_Viscosity_i/Density_i; - Ji = ScalarVar_i[0]/nu; - Ji_2 = Ji*Ji; - Ji_3 = Ji_2*Ji; - fv1 = Ji_3/(Ji_3+cv1_3); - fv2 = 1.0 - Ji/(1.0+Ji*fv1); - ft2 = ct3*exp(-ct4*Ji_2); - S = Omega; - inv_k2_d2 = 1.0/(k2*dist_i_2); - - Shat = S + ScalarVar_i[0]*fv2*inv_k2_d2; - Shat = max(Shat, 1.0e-10); - inv_Shat = 1.0/Shat; - - /*--- Production term ---*/; - - // Original SA model - // Production = cb1*(1.0-ft2)*Shat*ScalarVar_i[0]*Volume; - - Production = cb1*Shat*ScalarVar_i[0]*Volume; - - /*--- Destruction term ---*/ - - r = min(ScalarVar_i[0]*inv_Shat*inv_k2_d2,10.0); - g = r + cw2*(pow(r,6.0)-r); - g_6 = pow(g,6.0); - glim = pow((1.0+cw3_6)/(g_6+cw3_6),1.0/6.0); - fw = g*glim; - - Destruction = cw1*fw*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production - Destruction + CrossProduction; - - /*--- Implicit part, production term ---*/ - - dfv1 = 3.0*Ji_2*cv1_3/(nu*pow(Ji_3+cv1_3,2.)); - dfv2 = -(1/nu-Ji_2*dfv1)/pow(1.+Ji*fv1,2.); - if ( Shat <= 1.0e-10 ) dShat = 0.0; - else dShat = (fv2+ScalarVar_i[0]*dfv2)*inv_k2_d2; - Jacobian_i[0] += cb1*(ScalarVar_i[0]*dShat+Shat)*Volume; - - /*--- Implicit part, destruction term ---*/ - - dr = (Shat-ScalarVar_i[0]*dShat)*inv_Shat*inv_Shat*inv_k2_d2; - if (r == 10.0) dr = 0.0; - dg = dr*(1.+cw2*(6.0*pow(r,5.0)-1.0)); - dfw = dg*glim*(1.-g_6/(g_6+cw3_6)); - Jacobian_i[0] -= cw1*(dfw*ScalarVar_i[0] + 2.0*fw)*ScalarVar_i[0]/dist_i_2*Volume; - - } - - else { - - /*--- Production term ---*/ - - dist_i_2 = dist_i*dist_i; - - /*--- Production term ---*/; - - Production = cb1*(1.0-ct3)*Omega*ScalarVar_i[0]*Volume; - - /*--- Destruction term ---*/ - - Destruction = cw1*ScalarVar_i[0]*ScalarVar_i[0]/dist_i_2*Volume; - - /*--- Diffusion term ---*/ - - norm2_Grad = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - norm2_Grad += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim]; - - CrossProduction = cb2_sigma*norm2_Grad*Volume; - - Residual = Production + Destruction + CrossProduction; - - /*--- Implicit part, production term ---*/ - - Jacobian_i[0] += cb1*(1.0-ct3)*Omega*Volume; - - /*--- Implicit part, destruction term ---*/ - - Jacobian_i[0] += 2.0*cw1*ScalarVar_i[0]/dist_i_2*Volume; - - } - - } - -// AD::SetPreaccOut(Residual); -// AD::EndPreacc(); - - return ResidualType<>(&Residual, &Jacobian_i, nullptr); - -} - -template -CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, - unsigned short val_nVar, - const su2double *constants, - su2double val_kine_Inf, - su2double val_omega_Inf, - const CConfig* config) : - CNumerics(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()) { - - incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); - sustaining_terms = (config->GetKind_Turb_Model() == TURB_MODEL::SST_SUST); - axisymmetric = config->GetAxisymmetric(); - - /*--- Closure constants ---*/ - sigma_k_1 = constants[0]; - sigma_k_2 = constants[1]; - sigma_w_1 = constants[2]; - sigma_w_2 = constants[3]; - beta_1 = constants[4]; - beta_2 = constants[5]; - beta_star = constants[6]; - a1 = constants[7]; - alfa_1 = constants[8]; - alfa_2 = constants[9]; - - /*--- Set the ambient values of k and omega to the free stream values. ---*/ - kAmb = val_kine_Inf; - omegaAmb = val_omega_Inf; - - /*--- "Allocate" the Jacobian using the static buffer. ---*/ - Jacobian_i[0] = Jacobian_Buffer; - Jacobian_i[1] = Jacobian_Buffer+2; - -} - -template -CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfig* config) { - - AD::StartPreacc(); - AD::SetPreaccIn(StrainMag_i); - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); - AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i); - AD::SetPreaccIn(F1_i); AD::SetPreaccIn(F2_i); AD::SetPreaccIn(CDkw_i); - AD::SetPreaccIn(PrimVar_Grad_i, nDim+idx.Velocity(), nDim); - AD::SetPreaccIn(Vorticity_i, 3); - - unsigned short iDim; - su2double alfa_blended, beta_blended; - su2double diverg, pk, pw, zeta; - su2double VorticityMag = sqrt(Vorticity_i[0]*Vorticity_i[0] + - Vorticity_i[1]*Vorticity_i[1] + - Vorticity_i[2]*Vorticity_i[2]); - - AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); - - Density_i = V_i[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; - - Residual[0] = 0.0; Residual[1] = 0.0; - Jacobian_i[0][0] = 0.0; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = 0.0; - - /*--- Computation of blended constants for the source terms---*/ - - alfa_blended = F1_i*alfa_1 + (1.0 - F1_i)*alfa_2; - beta_blended = F1_i*beta_1 + (1.0 - F1_i)*beta_2; - - if (dist_i > 1e-10) { - - /*--- Production ---*/ - - diverg = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - diverg += PrimVar_Grad_i[iDim+idx.Velocity()][iDim]; - - /* if using UQ methodolgy, calculate production using perturbed Reynolds stress matrix */ - - if (using_uq){ - ComputePerturbedRSM(nDim, Eig_Val_Comp, uq_permute, uq_delta_b, uq_urlx, - PrimVar_Grad_i+idx.Velocity(), Density_i, Eddy_Viscosity_i, - ScalarVar_i[0], MeanPerturbedRSM); - SetPerturbedStrainMag(ScalarVar_i[0]); - pk = Eddy_Viscosity_i*PerturbedStrainMag*PerturbedStrainMag - - 2.0/3.0*Density_i*ScalarVar_i[0]*diverg; - } - else { - pk = Eddy_Viscosity_i*StrainMag_i*StrainMag_i - 2.0/3.0*Density_i*ScalarVar_i[0]*diverg; - } - - pk = min(pk,20.0*beta_star*Density_i*ScalarVar_i[1]*ScalarVar_i[0]); - pk = max(pk,0.0); - - zeta = max(ScalarVar_i[1], VorticityMag*F2_i/a1); - - /* if using UQ methodolgy, calculate production using perturbed Reynolds stress matrix */ - - if (using_uq){ - pw = PerturbedStrainMag * PerturbedStrainMag - 2.0/3.0*zeta*diverg; - } - else { - pw = StrainMag_i*StrainMag_i - 2.0/3.0*zeta*diverg; - } - pw = alfa_blended*Density_i*max(pw,0.0); - - /*--- Sustaining terms, if desired. Note that if the production terms are - larger equal than the sustaining terms, the original formulation is - obtained again. This is in contrast to the version in literature - where the sustaining terms are simply added. This latter approach could - lead to problems for very big values of the free-stream turbulence - intensity. ---*/ - - if ( sustaining_terms ) { - const su2double sust_k = beta_star*Density_i*kAmb*omegaAmb; - const su2double sust_w = beta_blended*Density_i*omegaAmb*omegaAmb; - - pk = max(pk, sust_k); - pw = max(pw, sust_w); - } - - /*--- Add the production terms to the residuals. ---*/ - - Residual[0] += pk*Volume; - Residual[1] += pw*Volume; - - /*--- Dissipation ---*/ - - Residual[0] -= beta_star*Density_i*ScalarVar_i[1]*ScalarVar_i[0]*Volume; - Residual[1] -= beta_blended*Density_i*ScalarVar_i[1]*ScalarVar_i[1]*Volume; - - /*--- Cross diffusion ---*/ - - Residual[1] += (1.0 - F1_i)*CDkw_i*Volume; - - /*--- Contribution due to 2D axisymmetric formulation ---*/ - - if (axisymmetric) ResidualAxisymmetric(alfa_blended,zeta); - - /*--- Implicit part ---*/ - - Jacobian_i[0][0] = -beta_star*ScalarVar_i[1]*Volume; - Jacobian_i[0][1] = -beta_star*ScalarVar_i[0]*Volume; - Jacobian_i[1][0] = 0.0; - Jacobian_i[1][1] = -2.0*beta_blended*ScalarVar_i[1]*Volume; - } - - AD::SetPreaccOut(Residual, nVar); - AD::EndPreacc(); - - return ResidualType<>(Residual, Jacobian_i, nullptr); - -} - -template -void CSourcePieceWise_TurbSST::SetPerturbedStrainMag(su2double turb_ke) { - - /*--- Compute norm of perturbed strain rate tensor. ---*/ - - PerturbedStrainMag = 0; - 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); - - PerturbedStrainMag += pow(StrainRate_ij, 2.0); - } - } - PerturbedStrainMag = sqrt(2.0*PerturbedStrainMag); - -} - -/*--- Explicit instantiations until we don't move this to the hpp. ---*/ -template class CSourcePieceWise_TurbSA >; -template class CSourcePieceWise_TurbSA >; -template class CSourcePieceWise_TurbSA >; - -template class CSourcePieceWise_TurbSA_COMP >; -template class CSourcePieceWise_TurbSA_COMP >; -template class CSourcePieceWise_TurbSA_COMP >; - -template class CSourcePieceWise_TurbSA_E >; -template class CSourcePieceWise_TurbSA_E >; -template class CSourcePieceWise_TurbSA_E >; - -template class CSourcePieceWise_TurbSA_E_COMP >; -template class CSourcePieceWise_TurbSA_E_COMP >; -template class CSourcePieceWise_TurbSA_E_COMP >; - -template class CSourcePieceWise_TurbSA_Neg >; -template class CSourcePieceWise_TurbSA_Neg >; -template class CSourcePieceWise_TurbSA_Neg >; - -template class CSourcePieceWise_TurbSST >; -template class CSourcePieceWise_TurbSST >; -template class CSourcePieceWise_TurbSST >; diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index f88a04b35abd..d28520ef186d 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -275,7 +275,6 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); - const bool transition = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM); const bool transition_BC = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); @@ -304,12 +303,6 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai numerics->SetStrainMag(flowNodes->GetStrainMag(iPoint), 0.0); - /*--- Set intermittency ---*/ - - if (transition) { - numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetIntermittency(iPoint)); - } - /*--- Turbulent variables w/o reconstruction, and its gradient ---*/ numerics->SetScalarVar(nodes->GetSolution(iPoint), nullptr); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 121e4a12fe18..dcf255af0bc6 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -300,7 +300,7 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Menter's second blending function ---*/ - numerics->SetF2blending(nodes->GetF2blending(iPoint),0.0); + numerics->SetF2blending(nodes->GetF2blending(iPoint)); /*--- Set vorticity and strain rate magnitude ---*/ @@ -310,7 +310,7 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Cross diffusion ---*/ - numerics->SetCrossDiff(nodes->GetCrossDiff(iPoint),0.0); + numerics->SetCrossDiff(nodes->GetCrossDiff(iPoint)); if (axisymmetric){ /*--- Set y coordinate ---*/ diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref index 90fc07b20e38..8c4cf0f978ca 100644 --- a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref +++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref @@ -1,2 +1,2 @@ "VARIABLE" , "AVG_DENSITY[0]", "AVG_ENTHALPY[0]", "AVG_NORMALVEL[0]", "DRAG[0]" , "EFFICIENCY[0]" , "FORCE_X[0]" , "FORCE_Y[0]" , "FORCE_Z[0]" , "LIFT[0]" , "MOMENT_X[0]" , "MOMENT_Y[0]" , "MOMENT_Z[0]" , "SIDEFORCE[0]" , "SURFACE_MACH[0]", "SURFACE_MASSFLOW[0]", "SURFACE_MOM_DISTORTION[0]", "SURFACE_PRESSURE_DROP[0]", "SURFACE_SECONDARY[0]", "SURFACE_SECOND_OVER_UNIFORM[0]", "SURFACE_STATIC_PRESSURE[0]", "SURFACE_STATIC_TEMPERATURE[0]", "SURFACE_TOTAL_PRESSURE[0]", "SURFACE_TOTAL_TEMPERATURE[0]", "SURFACE_UNIFORMITY[0]", "AVG_TEMPERATURE[1]", "MAXIMUM_HEATFLUX[1]", "TOTAL_HEATFLUX[1]", "FINDIFF_STEP" -0 , 0.0 , -100000.01639127731, -2.2204999999731917e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -0.04999999997368221, 1.1109999995104307e-08, -2.0599999983605954 , 0.0 , 2.12000000054946 , 3.709999998879887 , 330.00000030369847 , -39.999997625272954 , 315.0000011942211 , -39.999997625272954 , -1.400000004814217 , -129.99999512430804, 0.0 , -509.99999530176865, 1e-08 +0 , 0.0 , -100000.01639127731, 1.1102999999850092e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -0.04999999997368221, -4.4410000000756306e-08, -2.0599999983605954 , 0.0 , 2.12000000054946 , 3.709999998879887 , 330.00000030369847 , -39.999997625272954 , 315.0000011942211 , -39.999997625272954 , -1.400000004814217 , -129.99999512430804, 0.0 , -509.99999530176865, 1e-08 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index b0052c9645a2..e59579bd53df 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -710,7 +710,7 @@ def main(): turbmod_sa_comp_edw_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_comp_edw_rae2822.cfg_file = "turb_SA_COMP_EDW_RAE2822.cfg" turbmod_sa_comp_edw_rae2822.test_iter = 20 - turbmod_sa_comp_edw_rae2822.test_vals = [-2.004687, 0.742306, 0.497310, -5.290769, 0.809485, 0.062036] + turbmod_sa_comp_edw_rae2822.test_vals = [-2.004685, 0.742307, 0.497311, -5.290750, 0.809487, 0.062045] turbmod_sa_comp_edw_rae2822.su2_exec = "mpirun -n 2 SU2_CFD" turbmod_sa_comp_edw_rae2822.timeout = 1600 turbmod_sa_comp_edw_rae2822.new_output = True @@ -1428,7 +1428,7 @@ def main(): sp_pinArray_3d_cht_mf_hf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_3d" sp_pinArray_3d_cht_mf_hf_tp.cfg_file = "configMaster.cfg" sp_pinArray_3d_cht_mf_hf_tp.test_iter = 30 - sp_pinArray_3d_cht_mf_hf_tp.test_vals = [-13.380025, -7.476945, -7.025285, -0.009675, 99.879812, 4.1920e+02] + sp_pinArray_3d_cht_mf_hf_tp.test_vals = [-13.380430, -7.476945, -7.025285, -0.009675, 99.879812, 4.1920e+02] sp_pinArray_3d_cht_mf_hf_tp.su2_exec = "mpirun -n 2 SU2_CFD" sp_pinArray_3d_cht_mf_hf_tp.timeout = 1600 sp_pinArray_3d_cht_mf_hf_tp.tol = 0.00001 diff --git a/TestCases/vandv.py b/TestCases/vandv.py index a3f7daaa5b3f..35d6b4a514c0 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -45,7 +45,7 @@ def main(): p30n30.cfg_dir = "vandv/rans/30p30n" p30n30.cfg_file = "config.cfg" p30n30.test_iter = 20 - p30n30.test_vals = [-10.628370, -10.299097, -10.485327, -10.238060, -13.517229, 0.050962, 2.828563, 1.317849, -0.227233] + p30n30.test_vals = [-10.642687, -10.303442, -10.496806, -10.253229, -13.517216, 0.050962, 2.828563, 1.317849, -0.215792] test_list.append(p30n30) #################