From 1577281b4023c6d4e51a68b7a1d744c6a0ffec83 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 30 Dec 2025 08:55:36 -0800 Subject: [PATCH 1/2] add smoother to the options --- Common/include/linear_algebra/CSysSolve.hpp | 6 ++ Common/include/option_structure.hpp | 2 + Common/src/linear_algebra/CSysSolve.cpp | 87 +++++++++++-------- .../integration/CNewtonIntegration.hpp | 14 ++- config_template.cfg | 6 +- 5 files changed, 75 insertions(+), 40 deletions(-) diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 6d9fa44b38f..998daffac36 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/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. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 8375668778c..ed249c2d03e 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2389,10 +2389,12 @@ static const MapType 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 Inner_Linear_Solver_Map = { MakePair("NONE", LINEAR_SOLVER_INNER::NONE) MakePair("BCGSTAB", LINEAR_SOLVER_INNER::BCGSTAB) + MakePair("SMOOTHER", LINEAR_SOLVER_INNER::SMOOTHER) }; diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index c13df2d36a5..66977e9fcbc 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -312,8 +312,8 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector& /*--- 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 ---*/ @@ -1032,8 +1032,8 @@ unsigned long CSysSolve::BCGSTAB_LinSolver(const CSysVectorGetComm_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 +1205,8 @@ unsigned long CSysSolve::Smoother_LinSolver(const CSysVectorGetComm_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,7 +1246,9 @@ unsigned long CSysSolve::Smoother_LinSolver(const CSysVector::Smoother_LinSolver(const CSysVector1 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 +1286,33 @@ unsigned long CSysSolve::Smoother_LinSolver(const CSysVector +bool CSysSolve::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>(LINEAR_SOLVER_MODE::STANDARD); + inner_solver->SetxIsZero(true); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + return true; + } + return false; +} + template unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, const CSysVector& LinSysRes, CSysVector& LinSysSol, CGeometry* geometry, @@ -1334,16 +1365,7 @@ 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 - } + const bool nested = SetupInnerSolver(KindSolver, config); /*--- Stop the recording for the linear solver ---*/ bool TapeActive = NO; @@ -1389,11 +1411,15 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con auto f = [&](const CSysVector& u, CSysVector& 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(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(f); } @@ -1545,16 +1571,7 @@ 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 - } + const bool nested = SetupInnerSolver(KindSolver, config); /*--- Set up preconditioner and matrix-vector product ---*/ @@ -1572,13 +1589,15 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c CPreconditioner* nested_prec = nullptr; if (nested) { auto f = [&](const CSysVector& u, CSysVector& 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(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(f); } diff --git a/SU2_CFD/include/integration/CNewtonIntegration.hpp b/SU2_CFD/include/integration/CNewtonIntegration.hpp index c51cd4b29b9..52a6190d6b2 100644 --- a/SU2_CFD/include/integration/CNewtonIntegration.hpp +++ b/SU2_CFD/include/integration/CNewtonIntegration.hpp @@ -159,10 +159,16 @@ class CNewtonIntegration final : public CIntegration { auto product = CSysMatrixVectorProduct(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 (config->GetKind_Linear_Solver_Inner()) { + 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; diff --git a/config_template.cfg b/config_template.cfg index 04f8aa6d905..375d9f1cecb 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1553,8 +1553,10 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % BCGSTAB, FGMRES, FGCRODR, RESTARTED_FGMRES, CONJUGATE_GRADIENT (self-adjoint problems only), SMOOTHER. LINEAR_SOLVER= FGMRES -% Inner solver used for nested preconditioning (only possible with flexible linear solvers). -% Options: NONE, BCGSTAB +% Inner solver used for nested preconditioning (BCGSTAB is only possible with flexible linear solvers). +% Options: NONE, BCGSTAB, SMOOTHER +% Inner BCGSTAB does up to the number of iterations of the outer solver or sqrt(LINEAR_SOLVER_ERROR). +% Inner SMOOTHER does max(2, sqrt(LINEAR_SOLVER_ITER)) LINEAR_SOLVER_INNER= NONE % Same for discrete adjoint (smoothers not supported), replaces LINEAR_SOLVER in SU2_*_AD codes. From 1a781eaf0034791ac7e5a138042c04d191db3231 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 31 Dec 2025 12:16:40 -0800 Subject: [PATCH 2/2] fix --- SU2_CFD/include/integration/CNewtonIntegration.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/integration/CNewtonIntegration.hpp b/SU2_CFD/include/integration/CNewtonIntegration.hpp index 52a6190d6b2..7d59d30259a 100644 --- a/SU2_CFD/include/integration/CNewtonIntegration.hpp +++ b/SU2_CFD/include/integration/CNewtonIntegration.hpp @@ -152,14 +152,16 @@ class CNewtonIntegration final : public CIntegration { template::value> = 0> inline unsigned long Preconditioner_impl(const CSysVector& u, CSysVector& 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(solvers[FLOW_SOL]->Jacobian, geometry, config); v = MixedScalar(0.0); MixedScalar eps_t = eps; - switch (config->GetKind_Linear_Solver_Inner()) { + 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;