Merge pull request #1067 from su2code/fix_inc_rotatingframe

Fixes for incompressible solver - rotating frame and convergence rate for unsteady problems

Author: TobiKattmann

Common/include/CConfig.hpp

@@ -600,8 +600,9 @@ class CConfig {
   Kappa_4th_Flow,           /*!< \brief JST 4th order dissipation coefficient for flow equations. */
   Kappa_2nd_Heat,           /*!< \brief 2nd order dissipation coefficient for heat equation. */
   Kappa_4th_Heat,           /*!< \brief 4th order dissipation coefficient for heat equation. */
-  Cent_Jac_Fix_Factor;              /*!< \brief Multiply the dissipation contribution to the Jacobian of central schemes
+  Cent_Jac_Fix_Factor,              /*!< \brief Multiply the dissipation contribution to the Jacobian of central schemes
                                                 by this factor to make the global matrix more diagonal dominant. */
+  Cent_Inc_Jac_Fix_Factor;          /*!< \brief Multiply the dissipation contribution to the Jacobian of incompressible central schemes */
   su2double Geo_Waterline_Location; /*!< \brief Location of the waterline. */
 
   su2double Min_Beta_RoeTurkel,     /*!< \brief Minimum value of Beta for the Roe-Turkel low Mach preconditioner. */
@@ -4625,6 +4626,12 @@ class CConfig {
    */
   su2double GetCent_Jac_Fix_Factor(void) const { return Cent_Jac_Fix_Factor; }
 
+  /*!
+   * \brief Factor by which to multiply the dissipation contribution to Jacobians of incompressible central schemes.
+   * \return The factor.
+   */
+  su2double GetCent_Inc_Jac_Fix_Factor(void) const { return Cent_Inc_Jac_Fix_Factor; }
+
   /*!
    * \brief Get the kind of integration scheme (explicit or implicit)
    *        for the adjoint flow equations.

Common/src/CConfig.cpp

@@ -1819,6 +1819,8 @@ void CConfig::SetConfig_Options() {
   addBoolOption("USE_ACCURATE_FLUX_JACOBIANS", Use_Accurate_Jacobians, false);
   /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Improve the numerical properties (diagonal dominance) of the global Jacobian matrix, 3 to 4 is "optimum" (central schemes) \ingroup Config*/
   addDoubleOption("CENTRAL_JACOBIAN_FIX_FACTOR", Cent_Jac_Fix_Factor, 4.0);
+  /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Control numerical properties of the global Jacobian matrix using a multiplication factor for incompressible central schemes \ingroup Config*/
+  addDoubleOption("CENTRAL_INC_JACOBIAN_FIX_FACTOR", Cent_Inc_Jac_Fix_Factor, 1.0);
 
   /*!\brief CONV_NUM_METHOD_ADJFLOW
    *  \n DESCRIPTION: Convective numerical method for the adjoint solver.
@@ -4646,12 +4648,6 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
     }
   }
 
-  /*--- Rotating frame is not yet supported with the incompressible solver. ---*/
-
-  if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (Kind_GridMovement == ROTATING_FRAME)) {
-    SU2_MPI::Error("Support for rotating frame simulation not yet implemented for incompressible flows.", CURRENT_FUNCTION);
-  }
-
   /*--- Assert that there are two markers being analyzed if the
    pressure drop objective function is selected. ---*/
 

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

@@ -252,6 +252,8 @@ class CCentLaxInc_Flow final : public CNumerics {
   variable_density,                /*!< \brief Variable density incompressible flows. */
   energy;                          /*!< \brief computation with the energy equation. */
 
+  su2double fix_factor;            /*!< \brief Fix factor for Jacobians. */
+
   su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */
   su2double** Jacobian_j = nullptr; /*!< \brief The Jacobian w.r.t. point j after computation. */
 
@@ -312,6 +314,8 @@ class CCentJSTInc_Flow final : public CNumerics {
   variable_density,      /*!< \brief Variable density incompressible flows. */
   energy;                /*!< \brief computation with the energy equation. */
 
+  su2double fix_factor;  /*!< \brief Fix factor for Jacobians. */
+
   su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */
   su2double** Jacobian_j = nullptr; /*!< \brief The Jacobian w.r.t. point j after computation. */
 

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

@@ -344,6 +344,7 @@ CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_n
   variable_density = (config->GetKind_DensityModel() == VARIABLE);
   /* 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();
+  fix_factor = config->GetCent_Inc_Jac_Fix_Factor();
   energy           = config->GetEnergy_Equation();
 
   /*--- Artificial dissipation part ---*/
@@ -531,8 +532,8 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi
     for (jVar = 0; jVar < nVar; jVar++) {
       ProjFlux[iVar] += Precon[iVar][jVar]*Epsilon_0*Diff_V[jVar]*StretchingFactor*MeanLambda;
       if (implicit) {
-        Jacobian_i[iVar][jVar] += Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda;
-        Jacobian_j[iVar][jVar] -= Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda;
+        Jacobian_i[iVar][jVar] += fix_factor*Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda;
+        Jacobian_j[iVar][jVar] -= fix_factor*Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda;
       }
     }
   }
@@ -564,6 +565,7 @@ CCentJSTInc_Flow::CCentJSTInc_Flow(unsigned short val_nDim, unsigned short val_n
   energy           = config->GetEnergy_Equation();
   /* 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();
+  fix_factor = config->GetCent_Inc_Jac_Fix_Factor();
 
   /*--- Artifical dissipation part ---*/
 
@@ -759,8 +761,8 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi
     for (jVar = 0; jVar < nVar; jVar++) {
       ProjFlux[iVar] += Precon[iVar][jVar]*(Epsilon_2*Diff_V[jVar] - Epsilon_4*Diff_Lapl[jVar])*StretchingFactor*MeanLambda;
       if (implicit) {
-        Jacobian_i[iVar][jVar] += Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_i+1))*StretchingFactor*MeanLambda;
-        Jacobian_j[iVar][jVar] -= Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_j+1))*StretchingFactor*MeanLambda;
+        Jacobian_i[iVar][jVar] += fix_factor*Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_i+1))*StretchingFactor*MeanLambda;
+        Jacobian_j[iVar][jVar] -= fix_factor*Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_j+1))*StretchingFactor*MeanLambda;
       }
     }
   }

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1527,9 +1527,9 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
 
     for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
-      /*--- Load the conservative variables ---*/
+      /*--- Load the primitive variables ---*/
 
-      numerics->SetConservative(nodes->GetSolution(iPoint), nullptr);
+      numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr);
 
       /*--- Set incompressible density ---*/
 
@@ -3051,33 +3051,19 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       LinSysRes.AddBlock(iPoint, Residual);
 
       if (implicit) {
-
-        SetPreconditioner(config, iPoint);
-        for (iVar = 0; iVar < nVar; iVar++) {
-          for (jVar = 0; jVar < nVar; jVar++) {
-            Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar];
-          }
-        }
-
-        for (iVar = 0; iVar < nVar; iVar++) {
-          for (jVar = 0; jVar < nVar; jVar++) {
-            if (config->GetTime_Marching() == DT_STEPPING_1ST)
-              Jacobian_i[iVar][jVar] *= Volume_nP1 / TimeStep;
-            if (config->GetTime_Marching() == DT_STEPPING_2ND)
-              Jacobian_i[iVar][jVar] *= (Volume_nP1*3.0)/(2.0*TimeStep);
-          }
-        }
-
-        if (!energy) {
-            for (iVar = 0; iVar < nVar; iVar++) {
-              Jacobian_i[iVar][nDim+1] = 0.0;
-              Jacobian_i[nDim+1][iVar] = 0.0;
-            }
+        for (iVar = 1; iVar < nVar; iVar++) {
+          if (config->GetTime_Marching() == DT_STEPPING_1ST)
+            Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep;
+          if (config->GetTime_Marching() == DT_STEPPING_2ND)
+            Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep);
         }
+        for (iDim = 0; iDim < nDim; iDim++)
+          Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1];
+        if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1];
 
         Jacobian.AddBlock2Diag(iPoint, Jacobian_i);
-
       }
+
     }
 
   }
@@ -3276,31 +3262,21 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
        to the dual time source term. ---*/
       if (!energy) Residual[nDim+1] = 0.0;
       LinSysRes.AddBlock(iPoint, Residual);
-      if (implicit) {
-        SetPreconditioner(config, iPoint);
-        for (iVar = 0; iVar < nVar; iVar++) {
-          for (jVar = 0; jVar < nVar; jVar++) {
-            Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar];
-          }
-        }
 
-        for (iVar = 0; iVar < nVar; iVar++) {
-          for (jVar = 0; jVar < nVar; jVar++) {
-            if (config->GetTime_Marching() == DT_STEPPING_1ST)
-              Jacobian_i[iVar][jVar] *= Volume_nP1 / TimeStep;
-            if (config->GetTime_Marching() == DT_STEPPING_2ND)
-              Jacobian_i[iVar][jVar] *= (Volume_nP1*3.0)/(2.0*TimeStep);
-          }
+      if (implicit) {
+        for (iVar = 1; iVar < nVar; iVar++) {
+          if (config->GetTime_Marching() == DT_STEPPING_1ST)
+            Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep;
+          if (config->GetTime_Marching() == DT_STEPPING_2ND)
+            Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep);
         }
+        for (iDim = 0; iDim < nDim; iDim++)
+          Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1];
+        if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1];
 
-        if (!energy) {
-          for (iVar = 0; iVar < nVar; iVar++) {
-            Jacobian_i[iVar][nDim+1] = 0.0;
-            Jacobian_i[nDim+1][iVar] = 0.0;
-          }
-        }
         Jacobian.AddBlock2Diag(iPoint, Jacobian_i);
       }
+
     }
   }
 

SU2_PY/SU2/eval/gradients.py

@@ -838,6 +838,14 @@ def findiff( config, state=None ):
         name = su2io.expand_time(name,config)
         link.extend(name)
 
+    # files: restart solution for dual-time stepping first and second order
+    if 'RESTART_FILE_1' in files:
+        name = files['RESTART_FILE_1']
+        pull.append(name)
+    if 'RESTART_FILE_2' in files:
+        name = files['RESTART_FILE_2']
+        pull.append(name)
+
     # files: target equivarea distribution
     if 'EQUIV_AREA' in special_cases and 'TARGET_EA' in files:
         pull.append(files['TARGET_EA'])

SU2_PY/finite_differences.py

@@ -77,6 +77,21 @@ def finite_differences( filename           ,
     # State
     state = SU2.io.State()
     state.find_files(config)
+
+    # add restart files to state.FILES
+    if config.get('TIME_DOMAIN', 'NO') == 'YES' and config.get('RESTART_SOL', 'NO') == 'YES':
+        restart_name = config['RESTART_FILENAME'].split('.')[0]
+        restart_filename = restart_name + '_' + str(int(config['RESTART_ITER'])-1).zfill(5) + '.dat'
+        if not os.path.isfile(restart_filename): # throw, if restart files does not exist
+            sys.exit("Error: Restart file <" + restart_filename + "> not found.")
+        state['FILES']['RESTART_FILE_1'] = restart_filename
+
+        # use only, if time integration is second order
+        if config.get('TIME_MARCHING', 'NO') == 'DUAL_TIME_STEPPING-2ND_ORDER':
+            restart_filename = restart_name + '_' + str(int(config['RESTART_ITER'])-2).zfill(5) + '.dat'
+            if not os.path.isfile(restart_filename): # throw, if restart files does not exist
+                sys.exit("Error: Restart file <" + restart_filename + "> not found.")
+            state['FILES']['RESTART_FILE_2'] =restart_filename
 
     # Finite Difference Gradients
     SU2.eval.gradients.findiff(config,state)

TestCases/parallel_regression.py

@@ -866,7 +866,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-2.991060, -6.866340, 1.476415, 0.419712]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.004011, -6.876230, 1.487888, 0.421869]
     unst_inc_turb_naca0015_sa.su2_exec  = "parallel_computation.py -f"
     unst_inc_turb_naca0015_sa.timeout   = 1600
     unst_inc_turb_naca0015_sa.tol       = 0.00001
@@ -1195,8 +1195,8 @@ def main():
     cht_incompressible_unsteady.cfg_dir   = "coupled_cht/incomp_2d_unsteady"
     cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
     cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-1.348104, -0.080392, -0.080391, -0.080391] #last 4 columns
-    cht_incompressible_unsteady.su2_exec  = "SU2_CFD"
+    cht_incompressible_unsteady.test_vals = [-1.305471, -0.080372, -0.080376, -0.080372] #last 4 columns
+    cht_incompressible_unsteady.su2_exec  = "mpirun -n 2 SU2_CFD"
     cht_incompressible_unsteady.timeout   = 1600
     cht_incompressible_unsteady.multizone = True
     cht_incompressible_unsteady.unsteady  = True

TestCases/serial_regression.py

@@ -1001,7 +1001,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-2.994996, -6.869781, 1.434864, 0.416626] #last 4 columns
+    unst_inc_turb_naca0015_sa.test_vals = [-3.007635, -6.879789, 1.445300, 0.419281] #last 4 columns
     unst_inc_turb_naca0015_sa.su2_exec  = "SU2_CFD"
     unst_inc_turb_naca0015_sa.timeout   = 1600
     unst_inc_turb_naca0015_sa.tol       = 0.00001
@@ -1407,7 +1407,7 @@ def main():
     cht_incompressible_unsteady.cfg_dir   = "coupled_cht/incomp_2d_unsteady"
     cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
     cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-1.348104, -0.080392, -0.080391, -0.080391] #last 4 columns
+    cht_incompressible_unsteady.test_vals = [-1.303588, -0.080377, -0.080380, -0.080377] #last 4 columns
     cht_incompressible_unsteady.su2_exec  = "SU2_CFD"
     cht_incompressible_unsteady.timeout   = 1600
     cht_incompressible_unsteady.multizone = True