Updating SA-Diffusion term Discretisation

.github/workflows/regression.yml

@@ -209,7 +209,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_updated_SA_diffusion -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:
@@ -255,7 +255,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c feature_updated_SA_diffusion -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536
         with:
@@ -300,7 +300,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c feature_updated_SA_diffusion -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536
         with:

Common/include/CConfig.hpp

@@ -590,6 +590,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. */
@@ -4499,6 +4500,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/src/CConfig.cpp

@@ -1996,6 +1996,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,7 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::Jacobian_j;
 
   const su2double sigma = 2.0/3.0;
+  const su2double cb2 = 0.622;
 
   /*!
    * \brief Adds any extra variables to AD
@@ -67,17 +68,41 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
 
     /*--- Compute mean effective viscosity ---*/
 
-    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]);
+    /* First Term */
+    su2double nu_i = Laminar_Viscosity_i/Density_i;
+    su2double nu_j = Laminar_Viscosity_j/Density_j;
+    su2double nu_e = 0.5*(nu_i + nu_j + (1+cb2)*(ScalarVar_i[0] + ScalarVar_j[0])); 
+    su2double term_1 = nu_e; 
+
+    /* Second Term */
+    su2double nu_tilde_i = ScalarVar_i[0];
+    su2double term_2 = cb2 * nu_tilde_i;
+
+    Flux[0] = (term_1 - term_2)*Proj_Mean_GradScalarVar[0]/sigma;
+    su2double dTerm1_dnut_i = (1+cb2)*0.5;
+    su2double dTerm1_dnut_j = (1+cb2)*0.5;
+
+    su2double dTerm2_dnut_i = cb2;
+    su2double dTerm2_dnut_j = 0.0;
+
+    su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i;
+    su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j;
 
-    Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
+    su2double dGrad_dnut_i = -proj_vector_ij;
+    su2double dGrad_dnut_j =  proj_vector_ij;
 
     /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+    su2double diffusion_coefficient = term_1 - term_2;
+    bool UseExactJacobians = config->GetUse_Accurate_Turb_Jacobians();
 
     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;
+      if (UseExactJacobians) {
+        Jacobian_i[0][0] = (dDC_dnut_i*Proj_Mean_GradScalarVar[0] + dGrad_dnut_i*diffusion_coefficient)/sigma; //exact
+        Jacobian_j[0][0] = (dDC_dnut_j*Proj_Mean_GradScalarVar[0] + dGrad_dnut_j*diffusion_coefficient)/sigma; //exact
+      } else {
+        Jacobian_i[0][0] = diffusion_coefficient*-proj_vector_ij/sigma; //frozen diffusion coefficients
+        Jacobian_j[0][0] = diffusion_coefficient* proj_vector_ij/sigma; //frozen diffusion coefficients
+      }    
     }
   }
 
@@ -118,6 +143,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
@@ -137,27 +163,41 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
     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]);
+    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) ---*/
+    su2double zeta_i = ((1 + cb2)*nu_tilde_ij - cb2*nu_tilde_i)/nu_ij;
+    su2double fn_i;
+    if (zeta_i >= 0.0) {
+      fn_i = 1.0;  
+    } else {
+      fn_i = (cn1 + pow(zeta_i,3.0))/(cn1 - pow(zeta_i,3.0));
+    }
 
-    su2double nu_e;
+    su2double term_1 = (nu_ij + (1 + cb2)*nu_tilde_ij*fn_i);
+    su2double term_2 = cb2*nu_tilde_i*fn_i;
+    Flux[0] = (term_1 - term_2)*Proj_Mean_GradScalarVar[0]/sigma;
 
-    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;
-    }
+    /*--- 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. ---*/
 
-    Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
+    su2double diffusion_coefficient = (term_1 - term_2);
 
-    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+    su2double dGrad_dnut_i = -proj_vector_ij;
+    su2double dGrad_dnut_j =  proj_vector_ij;
 
-    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;
-    }
+    if (implicit) { 
+        Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i/sigma;  //frozen diffusion
+        Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j/sigma;  //frozen diffusion
+      }
   }
 
 public:

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -243,11 +243,11 @@ 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);
 
@@ -423,24 +423,22 @@ 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.
  * \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 +456,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 +485,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);
-  }
 };
 };
 

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> EdgeFluxes_Diff; /*!< \brief Flux difference across 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

@@ -348,8 +348,12 @@ void CScalarSolver<VariableType>::SumEdgeFluxes(CGeometry* geometry) {
     for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
       if (iPoint == geometry->edges->GetNode(iEdge, 0))
         LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
-      else
+      else {
         LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
+        if (EdgeFluxes_Diff.GetLocSize() > 0) {
+          LinSysRes.SubtractBlock(iPoint, EdgeFluxes_Diff.GetBlock(iEdge));
+        }
+      }
     }
   }
   END_SU2_OMP_FOR

SU2_CFD/src/solvers/CTurbSASolver.cpp

@@ -73,9 +73,11 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
     LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0);
     System.SetxIsZero(true);
 
-    if (ReducerStrategy)
+    if (ReducerStrategy) {
       EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr);
-
+      EdgeFluxes_Diff.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr);
+    }
+
     if (config->GetExtraOutput()) {
       if (nDim == 2) { nOutputVariables = 13; }
       else if (nDim == 3) { nOutputVariables = 19; }
@@ -301,8 +303,95 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry*
   };
 
   /*--- Now instantiate the generic implementation with the functor above. ---*/
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+    CFlowVariable* flowNodes = solver_container[FLOW_SOL] ?
+        su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes()) : nullptr;
+
+  /*--- Points in edge ---*/
+  auto iPoint = geometry->edges->GetNode(iEdge, 0);
+  auto jPoint = geometry->edges->GetNode(iEdge, 1);
+
+  /*--- Helper function to compute the flux ---*/
+  auto ComputeFlux = [&](unsigned long point_i, unsigned long point_j, const su2double* normal) {
+    numerics->SetCoord(geometry->nodes->GetCoord(point_i),geometry->nodes->GetCoord(point_j));
+    numerics->SetNormal(normal);
+
+    if (flowNodes) {
+      numerics->SetPrimitive(flowNodes->GetPrimitive(point_i), flowNodes->GetPrimitive(point_j));
+    }
+
+    numerics->SetScalarVar(nodes->GetSolution(point_i), nodes->GetSolution(point_j));
+    numerics->SetScalarVarGradient(nodes->GetGradient(point_i), nodes->GetGradient(point_j));
+
+    return numerics->ComputeResidual(config); 
+  };
+
+  SolverSpecificNumerics(iPoint, jPoint);
+
+  /*--- Compute fluxes and jacobians i->j ---*/
+  const su2double* normal = geometry->edges->GetNormal(iEdge);
+  auto residual_ij = ComputeFlux(iPoint, jPoint, normal);
+  if (ReducerStrategy) {
+    EdgeFluxes.SubtractBlock(iEdge, residual_ij);
+    EdgeFluxes_Diff.SetBlock(iEdge, residual_ij);
+    if (implicit) {
+      auto* Block_ij = Jacobian.GetBlock(iPoint, jPoint);
+      auto* Block_ji = Jacobian.GetBlock(jPoint, iPoint);
+
+      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) {
+      auto* Block_ii = Jacobian.GetBlock(iPoint, iPoint);
+      auto* Block_ij = Jacobian.GetBlock(iPoint, jPoint);
+
+      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]);
+        }
+    }
+  }
 
-  Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config);
+  /*--- Compute fluxes and jacobians j->i ---*/
+  su2double flipped_normal[MAXNDIM];
+  for (int iDim=0; iDim<nDim; iDim++)
+    flipped_normal[iDim] = -normal[iDim];
+
+  auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal);
+  if (ReducerStrategy) {
+    EdgeFluxes_Diff.AddBlock(iEdge, residual_ji);
+    if (implicit) {
+      auto* Block_ij = Jacobian.GetBlock(iPoint, jPoint);
+      auto* Block_ji = Jacobian.GetBlock(jPoint, iPoint);
+
+      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) {
+      auto* Block_ji = Jacobian.GetBlock(jPoint, iPoint);
+      auto* Block_jj = Jacobian.GetBlock(jPoint, jPoint);
+
+      //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]);
+        }
+    }
+  }
 }
 
 void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,