Newton Krylov improvements

Common/include/CConfig.hpp

@@ -431,7 +431,8 @@ class CConfig {
   bool UseVectorization;       /*!< \brief Whether to use vectorized numerics schemes. */
   bool NewtonKrylov;           /*!< \brief Use a coupled Newton method to solve the flow equations. */
   array<unsigned short,3> NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */
-  array<su2double,4> NK_DblParam{{-2.0, 0.1, -3.0, 1e-4}}; /*!< \brief Floating-point parameters for NK method. */
+  array<su2double,5> NK_DblParam{{-2.0, 0.1, -3.0, 1e-4, 1.0}}; /*!< \brief Floating-point parameters for NK method. */
+  su2double NK_Relaxation = 1.0;
 
   unsigned short nMGLevels;    /*!< \brief Number of multigrid levels (coarse levels). */
   unsigned short nCFL;         /*!< \brief Number of CFL, one for each multigrid level. */
@@ -1332,9 +1333,9 @@ class CConfig {
   template <class Tenum, class Tfield>
   void addEnumListOption(const string name, unsigned short& input_size, Tfield*& option_field, const map<string,Tenum>& enum_map);
 
-  void addDoubleArrayOption(const string& name, const int size, su2double* option_field);
+  void addDoubleArrayOption(const string& name, int size, bool allow_fewer, su2double* option_field);
 
-  void addUShortArrayOption(const string& name, const int size, unsigned short* option_field);
+  void addUShortArrayOption(const string& name, int size, bool allow_fewer, unsigned short* option_field);
 
   void addDoubleListOption(const string& name, unsigned short & size, su2double * & option_field);
 
@@ -4375,12 +4376,22 @@ class CConfig {
   /*!
    * \brief Get Newton-Krylov integer parameters.
    */
-  array<unsigned short,3> GetNewtonKrylovIntParam(void) const { return NK_IntParam; }
+  array<unsigned short,3> GetNewtonKrylovIntParam() const { return NK_IntParam; }
 
   /*!
    * \brief Get Newton-Krylov floating-point parameters.
    */
-  array<su2double,4> GetNewtonKrylovDblParam(void) const { return NK_DblParam; }
+  array<su2double,5> GetNewtonKrylovDblParam() const { return NK_DblParam; }
+
+  /*!
+   * \brief Get the Newton-Krylov relaxation.
+   */
+  su2double GetNewtonKrylovRelaxation() const { return NK_Relaxation; }
+
+  /*!
+   * \brief Set the Newton-Krylov relaxation.
+   */
+  void SetNewtonKrylovRelaxation(const su2double& relaxation) { NK_Relaxation = relaxation; }
 
   /*!
    * \brief Returns the Roe kappa (multipler of the dissipation term).

Common/include/option_structure.inl

@@ -233,18 +233,19 @@ class COptionEnumList final : public COptionBase {
 
 template <class Type>
 class COptionArray final : public COptionBase {
-  string name;     // Identifier for the option
-  const int size;  // Number of elements
-  Type* field;     // Reference to the field
+  string name;             // Identifier for the option
+  const int size;          // Number of elements
+  const bool allow_fewer;  // Allow smaller size
+  Type* field;             // Reference to the field
 
  public:
-  COptionArray(string option_field_name, const int list_size, Type* option_field)
-      : name(option_field_name), size(list_size), field(option_field) {}
+  COptionArray(string option_field_name, const int list_size, const bool allow_fewer, Type* option_field)
+      : name(std::move(option_field_name)), size(list_size), allow_fewer(allow_fewer), field(option_field) {}
 
   string SetValue(const vector<string>& option_value) override {
     COptionBase::SetValue(option_value);
     // Check that the size is correct
-    if (option_value.size() != (unsigned long)this->size) {
+    if ((option_value.size() < size_t(size) && !allow_fewer) || option_value.size() > size_t(size)) {
       string newstring;
       newstring.append(this->name);
       newstring.append(": wrong number of arguments: ");
@@ -258,7 +259,7 @@ class COptionArray final : public COptionBase {
       newstring.append(" found");
       return newstring;
     }
-    for (int i = 0; i < this->size; i++) {
+    for (size_t i = 0; i < option_value.size(); i++) {
       istringstream is(option_value[i]);
       if (!(is >> field[i])) {
         return badValue(" array", this->name);

Common/src/linear_algebra/CSysSolve.cpp

@@ -47,7 +47,7 @@ constexpr T linSolEpsilon() {
 }
 template <>
 constexpr float linSolEpsilon<float>() {
-  return 1e-12;
+  return 1e-14;
 }
 }  // namespace
 

SU2_CFD/include/integration/CNewtonIntegration.hpp

@@ -67,8 +67,10 @@ class CNewtonIntegration final : public CIntegration {
   enum class ResEvalType {EXPLICIT, DEFAULT};
 
   bool setup = false;
+  bool autoRelaxation = false;
   Scalar finDiffStepND = 0.0;
   Scalar finDiffStep = 0.0; /*!< \brief Based on RMS(solution), used in matrix-free products. */
+  Scalar nkRelaxation = 1.0;
   unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */
 
   /*--- Number of iterations and tolerance for the linear preconditioner,
@@ -81,7 +83,7 @@ class CNewtonIntegration final : public CIntegration {
    * criteria are zero, or the solver does not provide a linear
    * preconditioner, there is no startup phase. ---*/
   bool startupPeriod = false;
-  unsigned short startupIters = 0;
+  unsigned short startupIters = 0, iter = 0;
   su2double startupResidual = 0.0;
   su2double firstResidual = -20.0;
 
@@ -96,8 +98,9 @@ class CNewtonIntegration final : public CIntegration {
   CNumerics*** numerics = nullptr;
 
   /*--- Residual and linear solver. ---*/
-  CSysVector<Scalar> LinSysRes;
+  CSysVector<Scalar> LinSysRes, LinSysResRelax;
   CSysSolve<Scalar> LinSolver;
+  const CSysVector<Scalar>* LinSysRes0 = nullptr;
 
   /*--- If possible the solution vector of the solver is re-used, otherwise this temporary is used. ---*/
   CSysVector<Scalar> LinSysSol;

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -114,7 +114,7 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
                                                       const VariableType& solution) {
-  static_assert(ReconVarType::nVar <= PrimVarType::nVar,"");
+  static_assert(ReconVarType::nVar <= PrimVarType::nVar);
 
   const auto& gradients = solution.GetGradient_Reconstruction();
   const auto& limiters = solution.GetLimiter_Primitive();
@@ -129,7 +129,6 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
   if (muscl) {
     /*--- Recompute density and enthalpy instead of reconstructing. ---*/
     constexpr auto nVarGrad = ReconVarType::nVar - 2;
-
     switch (limiterType) {
     case LIMITER::NONE:
       musclUnlimited<nVarGrad>(iPoint, vector_ij, 0.5, gradients, V.i.all);

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -96,6 +96,7 @@ class CRoeBase : public Base {
     AD::StartPreacc();
 
     const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT);
+    const auto nkRelax = config.GetNewtonKrylovRelaxation();
     const auto& solution = static_cast<const CEulerVariable&>(solution_);
 
     const auto iPoint = geometry.edges->GetNode(iEdge,0);
@@ -118,8 +119,13 @@ class CRoeBase : public Base {
     V1st.i.all = gatherVariables<nPrimVar>(iPoint, solution.GetPrimitive());
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
+    VectorDbl<nDim> mod_vector_ij;
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      mod_vector_ij(iDim) = nkRelax * vector_ij(iDim);
+    }
+    /*--- Recompute density and enthalpy instead of reconstructing. ---*/
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution);
+        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, mod_vector_ij, solution);
 
     /*--- Compute conservative variables. ---*/
 

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -30,6 +30,11 @@
 #include <limits>
 #include "CFlowVariable.hpp"
 
+/*!
+ * \brief Returns the number of primitive variables for which to compute gradients.
+ */
+unsigned long EulerNPrimVarGrad(const CConfig *config, unsigned long ndim);
+
 /*!
  * \class CEulerVariable
  * \brief Class for defining the variables of the compressible Euler solver.

SU2_CFD/src/integration/CNewtonIntegration.cpp

@@ -66,13 +66,18 @@ void CNewtonIntegration::Setup() {
   auto iparam = config->GetNewtonKrylovIntParam();
   auto dparam = config->GetNewtonKrylovDblParam();
 
-  startupIters = iparam[0];
+  startupIters = iter = iparam[0];
   startupResidual = dparam[0];
   precondIters = iparam[1];
   precondTol = dparam[1];
   tolRelaxFactor = iparam[2];
   fullTolResidual = dparam[2];
   finDiffStepND = SU2_TYPE::GetValue(dparam[3]);
+  nkRelaxation = fmin(SU2_TYPE::GetValue(dparam[4]), 1);
+  if (nkRelaxation < 0) {
+    autoRelaxation = true;
+    nkRelaxation = 1;
+  }
 
   const auto nVar = solvers[FLOW_SOL]->GetnVar();
   const auto nPoint = geometry->GetnPoint();
@@ -83,6 +88,9 @@ void CNewtonIntegration::Setup() {
   LinSolver.SetxIsZero(true);
 
   LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0);
+  if (autoRelaxation || nkRelaxation < 1) {
+    LinSysResRelax.Initialize(nPoint, nPointDomain, nVar, 0.0);
+  }
 
   if (!std::is_same<Scalar,su2double>::value) {
     LinSysSol.Initialize(nPoint, nPointDomain, nVar, nullptr);
@@ -175,6 +183,17 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
 
   if (!setup) { Setup(); setup = true; }
 
+  // Ramp from 1st to 2nd order during the startup.
+  su2double baseNkRelaxation = 1;
+  if (startupPeriod && startupIters > 0 && !config->GetRestart()) {
+    baseNkRelaxation = su2double(startupIters - iter) / startupIters;
+  }
+  config->SetNewtonKrylovRelaxation(baseNkRelaxation);
+
+  // When using NK relaxation (not fully 2nd order Jacobian products) we need an additional
+  // residual evaluation that is used as the reference for finite differences.
+  LinSysRes0 = (!startupPeriod && nkRelaxation < 1) ? &LinSysResRelax : &LinSysRes;
+
   SU2_OMP_PARALLEL_(if(solvers[FLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Save the current solution to be able to perturb it. ---*/
@@ -191,10 +210,20 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
 
   if (preconditioner) preconditioner->Build();
 
-  SU2_OMP_FOR_STAT(omp_chunk_size)
-  for (auto i = 0ul; i < LinSysRes.GetNElmDomain(); ++i)
-    LinSysRes[i] = SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes[i]);
-  END_SU2_OMP_FOR
+  auto CopyLinSysRes = [&](int sign, auto& dst) {
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (auto i = 0ul; i < dst.GetNElmDomain(); ++i)
+      dst[i] = sign * SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes[i]);
+    END_SU2_OMP_FOR
+  };
+  CopyLinSysRes(1, LinSysRes);
+
+  if (!startupPeriod && nkRelaxation < 1) {
+    SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetNewtonKrylovRelaxation(nkRelaxation);)
+    ComputeResiduals(ResEvalType::EXPLICIT);
+    // Here the sign was not flipped by PrepareImplicitIteration.
+    CopyLinSysRes(-1, LinSysResRelax);
+  }
 
   su2double residual = 0.0;
   for (auto iVar = 0ul; iVar < LinSysRes.GetNVar(); ++iVar)
@@ -207,10 +236,10 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
   if (startupPeriod) {
     BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       firstResidual = max(firstResidual, residual);
-      if (startupIters) startupIters -= 1;
+      if (iter) iter -= 1;
     }
     END_SU2_OMP_SAFE_GLOBAL_ACCESS
-    endStartup = (startupIters == 0) && (residual - firstResidual < startupResidual);
+    endStartup = (iter == 0) && (residual - firstResidual < startupResidual);
   }
 
   /*--- The NK solves are expensive, the tolerance is relaxed while the residuals are high. ---*/
@@ -232,8 +261,7 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
 
   if (startupPeriod) {
     iter = Preconditioner_impl(LinSysRes, linSysSol, iter, eps);
-  }
-  else {
+  } else {
     ComputeFinDiffStep();
 
     eps *= toleranceFactor;
@@ -247,6 +275,15 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **
   BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
     solvers[FLOW_SOL]->SetIterLinSolver(iter);
     solvers[FLOW_SOL]->SetResLinSolver(eps);
+
+    if (!startupPeriod && autoRelaxation) {
+      const su2double adaptTol = config->GetCFL_Adapt() ? config->GetCFL_AdaptParam(4) : 0;
+      if (eps > fmax(config->GetLinear_Solver_Error(), adaptTol)) {
+        nkRelaxation *= 0.9;
+      } else if (eps < 0.9 * fmax(config->GetLinear_Solver_Error(), adaptTol)) {
+        nkRelaxation = fmin(nkRelaxation * 1.05, 1);
+      }
+    }
   }
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
@@ -303,11 +340,11 @@ void CNewtonIntegration::MatrixFreeProduct(const CSysVector<Scalar>& u, CSysVect
     su2double delta = (geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint)) /
                       max(EPS, solvers[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint));
     SU2_OMP_SIMD
-    for (auto iVar = 0ul; iVar < LinSysRes.GetNVar(); ++iVar) {
+    for (auto iVar = 0ul; iVar < LinSysRes0->GetNVar(); ++iVar) {
       Scalar perturbRes = SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes(iPoint,iVar));
 
       /*--- The global residual had its sign flipped, so we add to get the difference. ---*/
-      v(iPoint,iVar) = (perturbRes + LinSysRes(iPoint,iVar)) * factor;
+      v(iPoint,iVar) = (perturbRes + (*LinSysRes0)(iPoint,iVar)) * factor;
 
       /*--- Pseudotime term of the true Jacobian. ---*/
       v(iPoint,iVar) += SU2_TYPE::GetValue(delta) * u(iPoint,iVar);

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -277,6 +277,7 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){
   SetVolumeOutputFieldsScalarResidual(config);
 
   if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) {
+    AddVolumeOutput("LIMITER_TEMPERATURE", "Limiter_Temperature", "LIMITER", "Limiter value of the temperature");
     AddVolumeOutput("LIMITER_VELOCITY-X", "Limiter_Velocity_x", "LIMITER", "Limiter value of the x-velocity");
     AddVolumeOutput("LIMITER_VELOCITY-Y", "Limiter_Velocity_y", "LIMITER", "Limiter value of the y-velocity");
     if (nDim == 3) {
@@ -364,14 +365,17 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
   }
 
   if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) {
+    SetVolumeOutputValue("LIMITER_TEMPERATURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 0));
     SetVolumeOutputValue("LIMITER_VELOCITY-X", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 1));
     SetVolumeOutputValue("LIMITER_VELOCITY-Y", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 2));
     if (nDim == 3){
       SetVolumeOutputValue("LIMITER_VELOCITY-Z", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 3));
     }
     SetVolumeOutputValue("LIMITER_PRESSURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+1));
-    SetVolumeOutputValue("LIMITER_DENSITY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+2));
-    SetVolumeOutputValue("LIMITER_ENTHALPY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+3));
+    if (solver[FLOW_SOL]->GetnPrimVarGrad() > nDim + 2) {
+      SetVolumeOutputValue("LIMITER_DENSITY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+2));
+      SetVolumeOutputValue("LIMITER_ENTHALPY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+3));
+    }
   }
 
   if (config->GetKind_RoeLowDiss() != NO_ROELOWDISS){

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -64,7 +64,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
                          (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
   const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
   const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint();
-  const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED;
 
   int Unst_RestartIter = 0;
   unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0;
@@ -120,7 +119,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
   nVar = nDim + 2;
   nPrimVar = nDim + 9;
   /*--- Centered schemes only need gradients for viscous fluxes (T and v). ---*/
-  nPrimVarGrad = nDim + (centered && !config->GetContinuous_Adjoint() ? 1 : 4);
+  nPrimVarGrad = EulerNPrimVarGrad(config, nDim);
   nSecondaryVar = nSecVar;
   nSecondaryVarGrad = 2;
 
@@ -1796,6 +1795,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
   }
 
   const bool implicit         = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const su2double nkRelax     = config->GetNewtonKrylovRelaxation();
 
   const bool roe_turkel       = (config->GetKind_Upwind_Flow() == UPWIND::TURKEL);
   const auto kind_dissipation = config->GetKind_RoeLowDiss();
@@ -1875,7 +1875,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
 
       su2double Vector_ij[MAXNDIM] = {0.0};
       for (iDim = 0; iDim < nDim; iDim++) {
-        Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]);
+        Vector_ij[iDim] = nkRelax * 0.5 * (Coord_j[iDim] - Coord_i[iDim]);
       }
 
       auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint);

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1231,6 +1231,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
   const bool limiter    = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE);
   const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE);
   const bool bounded_scalar = config->GetBounded_Scalar();
+  const su2double nkRelax = config->GetNewtonKrylovRelaxation();
 
   /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of
   * variables, otherwise switch to the faster adjoint evaluation mode. ---*/
@@ -1271,7 +1272,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
 
       su2double Vector_ij[MAXNDIM] = {0.0};
       for (iDim = 0; iDim < nDim; iDim++) {
-        Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]);
+        Vector_ij[iDim] = nkRelax * 0.5 * (Coord_j[iDim] - Coord_i[iDim]);
       }
 
       auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint);