[WIP] defect correction of Nishikawa for skewed meshes

Common/include/CConfig.hpp

@@ -569,6 +569,7 @@ class CConfig {
   Kind_Centered_Template;       /*!< \brief Centered scheme for the template model. */
 
 
+  VISCOUS_GRAD_CORR Kind_Visc_Corr; /*!< \brief Gradient correction damping. */
   FEM_SHOCK_CAPTURING_DG Kind_FEM_Shock_Capturing_DG; /*!< \brief Shock capturing method for the FEM DG solver. */
   BGS_RELAXATION Kind_BGS_RelaxMethod; /*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */
   bool ReconstructionGradientRequired; /*!< \brief Enable or disable a second gradient calculation for upwind reconstruction only. */
@@ -4595,6 +4596,12 @@ class CConfig {
    */
   UPWIND GetKind_Upwind(void) const { return Kind_Upwind; }
 
+  /*!
+   * \brief Get kind of gradient correction damping term for the viscous terms.
+   * \return Kind of gradient correction damping term.
+   */
+  VISCOUS_GRAD_CORR GetKind_ViscousGradCorr(void) const { return Kind_Visc_Corr; }
+
   /*!
    * \brief Get if the upwind scheme used MUSCL or not.
    * \note This is the information that the code will use, the method will

Common/include/option_structure.hpp

@@ -878,6 +878,22 @@ static const MapType<std::string, UPWIND> Upwind_Map = {
   MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH)
 };
 
+/*!
+ * \brief Types of viscous gradient correction spatial discretizations
+ */
+enum class VISCOUS_GRAD_CORR {
+  NONE,           /*!< \brief No correction is used. */
+  EDGE_NORMAL,    /*!< \brief Edge Normal Correction. */
+  FACE_TANGENT,   /*!< \brief Face Tangent correction. */
+  ALPHA_DAMPING   /*!< \brief Alpha Damping with a=5/4. */
+};
+static const MapType<std::string, VISCOUS_GRAD_CORR> Viscous_Grad_Corr_Map = {
+    MakePair("NONE", VISCOUS_GRAD_CORR::NONE)
+    MakePair("EDGE_NORMAL", VISCOUS_GRAD_CORR::EDGE_NORMAL)
+    MakePair("FACE_TANGENT", VISCOUS_GRAD_CORR::FACE_TANGENT)
+    MakePair("ALPHA_DAMPING", VISCOUS_GRAD_CORR::ALPHA_DAMPING)
+};
+
 /*!
  * \brief Types of FEM spatial discretizations
  */

Common/src/CConfig.cpp

@@ -1693,6 +1693,7 @@ void CConfig::SetConfig_Options() {
   addStringDoubleListOption("MARKER_HEATFLUX", nMarker_HeatFlux, Marker_HeatFlux, Heat_Flux);
   /*!\brief INTEGRATED_HEATFLUX \n DESCRIPTION: Prescribe Heatflux in [W] instead of [W/m^2] \ingroup Config \default false */
   addBoolOption("INTEGRATED_HEATFLUX", Integrated_HeatFlux, false);
+  addEnumOption("VISC_GRAD_CORR", Kind_Visc_Corr, Viscous_Grad_Corr_Map, VISCOUS_GRAD_CORR::EDGE_NORMAL);
   /*!\brief MARKER_HEATTRANSFER DESCRIPTION: Heat flux with specified heat transfer coefficient boundary marker(s)\n
    * Format: ( Heat transfer marker, heat transfer coefficient, wall temperature (static), ... ) \ingroup Config  */
   addExhaustOption("MARKER_HEATTRANSFER", nMarker_HeatTransfer, Marker_HeatTransfer, HeatTransfer_Coeff, HeatTransfer_WallTemp);

SU2_CFD/include/numerics/CNumerics.hpp

@@ -654,20 +654,39 @@
   FORCEINLINE static su2double ComputeProjectedGradient(int nDim, int nVar, const Vec1& normal,
                                                         const Vec1& coord_i, const Vec1& coord_j,
                                                         const Mat& grad_i, const Mat& grad_j,
-                                                        bool correct,
+                                                        const VISCOUS_GRAD_CORR correct,
                                                         const Vec2& var_i, const Vec2& var_j,
                                                         su2double* projNormal,
                                                         su2double* projCorrected) {
     assert(nDim == 2 || nDim == 3);
     nDim = (nDim > 2)? 3 : 2;
     su2double edgeVec[MAXNDIM], dist_ij_2 = 0.0, proj_vector_ij = 0.0;
+    su2double diss = 0.0;
+    su2double Area2 = 0.0;
+    for (int iDim = 0; iDim < nDim; ++iDim) {
+      Area2 += pow(normal[iDim],2);
+    }
 
     for (int iDim = 0; iDim < nDim; iDim++) {
       edgeVec[iDim] = coord_j[iDim] - coord_i[iDim];
       dist_ij_2 += pow(edgeVec[iDim], 2);
       proj_vector_ij += edgeVec[iDim] * normal[iDim];
     }
-    proj_vector_ij /= max(dist_ij_2,EPS);
+
+    switch (correct) {
+      case VISCOUS_GRAD_CORR::EDGE_NORMAL:
+        diss = proj_vector_ij / max(dist_ij_2,EPS*EPS);
+        break;
+      case VISCOUS_GRAD_CORR::FACE_TANGENT:
+        diss = Area2 / max(proj_vector_ij,EPS);
+        break;
+      case VISCOUS_GRAD_CORR::ALPHA_DAMPING:
+        const su2double alpha = 4.0 / 3.0;        
+        diss = alpha * Area2 / max(abs(proj_vector_ij),EPS);
+        break;
+    }
+
+    //proj_vector_ij /= max(dist_ij_2,EPS);
 
     /*--- Mean gradient approximation. ---*/
     for (int iVar = 0; iVar < nVar; iVar++) {
@@ -677,11 +696,11 @@
       for (int iDim = 0; iDim < nDim; iDim++) {
         su2double meanGrad = 0.5 * (grad_i[iVar][iDim] + grad_j[iVar][iDim]);
         projNormal[iVar] += meanGrad * normal[iDim];
-        if (correct) edgeProj += meanGrad * edgeVec[iDim];
+        edgeProj += meanGrad * edgeVec[iDim];
       }
 
       projCorrected[iVar] = projNormal[iVar];
-      if (correct) projCorrected[iVar] -= (edgeProj - (var_j[iVar]-var_i[iVar])) * proj_vector_ij;
+      projCorrected[iVar] += ((var_j[iVar]-var_i[iVar]) - edgeProj) * diss;
     }
 
     return proj_vector_ij;

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -44,7 +44,7 @@
 class CAvgGrad_Base : public CNumerics {
 protected:
   const unsigned short nPrimVar;          /*!< \brief The size of the primitive variable array used in the numerics class. */
-  const bool correct_gradient;            /*!< \brief Apply a correction to the gradient term */
+  const VISCOUS_GRAD_CORR correct_gradient; /*!< \brief Apply a correction to the gradient term */
   bool implicit = false;                  /*!< \brief Implicit calculus. */
   su2double
   heat_flux_vector[MAXNDIM] = {0.0},      /*!< \brief Flux of total energy due to molecular and turbulent diffusion */
@@ -164,6 +164,15 @@ class CAvgGrad_Base : public CNumerics {
                        su2double val_dist_ij_2,
                        const unsigned short val_nPrimVar);
 
+  void GradientCorrection(su2double** GradPrimVar,
+                                    const su2double* val_PrimVar_i,
+                                    const su2double* val_PrimVar_j,
+                                    const su2double* val_edge_vector,
+                                    su2double val_dist_ij_2,
+                                    const su2double* diss,
+                                    const unsigned short val_nPrimVar);
+
+
 public:
 
   /*!
@@ -176,7 +185,7 @@ class CAvgGrad_Base : public CNumerics {
    */
   CAvgGrad_Base(unsigned short val_nDim, unsigned short val_nVar,
                 unsigned short val_nPrimVar,
-                bool val_correct_grad, const CConfig* config);
+               VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config);
 
   /*!
    * \brief Destructor of the class.
@@ -251,7 +260,7 @@ class CAvgGrad_Flow final : public CAvgGrad_Base {
    * \param[in] config - Definition of the particular problem.
    */
   CAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar,
-                bool val_correct_grad, const CConfig* config);
+               VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config);
 
   /*!
    * \brief Compute the viscous flow residual using an average of gradients.
@@ -330,7 +339,7 @@ class CAvgGradInc_Flow final : public CAvgGrad_Base {
    * \param[in] config - Definition of the particular problem.
    */
   CAvgGradInc_Flow(unsigned short val_nDim, unsigned short val_nVar,
-                   bool val_correct_grad, const CConfig* config);
+                  VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config);
 
   /*!
    * \brief Compute the viscous flow residual using an average of gradients.
@@ -383,7 +392,7 @@ class CGeneralAvgGrad_Flow final : public CAvgGrad_Base {
    * \param[in] val_correct_grad - Apply a correction to the gradient
    * \param[in] config - Definition of the particular problem.
    */
-  CGeneralAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, bool val_correct_grad, const CConfig* config);
+  CGeneralAvgGrad_Flow(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR val_correct_grad, const CConfig* config);
 
   /*!
    * \brief Compute the viscous flow residual using an average of gradients.

SU2_CFD/include/numerics/heat.hpp

@@ -80,7 +80,7 @@ class CAvgGrad_Heat final : public CAvgGrad_Scalar<CNoFlowIndices> {
    * \param[in] config - Definition of the particular problem.
    * \param[in] correct - Whether to correct the gradient.
    */
-  CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, bool correct)
+  CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, VISCOUS_GRAD_CORR correct)
     : CAvgGrad_Scalar<CNoFlowIndices>(val_nDim, 1, correct, config) {}
 
  private:

SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp

@@ -69,7 +69,8 @@ class CAvgGrad_Scalar : public CNumerics {
   su2double* Jacobian_j[MAXNVAR];             /*!< \brief Flux Jacobian w.r.t. node j. */
   su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */
 
-  const bool correct_gradient = false, incompressible = false;
+  const bool incompressible = false;
+  const VISCOUS_GRAD_CORR correct_gradient;
 
   /*!
    * \brief A pure virtual function; Adds any extra variables to AD
@@ -91,7 +92,7 @@ class CAvgGrad_Scalar : public CNumerics {
    * \param[in] correct_gradient - Whether to correct gradient for skewness.
    * \param[in] config - Definition of the particular problem.
    */
-  CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad,
+  CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad,
                   const CConfig* config)
     : CNumerics(val_nDim, val_nVar, config),
       idx(val_nDim, config->GetnSpecies()),
@@ -123,7 +124,7 @@ class CAvgGrad_Scalar : public CNumerics {
     AD::SetPreaccIn(Normal, nDim);
     AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim);
     AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim);
-    if (correct_gradient) {
+    if (correct_gradient != VISCOUS_GRAD_CORR::NONE) {
       AD::SetPreaccIn(ScalarVar_i, nVar);
       AD::SetPreaccIn(ScalarVar_j, nVar);
     }

SU2_CFD/include/numerics/species/species_diffusion.hpp

@@ -104,7 +104,7 @@ class CAvgGrad_Species final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] correct_grad - Whether to correct gradient for skewness.
    * \param[in] config - Definition of the particular problem.
    */
-  CAvgGrad_Species(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config)
+  CAvgGrad_Species(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config),
       turbulence(config->GetKind_Turb_Model() != TURB_MODEL::NONE) {}
 };

SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp

@@ -98,7 +98,7 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] correct_grad - Whether to correct gradient for skewness.
    * \param[in] config - Definition of the particular problem.
    */
-  CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config)
+  CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, VISCOUS_GRAD_CORR correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config){
   }
 

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -115,7 +115,7 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar,
-                  bool correct_grad, const CConfig* config)
+                  VISCOUS_GRAD_CORR correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config),
       use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {}
 };
@@ -208,7 +208,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   CAvgGrad_TurbSA_Neg(unsigned short val_nDim, unsigned short val_nVar,
-                      bool correct_grad, const CConfig* config)
+                     VISCOUS_GRAD_CORR correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config) {}
 };
 
@@ -275,7 +275,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
 
     /*--- 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). 
+    * 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;
@@ -284,7 +284,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
     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];
@@ -314,7 +314,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
         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];
-      }      
+      }
     }
   }
 
@@ -328,7 +328,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   CAvgGrad_TurbSST(unsigned short val_nDim, unsigned short val_nVar,
-                   const su2double* constants, bool correct_grad, const CConfig* config)
+                   const su2double* constants, VISCOUS_GRAD_CORR correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config),
       sigma_k1(constants[0]),
       sigma_k2(constants[1]),

SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp

@@ -49,18 +49,31 @@ FORCEINLINE MatrixDbl<nVar,nDim> averageGradient(Int iPoint, Int jPoint,
   return avgGrad;
 }
 
+template<size_t nSecVar, class SecondaryType>
+FORCEINLINE CCompressibleSecondary<nSecVar> averageSecondary(Int iPoint, Int jPoint,
+                                                             const SecondaryType& secondary) {
+
+    CCompressibleSecondary<nSecVar> out;
+    auto second_i = gatherVariables<nSecVar>(iPoint, secondary);
+    auto second_j = gatherVariables<nSecVar>(jPoint, secondary);
+    for (size_t iVar = 0; iVar < nSecVar; ++iVar) {
+        out.all(iVar) = 0.5 * (second_i(iVar) + second_j(iVar));
+    }
+    return out;
+}
+
 /*!
  * \brief Correct average gradient with the directional derivative to avoid decoupling.
  */
 template<size_t nVar, size_t nDim, class PrimitiveType>
 FORCEINLINE void correctGradient(const PrimitiveType& V,
                                  const VectorDbl<nDim>& vector_ij,
-                                 Double dist2_ij,
+                                 const VectorDbl<nDim>& diss,
                                  MatrixDbl<nVar,nDim>& avgGrad) {
   for (size_t iVar = 0; iVar < nVar; ++iVar) {
-    Double corr = (dot(avgGrad[iVar],vector_ij) - V.j.all(iVar) + V.i.all(iVar)) / dist2_ij;
+    Double corr = ( V.j.all(iVar) - V.i.all(iVar) - dot(avgGrad[iVar],vector_ij));
     for (size_t iDim = 0; iDim < nDim; ++iDim) {
-      avgGrad(iVar,iDim) -= corr * vector_ij(iDim);
+      avgGrad(iVar,iDim) += corr * diss(iDim);
     }
   }
 }

SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp

@@ -71,12 +71,16 @@
 class CCompressibleViscousFluxBase : public CNumericsSIMD {
 protected:
   static constexpr size_t nDim = NDIM;
-  static constexpr size_t nPrimVarGrad = nDim+1;
+  static constexpr size_t nPrimVarGrad = nDim+3;
+  static constexpr size_t nSeconVar = 8;
 
   const su2double gamma;
   const su2double gasConst;
   const su2double prandtlTurb;
-  const bool correct;
+  //const su2double cp;
+  const bool correct_EN;
+  const bool correct_FT;
+  const bool correct_AD;
   const bool useSA_QCR;
   const bool wallFun;
   const bool uq;
@@ -96,7 +100,9 @@
     gamma(config.GetGamma()),
     gasConst(config.GetGas_ConstantND()),
     prandtlTurb(config.GetPrandtl_Turb()),
-    correct(iMesh == MESH_0),
+    correct_EN(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::EDGE_NORMAL),
+    correct_FT(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::FACE_TANGENT),
+    correct_AD(iMesh == MESH_0 && config.GetKind_ViscousGradCorr() == VISCOUS_GRAD_CORR::ALPHA_DAMPING),
     useSA_QCR(config.GetSAParsedOptions().qcr2000),
     wallFun(config.GetWall_Functions()),
     uq(config.GetSSTParsedOptions().uq),
@@ -135,6 +141,7 @@
 
     const auto& solution = static_cast<const CNSVariable&>(solution_);
     const auto& gradient = solution.GetGradient_Primitive();
+    const auto& secondary = solution.GetSecondary();
 
     /*--- Compute distance and handle zero without "ifs" by making it large. ---*/
 
@@ -145,7 +152,22 @@
     /*--- Compute the corrected mean gradient. ---*/
 
     auto avgGrad = averageGradient<nPrimVarGrad,nDim>(iPoint, jPoint, gradient);
-    if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad);
+    auto avgSecond = averageSecondary<nSeconVar>(iPoint, jPoint, secondary);
+
+    Double eDn = dot(vector_ij,unitNormal);
+    mask = eDn < EPS;
+    eDn += mask / (EPS);
+    const Double eDotN = Double(1.0) / eDn;
+    const Double alpha = Double(4.0) / 3.0;
+
+    VectorDbl<nDim> diss;
+
+    if(correct_EN)      for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = vector_ij(iDim) / dist2_ij;
+    else if(correct_FT) for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = unitNormal(iDim) * eDotN;
+    else if(correct_AD) for (int iDim = 0; iDim < nDim; ++iDim) diss(iDim) = unitNormal(iDim) * alpha * abs(eDotN);
+
+    if (correct_EN || correct_FT || correct_AD)
+        correctGradient(V, vector_ij, diss, avgGrad);
 
     /*--- Stress and heat flux tensors. ---*/