Merge pull request #1556 from su2code/fix_newton_krylov_dual_time

Fix Newton-Krylov for unsteady problems

Author: pcarruscag

Common/include/CConfig.hpp

@@ -148,7 +148,6 @@ class CConfig {
   su2double CL_Target;         /*!< \brief Fixed Cl mode Target Cl. */
   TIME_MARCHING TimeMarching;        /*!< \brief Steady or unsteady (time stepping or dual time stepping) computation. */
   unsigned short Dynamic_Analysis;   /*!< \brief Static or dynamic structural analysis. */
-  unsigned short nStartUpIter;       /*!< \brief Start up iterations using the fine grid. */
   su2double FixAzimuthalLine;        /*!< \brief Fix an azimuthal line due to misalignments of the nearfield. */
   su2double **DV_Value;              /*!< \brief Previous value of the design variable. */
   su2double Venkat_LimiterCoeff;     /*!< \brief Limiter coefficient */
@@ -1977,12 +1976,6 @@ class CConfig {
    */
   su2double GetLength_Reynolds(void) const { return Length_Reynolds; }
 
-  /*!
-   * \brief Get the start up iterations using the fine grid, this works only for multigrid problems.
-   * \return Start up iterations using the fine grid.
-   */
-  unsigned short GetnStartUpIter(void) const { return nStartUpIter; }
-
   /*!
    * \brief Get the reference area for non dimensional coefficient computation. If the value from the
    *        is 0 then, the code will compute the reference area using the projection of the shape into

Common/src/CConfig.cpp

@@ -1810,11 +1810,8 @@ void CConfig::SetConfig_Options() {
   /*!\brief CONV_WINDOW_FIELD
    * \n DESCRIPTION: Output fields  for the Cauchy criterium for the TIME iteration. The criterium is applied to the windowed time average of the chosen funcion. */
   addStringListOption("CONV_WINDOW_FIELD",nWndConvField, WndConvField);
-  /*!\par CONFIG_CATEGORY: Multi-grid \ingroup Config*/
-  /*--- Options related to Multi-grid ---*/
 
-  /*!\brief START_UP_ITER \n DESCRIPTION: Start up iterations using the fine grid only. DEFAULT: 0 \ingroup Config*/
-  addUnsignedShortOption("START_UP_ITER", nStartUpIter, 0);
+  /*!\par CONFIG_CATEGORY: Multi-grid \ingroup Config*/
   /*!\brief MGLEVEL\n DESCRIPTION: Multi-grid Levels. DEFAULT: 0 \ingroup Config*/
   addUnsignedShortOption("MGLEVEL", nMGLevels, 0);
   /*!\brief MGCYCLE\n DESCRIPTION: Multi-grid cycle. OPTIONS: See \link MG_Cycle_Map \endlink. Defualt V_CYCLE \ingroup Config*/
@@ -5731,7 +5728,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) {
         break;
       }
     }
-    
+
     if(!found) {
       if (nZone==1)
         SU2_MPI::Error("DV_MARKER contains marker names that do not exist in the lists of BCs in the config file.", CURRENT_FUNCTION);
@@ -6818,7 +6815,6 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
 
     if (nMGLevels !=0) {
 
-      if (nStartUpIter != 0) cout << "A total of " << nStartUpIter << " start up iterations on the fine grid."<< endl;
       if (MGCycle == V_CYCLE) cout << "V Multigrid Cycle, with " << nMGLevels << " multigrid levels."<< endl;
       if (MGCycle == W_CYCLE) cout << "W Multigrid Cycle, with " << nMGLevels << " multigrid levels."<< endl;
       if (MGCycle == FULLMG_CYCLE) cout << "Full Multigrid Cycle, with " << nMGLevels << " multigrid levels."<< endl;

SU2_CFD/include/integration/CNewtonIntegration.hpp

@@ -62,7 +62,7 @@ class CNewtonIntegration final : public CIntegration {
 
 private:
   /*--- Residual evaluation modes, explicit for products, default to allow preconditioners to be built. ---*/
-  enum ResEvalType {EXPLICIT, DEFAULT};
+  enum class ResEvalType {EXPLICIT, DEFAULT};
 
   bool setup = false;
   Scalar finDiffStepND = 0.0;

SU2_CFD/src/integration/CNewtonIntegration.cpp

@@ -120,7 +120,7 @@ void CNewtonIntegration::ComputeResiduals(ResEvalType type) {
 
   /*--- Save the default integration scheme, and force to explicit if required. ---*/
   auto TimeIntScheme = config->GetKind_TimeIntScheme();
-  if (type == EXPLICIT) {
+  if (type == ResEvalType::EXPLICIT) {
     SU2_OMP_MASTER
     config->SetKind_TimeIntScheme(EULER_EXPLICIT);
     END_SU2_OMP_MASTER
@@ -129,10 +129,14 @@ void CNewtonIntegration::ComputeResiduals(ResEvalType type) {
 
   solvers[FLOW_SOL]->Preprocessing(geometry, solvers, config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, false);
 
+  if (type == ResEvalType::DEFAULT) {
+    solvers[FLOW_SOL]->SetTime_Step(geometry, solvers, config, MESH_0, config->GetTimeIter());
+  }
+
   Space_Integration(geometry, solvers, numerics[FLOW_SOL], config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS);
 
   /*--- Restore default. ---*/
-  if (type == EXPLICIT) {
+  if (type == ResEvalType::EXPLICIT) {
     SU2_OMP_MASTER
     config->SetKind_TimeIntScheme(TimeIntScheme);
     END_SU2_OMP_MASTER
@@ -187,12 +191,10 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
 
   /*--- Current residual. ---*/
 
-  ComputeResiduals(DEFAULT);
+  ComputeResiduals(ResEvalType::DEFAULT);
 
   /*--- Compute the approximate Jacobian for preconditioning. ---*/
 
-  solvers[FLOW_SOL]->SetTime_Step(geometry, solvers, config, MESH_0, config->GetTimeIter());
-
   solvers[FLOW_SOL]->PrepareImplicitIteration(geometry, solvers, config);
 
   if (preconditioner) preconditioner->Build();
@@ -304,7 +306,7 @@ void CNewtonIntegration::MatrixFreeProduct(const CSysVector<Scalar>& u, CSysVect
 
   PerturbSolution(u, factor);
 
-  ComputeResiduals(EXPLICIT);
+  ComputeResiduals(ResEvalType::EXPLICIT);
 
   /*--- Finalize product. ---*/
   factor = 1.0 / factor;
@@ -339,9 +341,13 @@ void CNewtonIntegration::Preconditioner(const CSysVector<Scalar>& u, CSysVector<
   else {
     /*--- Approximate diagonal preconditioner. ---*/
 
+    const bool dt1st = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
+    const bool dt2nd = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
+    const su2double dt = config->GetDelta_UnstTimeND() * (dt1st + 1.5 * dt2nd);
+
     SU2_OMP_FOR_STAT(omp_chunk_size)
     for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); ++iPoint) {
-      su2double delta = solvers[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint) /
+      su2double delta = (solvers[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint) + dt) /
                         (geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint));
       SU2_OMP_SIMD
       for (auto iVar = 0ul; iVar < u.GetNVar(); ++iVar)