Apply kappa-blending to Roe MUSCL reconstruction

Author: bmunguia

Common/include/CConfig.hpp

@@ -592,11 +592,12 @@ class CConfig {
   MUSCL_AdjFlow,           /*!< \brief MUSCL scheme for the adj flow equations.*/
   MUSCL_AdjTurb,           /*!< \brief MUSCL scheme for the adj turbulence equations.*/
   MUSCL_Species;           /*!< \brief MUSCL scheme for the species equations.*/
-  su2double MUSCL_Kappa_Flow,    /*!< \brief Blending coefficient for MUSCL scheme for the flow equations.*/
-  MUSCL_Kappa_Turb,              /*!< \brief Blending coefficient for MUSCL scheme for the turbulence equations.*/
-  MUSCL_Kappa_Heat,              /*!< \brief Blending coefficient for MUSCL scheme for the (fvm) heat equation.*/
-  MUSCL_Kappa_AdjFlow,           /*!< \brief Blending coefficient for MUSCL scheme for the adj flow equations.*/
-  MUSCL_Kappa_Species;      /*!< \brief MUSCL scheme for the species equations.*/
+  su2double MUSCL_Kappa,   /*!< \brief Blending coefficient for MUSCL scheme. */
+  MUSCL_Kappa_Flow,        /*!< \brief Blending coefficient for MUSCL scheme for the flow equations.*/
+  MUSCL_Kappa_Turb,        /*!< \brief Blending coefficient for MUSCL scheme for the turbulence equations.*/
+  MUSCL_Kappa_Heat,        /*!< \brief Blending coefficient for MUSCL scheme for the (fvm) heat equation.*/
+  MUSCL_Kappa_AdjFlow,     /*!< \brief Blending coefficient for MUSCL scheme for the adj flow equations.*/
+  MUSCL_Kappa_Species;     /*!< \brief Blending coefficient for 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. */
@@ -2461,11 +2462,12 @@ class CConfig {
    * \param[in] val_kind_upwind - If upwind scheme, kind of upwind scheme (Roe, etc.).
    * \param[in] val_kind_slopelimit - If upwind scheme, kind of slope limit.
    * \param[in] val_muscl - Define if we apply a MUSCL scheme or not.
+   * \param[in] val_muscl_kappa - Define the blending coefficient for the MUSCL scheme.
    * \param[in] val_kind_fem - If FEM, what kind of FEM discretization.
    */
   void SetKind_ConvNumScheme(unsigned short val_kind_convnumscheme, CENTERED val_kind_centered,
                              UPWIND val_kind_upwind, LIMITER val_kind_slopelimit,
-                             bool val_muscl,  unsigned short val_kind_fem);
+                             bool val_muscl, su2double val_muscl_kappa, unsigned short val_kind_fem);
 
   /*!
    * \brief Get the value of limiter coefficient.

Common/src/CConfig.cpp

@@ -8697,13 +8697,14 @@ unsigned short CConfig::GetContainerPosition(unsigned short val_eqsystem) {
 void CConfig::SetKind_ConvNumScheme(unsigned short val_kind_convnumscheme,
                                     CENTERED val_kind_centered, UPWIND val_kind_upwind,
                                     LIMITER val_kind_slopelimit, bool val_muscl,
-                                    unsigned short val_kind_fem) {
+                                    su2double val_muscl_kappa, unsigned short val_kind_fem) {
   Kind_ConvNumScheme = val_kind_convnumscheme;
   Kind_Centered = val_kind_centered;
   Kind_Upwind = val_kind_upwind;
   Kind_FEM = val_kind_fem;
   Kind_SlopeLimit = val_kind_slopelimit;
   MUSCL = val_muscl;
+  MUSCL_Kappa = val_muscl_kappa;
 
 }
 
@@ -8721,7 +8722,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
     if (val_system == RUNTIME_FLOW_SYS) {
       SetKind_ConvNumScheme(Kind_ConvNumScheme_Flow, Kind_Centered_Flow,
                             Kind_Upwind_Flow, Kind_SlopeLimit_Flow,
-                            MUSCL_Flow, NONE);
+                            MUSCL_Flow, MUSCL_Kappa_Flow, NONE);
       SetKind_TimeIntScheme(Kind_TimeIntScheme_Flow);
     }
   };
@@ -8730,15 +8731,16 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
     if (val_system == RUNTIME_TURB_SYS) {
       SetKind_ConvNumScheme(Kind_ConvNumScheme_Turb, Kind_Centered_Turb,
                             Kind_Upwind_Turb, Kind_SlopeLimit_Turb,
-                            MUSCL_Turb, NONE);
+                            MUSCL_Turb, MUSCL_Kappa_Turb, NONE);
       SetKind_TimeIntScheme(Kind_TimeIntScheme_Turb);
     }
   };
 
   auto SetHeatParam = [&]() {
     if (val_system == RUNTIME_HEAT_SYS) {
       SetKind_ConvNumScheme(Kind_ConvNumScheme_Heat, Kind_Centered_Heat,
-                            Kind_Upwind_Heat, Kind_SlopeLimit_Heat, MUSCL_Heat, NONE);
+                            Kind_Upwind_Heat, Kind_SlopeLimit_Heat, MUSCL_Heat,
+                            MUSCL_Kappa_Heat, NONE);
       SetKind_TimeIntScheme(Kind_TimeIntScheme_Heat);
     }
   };
@@ -8747,7 +8749,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
     if (val_system == RUNTIME_SPECIES_SYS) {
       SetKind_ConvNumScheme(Kind_ConvNumScheme_Species, Kind_Centered_Species,
                             Kind_Upwind_Species, Kind_SlopeLimit_Species,
-                            MUSCL_Species, NONE);
+                            MUSCL_Species, MUSCL_Kappa_Species, NONE);
       SetKind_TimeIntScheme(Kind_TimeIntScheme_Species);
     }
   };
@@ -8756,7 +8758,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
     if (val_system == RUNTIME_ADJFLOW_SYS) {
       SetKind_ConvNumScheme(Kind_ConvNumScheme_AdjFlow, Kind_Centered_AdjFlow,
                             Kind_Upwind_AdjFlow, Kind_SlopeLimit_AdjFlow,
-                            MUSCL_AdjFlow, NONE);
+                            MUSCL_AdjFlow, MUSCL_Kappa_AdjFlow, NONE);
       SetKind_TimeIntScheme(Kind_TimeIntScheme_AdjFlow);
     }
   };
@@ -8782,7 +8784,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
       if (val_system == RUNTIME_TRANS_SYS) {
         SetKind_ConvNumScheme(Kind_ConvNumScheme_Turb, Kind_Centered_Turb,
                               Kind_Upwind_Turb, Kind_SlopeLimit_Turb,
-                              MUSCL_Turb, NONE);
+                              MUSCL_Turb, MUSCL_Kappa_Turb, NONE);
         SetKind_TimeIntScheme(Kind_TimeIntScheme_Turb);
       }
       break;
@@ -8796,7 +8798,7 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
       if (val_system == RUNTIME_FLOW_SYS) {
         SetKind_ConvNumScheme(Kind_ConvNumScheme_FEM_Flow, Kind_Centered_Flow,
                               Kind_Upwind_Flow, Kind_SlopeLimit_Flow,
-                              MUSCL_Flow, Kind_FEM_Flow);
+                              MUSCL_Flow, MUSCL_Kappa_Flow, Kind_FEM_Flow);
         SetKind_TimeIntScheme(Kind_TimeIntScheme_FEM_Flow);
       }
       break;
@@ -8813,22 +8815,22 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
       if (val_system == RUNTIME_ADJTURB_SYS) {
         SetKind_ConvNumScheme(Kind_ConvNumScheme_AdjTurb, Kind_Centered_AdjTurb,
                               Kind_Upwind_AdjTurb, Kind_SlopeLimit_AdjTurb,
-                              MUSCL_AdjTurb, NONE);
+                              MUSCL_AdjTurb, 0.0, NONE);
         SetKind_TimeIntScheme(Kind_TimeIntScheme_AdjTurb);
       }
       break;
     case MAIN_SOLVER::HEAT_EQUATION:
     case MAIN_SOLVER::DISC_ADJ_HEAT:
       if (val_system == RUNTIME_HEAT_SYS) {
-        SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE);
+        SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, 0.0, NONE);
         SetKind_TimeIntScheme(Kind_TimeIntScheme_Heat);
       }
       break;
 
     case MAIN_SOLVER::FEM_ELASTICITY:
     case MAIN_SOLVER::DISC_ADJ_FEM:
       if (val_system == RUNTIME_FEA_SYS) {
-        SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE);
+        SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, 0.0, NONE);
         SetKind_TimeIntScheme(NONE);
       }
       break;

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

@@ -35,16 +35,47 @@
 /*!
  * \brief Unlimited reconstruction.
  */
-template<size_t nVarGrad_ = 0, size_t nVar, size_t nDim, class Gradient_t>
+// template<size_t nVarGrad_ = 0, size_t nVar, size_t nDim, class Gradient_t>
+// FORCEINLINE void musclUnlimited(Int iPoint,
+//                                 const VectorDbl<nDim>& vector_ij,
+//                                 Double kappa,
+//                                 const Gradient_t& gradient,
+//                                 VectorDbl<nVar>& vars) {
+//   constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : nVar;
+//   auto grad = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+//   for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
+//     vars(iVar) += scale * dot(grad[iVar], vector_ij);
+//   }
+// }
+template<size_t nVarGrad_ = 0, size_t nDim, class VarType, class Gradient_t>
 FORCEINLINE void musclUnlimited(Int iPoint,
+                                Int jPoint,
                                 const VectorDbl<nDim>& vector_ij,
-                                Double scale,
                                 const Gradient_t& gradient,
-                                VectorDbl<nVar>& vars) {
-  constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : nVar;
-  auto grad = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+                                CPair<VarType>& V,
+                                Double kappa) {
+  constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
+
+  auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+  auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
+
   for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
-    vars(iVar) += scale * dot(grad[iVar], vector_ij);
+    /*--- Centered difference: dV_ij^cent = Vj - Vi ---*/
+    const Double delta_cent = V.j.all(iVar) - V.i.all(iVar);
+
+    /*--- Upwind difference: dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cent ---*/
+    const Double grad_proj_i = 2.0 * dot(grad_i[iVar], vector_ij);
+    const Double grad_proj_j = 2.0 * dot(grad_j[iVar], vector_ij);
+    const Double delta_upwind_i = grad_proj_i - delta_cent;
+    const Double delta_upwind_j = grad_proj_j - delta_cent;
+
+    /*--- Blended difference: dV_ij^kappa = (1-kappa)ΔV_ij^upw + (1+kappa)ΔV_ij^cent ---*/
+    const Double delta_kappa_i = (1.0 - kappa) * delta_upwind_i + (1.0 + kappa) * delta_cent;
+    const Double delta_kappa_j = (1.0 - kappa) * delta_upwind_j + (1.0 + kappa) * delta_cent;
+
+    /*--- Apply reconstruction: V_ij = Vi + 0.25 * dV_ij^kappa ---*/
+    V.i.all(iVar) += 0.25 * delta_kappa_i;
+    V.j.all(iVar) -= 0.25 * delta_kappa_j;
   }
 }
 
@@ -100,6 +131,7 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
  * \param[in] gamma - Heat capacity ratio.
  * \param[in] gasConst - Specific gas constant.
  * \param[in] muscl - If true, reconstruct, else simply fetch.
+ * \param[in] musclKappa - Blending coefficient for MUSCL reconstruction.
  * \param[in] limiterType - Type of flux limiter.
  * \param[in] V1st - Pair of compressible flow primitives for nodes i,j.
  * \param[in] vector_ij - Distance vector from i to j.
@@ -110,7 +142,8 @@ template<class ReconVarType, class PrimVarType, size_t nDim, class VariableType>
 FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint,
                                                       const su2double& gamma,
                                                       const su2double& gasConst,
-                                                      bool muscl, LIMITER limiterType,
+                                                      bool muscl, const su2double& musclKappa,
+                                                      LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
                                                       const VariableType& solution) {
@@ -132,8 +165,9 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
 
     switch (limiterType) {
     case LIMITER::NONE:
-      musclUnlimited<nVarGrad>(iPoint, vector_ij, 0.5, gradients, V.i.all);
-      musclUnlimited<nVarGrad>(jPoint, vector_ij,-0.5, gradients, V.j.all);
+      // musclUnlimited<nVarGrad>(iPoint, vector_ij, 0.5, gradients, V.i.all);
+      // musclUnlimited<nVarGrad>(jPoint, vector_ij,-0.5, gradients, V.j.all);
+      musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, musclKappa);
       break;
     case LIMITER::VAN_ALBADA_EDGE:
       musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V);

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -58,6 +58,7 @@ class CRoeBase : public Base {
   const su2double gamma;
   const su2double gasConst;
   const su2double entropyFix;
+  const su2double muscl_kappa;
   const bool finestGrid;
   const bool dynamicGrid;
   const bool muscl;
@@ -75,6 +76,7 @@ class CRoeBase : public Base {
     finestGrid(iMesh == MESH_0),
     dynamicGrid(config.GetDynamic_Grid()),
     muscl(finestGrid && config.GetMUSCL_Flow()),
+    muscl_kappa(config.GetMUSCL_Kappa_Flow()),
     typeLimiter(config.GetKind_SlopeLimit_Flow()) {
   }
 
@@ -119,7 +121,7 @@ class CRoeBase : public Base {
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution);
+        iEdge, iPoint, jPoint, gamma, gasConst, muscl, muscl_kappa, typeLimiter, V1st, vector_ij, solution);
 
     /*--- Compute conservative variables. ---*/