[WIP] updating cross-diffusion discretisation in sst

Common/include/CConfig.hpp

@@ -635,6 +635,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 SST_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent Variables*/
+  su2double SA_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent Variables*/
+  su2double Flow_UnderRelaxation_Factor; /*!< \brief UnderRelaxation Factor for SST Turbulent 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). */
@@ -4207,6 +4210,24 @@ class CConfig {
    */
   su2double GetRelaxation_Factor_CHT(void) const { return Relaxation_Factor_CHT; }
 
+  /*!
+   * \brief Get the under-relaxation for flow variables density and energy.
+   * \return under-relaxation for flow variables.
+   */
+  su2double GetUnderRelax_Flow(void) const { return Flow_UnderRelaxation_Factor; }
+
+  /*!
+   * \brief Get the under-relaxation for SA variable, nu_tilde.
+   * \return under-relaxation for SA variables.
+   */
+  su2double GetUnderRelax_SA(void) const { return SA_UnderRelaxation_Factor; }
+
+  /*!
+   * \brief Get the under-relaxation for SST turbulence variables k and omega.
+   * \return under-relaxation for SST variables.
+   */
+  su2double GetUnderRelax_SST(void) const { return SST_UnderRelaxation_Factor; }
+
   /*!
    * \brief Get the number of samples used in quasi-Newton methods.
    */

Common/src/CConfig.cpp

@@ -1887,6 +1887,12 @@ 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 */
 
+  addDoubleOption("UR_FACTOR_FLOW", Flow_UnderRelaxation_Factor, 0.2);
+  /* DESCRIPTION: Under-Relaxation Factor for Density and Energy Variables */
+  addDoubleOption("UR_FACTOR_SA", SA_UnderRelaxation_Factor, 0.99);
+  /* DESCRIPTION: Under-Relaxation Factor for SA Turbulence Variables */
+  addDoubleOption("UR_FACTOR_SST", SST_UnderRelaxation_Factor, 1.0);
+  /* DESCRIPTION: Under-Relaxation Factor for SST Turbulence Variables */
   /*!\par CONFIG_CATEGORY: Convergence\ingroup Config*/
   /*--- Options related to convergence ---*/
 

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -223,28 +223,43 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
     const su2double sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2;
     const su2double sigma_omega_i = F1_i*sigma_om1 + (1.0 - F1_i)*sigma_om2;
     const su2double sigma_omega_j = F1_j*sigma_om1 + (1.0 - F1_j)*sigma_om2;
+    const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i / ScalarVar_i[1];
+    const su2double lambda_j = 2 * (1 - F1_j) * Density_j * sigma_omega_j / ScalarVar_j[1];
+    const su2double lambda_ij = 0.5 * (lambda_i + lambda_j);
 
     /*--- Compute mean effective dynamic viscosity ---*/
     const su2double diff_i_kine = Laminar_Viscosity_i + sigma_kine_i*Eddy_Viscosity_i;
     const su2double diff_j_kine = Laminar_Viscosity_j + sigma_kine_j*Eddy_Viscosity_j;
-    const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i;
-    const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j;
+    const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i + lambda_i * ScalarVar_i[0];
+    const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j + lambda_j * ScalarVar_j[0];
+
+    //non-conservative term:
+    const su2double diff_omega_t1 = 0.5 * (diff_i_omega + diff_j_omega);
+    const su2double diff_omega_t2 = -ScalarVar_i[0] * lambda_ij;
 
     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 = diff_omega_t1 + diff_omega_t2;
 
     Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0];
     Flux[1] = diff_omega*Proj_Mean_GradScalarVar[1];
 
-    /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/
+    /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients 
+    * Using chain rule for dFlux[1]/dx_{i.j} = d(diff_coef)/dx_{i,j}*Proj + d(Proj)/dx_{i,j}*diff_ceof---*/
+
     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 dlambda_domega_i = - (2 * (1 - F1_i) * Density_i * (sigma_omega_i))/pow(ScalarVar_i[1], 2);
+      Jacobian_i[0][0] = -diff_kine*proj_on_rho_i;  
+      Jacobian_i[0][1] = 0.0;
+      Jacobian_i[1][0] = 0.0; //dF_w/dk_i
+      Jacobian_i[1][1] = (diff_omega) * -proj_on_rho_i; //lambda not frozen
 
       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;
+      const su2double dlambda_domega_j = - (2 * (1 - F1_j) * Density_j * (sigma_omega_j))/pow(ScalarVar_j[1], 2);    
+      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; //frozen diff_ceoff
     }
   }
 

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -924,7 +924,7 @@
 
       /*--- Cross diffusion ---*/
 
-      Residual[1] += (1.0 - F1_i) * CDkw_i * Volume;
+      // Residual[1] += (1.0 - F1_i) * CDkw_i * Volume;
 
       /*--- 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->GetUnderRelax_Flow();
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {

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/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/src/solvers/CTurbSASolver.cpp

@@ -1560,7 +1560,7 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) {
    system for this nonlinear iteration. */
 
   su2double localUnderRelaxation =  1.00;
-  const su2double allowableRatio =  0.99;
+  const su2double allowableRatio =  config->GetUnderRelax_SA();
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -75,8 +75,10 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh
     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);
+    }
 
     /*--- Initialize the BGS residuals in multizone problems. ---*/
     if (multizone){
@@ -295,8 +297,95 @@ void CTurbSSTSolver::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 CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,
@@ -1019,3 +1108,38 @@ void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMark
   }
 
 }
+
+void CTurbSSTSolver::ComputeUnderRelaxationFactor(const CConfig *config) {
+
+  const su2double allowableRatio = config->GetUnderRelax_SST();
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  /* Loop over the solution update given by relaxing the linear
+   system for this nonlinear iteration. */
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    su2double localUnderRelaxation = 1.0;
+
+    for (unsigned short iVar =0; iVar < nVar; iVar++) {
+      const unsigned long index = iPoint * nVar + iVar;
+      su2double ratio = fabs(LinSysSol[index])/(fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
+      /* We impose a limit on the maximum percentage that the
+      turbulence variables can change over a nonlinear iteration. */
+      if (ratio > allowableRatio) {
+        localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
+      }
+
+    }
+
+    /* Threshold the relaxation factor in the event that there is
+     a very small value. This helps avoid catastrophic crashes due
+     to non-realizable states by canceling the update. */
+
+    if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0;
+
+    /* Store the under-relaxation factor for this point. */
+
+    nodes->SetUnderRelaxation(iPoint, localUnderRelaxation);
+  }
+
+  END_SU2_OMP_FOR
+}