Updating the Cross-Diffusion Term Discretisation in k-w SST Model

Common/include/CConfig.hpp

@@ -639,6 +639,9 @@ class CConfig {
   su2double Relaxation_Factor_Adjoint;  /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
   su2double Relaxation_Factor_CHT;      /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
   su2double EntropyFix_Coeff;           /*!< \brief Entropy fix coefficient. */
+  su2double MaxUpdateSST;             /*!< \brief Cap for the Under-Relaxation Factor for SST Turbulent Variables*/
+  su2double MaxUpdateSA;              /*!< \brief Cap for the Under-Relaxation Factor for SA Turbulent Variables*/
+  su2double MaxUpdateFlow;            /*!< \brief Cap for the Under-Relaxation Factor for Flow Density and Energy Variables*/
   unsigned short nLocationStations,     /*!< \brief Number of section cuts to make when outputting mesh and cp . */
   nWingStations;                        /*!< \brief Number of section cuts to make when calculating internal volume. */
   su2double Kappa_1st_AdjFlow,  /*!< \brief Lax 1st order dissipation coefficient for adjoint flow equations (coarse multigrid levels). */
@@ -4233,6 +4236,24 @@ class CConfig {
    */
   su2double GetRelaxation_Factor_CHT(void) const { return Relaxation_Factor_CHT; }
 
+  /*!
+   * \brief Get the maximum update ratio for flow variables- density and energy.
+   * \return Maximum allowable update ratio for flow variables.
+   */
+  su2double GetMaxUpdateFractionFlow(void) const { return MaxUpdateFlow; }
+
+  /*!
+   * \brief Get the maximum update ratio for SA variable nu_tilde.
+   * \return Maximum allowable update ratio for SA variables.
+   */
+  su2double GetMaxUpdateFractionSA(void) const { return MaxUpdateSA; }
+
+  /*!
+   * \brief Get the maximum update ratio for SST turbulence variables TKE and Omega.
+   * \return Maximum allowable update ratio for SST variables.
+   */
+  su2double GetMaxUpdateFractionSST(void) const { return MaxUpdateSST; }
+
   /*!
    * \brief Get the number of samples used in quasi-Newton methods.
    */

Common/src/CConfig.cpp

@@ -1900,6 +1900,13 @@ void CConfig::SetConfig_Options() {
   addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU);
   /* DESCRIPTION: Linear solver for the discete adjoint systems */
 
+  /* DESCRIPTION: Maximum update ratio value for flow density and energy variables */
+  addDoubleOption("MAX_UPDATE_FLOW", MaxUpdateFlow, 0.2);
+  /* DESCRIPTION: Maximum update ratio value for SA turbulence variable nu_tilde */
+  addDoubleOption("MAX_UPDATE_SA", MaxUpdateSA, 0.99);
+  /* DESCRIPTION: Maximum update ratio value for SST turbulence variables TKE and Omega */
+  addDoubleOption("MAX_UPDATE_SST", MaxUpdateSST, 1.0);
+
   /*!\par CONFIG_CATEGORY: Convergence\ingroup Config*/
   /*--- Options related to convergence ---*/
 

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -240,6 +240,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
   const su2double sigma_k2;
   const su2double sigma_om1;
   const su2double sigma_om2;
+  const bool use_accurate_jacobians;
 
   su2double F1_i, F1_j; /*!< \brief Menter's first blending function */
 
@@ -270,20 +271,50 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
     const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j;
 
     const su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine);
-    const su2double diff_omega = 0.5*(diff_i_omega + diff_j_omega);
+    const su2double diff_omega_T1 = 0.5*(diff_i_omega + diff_j_omega);
+
+    /*--- We aim to treat the cross-diffusion as a diffusion term rather than a source term.
+    * Re-writing the cross-diffusion contribution as λ/w ∇w ∇k, where λ = (2 (1- F1) ρ σ_ω2)
+    * and expanding using the product rule for divergence theorem gives: ∇(w λ/w ∇k) - w ∇(λ/w ∇k). 
+    * Discretising using FVM, gives: (λ)_ij ∇k - w_c (λ/w)_ij ∇k. where w_c is the cell centre value ---*/
+
+    const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i;
+    const su2double lambda_j = 2 * (1 - F1_j) * Density_j * sigma_omega_j;
+    const su2double lambda_ij = 0.5 * (lambda_i + lambda_j);
+    const su2double w_ij = 0.5 * (ScalarVar_i[1] + ScalarVar_j[1]);
+
+    const su2double diff_omega_T2 = lambda_ij;
+
+    const su2double diff_omega_T3 = -ScalarVar_i[1] * lambda_ij/w_ij;
 
     Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0];
-    Flux[1] = diff_omega*Proj_Mean_GradScalarVar[1];
+    Flux[1] = diff_omega_T1*Proj_Mean_GradScalarVar[1] + (diff_omega_T2 + diff_omega_T3)*Proj_Mean_GradScalarVar[0];
 
     /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/
     if (implicit) {
       const su2double proj_on_rho_i = proj_vector_ij/Density_i;
-      Jacobian_i[0][0] = -diff_kine*proj_on_rho_i;  Jacobian_i[0][1] = 0.0;
-      Jacobian_i[1][0] = 0.0;                       Jacobian_i[1][1] = -diff_omega*proj_on_rho_i;
-
       const su2double proj_on_rho_j = proj_vector_ij/Density_j;
-      Jacobian_j[0][0] = diff_kine*proj_on_rho_j;   Jacobian_j[0][1] = 0.0;
-      Jacobian_j[1][0] = 0.0;                       Jacobian_j[1][1] = diff_omega*proj_on_rho_j;
+      Jacobian_i[0][0] = -diff_kine*proj_on_rho_i;
+      Jacobian_i[0][1] = 0.0;
+      Jacobian_i[1][0] = (diff_omega_T2+diff_omega_T3)*-proj_on_rho_i;
+      Jacobian_i[1][1] = -diff_omega_T1*proj_on_rho_i;
+
+      Jacobian_j[0][0] = diff_kine*proj_on_rho_j;
+      Jacobian_j[0][1] = 0.0;
+      Jacobian_j[1][0] = (diff_omega_T2+diff_omega_T3)*proj_on_rho_j;
+      Jacobian_j[1][1] = diff_omega_T1*proj_on_rho_j;
+
+      if (use_accurate_jacobians) {
+        Jacobian_i[0][0] = -diff_kine*proj_on_rho_i;
+        Jacobian_i[0][1] = 0.0;
+        Jacobian_i[1][0] = (diff_omega_T2 + diff_omega_T3)*-proj_on_rho_i;
+        Jacobian_i[1][1] = -proj_on_rho_i * diff_omega_T1 - 2*lambda_ij*ScalarVar_j[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0];
+
+        Jacobian_j[0][0] = diff_kine*proj_on_rho_j;
+        Jacobian_j[0][1] = 0.0;
+        Jacobian_j[1][0] = (diff_omega_T2 + diff_omega_T3)*proj_on_rho_j;
+        Jacobian_j[1][1] = proj_on_rho_j * diff_omega_T1 + 2*lambda_ij*ScalarVar_i[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0];
+      }      
     }
   }
 
@@ -302,7 +333,8 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
       sigma_k1(constants[0]),
       sigma_k2(constants[1]),
       sigma_om1(constants[2]),
-      sigma_om2(constants[3]) {
+      sigma_om2(constants[3]),
+      use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {
   }
 
   /*!

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -772,7 +772,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
     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()]);
@@ -911,9 +910,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       Residual[0] -= dk * Volume;
       Residual[1] -= dw * Volume;
 
-      /*--- Cross diffusion ---*/
-
-      Residual[1] += (1.0 - F1_i) * CDkw_i * Volume;
+      /*--- Cross diffusion is included in the viscous fluxes, discretisation in turb_diffusion.hpp ---*/
 
       /*--- Contribution due to 2D axisymmetric formulation ---*/
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -558,7 +558,7 @@ void CFVMFlowSolverBase<V, R>::ComputeUnderRelaxationFactor(const CConfig* confi
   /* Loop over the solution update given by relaxing the linear
    system for this nonlinear iteration. */
 
-  const su2double allowableRatio = 0.2;
+  const su2double allowableRatio = config->GetMaxUpdateFractionFlow();
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -146,6 +146,102 @@ class CScalarSolver : public CSolver {
     }
   }
 
+  /*!
+   * \brief Compute the viscous flux for the turbulence equations at a particular edge for a non-conservative discretisation.
+   * \tparam SolverSpecificNumericsTemp - lambda-function, to implement solver specific contributions to numerics.
+   * \note The functor has to implement (iPoint, jPoint)
+   * \param[in] iEdge - Edge for which we want to compute the flux
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   */
+  template<typename SolverSpecificNumericsFunc>
+  void Viscous_Residual_NonCons(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container,
+                                        CNumerics* numerics, const CConfig* config, SolverSpecificNumericsFunc&& SolverSpecificNumerics) {
+    const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+      CFlowVariable* flowNodes = solver_container[FLOW_SOL] ?
+          su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes()) : nullptr;
+
+    const auto iPoint = geometry->edges->GetNode(iEdge, 0);
+    const auto jPoint = geometry->edges->GetNode(iEdge, 1);
+
+    /*--- Lambda function to compute the flux ---*/
+    auto ComputeFlux = [&](unsigned long iPoint, unsigned long jPoint, const su2double* normal) {
+      numerics->SetCoord(geometry->nodes->GetCoord(iPoint),geometry->nodes->GetCoord(jPoint));
+      numerics->SetNormal(normal);
+
+      if (flowNodes) {
+        numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), flowNodes->GetPrimitive(jPoint));
+      }
+
+    /*--- Solver specific numerics contribution. ---*/
+      SolverSpecificNumerics(iPoint, jPoint);
+
+      numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint));
+      numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(jPoint));
+
+      return numerics->ComputeResidual(config); 
+    };
+
+    /*--- Compute fluxes and jacobians i->j ---*/
+    const su2double* normal = geometry->edges->GetNormal(iEdge);
+    auto residual_ij = ComputeFlux(iPoint, jPoint, normal);
+
+    JacobianScalarType *Block_ii = nullptr, *Block_ij = nullptr, *Block_ji = nullptr, *Block_jj = nullptr;
+    if (implicit) {
+      Jacobian.GetBlocks(iEdge, iPoint, jPoint, Block_ii, Block_ij, Block_ji, Block_jj);
+    }
+    if (ReducerStrategy) {
+      EdgeFluxes.SubtractBlock(iEdge, residual_ij);
+      EdgeFluxesDiff.SetBlock(iEdge, residual_ij);
+      if (implicit) {
+        /*--- For the reducer strategy the Jacobians are averaged for simplicity. ---*/
+        for (int iVar=0; iVar<nVar; iVar++)
+          for (int jVar=0; jVar<nVar; jVar++) {
+            Block_ij[iVar*nVar + jVar] -= 0.5 * SU2_TYPE::GetValue(residual_ij.jacobian_j[iVar][jVar]);
+            Block_ji[iVar*nVar + jVar] += 0.5 * SU2_TYPE::GetValue(residual_ij.jacobian_i[iVar][jVar]);
+          }
+      }
+    } else {
+      LinSysRes.SubtractBlock(iPoint, residual_ij);
+      if (implicit) {
+        for (int iVar=0; iVar<nVar; iVar++)
+          for (int jVar=0; jVar<nVar; jVar++) {
+            Block_ii[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_i[iVar][jVar]);
+            Block_ij[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_j[iVar][jVar]);
+          }
+      }
+    }
+
+    /*--- Compute fluxes and jacobians j->i ---*/
+    su2double flipped_normal[MAXNDIM];
+    for (auto iDim = 0u; iDim < nDim; iDim++) flipped_normal[iDim] = -normal[iDim];
+
+    auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal);
+    if (ReducerStrategy) {
+      EdgeFluxesDiff.AddBlock(iEdge, residual_ji);
+      if (implicit) {
+        for (int iVar=0; iVar<nVar; iVar++)
+          for (int jVar=0; jVar<nVar; jVar++) {
+            Block_ij[iVar*nVar + jVar] += 0.5 * SU2_TYPE::GetValue(residual_ji.jacobian_i[iVar][jVar]);
+            Block_ji[iVar*nVar + jVar] -= 0.5 * SU2_TYPE::GetValue(residual_ji.jacobian_j[iVar][jVar]);
+          }
+      }
+    } else {
+      LinSysRes.SubtractBlock(jPoint, residual_ji);
+      if (implicit) {
+        /*--- The order of arguments were flipped in the evaluation of residual_ji, the Jacobian
+        * associated with point i is stored in jacobian_j and point j in jacobian_i. ---*/
+        for (int iVar=0; iVar<nVar; iVar++)
+          for (int jVar=0; jVar<nVar; jVar++) {
+            Block_ji[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ji.jacobian_j[iVar][jVar]);
+            Block_jj[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ji.jacobian_i[iVar][jVar]);
+          }
+      }
+    }
+  }
+
   /*!
    * \brief Generic implementation of the fluid interface boundary condition for scalar solvers.
    * \tparam SolverSpecificNumericsFunc - lambda that implements solver specific contributions to viscous numerics.

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -53,6 +53,14 @@ class CTurbSSTSolver final : public CTurbSolver {
                      CSolver **solver_container,
                      const CConfig *config,
                      unsigned short val_marker);
+
+  /*!
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
+   * a nonlinear iteration for stability.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeUnderRelaxationFactor(const CConfig *config);
+
 public:
   /*!
    * \brief Constructor.

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -147,4 +147,11 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
    * \returns The number of extra variables.
    */
   unsigned long RegisterSolutionExtra(bool input, const CConfig* config) final;
+
+  /*!
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
+   * a nonlinear iteration for stability.
+   * \param[in] allowableRatio - Maximum percentage update in variable per iteration.
+   */
+  void ComputeUnderRelaxationFactorHelper(su2double allowableRatio);
 };

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -100,6 +100,12 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \brief Set the Diffusion Coefficients of TKE and omega equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_DC_kw - diffusion coefficient value
+   */
+
   /*!
    * \brief Register eddy viscosity (muT) as Input or Output of an AD recording.
    * \param[in] input - Boolean whether In- or Output should be registered.