Merge pull request #1036 from su2code/fix_cfl_adapt

Update adaptive CFL

Author: pcarruscag

Common/include/CConfig.hpp

@@ -735,7 +735,8 @@ class CConfig {
   unsigned short Tab_FileFormat;      /*!< \brief Format of the output files. */
   unsigned short ActDisk_Jump;        /*!< \brief Format of the output files. */
   unsigned long StartWindowIteration; /*!< \brief Starting Iteration for long time Windowing apporach . */
-  bool CFL_Adapt;        /*!< \brief Adaptive CFL number. */
+  unsigned short nCFL_AdaptParam;     /*!< \brief Number of CFL parameters provided in config. */
+  bool CFL_Adapt;        /*!< \brief Use adaptive CFL number. */
   bool HB_Precondition;  /*!< \brief Flag to turn on harmonic balance source term preconditioning */
   su2double RefArea,     /*!< \brief Reference area for coefficient computation. */
   RefElemLength,         /*!< \brief Reference element length for computing the slope limiting epsilon. */
@@ -758,12 +759,10 @@ class CConfig {
   su2double* Mesh_Box_Length;    /*!< \brief Array containing the length in the x-, y-, and z-directions for the analytic RECTANGLE and BOX grid formats. */
   su2double* Mesh_Box_Offset;    /*!< \brief Array containing the offset from 0.0 in the x-, y-, and z-directions for the analytic RECTANGLE and BOX grid formats. */
   string Mesh_FileName,          /*!< \brief Mesh input file. */
-  Target_Mesh_FileName,          /*!< \brief Mesh name to be adapted. */
   Mesh_Out_FileName,             /*!< \brief Mesh output file. */
   Solution_FileName,             /*!< \brief Flow solution input file. */
   Solution_AdjFileName,          /*!< \brief Adjoint solution input file for drag functional. */
   Volume_FileName,               /*!< \brief Flow variables output file. */
-  Residual_FileName,             /*!< \brief Residual variables output file. */
   Conv_FileName,                 /*!< \brief Convergence history output file. */
   Breakdown_FileName,            /*!< \brief Breakdown output file. */
   Restart_FileName,              /*!< \brief Restart file for flow variables. */
@@ -1057,10 +1056,10 @@ class CConfig {
   bool Radiation;                      /*!< \brief Determines if a radiation model is incorporated. */
   su2double CFL_Rad;                   /*!< \brief CFL Number for the radiation solver. */
 
+  array<su2double,5> default_cfl_adapt;  /*!< \brief Default CFL adapt param array for the COption class. */
   su2double default_vel_inf[3],  /*!< \brief Default freestream velocity array for the COption class. */
   default_eng_cyl[7],            /*!< \brief Default engine box array for the COption class. */
   default_eng_val[5],            /*!< \brief Default engine box array values for the COption class. */
-  default_cfl_adapt[4],          /*!< \brief Default CFL adapt param array for the COption class. */
   default_jst_coeff[2],          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
   default_ffd_coeff[3],          /*!< \brief Default artificial dissipation (flow) array for the COption class. */
   default_mixedout_coeff[3],     /*!< \brief Default default mixedout algorithm coefficients for the COption class. */
@@ -1157,7 +1156,7 @@ class CConfig {
   monoatomic;                               /*!< \brief Flag for monoatomic mixture. */
   string GasModel,                          /*!< \brief Gas Model. */
   *Wall_Catalytic;                          /*!< \brief Pointer to catalytic walls. */
-  
+
   /*!
    * \brief Set the default values of config options not set in the config file using another config object.
    * \param config - Config object to use the default values from.
@@ -1539,13 +1538,6 @@ class CConfig {
    */
   su2double GetCFL_AdaptParam(unsigned short val_index) const { return CFL_AdaptParam[val_index]; }
 
-  /*!
-   * \brief Set the values of the CFL adaption parameters.
-   * \param[in] val_index     - Index of the particular CFL adaption parameter
-   * \param[in] val_cfl_param - Value of the CFL adaption parameter
-   */
-  void SetCFL_AdaptParam(unsigned short val_index, su2double val_cfl_param) { CFL_AdaptParam[val_index] = val_cfl_param; }
-
   /*!
    * \brief Get the value of the CFL adaption flag.
    * \return <code>TRUE</code> if CFL adaption is active; otherwise <code>FALSE</code>.
@@ -5112,11 +5104,6 @@ class CConfig {
    */
   bool GetRestart(void) const { return Restart; }
 
-  /*!
-   * \brief Sets the restart information.
-   */
-  void SetRestart(bool val_restart) { Restart = val_restart; }
-
   /*!
    * \brief Flag for whether binary SU2 native restart files are read.
    * \return Flag for whether binary SU2 native restart files are read, if <code>TRUE</code> then the code will load binary restart files.
@@ -5214,12 +5201,6 @@ class CConfig {
    */
   string GetSolution_AdjFileName(void) const { return Solution_AdjFileName; }
 
-  /*!
-   * \brief Get the name of the file with the residual of the problem.
-   * \return Name of the file with the residual of the problem.
-   */
-  string GetResidual_FileName(void);
-
   /*!
    * \brief Get the format of the input/output grid.
    * \return Format of the input/output grid.

Common/include/option_structure.hpp

@@ -94,7 +94,7 @@ const unsigned int INST_0 = 0;  /*!< \brief Definition of the first instance per
 
 const su2double STANDARD_GRAVITY = 9.80665;           /*!< \brief Acceleration due to gravity at surface of earth. */
 const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598;   /*!< \brief Universal gas constant in J/(mol*K) */
-const su2double BOLTZMANN_CONSTANT = 1.3806503E-23; /*! \brief Boltzmann's constant [J K^-1] */
+const su2double BOLTZMANN_CONSTANT = 1.3806503E-23;   /*! \brief Boltzmann's constant [J K^-1] */
 const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogardro's constant, number of particles in one kmole. */
 
 const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
@@ -541,9 +541,9 @@ enum ENUM_FLUIDMODEL {
   PR_GAS = 3,             /*!< \brief Perfect Real gas model. */
   CONSTANT_DENSITY = 4,   /*!< \brief Constant density gas model. */
   INC_IDEAL_GAS = 5,      /*!< \brief Incompressible ideal gas model. */
-  INC_IDEAL_GAS_POLY = 6,  /*!< \brief Inc. ideal gas, polynomial gas model. */
+  INC_IDEAL_GAS_POLY = 6, /*!< \brief Inc. ideal gas, polynomial gas model. */
   MUTATIONPP = 7,         /*!< \brief Mutation++ gas model for nonequilibrium flow. */
-  SU2_NONEQ = 8        /*!< \brief User defined gas model for nonequilibrium flow. */
+  SU2_NONEQ = 8           /*!< \brief User defined gas model for nonequilibrium flow. */
 };
 static const MapType<string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)

Common/src/CConfig.cpp

@@ -1618,11 +1618,12 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Activate The adaptive CFL number. */
   addBoolOption("CFL_ADAPT", CFL_Adapt, false);
   /* !\brief CFL_ADAPT_PARAM
-   * DESCRIPTION: Parameters of the adaptive CFL number (factor down, factor up, CFL limit (min and max) )
+   * DESCRIPTION: Parameters of the adaptive CFL number (factor down, factor up, CFL limit (min and max), acceptable linear residual )
    * Factor down generally <1.0, factor up generally > 1.0 to cause the CFL to increase when the under-relaxation parameter is 1.0
    * and to decrease when the under-relaxation parameter is less than 0.1. Factor is multiplicative. \ingroup Config*/
   default_cfl_adapt[0] = 1.0; default_cfl_adapt[1] = 1.0; default_cfl_adapt[2] = 10.0; default_cfl_adapt[3] = 100.0;
-  addDoubleArrayOption("CFL_ADAPT_PARAM", 4, CFL_AdaptParam, default_cfl_adapt);
+  default_cfl_adapt[4] = 0.001;
+  addDoubleListOption("CFL_ADAPT_PARAM", nCFL_AdaptParam, CFL_AdaptParam);
   /* DESCRIPTION: Reduction factor of the CFL coefficient in the adjoint problem */
   addDoubleOption("CFL_REDUCTION_ADJFLOW", CFLRedCoeff_AdjFlow, 0.8);
   /* DESCRIPTION: Reduction factor of the CFL coefficient in the level set problem */
@@ -4143,6 +4144,19 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
   CFL = new su2double[nCFL];
   CFL[0] = CFLFineGrid;
 
+  /*--- Handle optional CFL adapt parameter values ---*/
+
+  if (nCFL_AdaptParam < default_cfl_adapt.size()) {
+    auto newParam = new su2double [default_cfl_adapt.size()];
+    for (iCFL = 0; iCFL < default_cfl_adapt.size(); ++iCFL) {
+      if (iCFL < nCFL_AdaptParam) newParam[iCFL] = CFL_AdaptParam[iCFL];
+      else newParam[iCFL] = default_cfl_adapt[iCFL];
+    }
+    swap(newParam, CFL_AdaptParam);
+    delete [] newParam;
+    nCFL_AdaptParam = default_cfl_adapt.size();
+  }
+
   /*--- Evaluate when the Cl should be evaluated ---*/
 
   Iter_Fixed_CM        = SU2_TYPE::Int(nInnerIter / (su2double(Update_iH)+1));
@@ -6540,7 +6554,8 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) {
 
       if (!CFL_Adapt) cout << "No CFL adaptation." << endl;
       else cout << "CFL adaptation. Factor down: "<< CFL_AdaptParam[0] <<", factor up: "<< CFL_AdaptParam[1]
-        <<",\n                lower limit: "<< CFL_AdaptParam[2] <<", upper limit: " << CFL_AdaptParam[3] <<"."<< endl;
+        <<",\n                lower limit: "<< CFL_AdaptParam[2] <<", upper limit: " << CFL_AdaptParam[3]
+        <<",\n                acceptable linear residual: "<< CFL_AdaptParam[4] << "." << endl;
 
       if (nMGLevels !=0) {
         PrintingToolbox::CTablePrinter MGTable(&std::cout);

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2026,13 +2026,11 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
+  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0, dTn, dTven,
             RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
-            thermal_conductivity, thermal_conductivity_tr, thermal_conductivity_ve = 0.0,
-            MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
+            thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
-  su2double **Grad_PrimVar = nullptr, dTn, dTven;
 
   string Marker_Tag, Monitoring_Tag;
 
@@ -2170,12 +2168,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       }
       Density = nodes->GetDensity(iPoint);
 
-      if (nemo) {
-        thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
-        thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
-        Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
-      }
-
       Area = GeometryToolbox::Norm(nDim, Normal);
       for (iDim = 0; iDim < nDim; iDim++) {
         UnitNormal[iDim] = Normal[iDim] / Area;
@@ -2247,6 +2239,10 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       } else {
 
+        const auto thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
+        const auto thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
+        const auto Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
+
         dTn = 0.0; dTven = 0.0;
         for (iDim = 0; iDim < nDim; iDim++) {
           dTn   += Grad_PrimVar[T_INDEX][iDim]*UnitNormal[iDim];

SU2_CFD/include/solvers/CSolver.hpp

@@ -69,8 +69,6 @@ class CSolver {
   unsigned short MGLevel;        /*!< \brief Multigrid level of this solver object. */
   unsigned short IterLinSolver;  /*!< \brief Linear solver iterations. */
   su2double ResLinSolver;        /*!< \brief Final linear solver residual. */
-  su2double NonLinRes_Value,     /*!< \brief Summed value of the nonlinear residual indicator. */
-  NonLinRes_Func;                /*!< \brief Current value of the nonlinear residual indicator at one iteration. */
   unsigned short NonLinRes_Counter;   /*!< \brief Number of elements of the nonlinear residual indicator series. */
   vector<su2double> NonLinRes_Series; /*!< \brief Vector holding the nonlinear residual indicator series. */
   su2double Old_Func,  /*!< \brief Old value of the nonlinear residual indicator. */

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -58,10 +58,11 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe
                               CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) {
   unsigned long InnerIter, TimeIter;
 
-  bool unsteady = (config[val_iZone]->GetTime_Marching() == DT_STEPPING_1ST) ||
-                  (config[val_iZone]->GetTime_Marching() == DT_STEPPING_2ND);
-  bool frozen_visc = (config[val_iZone]->GetContinuous_Adjoint() && config[val_iZone]->GetFrozen_Visc_Cont()) ||
-                     (config[val_iZone]->GetDiscrete_Adjoint() && config[val_iZone]->GetFrozen_Visc_Disc());
+  const bool unsteady = (config[val_iZone]->GetTime_Marching() == DT_STEPPING_1ST) ||
+                        (config[val_iZone]->GetTime_Marching() == DT_STEPPING_2ND);
+  const bool frozen_visc = (config[val_iZone]->GetContinuous_Adjoint() && config[val_iZone]->GetFrozen_Visc_Cont()) ||
+                           (config[val_iZone]->GetDiscrete_Adjoint() && config[val_iZone]->GetFrozen_Visc_Disc());
+  const bool disc_adj = (config[val_iZone]->GetDiscrete_Adjoint());
   TimeIter = config[val_iZone]->GetTimeIter();
 
   /* --- Setting up iteration values depending on if this is a
@@ -134,7 +135,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe
 
   /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/
 
-  if (config[val_iZone]->GetCFL_Adapt() == YES) {
+  if ((config[val_iZone]->GetCFL_Adapt() == YES) && (!disc_adj)) {
     SU2_OMP_PARALLEL
     solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst],
                                                                    solver[val_iZone][val_iInst], config[val_iZone]);

SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp

@@ -102,7 +102,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeChemistry(const CConfig *config)
 
   ws = fluidmodel->ComputeNetProductionRates();
 
-  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) 
+  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     residual[iSpecies] = ws[iSpecies] *Volume;
 
   //if (implicit) {
@@ -230,7 +230,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeChemistry(const CConfig *config)
 //    } // implicit
 
   return ResidualType<>(residual, nullptr, nullptr);
- 
+
 }
 
 CNumerics::ResidualType<> CSource_NEMO::ComputeVibRelaxation(const CConfig *config) {
@@ -286,19 +286,19 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeVibRelaxation(const CConfig *conf
 //  }
 
   if(config->GetVTTransferResidualLimiting()){
-    if(residual[nSpecies+nDim+1]>res_max) residual[nSpecies+nDim+1]=res_max; 
+    if(residual[nSpecies+nDim+1]>res_max) residual[nSpecies+nDim+1]=res_max;
     if(residual[nSpecies+nDim+1]<res_min) residual[nSpecies+nDim+1]=res_min;
-  } 
+  }
 
   return ResidualType<>(residual, nullptr, nullptr);
 }
 
 CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *config) {
 
   unsigned short iDim, iSpecies, iVar;
-  su2double rho, rhov, vel2, H, yinv, T, Tve, qy_ve, Ru, RuSI;
+  su2double rho, rhov, vel2, H, yinv, T, Tve, Ru, RuSI;
   su2double *Ds, **GV, ktr, kve;
-  
+
   /*--- Rename for convenience ---*/
   Ds = Diffusion_Coeff_i;
   ktr = Thermal_Conductivity_i;
@@ -311,7 +311,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
   Ru  = 1000.0*RuSI;
 
   auto& Ms = fluidmodel->GetSpeciesMolarMass();
-  
+
   bool viscous = config->GetViscous();
   bool rans = (config->GetKind_Turb_Model() != NONE);
   hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, eve_i);
@@ -352,23 +352,23 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
 
     for (iSpecies=0; iSpecies<nSpecies; iSpecies++)
       Mass += V_i[iSpecies]*Ms[iSpecies];
-        
+
     su2double heat_capacity_cp_i   = V_i[RHOCVTR_INDEX]+Ru/Mass;
     su2double total_viscosity_i    = Laminar_Viscosity_i + Eddy_Viscosity_i;
-    su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;      
+    su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;
     su2double u                    = V_i[VEL_INDEX];
     su2double v                    = V_i[VEL_INDEX+1];
     su2double qy_ve                = kve*GV[TVE_INDEX][1];
 
     /*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
-    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { 
+    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
       sumJhs_y  += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
       sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
     }
 
     for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
       residual[iSpecies] -= 0.0;
-    residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0]) 
+    residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0])
                                                                     -TWO3*AuxVar_Grad_i[0][0]);
     residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
                                                                     -TWO3*AuxVar_Grad_i[0][1]);

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -3854,12 +3854,12 @@ void CEulerSolver::GetPower_Properties(CGeometry *geometry, CConfig *config, uns
               MassFlow -= Vector[iDim]*Velocity[iDim]*Density;
             }
 
+            Area = GeometryToolbox::Norm(nDim, Vector);
             Vn = 0.0; ReverseFlow = false;
             for (iDim = 0; iDim < nDim; iDim++) {  Vn -= Velocity[iDim]*Vector[iDim]/Area; }
             if (Vn < 0.0) { ReverseFlow = true; }
 
             Vel_Infty = sqrt (Vel_Infty2);
-            Area = GeometryToolbox::Norm(nDim, Vector);
             Mach = sqrt(Velocity2)/SoundSpeed;
             TotalPressure = Pressure * pow( 1.0 + Mach * Mach * 0.5 * (Gamma - 1.0), Gamma    / (Gamma - 1.0));
             TotalTemperature = Temperature * (1.0 + Mach * Mach * 0.5 * (Gamma - 1.0));