Improving the discretisation of the SA non-linear diffusion term $(\tilde{\nu})^2$

Common/include/CConfig.hpp

@@ -591,6 +591,7 @@ class CConfig {
   MUSCL_AdjTurb;           /*!< \brief MUSCL scheme for the adj turbulence equations.*/
   bool MUSCL_Species;      /*!< \brief MUSCL scheme for the species equations.*/
   bool Use_Accurate_Jacobians;  /*!< \brief Use numerically computed Jacobians for AUSM+up(2) and SLAU(2). */
+  bool Use_Accurate_Turb_Jacobians; /*!< \brief Use numerically computed Jacobians for standard SA turbulence model. */
   bool EulerPersson;       /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */
   bool FSI_Problem = false,/*!< \brief Boolean to determine whether the simulation is FSI or not. */
   Multizone_Problem;       /*!< \brief Boolean to determine whether we are solving a multizone problem. */
@@ -4516,6 +4517,12 @@ class CConfig {
    */
   bool GetUse_Accurate_Jacobians(void) const { return Use_Accurate_Jacobians; }
 
+  /*!
+   * \brief Get whether to "Use Accurate Jacobians" for Standard SA turbulence model.
+   * \return yes/no.
+   */
+  bool GetUse_Accurate_Turb_Jacobians(void) const { return Use_Accurate_Turb_Jacobians; }
+
   /*!
    * \brief Get the kind of integration scheme (explicit or implicit)
    *        for the flow equations.

Common/include/linear_algebra/CSysMatrix.hpp

@@ -553,6 +553,22 @@ class CSysMatrix {
     SetBlock<MatrixType, false>(block_i, block_j, val_block, alpha);
   }
 
+  /*!
+   * \brief Returns the 4 blocks ii, ij, ji, jj used by "UpdateBlocks".
+   * \note This method assumes an FVM-type sparse pattern.
+   * \param[in] edge - Index of edge that connects iPoint and jPoint.
+   * \param[in] iPoint - Row to which we add the blocks.
+   * \param[in] jPoint - Row from which we subtract the blocks.
+   * \param[out] bii, bij, bji, bjj - Blocks of the matrix.
+   */
+  inline void GetBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, ScalarType*& bii,
+                        ScalarType*& bij, ScalarType*& bji, ScalarType*& bjj) {
+    bii = &matrix[dia_ptr[iPoint] * nVar * nEqn];
+    bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn];
+    bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn];
+    bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn];
+  }
+
   /*!
    * \brief Update 4 blocks ii, ij, ji, jj (add to i* sub from j*).
    * \note This method assumes an FVM-type sparse pattern.
@@ -566,10 +582,8 @@ class CSysMatrix {
   template <class MatrixType, class OtherType = ScalarType>
   inline void UpdateBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& block_i,
                            const MatrixType& block_j, OtherType scale = 1) {
-    ScalarType* bii = &matrix[dia_ptr[iPoint] * nVar * nEqn];
-    ScalarType* bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn];
-    ScalarType* bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn];
-    ScalarType* bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn];
+    ScalarType *bii, *bij, *bji, *bjj;
+    GetBlocks(iEdge, iPoint, jPoint, bii, bij, bji, bjj);
 
     unsigned long iVar, jVar, offset = 0;
 

Common/src/CConfig.cpp

@@ -2006,6 +2006,8 @@ void CConfig::SetConfig_Options() {
   /*!\brief CONV_NUM_METHOD_TURB
    *  \n DESCRIPTION: Convective numerical method \ingroup Config*/
   addConvectOption("CONV_NUM_METHOD_TURB", Kind_ConvNumScheme_Turb, Kind_Centered_Turb, Kind_Upwind_Turb);
+  /*!\brief USE_ACCURATE_TURB_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for Standard SA model \ingroup Config*/
+  addBoolOption("USE_ACCURATE_TURB_JACOBIANS", Use_Accurate_Turb_Jacobians, false);
 
   /*!\brief MUSCL_ADJTURB \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
   addBoolOption("MUSCL_ADJTURB", MUSCL_AdjTurb, false);

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -52,6 +52,9 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::Jacobian_j;
 
   const su2double sigma = 2.0/3.0;
+  const su2double cb2 = 0.622;
+
+  const bool use_accurate_jacobians;
 
   /*!
    * \brief Adds any extra variables to AD
@@ -67,17 +70,39 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
 
     /*--- Compute mean effective viscosity ---*/
 
+    /*--- First Term. Normal diffusion, and conservative part of the quadratic diffusion.
+     * ||grad nu_t||^2 = div(nu_t grad nu_t) - nu_t div grad nu_t ---*/
     const su2double nu_i = Laminar_Viscosity_i/Density_i;
     const su2double nu_j = Laminar_Viscosity_j/Density_j;
-    const su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]);
+    const su2double nu_e = 0.5 * (nu_i + nu_j + (1 + cb2) * (ScalarVar_i[0] + ScalarVar_j[0]));
+    const su2double term_1 = nu_e;
 
-    Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
+    /* Second Term (quadratic diffusion, non conservative). */
+    const su2double nu_tilde_i = ScalarVar_i[0];
+    const su2double term_2 = cb2 * nu_tilde_i;
 
-    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+    const su2double diffusion_coefficient = term_1 - term_2;
+    Flux[0] = diffusion_coefficient * Proj_Mean_GradScalarVar[0] / sigma;
 
     if (implicit) {
-      Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma;
-      Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma;
+      /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+      Jacobian_i[0][0] = -diffusion_coefficient * proj_vector_ij / sigma;
+      Jacobian_j[0][0] = diffusion_coefficient * proj_vector_ij / sigma;
+
+      if (use_accurate_jacobians) {
+        /*--- The diffusion coefficient is also a function of nu_t. ---*/
+        const su2double dTerm1_dnut_i = (1 + cb2) * 0.5;
+        const su2double dTerm1_dnut_j = (1 + cb2) * 0.5;
+
+        const su2double dTerm2_dnut_i = cb2;
+        const su2double dTerm2_dnut_j = 0.0;
+
+        const su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i;
+        const su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j;
+
+        Jacobian_i[0][0] += dDC_dnut_i * Proj_Mean_GradScalarVar[0] / sigma;
+        Jacobian_j[0][0] += dDC_dnut_j * Proj_Mean_GradScalarVar[0] / sigma;
+      }
     }
   }
 
@@ -91,7 +116,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
    */
   CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar,
                   bool correct_grad, const CConfig* config)
-    : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config) {}
+    : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config),
+      use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {}
 };
 
 /*!
@@ -118,6 +144,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
 
   const su2double sigma = 2.0/3.0;
   const su2double cn1 = 16.0;
+  const su2double cb2 = 0.622;
 
   /*!
    * \brief Adds any extra variables to AD
@@ -136,27 +163,39 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
     const su2double nu_i = Laminar_Viscosity_i/Density_i;
     const su2double nu_j = Laminar_Viscosity_j/Density_j;
 
-    const su2double nu_ij = 0.5*(nu_i+nu_j);
-    const su2double nu_tilde_ij = 0.5*(ScalarVar_i[0] + ScalarVar_j[0]);
-
-    su2double nu_e;
-
-    if (nu_tilde_ij > 0.0) {
-      nu_e = nu_ij + nu_tilde_ij;
-    }
-    else {
-      const su2double Xi = nu_tilde_ij/nu_ij;
-      const su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi);
-      nu_e = nu_ij + fn*nu_tilde_ij;
+    const su2double nu_ij = 0.5 * (nu_i + nu_j);
+    const su2double nu_tilde_i = ScalarVar_i[0];
+    const su2double nu_tilde_j = ScalarVar_j[0];
+    const su2double nu_tilde_ij = 0.5 * (nu_tilde_i + nu_tilde_j);
+
+    /*--- Following Diskin's implementation from 10.2514/1.J064629, they propose a new fn function
+     * to be evaluated at the cell to maintain positivity in the diffusion coefficient, which is
+     * used in both terms. The new fn term averaged across the face reverts to the original fn
+     * function. ---*/
+
+    /*--- Second Term (LHS) ---*/
+    const su2double zeta_i = ((1 + cb2) * nu_tilde_ij - cb2 * nu_tilde_i) / nu_ij;
+    su2double fn_i = 1.0;
+    if (zeta_i < 0.0) {
+      fn_i = (cn1 + pow(zeta_i,3)) / (cn1 - pow(zeta_i,3));
     }
 
-    Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
+    const su2double term_1 = (nu_ij + (1 + cb2) * nu_tilde_ij * fn_i);
+    const su2double term_2 = cb2 * nu_tilde_i * fn_i;
+    Flux[0] = (term_1 - term_2) * Proj_Mean_GradScalarVar[0] / sigma;
 
-    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients
+    * Exact Jacobians were tested on multiple cases but resulted in divergence of all
+    * simulations, hence only frozen diffusion coefficient (approximate) Jacobians are used. ---*/
 
     if (implicit) {
-      Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma;
-      Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma;
+      const su2double diffusion_coefficient = (term_1 - term_2);
+
+      const su2double dGrad_dnut_i = -proj_vector_ij;
+      const su2double dGrad_dnut_j = proj_vector_ij;
+
+      Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i / sigma;
+      Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j / sigma;
     }
   }
 

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -86,7 +86,7 @@ class CSourceBase_TurbSA : public CNumerics {
    */
   inline void ResidualAxisymmetricDiffusion(su2double sigma) {
     if (Coord_i[1] < EPS) return;
-    
+
     const su2double yinv = 1.0 / Coord_i[1];
     const su2double& nue = ScalarVar_i[0];
 
@@ -243,14 +243,14 @@ class CSourceBase_TurbSA : public CNumerics {
         var.interDestrFactor = 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]);
+      /*--- Compute production, destruction and jacobian ---*/
+      su2double Production = 0.0, Destruction = 0.0;
+      SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, Jacobian_i[0]);
 
-      Residual = (Production - Destruction + CrossProduction) * Volume;
+      Residual = (Production - Destruction) * Volume;
 
       if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma);
-      
+
       Jacobian_i[0] *= Volume;
     }
 
@@ -423,24 +423,23 @@ struct Edw {
 };
 
 /*!
- * \brief SA source terms classes: production, destruction and cross-productions term and their derivative.
+ * \brief SA source terms classes: production and destruction term and their derivative.
+ * \note Quadratic diffusion is included in the viscous fluxes.
  * \ingroup SourceDiscr
  * \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.
  */
 struct 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) {
+                  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,
@@ -458,23 +457,17 @@ struct Bsl {
     jacobian -= var.interDestrFactor * ((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) {
+                  su2double& jacobian) {
     if (nue > 0.0) {
-      Bsl::get(nue, var, production, destruction, cross_production, jacobian);
+      Bsl::get(nue, var, production, destruction, jacobian);
     } else {
       ComputeProduction(nue, var, production, jacobian);
       ComputeDestruction(nue, var, destruction, jacobian);
-      ComputeCrossProduction(nue, var, cross_production, jacobian);
     }
   }
 
@@ -493,10 +486,6 @@ struct Neg {
     jacobian -= 2 * dD_dnu * var.interDestrFactor;
   }
 
-  static void ComputeCrossProduction(const su2double& nue, const CSAVariables& var, su2double& cross_production,
-                                     su2double& jacobian) {
-    Bsl::ComputeCrossProduction(nue, var, cross_production, jacobian);
-  }
 };
 };
 
@@ -562,7 +551,7 @@ class CCompressibilityCorrection final : public ParentClass {
       const su2double v = V_i[idx.Velocity() + 1];
 
       const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume;
-      const su2double axiCorrection = 0.5 * nue * d_axiCorrection; 
+      const su2double axiCorrection = 0.5 * nue * d_axiCorrection;
 
       this->Residual -= axiCorrection;
       this->Jacobian_i[0] -= d_axiCorrection;

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -77,6 +77,7 @@ class CScalarSolver : public CSolver {
 
   /*--- Edge fluxes for reducer strategy (see the notes in CEulerSolver.hpp). ---*/
   CSysVector<su2double> EdgeFluxes; /*!< \brief Flux across each edge. */
+  CSysVector<su2double> EdgeFluxesDiff; /*!< \brief Flux difference between ij and ji for non-conservative discretisation. */
 
   /*!
    * \brief The highest level in the variable hierarchy this solver can safely use.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -341,15 +341,21 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
 
 template <class VariableType>
 void CScalarSolver<VariableType>::SumEdgeFluxes(CGeometry* geometry) {
+  const bool nonConservative = EdgeFluxesDiff.GetLocSize() > 0;
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
     LinSysRes.SetBlock_Zero(iPoint);
 
     for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
-      if (iPoint == geometry->edges->GetNode(iEdge, 0))
+      if (iPoint == geometry->edges->GetNode(iEdge, 0)) {
         LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
-      else
+      } else {
         LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
+        if (nonConservative) {
+          LinSysRes.SubtractBlock(iPoint, EdgeFluxesDiff.GetBlock(iEdge));
+        }
+      }
     }
   }
   END_SU2_OMP_FOR

SU2_CFD/include/solvers/CSolver.hpp

@@ -191,12 +191,12 @@ class CSolver {
   CSysVector<su2double> LinSysSol;    /*!< \brief vector to store iterative solution of implicit linear system. */
   CSysVector<su2double> LinSysRes;    /*!< \brief vector to store iterative residual of implicit linear system. */
 #ifndef CODI_FORWARD_TYPE
-  CSysMatrix<su2mixedfloat> Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */
-  CSysSolve<su2mixedfloat>  System;   /*!< \brief Linear solver/smoother. */
+  using JacobianScalarType = su2mixedfloat;
 #else
-  CSysMatrix<su2double> Jacobian;
-  CSysSolve<su2double>  System;
+  using JacobianScalarType = su2double;
 #endif
+  CSysMatrix<JacobianScalarType> Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */
+  CSysSolve<JacobianScalarType> System;    /*!< \brief Linear solver/smoother. */
 
   CSysVector<su2double> OutputVariables;    /*!< \brief vector to store the extra variables to be written. */
   string* OutputHeadingNames;               /*!< \brief vector of strings to store the headings for the exra variables */