Added draft of Uzawa preconditioner

dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh

@@ -38,6 +38,9 @@ public:
   {
     using MatrixType = Dune::PDELab::Backend::Native<M>;
     MatrixType& mat = Dune::PDELab::Backend::native(A);
+//     Dune::writeMatrixToMatlab(mat, "matrix-M");
+//     auto& matA = mat[Dune::Indices::_0][Dune::Indices::_0];
+//     Dune::writeMatrixToMatlab(matA, "matrix-A");
     using VectorType = Dune::PDELab::Backend::Native<W>;
 #if HAVE_MPI
     typedef typename Dune::PDELab::ISTL::CommSelector<96,Dune::MPIHelper::isFake>::type Comm;

dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh

+#pragma once
+
+#include<dune/pdelab/backend/istl.hh>
+#include "SEQ_Uzawa_preconditioner.hh"
+
+template<template<class> class Solver>
+class ISTLBackend_SEQ_Richardson
+{
+public:
+  /*! \brief make a linear solver object
+      
+      \param[in] maxiter_ maximum number of iterations to do
+      \param[in] verbose_ print messages if true
+  */
+  explicit ISTLBackend_SEQ_Richardson(unsigned maxiter_=5000, int verbose_=1)
+    : maxiter(maxiter_), verbose(verbose_)
+  {}
+  
+  /*! \brief compute global norm of a vector
+
+      \param[in] v the given vector
+  */
+  template<class V>
+  typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(const V& v) const
+  {
+    return v.two_norm();
+  }
+
+  /*! \brief solve the given linear system
+
+      \param[in] A the given matrix
+      \param[out] z the solution vector to be computed
+      \param[in] r right hand side
+      \param[in] reduction to be achieved
+  */
+  template<class M, class V, class W>
+  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
+  {
+    using Dune::PDELab::Backend::Native;
+    using Dune::PDELab::Backend::native;
+
+    Dune::MatrixAdapter<Native<M>,
+                        Native<V>,
+                        Native<W>> opa(native(A));
+        
+//     using MatrixType = Native<M>;
+//     MatrixType& mat = native(A);
+//     Dune::writeMatrixToMatlab(mat, "matrix-M");
+//     auto& matA = mat[Dune::Indices::_0][Dune::Indices::_0];                    
+//     Dune::writeMatrixToMatlab(matA, "matrix-A");
+//     auto& matB = mat[Dune::Indices::_0][Dune::Indices::_1];                    
+//     Dune::writeMatrixToMatlab(matB, "matrix-B");
+//     auto& matC = mat[Dune::Indices::_1][Dune::Indices::_0];                    
+//     Dune::writeMatrixToMatlab(matC, "matrix-C");
+//     auto& matD = mat[Dune::Indices::_1][Dune::Indices::_1];                    
+//     Dune::writeMatrixToMatlab(matD, "matrix-D");
+        
+    Dune::Richardson<Native<V>,Native<W>> prec(0.7);
+    Solver<Native<V>> solver(opa, prec, reduction, maxiter, verbose);
+    
+    Dune::InverseOperatorResult stat;
+    solver.apply(native(z), native(r), stat);
+    res.converged  = stat.converged;
+    res.iterations = stat.iterations;
+    res.elapsed    = stat.elapsed;
+    res.reduction  = stat.reduction;
+    res.conv_rate  = stat.conv_rate;
+  }
+  
+  const Dune::PDELab::LinearSolverResult<double>& result() const
+  {
+    return res;
+  }
+
+private:
+  unsigned maxiter;
+  int verbose;
+  Dune::PDELab::LinearSolverResult<double> res;
+};
+    
+ 
+/**
+  * @brief Backend for sequential BiCGSTAB solver with Richardson
+  * preconditioner (equivalent to no preconditioner for Richards with
+  * parameter=1.0)
+  */
+class ISTLBackend_SEQ_BCGS_Richardson
+  : public ISTLBackend_SEQ_Richardson<Dune::BiCGSTABSolver>
+{
+public:
+  /*! \brief make a linear solver object
+      
+      \param[in] maxiter_ maximum number of iterations to do
+      \param[in] verbose_ print messages if true
+  */
+  explicit ISTLBackend_SEQ_BCGS_Richardson (unsigned maxiter_=5000, int verbose_=1)
+    : ISTLBackend_SEQ_Richardson<Dune::BiCGSTABSolver>(maxiter_, verbose_)
+  {}
+};
+
+/** \brief Linear solver backend for Restarted GMRes solver preconditioned with the Richardson preconditioner
+  * (equivalent to no preconditioner for Richards with parameter=1.0)
+  */
+class ISTLBackend_SEQ_GMRES_Richardson
+{
+public :
+  /*! \brief make a linear solver object
+
+      \param[in] restart_ number of iterations when GMRes has to be restarted
+      \param[in] maxiter_ maximum number of iterations to do
+      \param[in] verbose_ print messages if true
+  */
+  explicit ISTLBackend_SEQ_GMRES_Richardson(int restart_ = 200, int maxiter_ = 5000, int verbose_ = 1)
+    : restart(restart_), maxiter(maxiter_), verbose(verbose_)
+  {}
+  
+  /*! \brief compute global norm of a vector
+
+      \param[in] v the given vector
+  */
+  template<class V>
+  typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(const V& v) const
+  {
+    return v.two_norm();
+  }
+
+  /** \brief solve the given linear system
+
+      \param[in] A the given matrix
+      \param[out] z the solution vector to be computed
+      \param[in] r right hand side
+      \param[in] reduction to be achieved
+  */
+  template<class M, class V, class W>
+  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType>::real_type reduction)
+  {
+    using Dune::PDELab::Backend::Native;
+    using Dune::PDELab::Backend::native;
+    
+    Dune::MatrixAdapter<Native<M>,
+                        Native<V>,
+                        Native<W>> opa(native(A));
+                        
+    Dune::Richardson<Native<V>,Native<W>> prec(0.7);
+    Dune::RestartedGMResSolver<Native<V>> solver(opa, prec, reduction, restart, maxiter, verbose);
+    
+    Dune::InverseOperatorResult stat;
+    solver.apply(native(z), native(r), stat);
+    res.converged  = stat.converged;
+    res.iterations = stat.iterations;
+    res.elapsed    = stat.elapsed;
+    res.reduction  = stat.reduction;
+    res.conv_rate  = stat.conv_rate;
+  }
+  
+  const Dune::PDELab::LinearSolverResult<double>& result() const
+  {
+    return res;
+  }  
+
+private :
+  int restart, maxiter, verbose;
+  Dune::PDELab::LinearSolverResult<double> res;
+};
+
+/** \brief Linear solver backend for Restarted GMRes solver preconditioned with the Uzawa preconditioner
+  */
+class ISTLBackend_SEQ_GMRES_Uzawa
+{
+public :
+  /*! \brief make a linear solver object
+
+      \param[in] restart_ number of iterations when GMRes has to be restarted
+      \param[in] maxiter_ maximum number of iterations to do
+      \param[in] verbose_ print messages if true
+  */
+  explicit ISTLBackend_SEQ_GMRES_Uzawa(int restart_ = 200, int maxiter_ = 5000, int verbose_ = 1, const Dune::ParameterTree& params_ = {})
+    : restart(restart_), maxiter(maxiter_), verbose(verbose_), params(params_)
+  {}
+  
+  /*! \brief compute global norm of a vector
+
+      \param[in] v the given vector
+  */
+  template<class V>
+  typename Dune::template FieldTraits<typename V::ElementType >::real_type norm(const V& v) const
+  {
+    return v.two_norm();
+  }
+
+  /** \brief solve the given linear system
+
+      \param[in] A the given matrix
+      \param[out] z the solution vector to be computed
+      \param[in] r right hand side
+      \param[in] reduction to be achieved
+  */
+  template<class M, class V, class W>
+  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType>::real_type reduction)
+  {
+    using Dune::PDELab::Backend::Native;
+    using Dune::PDELab::Backend::native;
+    
+    Dune::MatrixAdapter<Native<M>,
+                        Native<V>,
+                        Native<W>> opa(native(A));
+                        
+    SeqUzawa<Native<M>,Native<V>,Native<W>> prec(native(A), params);
+    Dune::RestartedGMResSolver<Native<V>> solver(opa, prec, reduction, restart, maxiter, verbose);
+    
+    Dune::InverseOperatorResult stat;
+    solver.apply(native(z), native(r), stat);
+    res.converged  = stat.converged;
+    res.iterations = stat.iterations;
+    res.elapsed    = stat.elapsed;
+    res.reduction  = stat.reduction;
+    res.conv_rate  = stat.conv_rate;
+  }
+  
+  const Dune::PDELab::LinearSolverResult<double>& result() const
+  {
+    return res;
+  }  
+
+private :
+  int restart, maxiter, verbose;
+  Dune::PDELab::LinearSolverResult<double> res;
+  const Dune::ParameterTree& params;
+};
+
+
+// create Uzawa-, SIMPLE(R)-type preconditioners (istl)
+
+// create CH block preconditioner
+
+// create sequential GMRes solver with above preconditioners
+
+// test on partitioned case

dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh

+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/version.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/paamg/amg.hh>
+
+/*!
+ * \brief A preconditioner based on the Uzawa algorithm for saddle-point problems of the form
+ * \f$
+ \begin{pmatrix}
+    A & B \\
+    C & D
+ \end{pmatrix}
+
+ \begin{pmatrix}
+    u\\
+    p
+ \end{pmatrix}
+
+ =
+
+ \begin{pmatrix}
+    f\\
+    g
+ \end{pmatrix}
+ * \f$
+ *
+ * This preconditioner is especially suited for solving the incompressible (Navier-)Stokes equations.
+ * Here, \f$D = 0\f$ and \f$B = C^T\f$ if \f$\rho = 1\f$.
+ * We do not expect good convergence if energy or mass transport is considered.
+ *
+ * See: Benzi, M., Golub, G. H., & Liesen, J. (2005). Numerical solution of saddle point problems. Acta numerica, 14, 1-137
+ *      Ho, N., Olson, S. D., & Walker, H. F. (2017). Accelerating the Uzawa algorithm. SIAM Journal on Scientific Computing, 39(5), S461-S476
+ *
+ * \tparam M Type of the matrix.
+ * \tparam X Type of the update.
+ * \tparam Y Type of the defect.
+ * \tparam l Preconditioner block level (for compatibility reasons, unused).
+ */
+template<class M, class X, class Y, int l = 1>
+class SeqUzawa : public Dune::Preconditioner<X,Y>
+{
+  typedef std::decay_t<decltype(std::declval<M>()[Dune::Indices::_0][Dune::Indices::_0])> A;
+  typedef std::decay_t<decltype(std::declval<X>()[Dune::Indices::_0])> U;
+
+  typedef Dune::Amg::SequentialInformation Comm;
+  typedef Dune::MatrixAdapter<A, U, U> LinearOperator;
+  typedef Dune::SeqSSOR<A, U, U> Smoother;
+  typedef Dune::Amg::AMG<LinearOperator, U, Smoother, Comm> AMGSolverForA;
+  
+public:
+  //! \brief The matrix type the preconditioner is for.
+  typedef M matrix_type;
+  //! \brief The domain type of the preconditioner.
+  typedef X domain_type;
+  //! \brief The range type of the preconditioner.
+  typedef Y range_type;
+  //! \brief The field type of the preconditioner.
+  typedef typename X::field_type field_type;
+  //! \brief scalar type underlying the field_type
+  typedef Dune::Simd::Scalar<field_type> scalar_field_type;
+  //! \brief real scalar type underlying the field_type
+  typedef typename Dune::FieldTraits<scalar_field_type>::real_type real_field_type;
+  
+  /*!
+   * \brief Constructor
+   *
+   * \param mat The matrix to operate on.
+   * \param params Collection of paramters.
+   */
+  SeqUzawa(const M& op, const Dune::ParameterTree& params)
+    : matrix_(op)
+    , numIterations_(params.get<std::size_t>("NumIterations", 1))
+    , relaxationFactor_(params.get<scalar_field_type>("Relaxation", 1.0))
+    , verbosity_(params.get<int>("VerbosityLevel", 1))
+    , useDirectVelocitySolverForA_(params.get<bool>("DirectSolverForA", false))
+  {
+    const bool determineRelaxationFactor = params.get<bool>("DetermineRelaxationFactor", true);
+
+    // AMG is needed for determination of omega
+    if (determineRelaxationFactor || !useDirectVelocitySolverForA_)
+      initAMG_(params);
+
+    if (useDirectVelocitySolverForA_)
+      initUMFPack_();
+
+    if (determineRelaxationFactor)
+      relaxationFactor_ = estimateOmega_(params);
+  }
+  
+  /*!
+   * \brief Prepare the preconditioner.
+   */
+  virtual void pre(X& x, Y& b) {}
+
+  /*!
+   * \brief Apply the preconditioner
+   *
+   * \param update The update to be computed.
+   * \param currentDefect The current defect.
+   */
+  virtual void apply(X& update, const Y& currentDefect)
+  {
+    using namespace Dune::Indices;
+
+    auto& A = matrix_[_0][_0];
+    auto& B = matrix_[_0][_1];
+    auto& C = matrix_[_1][_0];
+    auto& D = matrix_[_1][_1];
+
+    const auto& f = currentDefect[_0];
+    const auto& g = currentDefect[_1];
+    auto& u = update[_0];
+    auto& p = update[_1];
+
+    // incorporate Dirichlet cell values
+    // TODO: pass Dirichlet constraint handler from outside
+    for (std::size_t i = 0; i < D.N(); ++i)
+    {
+      const auto& block = D[i][i];
+      for (auto rowIt = block.begin(); rowIt != block.end(); ++rowIt)
+        if (Dune::FloatCmp::eq<scalar_field_type>(rowIt->one_norm(), 1.0))
+          p[i][rowIt.index()] = g[i][rowIt.index()];
+    }
+    
+    // the actual Uzawa iteration
+    for (std::size_t k = 0; k < numIterations_; ++k)
+    {
+      // u_k+1 = u_k + Q_A^−1 * (f − (A*u_k + B*p_k))
+      auto uRhs = f;
+      A.mmv(u, uRhs);
+      B.mmv(p, uRhs);
+      auto uIncrement = u;
+      applySolverForA_(uIncrement, uRhs);
+      u += uIncrement;
+
+      // p_k+1 = p_k + omega * (g - C*u_k+1 - D*p_k)
+      auto pIncrement = g;
+      C.mmv(u, pIncrement);
+      D.mmv(p, pIncrement);
+      pIncrement *= relaxationFactor_;
+      p += pIncrement;
+      
+      if (verbosity_ > 1)
+      {
+        std::cout << "Uzawa iteration " << k << ", residual: " << uRhs.two_norm() + pIncrement.two_norm()/relaxationFactor_ << std::endl;
+      }
+    }
+  }
+
+  /*!
+   * \brief Clean up.
+   */
+  virtual void post(X& x) {}
+  
+  //! Category of the preconditioner (see SolverCategory::Category)
+  virtual Dune::SolverCategory::Category category() const
+  {
+    return Dune::SolverCategory::sequential;
+  }
+  
+private:
+  void initAMG_(const Dune::ParameterTree& params)
+  {
+    using namespace Dune::Indices;
+    auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
+    amgSolverForA_ = std::make_unique<AMGSolverForA>(linearOperator, params);
+  }
+
+  void initUMFPack_()
+  {
+#if HAVE_UMFPACK
+    using namespace Dune::Indices;
+    umfPackSolverForA_ = std::make_unique<Dune::UMFPack<A>>(matrix_[_0][_0]);
+#else
+    DUNE_THROW(Dune::InvalidStateException, "UMFPack not available. Use DirectSolverForA = false.");
+#endif
+  }
+  
+  /*!
+   * \brief Estimate the relaxation factor omega
+   *
+   * The optimal relaxation factor is omega = 2/(lambdaMin + lambdaMax), where lambdaMin and lambdaMax are the smallest and largest
+   * eigenvalues of the Schur complement -C*Ainv*B (assuming D = 0).
+   * lambdaMax can be easily determined using the power iteration algorithm (https://en.wikipedia.org/wiki/Power_iteration) and lambdaMin
+   * could be estimated in a similar manner. We do not consider lambdaMin because for certain cases, e.g., when C contains some rows of zeroes only,
+   * this estimate will fail.
+   *
+   * Instead we assume that lambdaMin is sufficiently close to lambdaMax such that omega = 1/lambdaMax.
+   * This seems to work rather well for various applications.
+   * We will underestimate omega by a factor of 2 in the worst case (i.e, lambdaMin = 0).
+   *
+   * When facing convergence issues, you may set LinearSolver.Preconditioner.Verbosity = 1 to see the estimate of lambdaMax.
+   * In a new simulation run, you can then set LinearSolver.Preconditioner.DetermineRelaxationFactor = false and set some other value
+   * for LinearSolver.Preconditioner.Relaxation based on the estimate of lambdaMax.
+   *
+   * See: Benzi, M., Golub, G. H., & Liesen, J. (2005). Numerical solution of saddle point problems. Acta numerica, 14, 1-137.
+   */
+  scalar_field_type estimateOmega_(const Dune::ParameterTree& params) const
+  {
+    using namespace Dune::Indices;
+    auto& A = matrix_[_0][_0];
+    auto& B = matrix_[_0][_1];
+    auto& C = matrix_[_1][_0];
+
+    U x(A.M());
+    x = 1.0;
+
+    scalar_field_type omega = 0.0;
+    scalar_field_type lambdaMax = 0.0;
+
+    static const auto iterations = params.get<std::size_t>("PowerLawIterations", 5);
+
+    // apply power iteration x_k+1 = M*x_k/|M*x_k| for the matrix M = -C*Ainv*B
+    for (std::size_t i = 0; i < iterations; ++i)
+    {
+      // bx = B*x
+      U bx(x.size());
+      B.mv(x, bx);
+
+      // ainvbx = Ainv*(B*x)
+      auto ainvbx = x;
+      applySolverForA_(ainvbx, bx);
+
+      // v = M*x = -C*(Ainv*B*x)
+      U v(x.size());
+      C.mv(ainvbx, v);
+      v *= -1.0;
+
+      // eigenvalue lambdaMax = xt*M*x/(xt*x) = xt*v/(xt*x);
+      lambdaMax = x.dot(v)/(x.dot(x));
+
+      // relaxation factor omega = 1/lambda;
+      omega = 1.0/lambdaMax;
+
+      // new iterate x = M*x/|M*x| = v/|v|
+      x = v;
+      x /= v.two_norm();
+    }
+
+    if (verbosity_ > 1)
+    {
+      std::cout << "\n*** Uzawa Preconditioner ***" << std::endl;
+      std::cout << "Estimating relaxation factor based on Schur complement" << std::endl;
+      std::cout << "Largest estimated eigenvalue lambdaMax = " << lambdaMax << std::endl;
+      std::cout << "Relaxation factor omega = 1/lambdaMax = " << omega << std::endl;
+    }
+
+    return omega;
+  }
+  
+  template<class Sol, class Rhs>
+  void applySolverForA_(Sol& sol, Rhs& rhs) const
+  {
+    if (useDirectVelocitySolverForA_)
+    {
+#if HAVE_UMFPACK
+      Dune::InverseOperatorResult res;
+      umfPackSolverForA_->apply(sol, rhs, res);
+#endif
+    }
+    else
+    {
+      amgSolverForA_->pre(sol, rhs);
+      amgSolverForA_->apply(sol, rhs);
+      amgSolverForA_->post(sol);
+    }
+  }
+  
+  //! \brief The matrix we operate on.
+  const M& matrix_;
+  //! \brief The number of steps to do in apply
+  const std::size_t numIterations_;
+  //! \brief The relaxation factor to use
+  scalar_field_type relaxationFactor_;
+  //! \brief The verbosity level
+  const int verbosity_;
+  
+  std::unique_ptr<AMGSolverForA> amgSolverForA_;
+#if HAVE_UMFPACK
+  std::unique_ptr<Dune::UMFPack<A>> umfPackSolverForA_;
+#endif
+  const bool useDirectVelocitySolverForA_;
+};
+// DUNE_REGISTER_PRECONDITIONER("uzawa", defaultPreconditionerBlockLevelCreator<SeqUzawa>());
+
+
+
+
+
+
+
+
+
+
+
+  
+

src/example_2p_fs/2p_driver.hh

@@ -14,6 +14,8 @@
 #include <dune/phasefield/localoperator/pf_diff_estimator.hh>
 
 #include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+#include <dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh>
+// #include <dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh>
 
 #include <dune/phasefield/chns_base.hh>
 #include <dune/phasefield/ff_chns.hh>
@@ -42,7 +44,8 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem,
   //_________________________________________________________________________________________________________________________
   // Make grid function spaces
 
-  typedef GFS_ff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+//   typedef GFS_ff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+  typedef GFS_ff_chns_blocked<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
   GFS gfs(es, vFem, pFem, phiFem);
   typedef typename GFS::NS NS_GFS;
   typedef typename GFS::CH CH_GFS;
@@ -156,16 +159,26 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem,
       typedef Dune::PDELab::GridOperator<CH_GFS, CH_GFS, CH_LOP, MBE, RF, RF, RF, typename GFS::CH_CC, typename GFS::CH_CC> CH_GO;
       CH_GO chGo(gfs.ch, gfs.chCC, gfs.ch, gfs.chCC, chLop, mbe);
 
+      // Linear solvers
+      //typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_GO, Dune::SeqILU, Dune::RestartedGMResSolver> CH_LS; 
+      //CH_LS ls_ch(chGo, 600, 200, 3, 1);
 
-      //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 CH_LS; CH_LS ls_ch(200, 600, 1);
-      //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 NS_LS; NS_LS ls_ns(200, 1000, 1);
+      //typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_GO, Dune::SeqILU, Dune::RestartedGMResSolver> NS_LS; 
+      //NS_LS ls_ns(nsGo, 1000, 200, 3, 1);
       
-      // Linear solver
-      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-                  <CH_GO, Dune::SeqILU, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo, 600, 200, 3, 1);
+      typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 CH_LS; 
+      CH_LS ls_ch(200, 600, 1);
       
-      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-                  <NS_GO, Dune::SeqILU, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo, 1000, 200, 3, 1);
+//       typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 NS_LS; 
+//       NS_LS ls_ns(200, 1000, 1);
+      
+//       typedef ISTLBackend_SEQ_BCGS_Richardson CH_LS; CH_LS ls_ch(1000, 1);
+
+//       typedef ISTLBackend_SEQ_BCGS_Richardson NS_LS; NS_LS ls_ns(1000, 1);
+      
+//       typedef ISTLBackend_SEQ_GMRES_Richardson NS_LS; NS_LS ls_ns(200, 10000, 1);
+      
+      typedef ISTLBackend_SEQ_GMRES_Uzawa NS_LS; NS_LS ls_ns(200, 10000, 1, ptree.sub("LinearSolver_Preconditioner"));
 
       // Nonlinear solver
       typedef Dune::PDELab::NewtonMethod<NS_GO, NS_LS> PDESolver1;