New Nonlinear Iteration Controller for Improved Fluid Solver Robustness & Convergence

SU2_CFD/include/output_structure.hpp

@@ -784,19 +784,6 @@ class COutput {
    */
   void SpecialOutput_ForcesBreakdown(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short val_iZone, bool output);
 
-  /*!
-   * \brief Write the history file and the convergence on the screen for serial computations.
-   * \param[in] ConvHist_file - Pointer to the convergence history file (which is defined in the main subroutine).
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] integration - Generic subroutines for space integration, time integration, and monitoring.
-   * \param[in] iExtIter - Current external (time) iteration.
-   * \param[in] timeused - Current number of clock tick in the computation (related with total time).
-   * \param[in] val_nZone - iZone index.
-   */
-  void SetCFL_Number(CSolver *****solver_container, CConfig **config, unsigned short val_iZone);
-
   /*!
    * \brief Write the sensitivity (including mesh sensitivity) computed with the discrete adjoint method
    *  on the surface and in the volume to a file.

SU2_CFD/include/solver_structure.hpp

@@ -155,7 +155,7 @@ class CSolver {
   su2double ***VertexTractionAdjoint;   /*- Also temporary -*/
 
 public:
-  
+
   CSysVector<su2double> LinSysSol;    /*!< \brief vector to store iterative solution of implicit linear system. */
   CSysVector<su2double> LinSysRes;    /*!< \brief vector to store iterative residual of implicit linear system. */
   CSysVector<su2double> LinSysAux;    /*!< \brief vector to store iterative residual of implicit linear system. */
@@ -1291,6 +1291,22 @@ class CSolver {
    */
   virtual void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config);
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] solver - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void ComputeUnderRelaxationFactor(CSolver **solver, CConfig *config);
+
+  /*!
+   * \brief Adapt the CFL number based on the local under-relaxation parameters
+   *        computed for each nonlinear iteration.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   */
+  void AdaptCFLNumber(CGeometry **geometry, CSolver ***solver_container, CConfig *config);
+
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -1427,7 +1443,7 @@ class CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  virtual void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config);
+  virtual void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction = false);
 
   /*!
    * \brief A virtual member.
@@ -5049,7 +5065,7 @@ class CEulerSolver : public CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config);
+  void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction = false);
 
   /*!
    * \brief Compute the limiter of the primitive variables.
@@ -5427,6 +5443,13 @@ class CEulerSolver : public CSolver {
    */
   void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config);
 
+  /*!
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear iteration for stability.
+   * \param[in] solver - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeUnderRelaxationFactor(CSolver **solver, CConfig *config);
+
   /*!
    * \brief Compute the pressure forces and all the adimensional coefficients.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -7317,7 +7340,7 @@ class CIncEulerSolver : public CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config);
+  void SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction = false);
 
   /*!
    * \brief Compute the limiter of the primitive variables.
@@ -7496,6 +7519,13 @@ class CIncEulerSolver : public CSolver {
    */
   void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config);
 
+  /*!
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear iteration for stability.
+   * \param[in] solver - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeUnderRelaxationFactor(CSolver **solver, CConfig *config);
+
   /*!
    * \brief Provide the non dimensional lift coefficient (inviscid contribution).
    * \param val_marker Surface where the coefficient is going to be computed.
@@ -9412,6 +9442,13 @@ class CTurbSolver : public CSolver {
   void SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config,
                             unsigned short iRKStep, unsigned short iMesh, unsigned short RunTime_EqSystem);
 
+  /*!
+   * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear iteration for stability.
+   * \param[in] solver - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeUnderRelaxationFactor(CSolver **solver, CConfig *config);
+
   /*!
    * \brief Load a solution from a restart file.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solver_structure.inl

@@ -104,7 +104,7 @@ inline su2double CSolver::GetPsiE_Inf(void) { return 0; }
 
 inline void CSolver::SetPrimitive_Gradient_GG(CGeometry *geometry, CConfig *config) { }
 
-inline void CSolver::SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config) { }
+inline void CSolver::SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction) { }
 
 inline void CSolver::SetPrimitive_Limiter_MPI(CGeometry *geometry, CConfig *config) { }
 
@@ -855,6 +855,8 @@ inline void CSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solv
 
 inline void CSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { }
 
+inline void CSolver::ComputeUnderRelaxationFactor(CSolver **solver_container, CConfig *config) { }
+
 inline void CSolver::ImplicitNewmark_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { }
 
 inline void CSolver::ImplicitNewmark_Update(CGeometry *geometry, CSolver **solver_container, CConfig *config) { }

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -57,6 +57,7 @@ class CEulerVariable : public CVariable {
 
   su2double *Primitive;  /*!< \brief Primitive variables (T, vx, vy, vz, P, rho, h, c) in compressible flows. */
   su2double **Gradient_Primitive;  /*!< \brief Gradient of the primitive variables (T, vx, vy, vz, P, rho). */
+  su2double **Gradient_Reconstruction;  /*!< \brief Gradient of the primitive variables for MUSCL reconstruction for the convective term (T, vx, vy, vz, P, rho). */
   su2double *Limiter_Primitive;    /*!< \brief Limiter of the primitive variables (T, vx, vy, vz, P, rho). */
 
   /*--- Secondary variable definition ---*/
@@ -255,6 +256,28 @@ class CEulerVariable : public CVariable {
    */
   inline su2double *GetLimiter_Secondary(void) {return Limiter_Secondary; }
 
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double GetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim) {return Gradient_Reconstruction[val_var][val_dim]; }
+
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \param[in] val_value - Value of the gradient.
+   */
+  inline void SetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim, su2double val_value) {Gradient_Reconstruction[val_var][val_dim] = val_value; }
+
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double **GetGradient_Reconstruction(void) {return Gradient_Reconstruction; }
+
   /*!
    * \brief A virtual member.
    */

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -51,9 +51,11 @@ class CIncEulerVariable : public CVariable {
 
   /*--- Primitive variable definition ---*/
 
-  su2double *Primitive;  /*!< \brief Primitive variables (T, vx, vy, vz, P, rho, h, c) in compressible flows. */
-  su2double **Gradient_Primitive;  /*!< \brief Gradient of the primitive variables (T, vx, vy, vz, P, rho). */
-  su2double *Limiter_Primitive;    /*!< \brief Limiter of the primitive variables (T, vx, vy, vz, P, rho). */
+  su2double *Primitive;  /*!< \brief Primitive variables (P, vx, vy, vz, T, rho ...) for incompressible flows. */
+  su2double **Gradient_Primitive;  /*!< \brief Gradient of the primitive variables */
+  su2double **Gradient_Reconstruction;  /*!< \brief Gradient of the primitive variables for MUSCL reconstruction for the convective term */
+
+  su2double *Limiter_Primitive;    /*!< \brief Limiter of the primitive variables */
 
   /*--- Old solution container for BGS iterations ---*/
 
@@ -159,6 +161,28 @@ class CIncEulerVariable : public CVariable {
    */
   inline su2double *GetLimiter_Primitive(void) {return Limiter_Primitive; }
 
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double GetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim) {return Gradient_Reconstruction[val_var][val_dim]; }
+
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \param[in] val_value - Value of the gradient.
+   */
+  inline void SetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim, su2double val_value) {Gradient_Reconstruction[val_var][val_dim] = val_value; }
+
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double **GetGradient_Reconstruction(void) {return Gradient_Reconstruction; }
+
   /*!
    * \brief Set the value of the pressure.
    */

SU2_CFD/include/variables/CVariable.hpp

@@ -62,6 +62,7 @@ class CVariable {
   su2double *Solution,    /*!< \brief Solution of the problem. */
   *Solution_Old;      /*!< \brief Old solution of the problem R-K. */
   bool Non_Physical;      /*!< \brief Non-physical points in the solution (force first order). */
+  unsigned short Non_Physical_Counter;      /*!< \brief Number of consecutive iterations that a point has been treated first-order. After a specified number of successful reconstructions, the point can be returned to second-order. */
   su2double *Solution_time_n,  /*!< \brief Solution of the problem at time n for dual-time stepping technique. */
   *Solution_time_n1;      /*!< \brief Solution of the problem at time n-1 for dual-time stepping technique. */
   su2double **Gradient;    /*!< \brief Gradient of the solution of the problem. */
@@ -81,6 +82,9 @@ class CVariable {
   su2double *Res_TruncError,  /*!< \brief Truncation error for multigrid cycle. */
   *Residual_Old,    /*!< \brief Auxiliar structure for residual smoothing. */
   *Residual_Sum;    /*!< \brief Auxiliar structure for residual smoothing. */
+  su2double UnderRelaxation;  /*!< \brief Value of the under-relxation parameter local to the control volume. */
+  su2double LocalCFL;  /*!< \brief Value of the CFL number local to the control volume. */
+
   static unsigned short nDim;    /*!< \brief Number of dimension of the problem. */
   unsigned short nVar;    /*!< \brief Number of variables of the problem,
                            note that this variable cannnot be static, it is possible to
@@ -147,13 +151,23 @@ class CVariable {
    * \brief Set the value of the non-physical point.
    * \param[in] val_value - identification of the non-physical point.
    */
-  inline void SetNon_Physical(bool val_value) { Non_Physical = !val_value; }
+  inline void SetNon_Physical(bool val_value) {
+    if (val_value) {
+      Non_Physical = val_value;
+      Non_Physical_Counter = 0;
+    } else {
+      Non_Physical_Counter++;
+      if (Non_Physical_Counter > 20) {
+        Non_Physical = val_value;
+      }
+    }
+  }
 
   /*!
    * \brief Get the value of the non-physical point.
    * \return Value of the Non-physical point.
    */
-  inline su2double GetNon_Physical(void) { return su2double(Non_Physical); }
+  inline bool GetNon_Physical(void) { return Non_Physical; }
 
   /*!
    * \brief Get the solution.
@@ -472,6 +486,30 @@ class CVariable {
       val_residual[iVar] = Residual_Sum[iVar];
   }
 
+  /*!
+   * \brief Set the value of the under-relaxation parameter for the current control volume (CV).
+   * \param[in] val_under_relaxation - the input value of the under-relaxation parameter for this CV.
+   */
+  inline void SetUnderRelaxation(su2double val_under_relaxation) { UnderRelaxation = val_under_relaxation; }
+
+  /*!
+   * \brief Get the value of the under-relaxation parameter for the current control volume (CV).
+   * \return Value of the under-relaxation parameter for this CV.
+   */
+  inline su2double GetUnderRelaxation(void) { return UnderRelaxation; }
+
+  /*!
+   * \brief Set the value of the local CFL number for the current control volume (CV).
+   * \param[in] val_cfl - the input value of the local CFL number for this CV.
+   */
+  inline void SetLocalCFL(su2double val_cfl) { LocalCFL = val_cfl; }
+
+  /*!
+   * \brief Get the value of the local CFL number for the current control volume (CV).
+   * \return Value of the local CFL number for this CV.
+   */
+  inline su2double GetLocalCFL(void) { return LocalCFL; }
+
   /*!
    * \brief Set auxiliar variables, we are looking for the gradient of that variable.
    * \param[in] val_auxvar - Value of the auxiliar variable.
@@ -1850,6 +1888,28 @@ class CVariable {
    */
   inline virtual su2double *GetLimiter_Secondary(void) {return NULL; }
 
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \return Value of the primitive variables gradient.
+   */
+  inline virtual su2double GetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim) { return 0; }
+
+  /*!
+   * \brief Set the value of the primitive gradient for MUSCL reconstruction.
+   * \param[in] val_var - Index of the variable.
+   * \param[in] val_dim - Index of the dimension.
+   * \param[in] val_value - Value of the gradient.
+   */
+  inline virtual void SetGradient_Reconstruction(unsigned short val_var, unsigned short val_dim, su2double val_value) { }
+
+  /*!
+   * \brief Get the value of the primitive gradient for MUSCL reconstruction.
+   * \return Value of the primitive gradient for MUSCL reconstruction.
+   */
+  inline virtual su2double **GetGradient_Reconstruction(void) { return NULL; }
+
   /*!
    * \brief Set the blending function for the blending of k-w and k-eps.
    * \param[in] val_viscosity - Value of the vicosity.

SU2_CFD/src/drivers/CDriver.cpp

@@ -3944,14 +3944,6 @@ bool CDriver::Monitor(unsigned long ExtIter) {
     }
   }
 
-  /*--- Evaluate the new CFL number (adaptive). ---*/
-  if (config_container[ZONE_0]->GetCFL_Adapt() == YES) {
-    for (iZone = 0; iZone < nZone; iZone++){
-      if (!(config_container[iZone]->GetMultizone_Problem())) // This needs to be changed everywhere in the code, in a future PR
-        output->SetCFL_Number(solver_container, config_container, iZone);
-    }
-  }
-
   /*--- Check whether the current simulation has reached the specified
    convergence criteria, and set StopCalc to true, if so. ---*/
 
@@ -4276,7 +4268,6 @@ void CTurbomachineryDriver::SetTurboPerformance(unsigned short targetZone){
 
 bool CTurbomachineryDriver::Monitor(unsigned long ExtIter) {
 
-  su2double CFL;
   su2double rot_z_ini, rot_z_final ,rot_z;
   su2double outPres_ini, outPres_final, outPres;
   unsigned long rampFreq, finalRamp_Iter;
@@ -4312,26 +4303,6 @@ bool CTurbomachineryDriver::Monitor(unsigned long ExtIter) {
           config_container, integration_container, false, UsedTime, iZone, iInst);
   }
 
-
-  /*--- Evaluate the new CFL number (adaptive). ---*/
-  if (config_container[ZONE_0]->GetCFL_Adapt() == YES) {
-    if(mixingplane){
-      CFL = 0;
-      for (iZone = 0; iZone < nZone; iZone++){
-        output->SetCFL_Number(solver_container, config_container, iZone);
-        CFL += config_container[iZone]->GetCFL(MESH_0);
-      }
-      /*--- For fluid-multizone the new CFL number is the same for all the zones and it is equal to the zones' minimum value. ---*/
-      for (iZone = 0; iZone < nZone; iZone++){
-        config_container[iZone]->SetCFL(MESH_0, CFL/nZone);
-      }
-    }
-    else{
-      output->SetCFL_Number(solver_container, config_container, ZONE_0);
-    }
-  }
-
-
   /*--- ROTATING FRAME Ramp: Compute the updated rotational velocity. ---*/
   if (config_container[ZONE_0]->GetGrid_Movement() && config_container[ZONE_0]->GetRampRotatingFrame()) {
     rampFreq       = SU2_TYPE::Int(config_container[ZONE_0]->GetRampRotatingFrame_Coeff(1));