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/include/CConfig.hpp b/Common/include/CConfig.hpp
index 8abcc77e4fbe..48aecba36956 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -515,6 +515,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. */
@@ -4288,6 +4289,10 @@ class CConfig {
    */
   unsigned short GetKind_Linear_Solver(void) const { return Kind_Linear_Solver; }
 
+  /*!
+   * \brief Get the inner linear solver used in nested Krylov linear solvers.
+   */
+  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 <functional>
 #include "../CConfig.hpp"
 #include "../geometry/CGeometry.hpp"
 #include "CSysVector.hpp"
@@ -305,6 +306,25 @@ class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
   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 ScalarType>
+class CAbstractPreconditioner final : public CPreconditioner<ScalarType> {
+ private:
+  std::function<void(const CSysVector<ScalarType>&, CSysVector<ScalarType>&)> impl;
+
+ public:
+  CAbstractPreconditioner() = delete;
+
+  template <class F>
+  explicit CAbstractPreconditioner(const F& function) : impl(function) {}
+
+  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { impl(u, v); }
+};
+
 template <class ScalarType>
 CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOLVER_PREC kind,
                                                                  CSysMatrix<ScalarType>& 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 <cmath>
+#include <memory>
 #include <vector>
 #include <iostream>
 #include <cstdlib>
@@ -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<CSysSolve<ScalarType>> 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<std::string, ENUM_LINEAR_SOLVER> 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<std::string, LINEAR_SOLVER_INNER> 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 8a84e679ab69..4d1a262d4b1f 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1887,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 */
@@ -1907,7 +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 */
 
   /* DESCRIPTION: Maximum update ratio value for flow density and energy variables */
   addDoubleOption("MAX_UPDATE_FLOW", MaxUpdateFlow, 0.2);
@@ -7260,10 +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
-                cout << "FGMRES 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;
+                }
+              }
               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;
@@ -7310,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 1a305c123dd3..889f5aa98658 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 <limits>
+#include <memory>
 
 /*!
  * \brief Epsilon used in CSysSolve depending on datatype to
@@ -110,7 +111,7 @@ void CSysSolve<ScalarType>::SolveReduced(int n, const su2matrix<ScalarType>& Hsb
 
 template <class ScalarType>
 void CSysSolve<ScalarType>::ModGramSchmidt(bool shared_hsbg, int i, su2matrix<ScalarType>& Hsbg,
-                                           vector<CSysVector<ScalarType> >& w) const {
+                                           vector<CSysVector<ScalarType>>& w) const {
   const auto thread = omp_get_thread_num();
 
   /*--- If Hsbg is shared by multiple threads calling this function, only one
@@ -902,9 +903,8 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& 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();
@@ -915,6 +915,17 @@ 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
+  }
+
   /*--- Stop the recording for the linear solver ---*/
   bool TapeActive = NO;
 
@@ -948,13 +959,25 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
 
     auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
 
-    const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
-
-    auto precond = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
-
     /*--- Build preconditioner. ---*/
 
-    precond->Build();
+    const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+    auto* normal_prec = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
+    normal_prec->Build();
+
+    CPreconditioner<ScalarType>* nested_prec = nullptr;
+    if (nested) {
+      nested_prec = new CAbstractPreconditioner<ScalarType>([&](const CSysVector<ScalarType>& u,
+                                                                CSysVector<ScalarType>& 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. ---*/
 
@@ -1000,7 +1023,8 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& 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(?). ---*/
@@ -1085,9 +1109,8 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& 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();
@@ -1098,19 +1121,43 @@ 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
+  }
+
   /*--- Set up preconditioner and matrix-vector product ---*/
 
-  const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+  auto mat_vec = CSysMatrixVectorProduct<ScalarType>(Jacobian, geometry, config);
 
-  auto precond = CPreconditioner<ScalarType>::Create(kindPrec, Jacobian, geometry, config);
+  const auto kindPrec = static_cast<ENUM_LINEAR_SOLVER_PREC>(KindPrecond);
+  auto* normal_prec = CPreconditioner<ScalarType>::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<ScalarType>(Jacobian, geometry, config);
+  CPreconditioner<ScalarType>* nested_prec = nullptr;
+  if (nested) {
+    nested_prec =
+        new CAbstractPreconditioner<ScalarType>([&](const CSysVector<ScalarType>& u, CSysVector<ScalarType>& 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 ---*/
 
@@ -1154,7 +1201,8 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
 
   HandleTemporariesOut(LinSysSol);
 
-  delete precond;
+  delete normal_prec;
+  delete nested_prec;
 
   SU2_OMP_MASTER {
     Residual = residual;
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<MixedScalar>(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;
   }
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/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)
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
old mode 100644
new mode 100755
index c5a3a055ad5a..9227d447d597
--- 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
old mode 100644
new mode 100755
index 1d06f1aa8d5d..cfd25784fa03
--- 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
diff --git a/config_template.cfg b/config_template.cfg
index aa7a1a97fc97..93dd76e84638 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1543,7 +1543,11 @@ 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
-%
+
+% 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
 %