Merge pull request #1784 from su2code/fix_sst_axisymmetry

Fix axisymmetric divergence for SST production

Author: pcarruscag

.github/workflows/regression.yml

@@ -128,7 +128,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:220614-1237
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop_tmp -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:220614-1237
         with:

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -610,14 +610,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   }
 
   /*!
-   * \brief Add contribution due to axisymmetric formulation to 2D residual
+   * \brief Add contribution from convection and diffusion due to axisymmetric formulation to 2D residual
    */
-  inline void ResidualAxisymmetric(su2double alfa_blended, su2double zeta) {
+  inline void ResidualAxisymmetricConvectionDiffusion(su2double alfa_blended, su2double zeta) {
     if (Coord_i[1] < EPS) return;
 
-    AD::SetPreaccIn(Coord_i[1]);
-    AD::SetPreaccIn(V_i[idx.Velocity() + 1]);
-
     const su2double yinv = 1.0 / Coord_i[1];
     const su2double rhov = Density_i * V_i[idx.Velocity() + 1];
     const su2double& k = ScalarVar_i[0];
@@ -627,18 +624,20 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
     const su2double sigma_k_i = F1_i * sigma_k_1 + (1.0 - F1_i) * sigma_k_2;
     const su2double sigma_w_i = F1_i * sigma_w_1 + (1.0 - F1_i) * sigma_w_2;
 
-    /*--- Production ---*/
-    const su2double pk_axi = max(
-        0.0, 2.0 / 3.0 * rhov * k * ((2.0 * yinv * V_i[idx.Velocity() + 1] - PrimVar_Grad_i[idx.Velocity()+1][1] - PrimVar_Grad_i[idx.Velocity()][0]) / zeta - 1.0));
-    const su2double pw_axi = alfa_blended * zeta / k * pk_axi;
-
     /*--- Convection-Diffusion ---*/
     const su2double cdk_axi = rhov * k - (Laminar_Viscosity_i + sigma_k_i * Eddy_Viscosity_i) * ScalarVar_Grad_i[0][1];
     const su2double cdw_axi = rhov * w - (Laminar_Viscosity_i + sigma_w_i * Eddy_Viscosity_i) * ScalarVar_Grad_i[1][1];
 
     /*--- Add terms to the residuals ---*/
-    Residual[0] += yinv * Volume * (pk_axi - cdk_axi);
-    Residual[1] += yinv * Volume * (pw_axi - cdw_axi);
+ 
+    Residual[0] -= yinv * Volume * cdk_axi;
+    Residual[1] -= yinv * Volume * cdw_axi;
+
+    Jacobian_i[0][0] -= yinv * Volume * rhov;
+    Jacobian_i[0][1] -= 0.0; 
+    Jacobian_i[1][0] -= 0.0; 
+    Jacobian_i[1][1] -= yinv * Volume * rhov;
+
   }
 
  public:
@@ -716,6 +715,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
     AD::SetPreaccIn(PrimVar_Grad_i, nDim + idx.Velocity(), nDim);
     AD::SetPreaccIn(Vorticity_i, 3);
     AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
+    AD::SetPreaccIn(V_i[idx.Velocity() + 1]);
 
     Density_i = V_i[idx.Density()];
     Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
@@ -745,6 +745,10 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       su2double diverg = 0.0;
       for (unsigned short iDim = 0; iDim < nDim; iDim++)
         diverg += PrimVar_Grad_i[iDim + idx.Velocity()][iDim];
+      if (axisymmetric && Coord_i[1] > EPS) {
+        AD::SetPreaccIn(Coord_i[1]);
+        diverg += V_i[idx.Velocity() + 1] / Coord_i[1];
+      }
 
       /*--- If using UQ methodolgy, calculate production using perturbed Reynolds stress matrix ---*/
 
@@ -836,7 +840,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       /*--- Contribution due to 2D axisymmetric formulation ---*/
 
-      if (axisymmetric) ResidualAxisymmetric(alfa_blended, zeta);
+      if (axisymmetric) ResidualAxisymmetricConvectionDiffusion(alfa_blended, zeta);
 
       /*--- Implicit part ---*/
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -1004,7 +1004,7 @@ class CFVMFlowSolverBase : public CSolver {
   /*!
    * \brief Evaluate the vorticity and strain rate magnitude.
    */
-  void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh);
+  void ComputeVorticityAndStrainMag(const CConfig& config, const CGeometry *geometry, unsigned short iMesh);
 
   /*!
    * \brief Destructor.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -584,7 +584,7 @@ void CFVMFlowSolverBase<V, R>::ImplicitEuler_Iteration(CGeometry *geometry, CSol
 }
 
 template <class V, ENUM_REGIME R>
-void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) {
+void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& config, const CGeometry *geometry, unsigned short iMesh) {
 
   auto& StrainMag = nodes->GetStrainMag();
 
@@ -611,23 +611,21 @@ void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& confi
 
     /*--- Strain Magnitude ---*/
 
+    const su2double vy = nodes->GetVelocity(iPoint, 1);
+    const su2double y = geometry->nodes->GetCoord(iPoint, 1);
     AD::StartPreacc();
     AD::SetPreaccIn(VelocityGradient, nDim, nDim);
-
-    su2double Div = 0.0;
-    for (unsigned long iDim = 0; iDim < nDim; iDim++)
-      Div += VelocityGradient(iDim, iDim);
-    Div /= 3.0;
+    AD::SetPreaccIn(vy, y);
 
     StrainMag(iPoint) = 0.0;
 
     /*--- Add diagonal part ---*/
 
     for (unsigned long iDim = 0; iDim < nDim; iDim++) {
-      StrainMag(iPoint) += pow(VelocityGradient(iDim, iDim) - Div, 2);
+      StrainMag(iPoint) += pow(VelocityGradient(iDim, iDim), 2);
     }
-    if (nDim == 2) {
-      StrainMag(iPoint) += pow(Div, 2);
+    if (config.GetAxisymmetric() && y > EPS) {
+      StrainMag(iPoint) += pow(vy / y, 2);
     }
 
     /*--- Add off diagonals ---*/

SU2_CFD/include/variables/CVariable.hpp

@@ -1133,20 +1133,6 @@ class CVariable {
    */
   inline virtual su2double GetThermalConductivity_ve(unsigned long iPoint) const { return 0.0; }
 
-  /*!
-   * \brief A virtual member.
-   * \param[in] iPoint - Point index.
-   * \return Sets separation intermittency
-   */
-  inline virtual void SetGammaSep(unsigned long iPoint, su2double gamma_sep) {}
-
-  /*!
-   * \brief A virtual member.
-   * \param[in] iPoint - Point index.
-   * \return Sets separation intermittency
-   */
-  inline virtual void SetGammaEff(unsigned long iPoint) {}
-
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.
@@ -1408,29 +1394,10 @@ class CVariable {
    */
   inline virtual su2double *GetdTvedU(unsigned long iPoint) { return nullptr; }
 
-  /*!
-   * \brief A virtual member.
-   * \param[in] val_velocity - Value of the velocity.
-   * \param[in] Gamma - Ratio of Specific heats
-   */
-  inline virtual void SetDeltaPressure(unsigned long iPoint, const su2double *val_velocity, su2double Gamma) {}
-
-  /*!
-   * \brief A virtual member.
-   * \param[in] Gamma - Ratio of specific heats.
-   */
-  inline virtual bool SetSoundSpeed(unsigned long iPoint, su2double Gamma) { return false; }
-
-  /*!
-   * \brief A virtual member.
-   * \param[in] config - Configuration parameters.
-   */
-  inline virtual bool SetSoundSpeed(unsigned long iPoint, CConfig *config) { return false; }
-
   /*!
    * \brief A virtual member.
    */
-  inline virtual bool SetSoundSpeed(unsigned long iPoint) { return false; }
+  inline virtual bool SetSoundSpeed(unsigned long iPoint, su2double soundspeed2) { return false; }
 
   /*!
    * \brief A virtual member.

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -102,7 +102,7 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container
     SetPrimitive_Limiter(geometry, config);
   }
 
-  ComputeVorticityAndStrainMag(*config, iMesh);
+  ComputeVorticityAndStrainMag(*config, geometry, iMesh);
 
   /*--- Compute the TauWall from the wall functions ---*/
 

SU2_CFD/src/solvers/CNEMONSSolver.cpp

@@ -98,7 +98,7 @@ void CNEMONSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
 
   /*--- Compute vorticity and strain mag. ---*/
 
-  ComputeVorticityAndStrainMag(*config, iMesh);
+  ComputeVorticityAndStrainMag(*config, geometry, iMesh);
 
   /*--- Compute the TauWall from the wall functions ---*/
 

SU2_CFD/src/solvers/CNSSolver.cpp

@@ -116,7 +116,7 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C
     SetPrimitive_Limiter(geometry, config);
   }
 
-  ComputeVorticityAndStrainMag(*config, iMesh);
+  ComputeVorticityAndStrainMag(*config, geometry, iMesh);
 
   /*--- Compute the TauWall from the wall functions ---*/
 

…symmetric_rans/air_nozzle/air_nozzle.cfg …c_rans/air_nozzle/air_nozzle_restart.cfgTestCases/axisymmetric_rans/air_nozzle/air_nozzle.cfg renamed to TestCases/axisymmetric_rans/air_nozzle/air_nozzle_restart.cfg TestCases/axisymmetric_rans/air_nozzle/air_nozzle.cfg renamed to TestCases/axisymmetric_rans/air_nozzle/air_nozzle_restart.cfg

@@ -86,7 +86,7 @@ CFL_REDUCTION_TURB= 1.0
 %
 % --------------------------- CONVERGENCE PARAMETERS --------------------------%
 %
-ITER= 1000
+ITER= 3000
 CONV_RESIDUAL_MINVAL= -12
 CONV_STARTITER= 10
 %
@@ -95,6 +95,7 @@ CONV_STARTITER= 10
 MESH_FILENAME= nozzle.su2
 %
 SOLUTION_FILENAME= solution_flow.dat
+RESTART_FILENAME= restart_flow.dat
 OUTPUT_WRT_FREQ= 1000
 %
 % Note: This cfg is used for a primal and adjoint Testcase, therefore both residuals are present here.

TestCases/hybrid_regression.py

@@ -166,15 +166,15 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-0.510635, 4.871104, 0.811904, 0.061614, -82395.000000]
+    rae2822_sst.test_vals = [-0.510634, 4.869669, 0.811906, 0.061614, -82395.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
     rae2822_sst_sust           = TestCase('rae2822_sst_sust')
     rae2822_sst_sust.cfg_dir   = "rans/rae2822"
     rae2822_sst_sust.cfg_file  = "turb_SST_SUST_RAE2822.cfg"
     rae2822_sst_sust.test_iter = 20
-    rae2822_sst_sust.test_vals = [-2.430589, 4.871104, 0.811903, 0.061614]
+    rae2822_sst_sust.test_vals = [-2.420843, 4.869669, 0.811906, 0.061614]
     test_list.append(rae2822_sst_sust)
 
     # Flat plate
@@ -207,7 +207,7 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals         = [-11.450472, -12.797872, -5.863655, 1.049989, 0.019163, -1.856265, -38.694000]
+    turb_naca0012_sst.test_vals         = [-11.418415, -12.801780, -5.864804, 1.049990, 0.019163, -1.669899, -38.695000]
     turb_naca0012_sst.test_vals_aarch64 = [-11.450473, -12.797872, -5.863655, 1.049989, 0.019163, -1.856266, -38.694000]
     test_list.append(turb_naca0012_sst)
 
@@ -216,7 +216,7 @@ def main():
     turb_naca0012_sst_sust.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_sust.cfg_file  = "turb_NACA0012_sst_sust.cfg"
     turb_naca0012_sst_sust.test_iter = 10
-    turb_naca0012_sst_sust.test_vals         = [-11.367052, -12.640670, -5.746919, 1.005233, 0.019017, -1.913902]
+    turb_naca0012_sst_sust.test_vals         = [-11.361657, -12.642245, -5.747099, 1.005234, 0.019017, -1.573783]
     turb_naca0012_sst_sust.test_vals_aarch64 = [-11.367052, -12.640670, -5.746919, 1.005233, 0.019017, -1.913907]
     test_list.append(turb_naca0012_sst_sust)
 
@@ -225,7 +225,7 @@ def main():
     turb_naca0012_sst_fixedvalues.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_fixedvalues.cfg_file  = "turb_NACA0012_sst_fixedvalues.cfg"
     turb_naca0012_sst_fixedvalues.test_iter = 10
-    turb_naca0012_sst_fixedvalues.test_vals         = [-5.192502, -9.575898, -1.568269, 1.022571, 0.040527, -2.384329]
+    turb_naca0012_sst_fixedvalues.test_vals         = [-5.192502, -9.575351, -1.568269, 1.022571, 0.040527, -2.384334]
     turb_naca0012_sst_fixedvalues.test_vals_aarch64 = [-5.192501, -9.575897, -1.568269, 1.022571, 0.040527, -2.384327]
     test_list.append(turb_naca0012_sst_fixedvalues)
 
@@ -250,13 +250,13 @@ def main():
     #######################################
 
     # Axisymmetric air nozzle (transonic)
-    axi_rans_air_nozzle           = TestCase('axi_rans_air_nozzle')
-    axi_rans_air_nozzle.cfg_dir   = "axisymmetric_rans/air_nozzle"
-    axi_rans_air_nozzle.cfg_file  = "air_nozzle.cfg"
-    axi_rans_air_nozzle.test_iter = 10
-    axi_rans_air_nozzle.test_vals         = [-12.093553, -6.630383, -8.798737, -2.399130, -1938.200000]
-    axi_rans_air_nozzle.test_vals_aarch64 = [-12.093539, -6.630357, -8.798732, -2.399130, -1938.200000]
-    test_list.append(axi_rans_air_nozzle)
+    axi_rans_air_nozzle_restart           = TestCase('axi_rans_air_nozzle_restart')
+    axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
+    axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
+    axi_rans_air_nozzle_restart.test_iter = 10
+    axi_rans_air_nozzle_restart.test_vals         = [-12.094405, -6.649212, -8.899548, -2.491859, -1924.900000]
+    axi_rans_air_nozzle_restart.test_vals_aarch64 = [-12.093539, -6.630357, -8.798732, -2.399130, -1938.200000]
+    test_list.append(axi_rans_air_nozzle_restart)
 
     #################################
     ## Compressible RANS Restart  ###
@@ -424,7 +424,7 @@ def main():
     inc_turb_naca0012_sst_sust.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012_sst_sust.cfg_file  = "naca0012_SST_SUST.cfg"
     inc_turb_naca0012_sst_sust.test_iter = 20
-    inc_turb_naca0012_sst_sust.test_vals = [-7.276424, 0.145860, 0.000003, 0.312011]
+    inc_turb_naca0012_sst_sust.test_vals = [-7.274028, 0.145888, 0.000002, 0.312014]
     test_list.append(inc_turb_naca0012_sst_sust)
 
     ######################################
@@ -457,7 +457,7 @@ def main():
     square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     square_cylinder.cfg_file  = "turb_square.cfg"
     square_cylinder.test_iter = 3
-    square_cylinder.test_vals = [-2.558001, -1.162564, 0.066372, 1.399788, 2.220402, 1.399743, 2.218603, -0.453110]
+    square_cylinder.test_vals = [-2.558128, -1.158062, 0.067932, 1.399788, 2.220402, 1.399743, 2.218603, -0.453110]
     square_cylinder.test_vals_aarch64 = [-2.558106, -1.162563, 0.066386, 1.399788, 2.220402, 1.399743, 2.218603, -0.453110]
     square_cylinder.unsteady  = True
     test_list.append(square_cylinder)
@@ -534,53 +534,34 @@ def main():
     ### Turbomachinery                 ###
     ######################################
 
-    # Jones APU Turbocharger
-    Jones_tc           = TestCase('jones_turbocharger')
-    Jones_tc.cfg_dir   = "turbomachinery/APU_turbocharger"
-    Jones_tc.cfg_file  = "Jones.cfg"
-    Jones_tc.test_iter = 5
-    Jones_tc.test_vals = [-5.279930, 0.379651, 72.212100, 1.277439]
-    Jones_tc.new_output = False
-    test_list.append(Jones_tc)
-
 	# Jones APU Turbocharger restart
-    Jones_tc_rst           = TestCase('jones_turbocharger_restart')
-    Jones_tc_rst.cfg_dir   = "turbomachinery/APU_turbocharger"
-    Jones_tc_rst.cfg_file  = "Jones_rst.cfg"
-    Jones_tc_rst.test_iter = 5
-    Jones_tc_rst.test_vals = [-4.625251, -1.568824, 33.995140, 10.181940]
-    Jones_tc_rst.new_output = False
-    test_list.append(Jones_tc_rst)
+    Jones_tc_restart           = TestCase('jones_turbocharger_restart')
+    Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
+    Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
+    Jones_tc_restart.test_iter = 5
+    Jones_tc_restart.test_vals = [-10.691521, -7.643723, 85.827890, 2.277151]
+    Jones_tc_restart.new_output = False
+    test_list.append(Jones_tc_restart)
 
     # 2D axial stage
     axial_stage2D           = TestCase('axial_stage2D')
     axial_stage2D.cfg_dir   = "turbomachinery/axial_stage_2D"
     axial_stage2D.cfg_file  = "Axial_stage2D.cfg"
     axial_stage2D.test_iter = 20
-    axial_stage2D.test_vals         = [-1.933115, 5.365584, 73.354510, 0.925901]
+    axial_stage2D.test_vals         = [-1.937003, 5.338805, 73.357170, 0.915734]
     axial_stage2D.test_vals_aarch64 = [-1.933115, 5.365583, 73.354510, 0.925902]
     axial_stage2D.new_output = False
     test_list.append(axial_stage2D)
 
-    # 2D transonic stator
-    transonic_stator           = TestCase('transonic_stator')
-    transonic_stator.cfg_dir   = "turbomachinery/transonic_stator_2D"
-    transonic_stator.cfg_file  = "transonic_stator.cfg"
-    transonic_stator.test_iter = 20
-    transonic_stator.test_vals         = [-0.565608, 5.833408, 96.476150, 0.062517]
-    transonic_stator.test_vals_aarch64 = [-0.565610, 5.833410, 96.476160, 0.062517]
-    transonic_stator.new_output = False
-    test_list.append(transonic_stator)
-
     # 2D transonic stator restart
-    transonic_stator_rst           = TestCase('transonic_stator_restart')
-    transonic_stator_rst.cfg_dir   = "turbomachinery/transonic_stator_2D"
-    transonic_stator_rst.cfg_file  = "transonic_stator_rst.cfg"
-    transonic_stator_rst.test_iter = 20
-    transonic_stator_rst.test_vals         = [-6.619122, -0.615690, 5.002986, 0.002951]
-    transonic_stator_rst.test_vals_aarch64 = [-6.619122, -0.615705, 5.002986, 0.002951]
-    transonic_stator_rst.new_output = False
-    test_list.append(transonic_stator_rst)
+    transonic_stator_restart           = TestCase('transonic_stator_restart')
+    transonic_stator_restart.cfg_dir   = "turbomachinery/transonic_stator_2D"
+    transonic_stator_restart.cfg_file  = "transonic_stator_restart.cfg"
+    transonic_stator_restart.test_iter = 20
+    transonic_stator_restart.test_vals         = [-6.796363, -0.743837, 5.003453, 0.002946]
+    transonic_stator_restart.test_vals_aarch64 = [-6.619122, -0.615705, 5.002986, 0.002951]
+    transonic_stator_restart.new_output = False
+    test_list.append(transonic_stator_restart)
 
     ######################################
     ### Sliding Mesh                   ###
@@ -652,7 +633,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.619686, -1.564595]
+    bars_SST_2D.test_vals = [13.000000, -0.609170, -1.523885]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)