Merge branch 'develop' into feature_SST_RoughBCs

Author: pcarruscag

AUTHORS.md

@@ -56,6 +56,7 @@ Aniket C. Aranake
 Antonio Rubino
 Arne Bachmann
 Arne Voß
+Ayush Kumar
 Beckett Y. Zhou
 Benjamin S. Kirk
 Brendan Tracey

Common/include/linear_algebra/CSysSolve.hpp

@@ -295,6 +295,12 @@ class CSysSolve {
                                       ScalarType& residual, bool monitoring, const CConfig* config,
                                       FgcrodrMode mode) const;
 
+  /*!
+   * \brief Creates the inner solver for nested preconditioning if the settings allow it.
+   * \returns True if the inner solver can be used.
+   */
+  bool SetupInnerSolver(unsigned short kind_solver, const CConfig* config);
+
  public:
   /*!
    * \brief default constructor of the class.

Common/include/option_structure.hpp

@@ -2441,10 +2441,12 @@ static const MapType<std::string, ENUM_LINEAR_SOLVER> Linear_Solver_Map = {
 enum class LINEAR_SOLVER_INNER {
   NONE,     /*!< \brief Do not use a nested linear solver. */
   BCGSTAB,  /*!< \brief Use BCGSTAB as the preconditioning linear solver. */
+  SMOOTHER, /*!< \brief Iterative smoother. */
 };
 static const MapType<std::string, LINEAR_SOLVER_INNER> Inner_Linear_Solver_Map = {
   MakePair("NONE", LINEAR_SOLVER_INNER::NONE)
   MakePair("BCGSTAB", LINEAR_SOLVER_INNER::BCGSTAB)
+  MakePair("SMOOTHER", LINEAR_SOLVER_INNER::SMOOTHER)
 };
 
 

Common/src/linear_algebra/CSysSolve.cpp

@@ -312,8 +312,8 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -626,14 +626,15 @@ unsigned long CSysSolve<ScalarType>::RFGMRES_LinSolver(const CSysVector<ScalarTy
   }
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
-  for (auto totalIter = 0ul; totalIter < MaxIter;) {
+  auto totalIter = 0ul;
+  while (totalIter < MaxIter) {
     /*--- Enforce a hard limit on total number of iterations ---*/
     auto iterLimit = min(restartIter, MaxIter - totalIter);
     auto iter = FGMRES_LinSolver(b, x, mat_vec, precond, tol, iterLimit, residual, monitoring, config);
     totalIter += iter;
-    if (residual <= tol || iter < iterLimit) return totalIter;
+    if (residual <= tol || iter < iterLimit) break;
   }
-  return 0;
+  return totalIter;
 }
 
 template <class ScalarType>
@@ -1032,8 +1033,8 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -1205,8 +1206,8 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
   /*--- Only compute the residuals in full communication mode. ---*/
 
   if (config->GetComm_Level() == COMM_FULL) {
-    norm_r = r.norm();
     norm0 = b.norm();
+    norm_r = xIsZero ? norm0 : r.norm();
 
     /*--- Set the norm to the initial initial residual value ---*/
 
@@ -1246,15 +1247,19 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
      current residual, the system is linear so this saves some computation
      compared to re-evaluating r = b-A*x. ---*/
 
-    mat_vec(z, A_x);
+    if (!fix_iter_mode || i != m - 1) {
+      mat_vec(z, A_x);
+    }
 
     /*--- Update solution and residual with relaxation omega. Mathematically this
      is a modified Richardson iteration for the left-preconditioned system
      M^{-1}(b-A*x) which converges if ||I-w*M^{-1}*A|| < 1. Combining this method
      with a Gauss-Seidel preconditioner and w>1 is NOT equivalent to SOR. ---*/
 
     x += omega * z;
-    r -= omega * A_x;
+    if (!fix_iter_mode || i != m - 1) {
+      r -= omega * A_x;
+    }
 
     /*--- Only compute the residuals in full communication mode. ---*/
     /*--- Check if solution has converged, else output the relative residual if necessary. ---*/
@@ -1282,6 +1287,33 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
   return i;
 }
 
+template <class ScalarType>
+bool CSysSolve<ScalarType>::SetupInnerSolver(unsigned short kind_solver, const CConfig* config) {
+  bool flexible = false;
+  switch (kind_solver) {
+    case FGMRES:
+    case FGCRODR:
+    case RESTARTED_FGMRES:
+    case SMOOTHER:
+      flexible = true;
+      break;
+    default:
+      flexible = false;
+  }
+  const bool is_linear = config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::SMOOTHER;
+
+  if (config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE && (flexible || is_linear)) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+    if (!inner_solver) {
+      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
+      inner_solver->SetxIsZero(true);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+    return true;
+  }
+  return false;
+}
+
 template <class ScalarType>
 unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, const CSysVector<su2double>& LinSysRes,
                                            CSysVector<su2double>& LinSysSol, CGeometry* geometry,
@@ -1334,16 +1366,7 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
     }
   }
 
-  const bool nested = (KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES || KindSolver == SMOOTHER) &&
-                      config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE;
-
-  if (nested && !inner_solver) {
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
-      inner_solver->SetxIsZero(true);
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
+  const bool nested = SetupInnerSolver(KindSolver, config);
 
   /*--- Stop the recording for the linear solver ---*/
   bool TapeActive = NO;
@@ -1389,11 +1412,15 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
       auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
         /*--- Initialize to 0 to be safe. ---*/
         v = ScalarType{};
-        ScalarType unused{};
+        ScalarType res{};
         /*--- Handle other types here if desired but do not call Solve because
          * that will create issues with the AD external function. ---*/
-        (void)inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, unused, false,
-                                              config);
+        if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
+          inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
+        } else {
+          const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
+          inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
+        }
       };
       nested_prec = new CAbstractPreconditioner<ScalarType>(f);
     }
@@ -1545,16 +1572,7 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
     }
   }
 
-  const bool nested = (KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES || KindSolver == SMOOTHER) &&
-                      config->GetKind_Linear_Solver_Inner() != LINEAR_SOLVER_INNER::NONE;
-
-  if (nested && !inner_solver) {
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
-      inner_solver = std::make_unique<CSysSolve<ScalarType>>(LINEAR_SOLVER_MODE::STANDARD);
-      inner_solver->SetxIsZero(true);
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
+  const bool nested = SetupInnerSolver(KindSolver, config);
 
   /*--- Set up preconditioner and matrix-vector product ---*/
 
@@ -1572,13 +1590,15 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
   CPreconditioner<ScalarType>* nested_prec = nullptr;
   if (nested) {
     auto f = [&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) {
-      /*--- Initialize to 0 to be safe. ---*/
+      /*--- See "Solve". ---*/
       v = ScalarType{};
-      ScalarType unused{};
-      /*--- Handle other types here if desired but do not call Solve because
-       * that will create issues with the AD external function. ---*/
-      (void)inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, unused, false,
-                                            config);
+      ScalarType res{};
+      if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) {
+        inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config);
+      } else {
+        const auto smooth_iter = static_cast<unsigned long>(std::round(fmax(2, sqrt(MaxIter))));
+        inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config);
+      }
     };
     nested_prec = new CAbstractPreconditioner<ScalarType>(f);
   }

SU2_CFD/include/integration/CNewtonIntegration.hpp

@@ -152,17 +152,25 @@ class CNewtonIntegration final : public CIntegration {
   template<class T, su2enable_if<std::is_same<T,MixedScalar>::value> = 0>
   inline unsigned long Preconditioner_impl(const CSysVector<T>& u, CSysVector<T>& v,
                                            unsigned long iters, Scalar& eps) const {
-    if (iters == 0) {
+    const auto inner_solver = config->GetKind_Linear_Solver_Inner();
+
+    if (iters == 0 || (iters == 1 && inner_solver == LINEAR_SOLVER_INNER::SMOOTHER)) {
       (*preconditioner)(u, v);
       return 0;
     }
     auto product = CSysMatrixVectorProduct<MixedScalar>(solvers[FLOW_SOL]->Jacobian, geometry, config);
     v = MixedScalar(0.0);
     MixedScalar eps_t = eps;
-    if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::NONE) {
-      iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
-    } else {
-      iters = solvers[FLOW_SOL]->System.BCGSTAB_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+    switch (inner_solver) {
+      case LINEAR_SOLVER_INNER::NONE:
+        iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+        break;
+      case LINEAR_SOLVER_INNER::BCGSTAB:
+        iters = solvers[FLOW_SOL]->System.BCGSTAB_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config);
+        break;
+      case LINEAR_SOLVER_INNER::SMOOTHER:
+        iters = solvers[FLOW_SOL]->System.Smoother_LinSolver(u, v, product, *preconditioner, 0, iters, eps_t, false, config);
+        break;
     }
     eps = eps_t;
     return iters;

SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

@@ -80,7 +80,8 @@ class CFEAElasticity : public CNumerics {
   su2double *DV_Val = nullptr;            /*!< \brief For optimization cases, value of the design variables. */
   unsigned short n_DV = 0;                /*!< \brief For optimization cases, number of design variables. */
 
-  bool plane_stress = false;              /*!< \brief Checks if we are solving a plane stress case */
+  bool plane_stress = false;              /*!< \brief Checks if we are solving a plane stress case. */
+  bool linear = false;                    /*!< \brief Checks if we are solving a linear elasticity case. */
 
 public:
   /*!
@@ -181,28 +182,21 @@ class CFEAElasticity : public CNumerics {
 
   /*!
    * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz.
+   * \note Szz is required in 2D whereas Sxz and Syz are assumed to be 0.
    */
-  template<class T>
+  template <class T>
   static su2double VonMisesStress(unsigned short nDim, const T& stress) {
-    if (nDim == 2) {
-      su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2];
-
-      su2double S1, S2; S1 = S2 = (Sxx+Syy)/2;
-      su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2));
-      S1 += tauMax;
-      S2 -= tauMax;
-
-      return sqrt(S1*S1+S2*S2-2*S1*S2);
-    }
-    else {
-      su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3];
-      su2double Sxy = stress[2], Sxz = stress[4], Syz = stress[5];
-
-      return sqrt(0.5*(pow(Sxx - Syy, 2) +
-                       pow(Syy - Szz, 2) +
-                       pow(Szz - Sxx, 2) +
-                       6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz)));
+    /*--- In 2D, we only have 4 components: Sxx, Syy, Sxy, Szz. ---*/
+    const auto& Sxx = stress[0];
+    const auto& Syy = stress[1];
+    const auto& Sxy = stress[2];
+    const auto& Szz = stress[3];
+    su2double Sxz = 0, Syz = 0;
+    if (nDim == 3) {
+      Sxz = stress[4];
+      Syz = stress[5];
     }
+    return sqrt(0.5 * (pow(Sxx - Syy, 2) + pow(Syy - Szz, 2) + pow(Szz - Sxx, 2) + 6 * (Sxy * Sxy + Sxz * Sxz + Syz * Syz)));
   }
 
 protected:
@@ -233,8 +227,12 @@ class CFEAElasticity : public CNumerics {
     Mu     = E / (2.0*(1.0 + Nu));
     Lambda = Nu*E/((1.0+Nu)*(1.0-2.0*Nu));
     Kappa  = Lambda + (2/3)*Mu;
-    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php ---*/
-    ThermalStressTerm = -Alpha * E / (1 - (plane_stress ? 1 : 2) * Nu);
+    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php
+     * The stress tensor for nonlinear problems is still 3x3 and plane stress is imposed
+     * by determining the deformation that makes sigma_33 = 0, hence this denominator is
+     * only changed for linear elasticity. ---*/
+    const auto nu_mult = (linear && plane_stress) ? 1 : 2;
+    ThermalStressTerm = -Alpha * E / (1 - nu_mult * Nu);
   }
 
   /*!

SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp

@@ -134,8 +134,9 @@ class CFEANonlinearElasticity : public CFEAElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  virtual void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) = 0;
+  virtual void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) = 0;
 
   /*!
    * \brief Compute the stress tensor.