From 848b0978844afed5aecbefb95fea58f6eef18757 Mon Sep 17 00:00:00 2001 From: Ezgi Orbay Akcengiz <68737959+eorbay-metu@users.noreply.github.com> Date: Thu, 25 Dec 2025 03:45:32 +0300 Subject: [PATCH 1/6] Add nested preconditioner: FGMRES with BiCGSTAB (#2566) * Add nested preconditioner: FGMRES with BiCGSTAB * Fix style: replace tab with spaces in CConfig.cpp * Update CConfig.cpp * Add files via upload * Nested preconditioner test cases * Remove nested .git from TestCases/TestCases and add nested preconditioner test scripts * Apply clang-format style to CSysSolve nested solver * Add nested FGMRES inner BiCGSTAB options and logging * Add nested FGMRES inner BiCGSTAB test case * Fix formatting for nested FGMRES changes * Remove duplicated TestCases/TestCases directory * Add nested FGMRES with BCGSTAB test case cfg --------- Co-authored-by: Ezgi Orbay Akcengiz Co-authored-by: Nijso Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- Common/include/CConfig.hpp | 18 +++ Common/src/CConfig.cpp | 33 ++++- Common/src/linear_algebra/CSysSolve.cpp | 125 +++++++++++++++++- TestCases/parallel_regression.py | 0 .../rans/naca0012/turb_NACA0012_sa_nested.cfg | 117 ++++++++++++++++ .../turb_NACA0012_sa_nested.cfg | 119 +++++++++++++++++ TestCases/rans/turb_NACA0012_sa_nested.cfg | 108 +++++++++++++++ TestCases/serial_regression.py | 0 config_template.cfg | 8 ++ 9 files changed, 521 insertions(+), 7 deletions(-) mode change 100644 => 100755 TestCases/parallel_regression.py create mode 100644 TestCases/rans/naca0012/turb_NACA0012_sa_nested.cfg create mode 100644 TestCases/rans/naca0012_nested_fgmres/turb_NACA0012_sa_nested.cfg create mode 100644 TestCases/rans/turb_NACA0012_sa_nested.cfg mode change 100644 => 100755 TestCases/serial_regression.py diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 8abcc77e4fbe..cef7587f73c4 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -53,6 +53,12 @@ using namespace std; +/*Inner solver for nested linear solver:*/ +enum ENUM_LINEAR_SOLVER_INNER { + INNER_SOLVER_NONE = 0, + INNER_SOLVER_BCGSTAB +}; + /*! * \class CConfig * \brief Main class for defining the problem; basically this class reads the configuration file, and @@ -430,6 +436,7 @@ class CConfig { unsigned short nQuasiNewtonSamples; /*!< \brief Number of samples used in quasi-Newton solution methods. */ bool UseVectorization; /*!< \brief Whether to use vectorized numerics schemes. */ bool NewtonKrylov; /*!< \brief Use a coupled Newton method to solve the flow equations. */ + bool Nested_Linear_Solver; /*!< \brief Enable nested Krylov linear solver (e.g. FGMRES + inner solver). */ array NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */ array NK_DblParam{{-2.0, 0.1, -3.0, 1e-4, 1.0}}; /*!< \brief Floating-point parameters for NK method. */ su2double NK_Relaxation = 1.0; @@ -525,6 +532,7 @@ class CConfig { Kind_Deform_Linear_Solver, /*!< Numerical method to deform the grid */ Kind_Deform_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ Kind_Linear_Solver, /*!< \brief Numerical solver for the implicit scheme. */ + Kind_Linear_Solver_Inner, /*!< \brief Inner solver used in nested Krylov schemes. */ Kind_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ Kind_DiscAdj_Linear_Solver, /*!< \brief Linear solver for the discrete adjoint system. */ Kind_DiscAdj_Linear_Prec, /*!< \brief Preconditioner of the discrete adjoint linear solver. */ @@ -4288,6 +4296,16 @@ class CConfig { */ unsigned short GetKind_Linear_Solver(void) const { return Kind_Linear_Solver; } + /*! + * \brief Check whether nested Krylov linear solver is enabled. + */ + bool GetNested_Linear_Solver(void) const { return Nested_Linear_Solver; } + + /*! + * \brief Get the inner linear solver used in nested Krylov schemes. + */ + unsigned short GetKind_Linear_Solver_Inner(void) const { return Kind_Linear_Solver_Inner; } + /*! * \brief Get the kind of preconditioner for the implicit solver. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8a84e679ab69..31ec45ebeced 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -62,6 +62,12 @@ vector GEMM_Profile_MaxTime; /*!< \brief Maximum time spent for thi //#pragma omp threadprivate(Profile_Function_tp, Profile_Time_tp, Profile_ID_tp, Profile_Map_tp) +/* --- Map from config string to inner solver enum --- */ +std::map Inner_Linear_Solver_Map = { + {"NONE", INNER_SOLVER_NONE}, + {"BCGSTAB", INNER_SOLVER_BCGSTAB} +}; + CConfig::CConfig(char case_filename[MAX_STRING_SIZE], SU2_COMPONENT val_software, bool verb_high) { /*--- Set the case name to the base config file name without extension ---*/ @@ -1909,6 +1915,20 @@ void CConfig::SetConfig_Options() { addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU); /* DESCRIPTION: Linear solver for the discete adjoint systems */ + /*!\brief LINEAR_SOLVER_NESTED + * \n DESCRIPTION: Enable nested Krylov linear solver (e.g. FGMRES with inner solver). + * \n OPTIONS: YES, NO \ingroup Config */ + addBoolOption("LINEAR_SOLVER_NESTED", Nested_Linear_Solver, false); + + /*!\brief LINEAR_SOLVER_INNER + * \n DESCRIPTION: Inner linear solver used when LINEAR_SOLVER_NESTED = YES. + * \n OPTIONS: see \link Inner_Linear_Solver_Map \endlink + * \n DEFAULT: BCGSTAB \ingroup Config */ + addEnumOption("LINEAR_SOLVER_INNER", + Kind_Linear_Solver_Inner, + Inner_Linear_Solver_Map, + INNER_SOLVER_BCGSTAB); + /* DESCRIPTION: Maximum update ratio value for flow density and energy variables */ addDoubleOption("MAX_UPDATE_FLOW", MaxUpdateFlow, 0.2); /* DESCRIPTION: Maximum update ratio value for SA turbulence variable nu_tilde */ @@ -7260,10 +7280,17 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case BCGSTAB: case FGMRES: case RESTARTED_FGMRES: - if (Kind_Linear_Solver == BCGSTAB) + if (Kind_Linear_Solver == BCGSTAB){ cout << "BCGSTAB is used for solving the linear system." << endl; - else - cout << "FGMRES is used for solving the linear system." << endl; + } + else { + if (GetNested_Linear_Solver()){ + cout << "Nested FGMRES (FGMRES with inner BiCGSTAB) is used for solving the linear system." << endl; + } + else { + cout << "FGMRES is used for solving the linear system." << endl; + } + } switch (Kind_Linear_Solver_Prec) { case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; case LINELET: cout << "Using a linelet preconditioning."<< endl; break; diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 1a305c123dd3..0f99bc0d900a 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -172,6 +172,112 @@ void CSysSolve::ModGramSchmidt(bool shared_hsbg, int i, su2matrix +void BCGSTABpre_parallel(const CSysVector& a, CSysVector& b_in, + const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, + const CConfig* config) { + // NOTE: norm0_in is currently unused. To avoid -Werror unused-variable warnings, + // we mark it [[maybe_unused]] so that future changes can reuse it without + // triggering a build failure. + ScalarType norm_r_in = 0.0; + [[maybe_unused]] ScalarType norm0_in = 0.0; + unsigned long ii = 0; + + CSysVector A_z_i; + CSysVector r_0_in; + CSysVector r_in; + CSysVector p; + CSysVector v; + CSysVector z_i; + + /*--- Allocate if not allocated yet ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + auto nVar = a.GetNVar(); + auto nBlk = a.GetNBlk(); + auto nBlkDomain = a.GetNBlkDomain(); + + A_z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr); + r_0_in.Initialize(nBlk, nBlkDomain, nVar, nullptr); + r_in.Initialize(nBlk, nBlkDomain, nVar, nullptr); + p.Initialize(nBlk, nBlkDomain, nVar, nullptr); + v.Initialize(nBlk, nBlkDomain, nVar, nullptr); + z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + /*--- Calculate the initial residual, compute norm, and check if system is already solved ---*/ + mat_vec(b_in, A_z_i); + r_in = a - A_z_i; + + /*--- Only compute the residuals in full communication mode. ---*/ + if (config->GetComm_Level() == COMM_FULL) { + norm_r_in = r_in.norm(); + norm0_in = a.norm(); + /*--- Set the norm to the initial residual value ---*/ + /*--- if (tol_type == LinearToleranceType::RELATIVE) norm0_in = norm_r_in; ---*/ + } + + /*--- Initialization ---*/ + ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0; + p = ScalarType(0.0); + v = ScalarType(0.0); + r_0_in = r_in; + + ScalarType tolerance = 1e-5; // Tolerance for the residual norm + + /*--- Loop over all search directions ---*/ + while (ii < 1000) { // Arbitrary high iteration limit for safety + /*--- Compute rho_prime ---*/ + rho_prime = rho; + + /*--- Compute rho_i ---*/ + rho = r_in.dot(r_0_in); + + /*--- Compute beta ---*/ + ScalarType beta_in = (rho / rho_prime) * (alpha / omega); + + /*--- Update p ---*/ + p = beta_in * (p - omega * v) + r_in; + + /*--- Preconditioning step ---*/ + precond(p, z_i); + mat_vec(z_i, v); + + /*--- Calculate step-length alpha ---*/ + ScalarType r_0_v = r_0_in.dot(v); + alpha = rho / r_0_v; + + /*--- Update solution and residual ---*/ + b_in += alpha * z_i; + r_in -= alpha * v; + + /*--- Preconditioning step ---*/ + precond(r_in, z_i); + mat_vec(z_i, A_z_i); + + /*--- Calculate step-length omega, avoid division by 0. ---*/ + omega = A_z_i.squaredNorm(); + if (omega == ScalarType(0)) break; + omega = A_z_i.dot(r_in) / omega; + + /*--- Update solution and residual ---*/ + b_in += omega * z_i; + r_in -= omega * A_z_i; + + /*--- Update the residual norm ---*/ + norm_r_in = r_in.norm(); + + /*--- Check if residual norm is below tolerance ---*/ + if (norm_r_in < tolerance) { + break; // Stop if the residual norm is below the desired tolerance + } + + ii++; // Increment iteration counter + } +} + template void CSysSolve::WriteHeader(const string& solver, ScalarType restol, ScalarType resinit) const { cout << "\n# " << solver << " residual history\n"; @@ -359,6 +465,11 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector 1; + /*--- Nested inner solver (BCGSTAB) option ---*/ + const bool useNestedBiCGSTAB = + config->GetNested_Linear_Solver() && + (config->GetKind_Linear_Solver_Inner() == static_cast(INNER_SOLVER_BCGSTAB)); + /*--- Check the subspace size ---*/ if (m < 1) { @@ -454,13 +565,19 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector Date: Wed, 24 Dec 2025 17:53:10 -0800 Subject: [PATCH 2/6] cleanup --- Common/include/CConfig.hpp | 18 +- .../linear_algebra/CPreconditioner.hpp | 20 ++ Common/include/linear_algebra/CSysSolve.hpp | 4 + Common/include/option_structure.hpp | 13 ++ Common/src/CConfig.cpp | 41 +--- Common/src/linear_algebra/CSysSolve.cpp | 210 ++++++------------ .../rans/naca0012/turb_NACA0012_sa_nested.cfg | 117 ---------- .../turb_NACA0012_sa_nested.cfg | 119 ---------- TestCases/rans/turb_NACA0012_sa_nested.cfg | 108 --------- config_template.cfg | 10 +- 10 files changed, 122 insertions(+), 538 deletions(-) delete mode 100644 TestCases/rans/naca0012/turb_NACA0012_sa_nested.cfg delete mode 100644 TestCases/rans/naca0012_nested_fgmres/turb_NACA0012_sa_nested.cfg delete mode 100644 TestCases/rans/turb_NACA0012_sa_nested.cfg diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index cef7587f73c4..aa454dc35bc2 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -53,12 +53,6 @@ using namespace std; -/*Inner solver for nested linear solver:*/ -enum ENUM_LINEAR_SOLVER_INNER { - INNER_SOLVER_NONE = 0, - INNER_SOLVER_BCGSTAB -}; - /*! * \class CConfig * \brief Main class for defining the problem; basically this class reads the configuration file, and @@ -522,6 +516,7 @@ class CConfig { Kind_SlopeLimit_AdjFlow, /*!< \brief Slope limiter for the adjoint equation.*/ Kind_SlopeLimit_Heat, /*!< \brief Slope limiter for the adjoint equation.*/ Kind_SlopeLimit_Species; /*!< \brief Slope limiter for the species equation.*/ + LINEAR_SOLVER_INNER Kind_Linear_Solver_Inner; /*!< \brief Inner solver used in nested Krylov schemes. */ unsigned short Kind_FluidModel, /*!< \brief Kind of the Fluid Model: Ideal, van der Waals, etc. */ Kind_InitOption, /*!< \brief Kind of Init option to choose if initializing with Reynolds number or with thermodynamic conditions */ Kind_GridMovement, /*!< \brief Kind of the static mesh movement. */ @@ -532,7 +527,6 @@ class CConfig { Kind_Deform_Linear_Solver, /*!< Numerical method to deform the grid */ Kind_Deform_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ Kind_Linear_Solver, /*!< \brief Numerical solver for the implicit scheme. */ - Kind_Linear_Solver_Inner, /*!< \brief Inner solver used in nested Krylov schemes. */ Kind_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ Kind_DiscAdj_Linear_Solver, /*!< \brief Linear solver for the discrete adjoint system. */ Kind_DiscAdj_Linear_Prec, /*!< \brief Preconditioner of the discrete adjoint linear solver. */ @@ -4297,15 +4291,9 @@ class CConfig { unsigned short GetKind_Linear_Solver(void) const { return Kind_Linear_Solver; } /*! - * \brief Check whether nested Krylov linear solver is enabled. - */ - bool GetNested_Linear_Solver(void) const { return Nested_Linear_Solver; } - - /*! - * \brief Get the inner linear solver used in nested Krylov schemes. + * \brief Get the inner linear solver used in nested Krylov linear solvers. */ - unsigned short GetKind_Linear_Solver_Inner(void) const { return Kind_Linear_Solver_Inner; } - + LINEAR_SOLVER_INNER GetKind_Linear_Solver_Inner(void) const { return Kind_Linear_Solver_Inner; } /*! * \brief Get the kind of preconditioner for the implicit solver. diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 077d0205ec4a..1a07d38813e8 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -28,6 +28,7 @@ #pragma once +#include #include "../CConfig.hpp" #include "../geometry/CGeometry.hpp" #include "CSysVector.hpp" @@ -305,6 +306,25 @@ class CPastixPreconditioner final : public CPreconditioner { inline void Build() override { sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact); } }; +/*! + * \class CAbstractPreconditioner + * \brief Applies a std::function as the preconditioning operation. + * \note This can be used to treat almost anything as a preconditioner. + */ +template +class CAbstractPreconditioner final : public CPreconditioner { + private: + std::function&, CSysVector&)> impl; + + public: + CAbstractPreconditioner() = delete; + + template + explicit CAbstractPreconditioner(const F& function) : impl(function) {} + + inline void operator()(const CSysVector& u, CSysVector& v) const override { impl(u, v); } +}; + template CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOLVER_PREC kind, CSysMatrix& jacobian, CGeometry* geometry, diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index efef9d65b74c..5db5c2a1cc3b 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -31,6 +31,7 @@ #include "../containers/C2DContainer.hpp" #include +#include #include #include #include @@ -115,6 +116,9 @@ class CSysSolve { bool recomputeRes = false; /*!< \brief Recompute the residual after inner iterations, if monitoring. */ unsigned long monitorFreq = 10; /*!< \brief Monitoring frequency. */ + /*!< \brief Inner solver for nested preconditioning. */ + std::unique_ptr> inner_solver; + /*! * \brief sign transfer function * \param[in] x - value having sign prescribed diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b1387..88bc7d4f8743 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2369,6 +2369,19 @@ static const MapType Linear_Solver_Map = { MakePair("PASTIX_LU", PASTIX_LU) }; +/*! + * \brief Inner solver for nested linear solver, only compatible with "flexible" linear solvers. + */ +enum class LINEAR_SOLVER_INNER { + NONE, /*!< \brief Do not use a nested linear solver. */ + BCGSTAB, /*!< \brief Use BCGSTAB as the preconditioning linear solver. */ +}; +static const MapType Inner_Linear_Solver_Map = { + MakePair("NONE", LINEAR_SOLVER_INNER::NONE) + MakePair("BCGSTAB", LINEAR_SOLVER_INNER::BCGSTAB) +}; + + /*! * \brief Types surface continuity at the intersection with the FFD */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 31ec45ebeced..4d1a262d4b1f 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -62,12 +62,6 @@ vector GEMM_Profile_MaxTime; /*!< \brief Maximum time spent for thi //#pragma omp threadprivate(Profile_Function_tp, Profile_Time_tp, Profile_ID_tp, Profile_Map_tp) -/* --- Map from config string to inner solver enum --- */ -std::map Inner_Linear_Solver_Map = { - {"NONE", INNER_SOLVER_NONE}, - {"BCGSTAB", INNER_SOLVER_BCGSTAB} -}; - CConfig::CConfig(char case_filename[MAX_STRING_SIZE], SU2_COMPONENT val_software, bool verb_high) { /*--- Set the case name to the base config file name without extension ---*/ @@ -1893,6 +1887,8 @@ void CConfig::SetConfig_Options() { addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0); /* DESCRIPTION: Custom number of threads used for additive domain decomposition for ILU and LU_SGS (0 is "auto"). */ addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0); + /* DESCRIPTION: Use an inner linear solver. */ + addEnumOption("LINEAR_SOLVER_INNER", Kind_Linear_Solver_Inner, Inner_Linear_Solver_Map, LINEAR_SOLVER_INNER::NONE); /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */ addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0); /* DESCRIPTION: Relaxation of the CHT coupling */ @@ -1913,21 +1909,6 @@ void CConfig::SetConfig_Options() { addEnumOption("DISCADJ_LIN_SOLVER", Kind_DiscAdj_Linear_Solver, Linear_Solver_Map, FGMRES); /* DESCRIPTION: Preconditioner for the discrete adjoint Krylov linear solvers */ addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU); - /* DESCRIPTION: Linear solver for the discete adjoint systems */ - - /*!\brief LINEAR_SOLVER_NESTED - * \n DESCRIPTION: Enable nested Krylov linear solver (e.g. FGMRES with inner solver). - * \n OPTIONS: YES, NO \ingroup Config */ - addBoolOption("LINEAR_SOLVER_NESTED", Nested_Linear_Solver, false); - - /*!\brief LINEAR_SOLVER_INNER - * \n DESCRIPTION: Inner linear solver used when LINEAR_SOLVER_NESTED = YES. - * \n OPTIONS: see \link Inner_Linear_Solver_Map \endlink - * \n DEFAULT: BCGSTAB \ingroup Config */ - addEnumOption("LINEAR_SOLVER_INNER", - Kind_Linear_Solver_Inner, - Inner_Linear_Solver_Map, - INNER_SOLVER_BCGSTAB); /* DESCRIPTION: Maximum update ratio value for flow density and energy variables */ addDoubleOption("MAX_UPDATE_FLOW", MaxUpdateFlow, 0.2); @@ -7280,17 +7261,15 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case BCGSTAB: case FGMRES: case RESTARTED_FGMRES: - if (Kind_Linear_Solver == BCGSTAB){ + if (Kind_Linear_Solver == BCGSTAB) { cout << "BCGSTAB is used for solving the linear system." << endl; + } else { + if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB){ + cout << "Nested FGMRES (FGMRES with inner BiCGSTAB) is used for solving the linear system." << endl; + } else { + cout << "FGMRES is used for solving the linear system." << endl; + } } - else { - if (GetNested_Linear_Solver()){ - cout << "Nested FGMRES (FGMRES with inner BiCGSTAB) is used for solving the linear system." << endl; - } - else { - cout << "FGMRES is used for solving the linear system." << endl; - } - } switch (Kind_Linear_Solver_Prec) { case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; case LINELET: cout << "Using a linelet preconditioning."<< endl; break; @@ -7337,6 +7316,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "FGMRES is used for solving the linear system." << endl; cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl; cout << "Max number of iterations: "<< Linear_Solver_Iter <<"."<< endl; + if (Kind_Linear_Solver_Inner == LINEAR_SOLVER_INNER::BCGSTAB) + cout << "Nested BiCGSTAB is used as the inner solver." << endl; break; case CONJUGATE_GRADIENT: cout << "A Conjugate Gradient method is used for solving the linear system." << endl; diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 0f99bc0d900a..69702c51507a 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -35,6 +35,7 @@ #include "../../include/linear_algebra/CPreconditioner.hpp" #include +#include /*! * \brief Epsilon used in CSysSolve depending on datatype to @@ -110,7 +111,7 @@ void CSysSolve::SolveReduced(int n, const su2matrix& Hsb template void CSysSolve::ModGramSchmidt(bool shared_hsbg, int i, su2matrix& Hsbg, - vector >& w) const { + vector>& w) const { const auto thread = omp_get_thread_num(); /*--- If Hsbg is shared by multiple threads calling this function, only one @@ -172,112 +173,6 @@ void CSysSolve::ModGramSchmidt(bool shared_hsbg, int i, su2matrix -void BCGSTABpre_parallel(const CSysVector& a, CSysVector& b_in, - const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, - const CConfig* config) { - // NOTE: norm0_in is currently unused. To avoid -Werror unused-variable warnings, - // we mark it [[maybe_unused]] so that future changes can reuse it without - // triggering a build failure. - ScalarType norm_r_in = 0.0; - [[maybe_unused]] ScalarType norm0_in = 0.0; - unsigned long ii = 0; - - CSysVector A_z_i; - CSysVector r_0_in; - CSysVector r_in; - CSysVector p; - CSysVector v; - CSysVector z_i; - - /*--- Allocate if not allocated yet ---*/ - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - auto nVar = a.GetNVar(); - auto nBlk = a.GetNBlk(); - auto nBlkDomain = a.GetNBlkDomain(); - - A_z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr); - r_0_in.Initialize(nBlk, nBlkDomain, nVar, nullptr); - r_in.Initialize(nBlk, nBlkDomain, nVar, nullptr); - p.Initialize(nBlk, nBlkDomain, nVar, nullptr); - v.Initialize(nBlk, nBlkDomain, nVar, nullptr); - z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr); - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - - /*--- Calculate the initial residual, compute norm, and check if system is already solved ---*/ - mat_vec(b_in, A_z_i); - r_in = a - A_z_i; - - /*--- Only compute the residuals in full communication mode. ---*/ - if (config->GetComm_Level() == COMM_FULL) { - norm_r_in = r_in.norm(); - norm0_in = a.norm(); - /*--- Set the norm to the initial residual value ---*/ - /*--- if (tol_type == LinearToleranceType::RELATIVE) norm0_in = norm_r_in; ---*/ - } - - /*--- Initialization ---*/ - ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0; - p = ScalarType(0.0); - v = ScalarType(0.0); - r_0_in = r_in; - - ScalarType tolerance = 1e-5; // Tolerance for the residual norm - - /*--- Loop over all search directions ---*/ - while (ii < 1000) { // Arbitrary high iteration limit for safety - /*--- Compute rho_prime ---*/ - rho_prime = rho; - - /*--- Compute rho_i ---*/ - rho = r_in.dot(r_0_in); - - /*--- Compute beta ---*/ - ScalarType beta_in = (rho / rho_prime) * (alpha / omega); - - /*--- Update p ---*/ - p = beta_in * (p - omega * v) + r_in; - - /*--- Preconditioning step ---*/ - precond(p, z_i); - mat_vec(z_i, v); - - /*--- Calculate step-length alpha ---*/ - ScalarType r_0_v = r_0_in.dot(v); - alpha = rho / r_0_v; - - /*--- Update solution and residual ---*/ - b_in += alpha * z_i; - r_in -= alpha * v; - - /*--- Preconditioning step ---*/ - precond(r_in, z_i); - mat_vec(z_i, A_z_i); - - /*--- Calculate step-length omega, avoid division by 0. ---*/ - omega = A_z_i.squaredNorm(); - if (omega == ScalarType(0)) break; - omega = A_z_i.dot(r_in) / omega; - - /*--- Update solution and residual ---*/ - b_in += omega * z_i; - r_in -= omega * A_z_i; - - /*--- Update the residual norm ---*/ - norm_r_in = r_in.norm(); - - /*--- Check if residual norm is below tolerance ---*/ - if (norm_r_in < tolerance) { - break; // Stop if the residual norm is below the desired tolerance - } - - ii++; // Increment iteration counter - } -} - template void CSysSolve::WriteHeader(const string& solver, ScalarType restol, ScalarType resinit) const { cout << "\n# " << solver << " residual history\n"; @@ -465,11 +360,6 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector 1; - /*--- Nested inner solver (BCGSTAB) option ---*/ - const bool useNestedBiCGSTAB = - config->GetNested_Linear_Solver() && - (config->GetKind_Linear_Solver_Inner() == static_cast(INNER_SOLVER_BCGSTAB)); - /*--- Check the subspace size ---*/ if (m < 1) { @@ -565,19 +455,8 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector::Solve(CSysMatrix& Jacobian, con break; } - /*--- Normal mode - * assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', but does not enforce it to avoid compiler warning. - * ---*/ + /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', + * but does not enforce it to avoid compiler warning. ---*/ default: { KindSolver = config->GetKind_Linear_Solver(); KindPrecond = config->GetKind_Linear_Solver_Prec(); @@ -1032,6 +910,17 @@ unsigned long CSysSolve::Solve(CSysMatrix& 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>(LINEAR_SOLVER_MODE::STANDARD); + inner_solver->SetxIsZero(true); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- Stop the recording for the linear solver ---*/ bool TapeActive = NO; @@ -1065,13 +954,25 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); - const auto kindPrec = static_cast(KindPrecond); - - auto precond = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); - /*--- Build preconditioner. ---*/ - precond->Build(); + const auto kindPrec = static_cast(KindPrecond); + auto* normal_prec = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); + normal_prec->Build(); + + CPreconditioner* nested_prec = nullptr; + if (nested) { + nested_prec = new CAbstractPreconditioner([&](const CSysVector& u, + CSysVector& v) { + /*--- Initialize to 0 to be safe. ---*/ + 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, SolverTol, MaxIter, unused, false, config); + }); + } + const auto* precond = nested ? nested_prec : normal_prec; /*--- Solve system. ---*/ @@ -1117,7 +1018,8 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con HandleTemporariesOut(LinSysSol); - delete precond; + delete normal_prec; + delete nested_prec; if (TapeActive) { /*--- To keep the behavior of SU2_DOT, but not strictly required since jacobian is symmetric(?). ---*/ @@ -1202,9 +1104,8 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c break; } - /*--- Normal mode - * assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', but does not enforce it to avoid compiler warning. - * ---*/ + /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', + * but does not enforce it to avoid compiler warning. ---*/ default: { KindSolver = config->GetKind_Linear_Solver(); KindPrecond = config->GetKind_Linear_Solver_Prec(); @@ -1215,19 +1116,43 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& 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>(LINEAR_SOLVER_MODE::STANDARD); + inner_solver->SetxIsZero(true); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + /*--- Set up preconditioner and matrix-vector product ---*/ - const auto kindPrec = static_cast(KindPrecond); + auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); - auto precond = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); + const auto kindPrec = static_cast(KindPrecond); + auto* normal_prec = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); /*--- If there was no call to solve first the preconditioner needs to be built here. ---*/ if (directCall) { Jacobian.TransposeInPlace(); - precond->Build(); + normal_prec->Build(); } - auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); + CPreconditioner* nested_prec = nullptr; + if (nested) { + nested_prec = + new CAbstractPreconditioner([&](const CSysVector& u, CSysVector& v) { + /*--- Initialize to 0 to be safe. ---*/ + 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, SolverTol, MaxIter, unused, false, config); + }); + } + const auto* precond = nested ? nested_prec : normal_prec; /*--- Solve the system ---*/ @@ -1271,7 +1196,8 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c HandleTemporariesOut(LinSysSol); - delete precond; + delete normal_prec; + delete nested_prec; SU2_OMP_MASTER { Residual = residual; diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa_nested.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa_nested.cfg deleted file mode 100644 index 29695e4ad4f6..000000000000 --- a/TestCases/rans/naca0012/turb_NACA0012_sa_nested.cfg +++ /dev/null @@ -1,117 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% SU2 configuration file % -% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % -% http://turbmodels.larc.nasa.gov/naca0012_val_sa.html % -% Author: Francisco Palacios % -% Institution: Stanford University % -% Date: Feb 18th, 2013 % -% File Version 8.3.0 "Harrier" % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% -% -SOLVER= RANS -KIND_TURB_MODEL= SA -SA_OPTIONS= NEGATIVE, EXPERIMENTAL -MATH_PROBLEM= DIRECT -RESTART_SOL= NO - -% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% -% -MACH_NUMBER= 0.15 -AOA= 0.0 -FREESTREAM_TEMPERATURE= 300.0 -REYNOLDS_NUMBER= 6.0E6 -REYNOLDS_LENGTH= 1.0 - -% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% -% -REF_ORIGIN_MOMENT_X = 0.25 -REF_ORIGIN_MOMENT_Y = 0.00 -REF_ORIGIN_MOMENT_Z = 0.00 -REF_LENGTH= 1.0 -REF_AREA= 1.0 -REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE - -% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% -% -MARKER_HEATFLUX= ( airfoil, 0.0 ) -MARKER_FAR= ( farfield ) -MARKER_PLOTTING= ( airfoil ) -MARKER_MONITORING= ( airfoil ) - -% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% -% -NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -NUM_METHOD_GRAD_RECON= LEAST_SQUARES -CFL_NUMBER= 1000.0 -MAX_DELTA_TIME= 1E10 -CFL_ADAPT= NO -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) -ITER= 99999 - -% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% -VENKAT_LIMITER_COEFF= 0.03 -LIMITER_ITER= 99999 - -% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% -% -LINEAR_SOLVER= FGMRES -% -% Use nested Krylov solver (YES, NO). -% When YES: FGMRES uses an inner Krylov solver (e.g. BCGSTAB). -LINEAR_SOLVER_NESTED= YES - -% Inner solver used only when LINEAR_SOLVER_NESTED = YES. -% Options: BCGSTAB -LINEAR_SOLVER_INNER= BCGSTAB - -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-5 -LINEAR_SOLVER_ITER= 5 - -% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% -% -CONV_NUM_METHOD_FLOW= ROE -MUSCL_FLOW= YES -JST_SENSOR_COEFF= ( 0.5, 0.02 ) -TIME_DISCRE_FLOW= EULER_IMPLICIT - -% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% -% -CONV_NUM_METHOD_TURB= SCALAR_UPWIND -MUSCL_TURB= YES -SLOPE_LIMITER_TURB= NONE -TIME_DISCRE_TURB= EULER_IMPLICIT -CFL_REDUCTION_TURB= 1.0 - -% --------------------------- CONVERGENCE PARAMETERS --------------------------% -CONV_RESIDUAL_MINVAL= -12 -CONV_STARTITER= 10 -CONV_CAUCHY_ELEMS= 100 -CONV_CAUCHY_EPS= 1E-6 - -% Output the performance summary to the console at the end of SU2_CFD -WRT_PERFORMANCE= YES - - -% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% -% -MESH_FILENAME= n0012_897-257.su2 -MESH_FORMAT= SU2 -MESH_OUT_FILENAME= mesh_out -SOLUTION_FILENAME= solution_flow_sa -SOLUTION_ADJ_FILENAME= solution_adj -TABULAR_FORMAT= CSV -CONV_FILENAME= history -RESTART_FILENAME= restart_flow -RESTART_ADJ_FILENAME= restart_adj -VOLUME_FILENAME= flow -VOLUME_ADJ_FILENAME= adjoint -GRAD_OBJFUNC_FILENAME= of_grad -SURFACE_FILENAME= surface_flow -SURFACE_ADJ_FILENAME= surface_adjoint -OUTPUT_WRT_FREQ= 10000 -SCREEN_OUTPUT=(WALL_TIME, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, LINSOL_ITER, LINSOL_RESIDUAL, LINSOL_ITER_TURB, LINSOL_RESIDUAL_TURB, TOTAL_HEATFLUX) diff --git a/TestCases/rans/naca0012_nested_fgmres/turb_NACA0012_sa_nested.cfg b/TestCases/rans/naca0012_nested_fgmres/turb_NACA0012_sa_nested.cfg deleted file mode 100644 index d9991fcf2d24..000000000000 --- a/TestCases/rans/naca0012_nested_fgmres/turb_NACA0012_sa_nested.cfg +++ /dev/null @@ -1,119 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% SU2 configuration file % -% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % -% http://turbmodels.larc.nasa.gov/naca0012_val_sa.html % -% Author: Francisco Palacios % -% Institution: Stanford University % -% Date: Feb 18th, 2013 % -% File Version 8.3.0 "Harrier" % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% -% -SOLVER= RANS -KIND_TURB_MODEL= SA -SA_OPTIONS= NEGATIVE, EXPERIMENTAL -MATH_PROBLEM= DIRECT -RESTART_SOL= NO - -% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% -% -MACH_NUMBER= 0.05 -AOA= 0.0 -FREESTREAM_TEMPERATURE= 300.0 -REYNOLDS_NUMBER= 6.0E6 -REYNOLDS_LENGTH= 1.0 - -% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% -% -REF_ORIGIN_MOMENT_X = 0.25 -REF_ORIGIN_MOMENT_Y = 0.00 -REF_ORIGIN_MOMENT_Z = 0.00 -REF_LENGTH= 1.0 -REF_AREA= 1.0 -REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE - -% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% -% -MARKER_HEATFLUX= ( airfoil, 0.0 ) -MARKER_FAR= ( farfield ) -MARKER_PLOTTING= ( airfoil ) -MARKER_MONITORING= ( airfoil ) - -% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% -% -NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -NUM_METHOD_GRAD_RECON= LEAST_SQUARES -CFL_NUMBER= 200.0 -MAX_DELTA_TIME= 1E10 -CFL_ADAPT= NO -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) -ITER= 99999 - -% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% -VENKAT_LIMITER_COEFF= 0.03 -LIMITER_ITER= 99999 - -% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% -% -% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% -% -LINEAR_SOLVER= FGMRES -% -% Use nested Krylov solver (YES, NO). -% When YES: FGMRES uses an inner Krylov solver (e.g. BCGSTAB). -LINEAR_SOLVER_NESTED= YES - -% Inner solver used only when LINEAR_SOLVER_NESTED = YES. -% Options: BCGSTAB -LINEAR_SOLVER_INNER= BCGSTAB - -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-5 -LINEAR_SOLVER_ITER= 5 - -% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% -% -CONV_NUM_METHOD_FLOW= ROE -MUSCL_FLOW= YES -JST_SENSOR_COEFF= ( 0.5, 0.02 ) -TIME_DISCRE_FLOW= EULER_IMPLICIT - -% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% -% -CONV_NUM_METHOD_TURB= SCALAR_UPWIND -MUSCL_TURB= YES -SLOPE_LIMITER_TURB= NONE -TIME_DISCRE_TURB= EULER_IMPLICIT -CFL_REDUCTION_TURB= 1.0 - -% --------------------------- CONVERGENCE PARAMETERS --------------------------% -CONV_RESIDUAL_MINVAL= -12 -CONV_STARTITER= 10 -CONV_CAUCHY_ELEMS= 100 -CONV_CAUCHY_EPS= 1E-6 - -% Output the performance summary to the console at the end of SU2_CFD -WRT_PERFORMANCE= YES - - -% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% -% -MESH_FILENAME= n0012_897-257.su2 -MESH_FORMAT= SU2 -MESH_OUT_FILENAME= mesh_out -SOLUTION_FILENAME= solution_flow_sa -SOLUTION_ADJ_FILENAME= solution_adj -TABULAR_FORMAT= CSV -CONV_FILENAME= history -RESTART_FILENAME= restart_flow -RESTART_ADJ_FILENAME= restart_adj -VOLUME_FILENAME= flow -VOLUME_ADJ_FILENAME= adjoint -GRAD_OBJFUNC_FILENAME= of_grad -SURFACE_FILENAME= surface_flow -SURFACE_ADJ_FILENAME= surface_adjoint -OUTPUT_WRT_FREQ= 10000 -SCREEN_OUTPUT=(WALL_TIME, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, LINSOL_ITER, LINSOL_RESIDUAL, LINSOL_ITER_TURB, LINSOL_RESIDUAL_TURB, TOTAL_HEATFLUX) diff --git a/TestCases/rans/turb_NACA0012_sa_nested.cfg b/TestCases/rans/turb_NACA0012_sa_nested.cfg deleted file mode 100644 index 9d6d3f973b9c..000000000000 --- a/TestCases/rans/turb_NACA0012_sa_nested.cfg +++ /dev/null @@ -1,108 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% SU2 configuration file % -% Case description: 2D NACA 0012 Airfoil Validation Case (compressible) % -% http://turbmodels.larc.nasa.gov/naca0012_val_sa.html % -% Author: Francisco Palacios % -% Institution: Stanford University % -% Date: Feb 18th, 2013 % -% File Version 8.3.0 "Harrier" % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% -% -SOLVER= RANS -KIND_TURB_MODEL= SA -SA_OPTIONS= NEGATIVE, EXPERIMENTAL -MATH_PROBLEM= DIRECT -RESTART_SOL= NO - -% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% -% -MACH_NUMBER= 0.05 -AOA= 0.0 -FREESTREAM_TEMPERATURE= 300.0 -REYNOLDS_NUMBER= 6.0E6 -REYNOLDS_LENGTH= 1.0 - -% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% -% -REF_ORIGIN_MOMENT_X = 0.25 -REF_ORIGIN_MOMENT_Y = 0.00 -REF_ORIGIN_MOMENT_Z = 0.00 -REF_LENGTH= 1.0 -REF_AREA= 1.0 -REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE - -% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% -% -MARKER_HEATFLUX= ( airfoil, 0.0 ) -MARKER_FAR= ( farfield ) -MARKER_PLOTTING= ( airfoil ) -MARKER_MONITORING= ( airfoil ) - -% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% -% -NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -NUM_METHOD_GRAD_RECON= LEAST_SQUARES -CFL_NUMBER= 200.0 -MAX_DELTA_TIME= 1E10 -CFL_ADAPT= NO -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) -ITER= 99999 - -% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% -VENKAT_LIMITER_COEFF= 0.03 -LIMITER_ITER= 99999 - -% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% -% -LINEAR_SOLVER= FGMRESandBCGSTAB2 -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-5 -LINEAR_SOLVER_ITER= 5 - -% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% -% -CONV_NUM_METHOD_FLOW= ROE -MUSCL_FLOW= YES -JST_SENSOR_COEFF= ( 0.5, 0.02 ) -TIME_DISCRE_FLOW= EULER_IMPLICIT - -% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% -% -CONV_NUM_METHOD_TURB= SCALAR_UPWIND -MUSCL_TURB= YES -SLOPE_LIMITER_TURB= NONE -TIME_DISCRE_TURB= EULER_IMPLICIT -CFL_REDUCTION_TURB= 1.0 - -% --------------------------- CONVERGENCE PARAMETERS --------------------------% -CONV_RESIDUAL_MINVAL= -12 -CONV_STARTITER= 10 -CONV_CAUCHY_ELEMS= 100 -CONV_CAUCHY_EPS= 1E-6 - -% Output the performance summary to the console at the end of SU2_CFD -WRT_PERFORMANCE= YES - - -% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% -% -MESH_FILENAME= n0012_897-257.su2 -MESH_FORMAT= SU2 -MESH_OUT_FILENAME= mesh_out -SOLUTION_FILENAME= solution_flow_sa -SOLUTION_ADJ_FILENAME= solution_adj -TABULAR_FORMAT= CSV -CONV_FILENAME= history -RESTART_FILENAME= restart_flow -RESTART_ADJ_FILENAME= restart_adj -VOLUME_FILENAME= flow -VOLUME_ADJ_FILENAME= adjoint -GRAD_OBJFUNC_FILENAME= of_grad -SURFACE_FILENAME= surface_flow -SURFACE_ADJ_FILENAME= surface_adjoint -OUTPUT_WRT_FREQ= 10000 -SCREEN_OUTPUT=(WALL_TIME, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, LINSOL_ITER, LINSOL_RESIDUAL, LINSOL_ITER_TURB, LINSOL_RESIDUAL_TURB, TOTAL_HEATFLUX) diff --git a/config_template.cfg b/config_template.cfg index 9dd297cd90c5..93dd76e84638 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1543,14 +1543,10 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % Linear solver or smoother for implicit formulations: % BCGSTAB, FGMRES, RESTARTED_FGMRES, CONJUGATE_GRADIENT (self-adjoint problems only), SMOOTHER. LINEAR_SOLVER= FGMRES -% -% Use nested Krylov solver (YES, NO). -% When YES: FGMRES uses an inner Krylov solver (e.g. BCGSTAB). -LINEAR_SOLVER_NESTED= NO -% Inner solver used only when LINEAR_SOLVER_NESTED = YES. -% Options: BCGSTAB -LINEAR_SOLVER_INNER= BCGSTAB +% Inner solver used for nested preconditioning (only possible with flexible linear solvers). +% Options: NONE, BCGSTAB +LINEAR_SOLVER_INNER= NONE % Same for discrete adjoint (smoothers not supported), replaces LINEAR_SOLVER in SU2_*_AD codes. DISCADJ_LIN_SOLVER= FGMRES From 331071684349d26591bc39bdc207966e6b7c824a Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 24 Dec 2025 18:10:55 -0800 Subject: [PATCH 3/6] apply to NK --- Common/include/CConfig.hpp | 1 - SU2_CFD/include/integration/CNewtonIntegration.hpp | 6 +++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index aa454dc35bc2..48aecba36956 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -430,7 +430,6 @@ class CConfig { unsigned short nQuasiNewtonSamples; /*!< \brief Number of samples used in quasi-Newton solution methods. */ bool UseVectorization; /*!< \brief Whether to use vectorized numerics schemes. */ bool NewtonKrylov; /*!< \brief Use a coupled Newton method to solve the flow equations. */ - bool Nested_Linear_Solver; /*!< \brief Enable nested Krylov linear solver (e.g. FGMRES + inner solver). */ array NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */ array NK_DblParam{{-2.0, 0.1, -3.0, 1e-4, 1.0}}; /*!< \brief Floating-point parameters for NK method. */ su2double NK_Relaxation = 1.0; diff --git a/SU2_CFD/include/integration/CNewtonIntegration.hpp b/SU2_CFD/include/integration/CNewtonIntegration.hpp index f11d4dadff47..7198c34a4ed9 100644 --- a/SU2_CFD/include/integration/CNewtonIntegration.hpp +++ b/SU2_CFD/include/integration/CNewtonIntegration.hpp @@ -158,7 +158,11 @@ class CNewtonIntegration final : public CIntegration { auto product = CSysMatrixVectorProduct(solvers[FLOW_SOL]->Jacobian, geometry, config); v = MixedScalar(0.0); MixedScalar eps_t = eps; - iters = solvers[FLOW_SOL]->System.FGMRES_LinSolver(u, v, product, *preconditioner, eps, iters, eps_t, false, config); + 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); + } eps = eps_t; return iters; } From a10037243947ba9638d8ad647665c3d9d8a54e2f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 24 Dec 2025 18:20:24 -0800 Subject: [PATCH 4/6] modify a test --- TestCases/incomp_rans/naca0012/naca0012.cfg | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/TestCases/incomp_rans/naca0012/naca0012.cfg b/TestCases/incomp_rans/naca0012/naca0012.cfg index a3a0b25229b3..1ef87a5a4a48 100644 --- a/TestCases/incomp_rans/naca0012/naca0012.cfg +++ b/TestCases/incomp_rans/naca0012/naca0012.cfg @@ -50,8 +50,8 @@ MARKER_MONITORING= ( airfoil ) % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES CFL_NUMBER= 10.0 -CFL_ADAPT= NO -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +CFL_ADAPT= YES +CFL_ADAPT_PARAM= ( 0.8, 1.1, 1.0, 100.0 ) ITER= 2500 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% @@ -64,6 +64,7 @@ SENS_REMOVE_SHARP= NO % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_INNER= BCGSTAB LINEAR_SOLVER_PREC= LU_SGS LINEAR_SOLVER_ERROR= 1E-4 LINEAR_SOLVER_ITER= 5 @@ -107,4 +108,4 @@ GRAD_OBJFUNC_FILENAME= of_grad SURFACE_FILENAME= surface_flow SURFACE_ADJ_FILENAME= surface_adjoint OUTPUT_WRT_FREQ= 100 -SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_NU_TILDE, LIFT, DRAG, TOTAL_HEATFLUX) +SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_NU_TILDE, LIFT, DRAG, LINSOL) From 7e024ed3a14f7aba76f2f0493f48e8681ad4d484 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 24 Dec 2025 18:27:47 -0800 Subject: [PATCH 5/6] authors and restore some comments --- AUTHORS.md | 1 + Common/src/linear_algebra/CSysSolve.cpp | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/AUTHORS.md b/AUTHORS.md index 4c33c55a5773..7f8afb021de6 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -71,6 +71,7 @@ Eduardo Molina Edwin van der Weide Eitan Aberman Ethan Alan Hereth +Ezgi Orbay Akcengiz Florian Dittmann Filip Hahs Francesco Poli diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 69702c51507a..889f5aa98658 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -455,7 +455,12 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector Date: Wed, 24 Dec 2025 21:19:42 -0800 Subject: [PATCH 6/6] update tests --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 4682755364ab..fc6ee002821f 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -404,7 +404,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788405, -11.040877, 0.000008, 0.309505] + inc_turb_naca0012.test_vals = [-4.758062, -10.974496, -0.000005, -0.028654, 4, -5.397415, 2, -6.426845] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index c5a3a055ad5a..9227d447d597 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -647,7 +647,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788595, -11.040942, -0.000002, 0.309519] + inc_turb_naca0012.test_vals = [-4.758062, -10.974496, -0.000006, -0.028654, 3, -5.740142, 2, -6.615162] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 1d06f1aa8d5d..cfd25784fa03 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -442,7 +442,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788495, -11.040895, 0.000023, 0.309502] + inc_turb_naca0012.test_vals = [-4.758063, -10.974497, -0.000005, -0.028655, 4, -6.041638, 2, -6.473951] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST