diff --git a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh index f87547532e537daaaab63a69d33754913ccd4eee..9792b3442c1a33f0d7e80edb2db8bb99acbc3c36 100644 --- a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh +++ b/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; diff --git a/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh b/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh new file mode 100644 index 0000000000000000000000000000000000000000..af1324b37d6ba965959ae3a095ce151c0126fa7e --- /dev/null +++ b/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh @@ -0,0 +1,238 @@ +#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 diff --git a/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh b/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh new file mode 100644 index 0000000000000000000000000000000000000000..78a22d0c60f1e2e24cef2a597757176e1fec2e03 --- /dev/null +++ b/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh @@ -0,0 +1,299 @@ +#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>()); + + + + + + + + + + + + + diff --git a/src/example_2p_fs/2p_driver.hh b/src/example_2p_fs/2p_driver.hh index 25308abc858ba250e70f2ab1e68fd1e2aee32d63..2b0bf7c9efe1dddb43af9ec6b88826a173fe4fcb 100644 --- a/src/example_2p_fs/2p_driver.hh +++ b/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); - - //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); - - // 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 ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES - <NS_GO, Dune::SeqILU, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo, 1000, 200, 3, 1); + // 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 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 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_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;