Merge pull request #2339 from su2code/fix_roe_kappa

Fix hllc jacobians

Common/include/CConfig.hpp

@@ -4335,8 +4335,7 @@ class CConfig {
   array<su2double,4> GetNewtonKrylovDblParam(void) const { return NK_DblParam; }
 
   /*!
-   * \brief Get the relaxation coefficient of the linear solver for the implicit formulation.
-   * \return relaxation coefficient of the linear solver for the implicit formulation.
+   * \brief Returns the Roe kappa (multipler of the dissipation term).
    */
   su2double GetRoe_Kappa(void) const { return Roe_Kappa; }
 

SU2_CFD/include/numerics/flow/convection/hllc.hpp

@@ -49,7 +49,7 @@ class CUpwHLLC_Flow final : public CNumerics {
 
   su2double sq_velRoe, RoeDensity, RoeEnthalpy, RoeSoundSpeed, RoeProjVelocity, ProjInterfaceVel;
 
-  su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho, kappa;
+  su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho;
 
   su2double Omega, RHO, OmegaSM;
   su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU;
@@ -102,7 +102,7 @@ class CUpwGeneralHLLC_Flow final : public CNumerics {
   su2double sq_velRoe, RoeDensity, RoeEnthalpy, RoeSoundSpeed, RoeProjVelocity, ProjInterfaceVel;
   su2double Kappa_i, Kappa_j, Chi_i, Chi_j, RoeKappa, RoeChi, RoeKappaStaticEnthalpy;
 
-  su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho, kappa;
+  su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho;
 
   su2double Omega, RHO, OmegaSM;
   su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU;

SU2_CFD/src/numerics/flow/convection/hllc.cpp

@@ -31,7 +31,7 @@
 CUpwHLLC_Flow::CUpwHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) {
 
   implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-  kappa = config->GetRoe_Kappa();
+
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
 
@@ -114,30 +114,19 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
   Density_j  = V_j[nDim+2];
   Enthalpy_j = V_j[nDim+3];
 
-
-  sq_vel_i = 0.0;
-  sq_vel_j = 0.0;
-
-  for (iDim = 0; iDim < nDim; iDim++) {
-    sq_vel_i += Velocity_i[iDim] * Velocity_i[iDim];
-    sq_vel_j += Velocity_j[iDim] * Velocity_j[iDim];
-  }
+  sq_vel_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i);
+  sq_vel_j = GeometryToolbox::SquaredNorm(nDim, Velocity_j);
 
   Energy_i = Enthalpy_i - Pressure_i / Density_i;
   Energy_j = Enthalpy_j - Pressure_j / Density_j;
 
-  SoundSpeed_i = sqrt( (Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One );
-  SoundSpeed_j = sqrt( (Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One );
+  SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One);
+  SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One);
 
   /*--- Projected velocities ---*/
 
-  ProjVelocity_i = 0;
-  ProjVelocity_j = 0;
-
-  for (iDim = 0; iDim < nDim; iDim++) {
-    ProjVelocity_i += Velocity_i[iDim] * UnitNormal[iDim];
-    ProjVelocity_j += Velocity_j[iDim] * UnitNormal[iDim];
-  }
+  ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal);
+  ProjVelocity_j = GeometryToolbox::DotProduct(nDim, Velocity_j, UnitNormal);
 
   /*--- Projected Grid Velocity ---*/
 
@@ -146,7 +135,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
   if (dynamic_grid) {
 
     for (iDim = 0; iDim < nDim; iDim++)
-      ProjInterfaceVel += 0.5 * ( GridVel_i[iDim] + GridVel_j[iDim] )*UnitNormal[iDim];
+      ProjInterfaceVel += 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * UnitNormal[iDim];
 
     SoundSpeed_i -= ProjInterfaceVel;
     SoundSpeed_j += ProjInterfaceVel;
@@ -159,14 +148,11 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
 
   Rrho = ( sqrt(Density_i) + sqrt(Density_j) );
 
-  sq_velRoe        = 0.0;
-  RoeProjVelocity  = - ProjInterfaceVel;
-
   for (iDim = 0; iDim < nDim; iDim++) {
     RoeVelocity[iDim] = ( Velocity_i[iDim] * sqrt(Density_i) + Velocity_j[iDim] * sqrt(Density_j) ) / Rrho;
-    sq_velRoe        +=  RoeVelocity[iDim] * RoeVelocity[iDim];
-    RoeProjVelocity  +=  RoeVelocity[iDim] * UnitNormal[iDim];
   }
+  sq_velRoe = GeometryToolbox::SquaredNorm(nDim, RoeVelocity);
+  RoeProjVelocity = GeometryToolbox::DotProduct(nDim, RoeVelocity, UnitNormal) - ProjInterfaceVel;
 
   /*--- Mean Roe variables iPoint and jPoint ---*/
 
@@ -175,10 +161,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
 
   /*--- Roe-averaged speed of sound ---*/
 
-  //RoeSoundSpeed2 = Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe );
   RoeSoundSpeed  = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe  ) ) - ProjInterfaceVel;
 
-
   /*--- Speed of sound at L and R ---*/
 
   sL = min( RoeProjVelocity - RoeSoundSpeed, ProjVelocity_i - SoundSpeed_i);
@@ -214,8 +198,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
 
       IntermediateState[0] = rhoSL * Density_i;
       for (iDim = 0; iDim < nDim; iDim++)
-        IntermediateState[iDim+1] = rhoSL * ( Density_i * Velocity_i[iDim] + ( pStar - Pressure_i ) / ( sL - ProjVelocity_i ) * UnitNormal[iDim] ) ;
-      IntermediateState[nVar-1] = rhoSL * ( Density_i * Energy_i - ( Pressure_i * ProjVelocity_i - pStar * sM) / ( sL - ProjVelocity_i ) );
+        IntermediateState[iDim+1] = rhoSL * Density_i * Velocity_i[iDim] + (pStar - Pressure_i) * UnitNormal[iDim] / (sL - sM);
+      IntermediateState[nVar-1] = rhoSL * Density_i * Energy_i - (Pressure_i * ProjVelocity_i - pStar * sM) / (sL - sM);
 
 
       Flux[0] = sM * IntermediateState[0];
@@ -243,8 +227,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
 
       IntermediateState[0] = rhoSR * Density_j;
       for (iDim = 0; iDim < nDim; iDim++)
-        IntermediateState[iDim+1] = rhoSR * ( Density_j * Velocity_j[iDim] + ( pStar - Pressure_j ) / ( sR - ProjVelocity_j ) * UnitNormal[iDim] ) ;
-      IntermediateState[nVar-1] = rhoSR * ( Density_j * Energy_j - ( Pressure_j * ProjVelocity_j - pStar * sM ) / ( sR - ProjVelocity_j ) );
+        IntermediateState[iDim+1] = rhoSR * Density_j * Velocity_j[iDim] + (pStar - Pressure_j) * UnitNormal[iDim] / (sR - sM);
+      IntermediateState[nVar-1] = rhoSR * Density_j * Energy_j - (Pressure_j * ProjVelocity_j - pStar * sM) / (sR - sM);
 
 
       Flux[0] = sM * IntermediateState[0];
@@ -254,14 +238,12 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
     }
   }
 
-
-  for (iVar = 0; iVar < nVar; iVar++)
-    Flux[iVar] *= Area;
+  for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area;
 
   /*--- Return early if the Jacobians do not need to be computed. ---*/
 
-  if (implicit)
-  {
+  if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
+
   if (sM > 0.0) {
 
     if (sL > 0.0) {
@@ -537,27 +519,20 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config)
     }
   }
 
-
-  /*--- Jacobians of the inviscid flux, scale = k because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/
-
-  Area *= kappa;
-
+  /*--- Scale Jacobians by area (from Flux *= Area). ---*/
   for (iVar = 0; iVar < nVar; iVar++) {
     for (jVar = 0; jVar < nVar; jVar++) {
-      Jacobian_i[iVar][jVar] *=   Area;
-      Jacobian_j[iVar][jVar] *=   Area;
+      Jacobian_i[iVar][jVar] *= Area;
+      Jacobian_j[iVar][jVar] *= Area;
     }
   }
-  } // end if implicit
-
   return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
-
 }
 
 CUpwGeneralHLLC_Flow::CUpwGeneralHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) {
 
   implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-  kappa = config->GetRoe_Kappa();
+
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
 
@@ -792,13 +767,12 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
     }
   }
 
-  for (iVar = 0; iVar < nVar; iVar++)
-    Flux[iVar] *= Area;
+  for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area;
 
   /*--- Return early if the Jacobians do not need to be computed. ---*/
 
-  if (implicit)
-  {
+  if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
+
   if (sM > 0.0) {
 
     if (sL > 0.0) {
@@ -1090,21 +1064,14 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c
     }
   }
 
-
-  /*--- Jacobians of the inviscid flux, scale = kappa because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/
-
-  Area *= kappa;
-
+  /*--- Scale Jacobians by area (from Flux *= Area). ---*/
   for (iVar = 0; iVar < nVar; iVar++) {
     for (jVar = 0; jVar < nVar; jVar++) {
       Jacobian_i[iVar][jVar] *= Area;
       Jacobian_j[iVar][jVar] *= Area;
     }
   }
-  } // end if implicit
-
   return ResidualType<>(Flux, Jacobian_i, Jacobian_j);
-
 }
 
 void CUpwGeneralHLLC_Flow::VinokurMontagne() {

SU2_CFD/src/numerics/flow/convection/roe.cpp

@@ -34,7 +34,7 @@ CUpwRoeBase_Flow::CUpwRoeBase_Flow(unsigned short val_nDim, unsigned short val_n
   implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
-  kappa = config->GetRoe_Kappa(); // 1 is unstable
+  kappa = config->GetRoe_Kappa();
 
   Gamma = config->GetGamma();
   Gamma_Minus_One = Gamma - 1.0;
@@ -683,7 +683,7 @@ CUpwGeneralRoe_Flow::CUpwGeneralRoe_Flow(unsigned short val_nDim, unsigned short
   implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
-  kappa = config->GetRoe_Kappa(); // 1 is unstable
+  kappa = config->GetRoe_Kappa();
 
 
   Flux = new su2double [nVar];

TestCases/hybrid_regression.py

@@ -67,7 +67,7 @@ def main():
     wedge.cfg_dir   = "euler/wedge"
     wedge.cfg_file  = "inv_wedge_HLLC.cfg"
     wedge.test_iter = 20
-    wedge.test_vals = [-1.396962, 4.262003, -0.244219, 0.043052]
+    wedge.test_vals = [-1.368091, 4.302736, -0.243433, 0.042906]
     test_list.append(wedge)
 
     # ONERA M6 Wing
@@ -121,8 +121,8 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals         = [-6.850130, -1.388096, -0.056036, 108.140809, 0.007988]
-    cylinder_lowmach.test_vals_aarch64 = [-6.850130, -1.388096, -0.056036, 108.140813, 0.007988]
+    cylinder_lowmach.test_vals         = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
+    cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)

TestCases/parallel_regression.py

@@ -236,7 +236,7 @@ def main():
     wedge.cfg_dir   = "euler/wedge"
     wedge.cfg_file  = "inv_wedge_HLLC.cfg"
     wedge.test_iter = 20
-    wedge.test_vals = [-1.406716, 4.253025, -0.244411, 0.043089]
+    wedge.test_vals = [-1.377543, 4.293870, -0.243566, 0.042930]
     test_list.append(wedge)
 
     # ONERA M6 Wing
@@ -282,7 +282,7 @@ def main():
     ea_naca64206.cfg_dir   = "optimization_euler/equivalentarea_naca64206"
     ea_naca64206.cfg_file  = "NACA64206.cfg"
     ea_naca64206.test_iter = 10
-    ea_naca64206.test_vals = [-1.151334, -0.455927, -0.003879, 67775.000000]
+    ea_naca64206.test_vals = [-1.188459, -0.522783, -0.003147, 67775.000000]
     test_list.append(ea_naca64206)
 
     # SUPERSONIC FLOW PAST A RAMP IN A CHANNEL
@@ -327,7 +327,7 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals = [-6.858484, -1.396528, -1.854558, 110.033249, 0.001951]
+    cylinder_lowmach.test_vals = [-6.841604, -1.379532, -1.266739, 76.118218, 0.000274]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)
@@ -392,15 +392,15 @@ def main():
     turb_flatplate_species.test_vals = [-4.120225, -0.634021, -1.706720, 1.363240, -3.250204, 9.000000, -6.697079, 5.000000, -6.978731, 10.000000, -6.013196, 0.996237, 0.996237]
     test_list.append(turb_flatplate_species)
 
-    # Flat plate SST compressibility  correction Wilcox
+    # Flat plate SST compressibility correction Wilcox
     turb_flatplate_CC_Wilcox = TestCase('turb_flatplate_CC_Wilcox')
     turb_flatplate_CC_Wilcox.cfg_dir   = "rans/flatplate"
     turb_flatplate_CC_Wilcox.cfg_file  = "turb_SST_flatplate_compressibility_Wilcox.cfg"
     turb_flatplate_CC_Wilcox.test_iter = 20
     turb_flatplate_CC_Wilcox.test_vals = [-1.280875, 1.974212, 1.440458, 5.038402, -4.051125, 8.521136]
     test_list.append(turb_flatplate_CC_Wilcox)
 
-    # Flat plate SST compressibility  correction Sarkar
+    # Flat plate SST compressibility correction Sarkar
     turb_flatplate_CC_Sarkar = TestCase('turb_flatplate_CC_Sarkar')
     turb_flatplate_CC_Sarkar.cfg_dir   = "rans/flatplate"
     turb_flatplate_CC_Sarkar.cfg_file  = "turb_SST_flatplate_compressibility_Sarkar.cfg"

TestCases/parallel_regression_AD.py

@@ -71,7 +71,7 @@ def main():
     ea_naca64206.cfg_dir      = "optimization_euler/equivalentarea_naca64206"
     ea_naca64206.cfg_file     = "NACA64206.cfg"
     ea_naca64206.test_iter    = 10
-    ea_naca64206.test_vals    = [3.182170, 2.473052, -5509000.000000, 5.551800]
+    ea_naca64206.test_vals    = [3.127605, 2.411805, -5505700.000000, 10.591000]
     test_list.append(ea_naca64206)
 
     ####################################

TestCases/serial_regression.py

@@ -109,7 +109,7 @@ def main():
     wedge.cfg_dir   = "euler/wedge"
     wedge.cfg_file  = "inv_wedge_HLLC.cfg"
     wedge.test_iter = 20
-    wedge.test_vals = [-1.396962, 4.262003, -0.244219, 0.043052]
+    wedge.test_vals = [-1.368091, 4.302736, -0.243433, 0.042906]
     test_list.append(wedge)
 
     # ONERA M6 Wing
@@ -182,7 +182,7 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals = [-6.850123, -1.388088, -0.056090, 108.140176, 0.007983]
+    cylinder_lowmach.test_vals = [-6.830989, -1.368842, -0.143838, 73.962440, 0.002454]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)

TestCases/tutorials.py

@@ -145,7 +145,7 @@ def main():
     tutorial_inv_wedge.cfg_dir    = "../Tutorials/compressible_flow/Inviscid_Wedge"
     tutorial_inv_wedge.cfg_file   = "inv_wedge_HLLC.cfg"
     tutorial_inv_wedge.test_iter  = 0
-    tutorial_inv_wedge.test_vals  = [-0.864206, 4.850246, -0.259185, 0.045567]
+    tutorial_inv_wedge.test_vals  = [-0.864206, 4.850246, -0.245674, 0.043209]
     tutorial_inv_wedge.no_restart = True
     test_list.append(tutorial_inv_wedge)
 

config_template.cfg

@@ -1567,7 +1567,7 @@ CONV_NUM_METHOD_FLOW= ROE
 % Roe Low Dissipation function for Hybrid RANS/LES simulations (FD, NTS, NTS_DUCROS)
 ROE_LOW_DISSIPATION= FD
 %
-% Roe coefficient
+% Roe dissipation coefficient
 ROE_KAPPA= 0.5
 %
 % Minimum value for beta for the Roe-Turkel preconditioner