diff --git a/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh b/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh
index af1324b37d6ba965959ae3a095ce151c0126fa7e..509fb9cab1f1f3c339370d3b2e4a88e6e50bb05a 100644
--- a/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh
+++ b/dune/phasefield/utility/ISTLBackend_SEQ_Solvers.hh
@@ -2,6 +2,7 @@
#include<dune/pdelab/backend/istl.hh>
#include "SEQ_Uzawa_preconditioner.hh"
+#include "SEQ_SIMPLE_preconditioner.hh"
template<template<class> class Solver>
class ISTLBackend_SEQ_Richardson
@@ -228,8 +229,268 @@ private :
const Dune::ParameterTree& params;
};
+/** \brief Linear solver backend for Restarted GMRes solver preconditioned with the SIMPLE preconditioner
+ */
+class ISTLBackend_SEQ_GMRES_SIMPLE
+{
+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_SIMPLE(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));
+
+ SeqSimple<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;
+};
+
+/** \brief Linear solver backend for Restarted Flexible GMRes preconditioned with ILU(0)
+ */
+class ISTLBackend_SEQ_FGMRES_ILU0
+{
+public :
+ /** \brief make 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_FGMRES_ILU0(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::SeqILU<Native<M>,
+ Native<V>,
+ Native<W>> ilu0(native(A), 1.0);
+
+ Dune::RestartedFlexibleGMResSolver<Native<V>> solver(opa,ilu0,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 Flexible GMRes solver preconditioned with the Uzawa preconditioner
+ */
+class ISTLBackend_SEQ_FGMRES_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_FGMRES_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::RestartedFlexibleGMResSolver<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;
+};
+
+/** \brief Linear solver backend for Restarted Flexible GMRes solver preconditioned with the SIMPLE preconditioner
+ */
+class ISTLBackend_SEQ_FGMRES_SIMPLE
+{
+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_FGMRES_SIMPLE(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));
+
+ SeqSimple<Native<M>,Native<V>,Native<W>> prec(native(A), params);
+ Dune::RestartedFlexibleGMResSolver<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
diff --git a/dune/phasefield/utility/SEQ_SIMPLE_preconditioner.hh b/dune/phasefield/utility/SEQ_SIMPLE_preconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a1a7c2f502c8ff018f1aedcb8026d7abf97d886e
--- /dev/null
+++ b/dune/phasefield/utility/SEQ_SIMPLE_preconditioner.hh
@@ -0,0 +1,668 @@
+#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>
+#include <dune/istl/matrixindexset.hh>
+
+#if HAVE_UMFPACK
+#include <dune/istl/umfpack.hh>
+#endif
+
+/*!
+ * \brief A preconditioner based on the SIMPLE 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$
+ *
+ * See: https://www.cs.umd.edu/~elman/papers/tax.pdf
+ *
+ * \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 SeqSimple : 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 std::decay_t<decltype(std::declval<M>()[Dune::Indices::_1][Dune::Indices::_1])> D;
+ typedef std::decay_t<decltype(std::declval<X>()[Dune::Indices::_1])> P;
+
+ 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> AMGForA;
+
+ typedef Dune::Amg::SequentialInformation CommSchur;
+ typedef Dune::MatrixAdapter<D, P, P> LinearOperatorSchur;
+ typedef Dune::SeqSSOR<D, P, P> SmootherSchur;
+ typedef Dune::Amg::AMG<LinearOperatorSchur, P, SmootherSchur, CommSchur> AMGForSchur;
+
+ class SimpleOperator : public Dune::LinearOperator<P, P>
+ {
+ public:
+ SimpleOperator(const M& m) : m_(m) {}
+
+ template<class S>
+ void setSolver(const S& solver)
+ { solver_ = solver; }
+
+ /*! \brief apply operator to x: \f$ y = A(x) \f$
+ The input vector is consistent and the output must also be
+ consistent on the interior+border partition.
+ */
+ void apply (const P& x, P& y) const final
+ {
+ // convenience variables
+ const auto& deltaP = x;
+ auto& rhs = y;
+
+ using namespace Dune::Indices;
+ auto& B = m_[_0][_1];
+ auto& C = m_[_1][_0];
+ auto& D = m_[_1][_1];
+
+ // Apply the (approximate) Schur complement S * deltaP = (D - C * Q_A^-1 * B) * deltaP:
+
+ // B * deltaP
+ // rhsTmp has different size than rhs
+ auto rhsTmp = U(B.N());
+ B.mv(deltaP, rhsTmp);
+
+ // Q_A^-1 * B * deltaP
+ auto rhsTmp2 = rhsTmp;
+ Dune::InverseOperatorResult result;
+ solver_(rhsTmp, rhsTmp2, result);
+
+ // C * Q_A^-1 * B * deltaP
+ C.mv(rhsTmp, rhs);
+
+ // - C * Q_A^-1 * B * deltaP
+ rhs *= -1.0;
+
+ // (D - C * Q_A^-1 * B) * deltaP
+ D.umv(deltaP, rhs);
+ }
+
+ //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
+ void applyscaleadd (typename P::field_type alpha, const P& x, P& y) const final
+ {
+ auto tmp = P(y.size());
+ apply(x, tmp);
+ y.axpy(alpha, tmp);
+ }
+
+ //! Category of the linear operator (see SolverCategory::Category)
+ Dune::SolverCategory::Category category() const final
+ {
+ return Dune::SolverCategory::sequential;
+ };
+
+ private:
+ //! The system matrix
+ const M& m_;
+
+ //! Applies the effect of Q_A^-1 [A^-1 or diag(A)^-1 (depending on what is chosen in setSolver)]
+ std::function<void(U&, U&, Dune::InverseOperatorResult&)> solver_;
+ };
+
+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.
+ */
+ SeqSimple(const M& mat, const Dune::ParameterTree& params)
+ : matrix_(mat)
+ , params_(params)
+ , numIterations_(params.get<std::size_t>("NumIterations", 1))
+ , relaxationFactorP_(params.get<scalar_field_type>("Relaxation", 1.0))
+ , relaxationFactorU_(relaxationFactorP_)
+ , verbosity_(params.get<int>("VerbosityLevel", 1))
+ , useStandardSolverForA_(params.get<bool>("UseStandardSolverForA", false))
+ , useDirectSolverForA_(params.get<bool>("UseDirectSolverForA", false))
+ , useDiagonalInSchurComplement_(params.get<bool>("UseDiagonalInSchurComplement", false))
+ , useDiagonalInUpdate_(params.get<bool>("UseDiagonalInUpdate", false))
+ , useStandardSolverForSchurComplement_(params.get<bool>("UseStandardSolverForSchurComplement", false))
+ , useDirectSolverForSchurComplement_(params.get<bool>("UseDirectSolverForSchurComplement", false))
+ , useStandardSolverForAInSchurComplement_(params.get<bool>("UseStandardSolverForAInSchurComplement", false))
+ {
+ if(useStandardSolverForA_){
+ initAMGPreconditionerForA_(params);
+ initStandardSolverForA_();
+ if(!useStandardSolverForAInSchurComplement_){
+ if(!useDirectSolverForA_)
+ initAMGSolverForA_(params);
+ else
+ initUMFPackForA_();
+ }
+ }
+ else if(!useDirectSolverForA_)
+ initAMGSolverForA_(params);
+ else
+ initUMFPackForA_();
+
+ relaxationFactorP_ = params.get<scalar_field_type>("RelaxationP", relaxationFactorP_);
+ relaxationFactorU_ = params.get<scalar_field_type>("RelaxationU", relaxationFactorU_);
+
+ if(useStandardSolverForSchurComplement_){
+ initAMGPreconditionerForSchur_(params);
+ initStandardSchurComplementSolver_();
+ }
+ else if(!useDirectSolverForSchurComplement_)
+ initAMGSolverForSchur_(params);
+ else
+ initUMFPackForSchur_();
+ }
+
+ /*!
+ * \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];
+
+ static const bool dirichletHack = params_.get<bool>("DirichletHack", false);
+ // incorporate Dirichlet cell values
+ // TODO: pass Dirichlet constraint handler from outside
+ if (dirichletHack)
+ 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()];
+ }
+
+ static const std::string SIMPLEType = params_.get<std::string>("SIMPLEType", "SIMPLE");
+ if (SIMPLEType == "SIMPLE") {
+ // the SIMPLE algorithm
+ for (std::size_t k = 0; k < numIterations_; ++k)
+ {
+ // 1. Solve A*uNew + B*p = f for uNew
+ auto uRhs = f;
+ B.mmv(p, uRhs);
+ auto uNew = u;
+ if(useStandardSolverForA_){
+ Dune::InverseOperatorResult result;
+ solverForA_->apply(uNew, uRhs, result);
+ }
+ else
+ applySolverForA_(uNew, uRhs);
+
+ // 2. Solve S * deltaP = g - (C * uNew + D*p) for deltaP
+ auto pRhs = g;
+ C.mmv(uNew, pRhs);
+ D.mmv(p, pRhs);
+ auto deltaP = p;
+ if(useStandardSolverForSchurComplement_){
+ Dune::InverseOperatorResult result1;
+ schurComplementSolver_->apply(deltaP, pRhs, result1);
+ }
+ else
+ applySolverForSchurComplement_(deltaP, pRhs);
+
+ // 3. Calculate correction of u:
+ // deltaU = (-Q_A^-1 * B) * deltaP
+ auto deltaU = U(B.N());
+ B.mv(deltaP, deltaU);
+ if(useDiagonalInUpdate_)
+ applyInverseOfDiagonalOfA_(deltaU);
+ else if(useStandardSolverForA_){
+ auto deltaU2 = deltaU;
+ Dune::InverseOperatorResult result2;
+ solverForA_->apply(deltaU, deltaU2, result2);
+ }
+ else
+ applySolverForA_(deltaU, deltaU);
+
+ deltaU *= -1.0;
+
+ // 4. Update p
+ deltaP *= relaxationFactorP_;
+ p += deltaP;
+
+ // 5. Update u
+ deltaU += uNew;
+ deltaU *= relaxationFactorU_;
+ u *= (1 - relaxationFactorU_);
+ u += deltaU;
+ }
+ }
+ else if (SIMPLEType == "SIMPLER") {
+ // the SIMPLER algorithm
+ for (std::size_t k = 0; k < numIterations_; ++k)
+ {
+ // 1. Solve S * pNew = g - C*u - C * diag(A)^-1 * (f - A*u)
+ auto uRhs = f;
+ A.mmv(u, uRhs);
+
+ applyInverseOfDiagonalOfA_(uRhs);
+
+ auto pRhs = p;
+ C.mv(uRhs, pRhs);
+
+ C.umv(u, pRhs);
+ pRhs *= -1.0;
+ pRhs += g;
+
+ auto pNew = p;
+ if(useStandardSolverForSchurComplement_){
+ Dune::InverseOperatorResult result;
+ schurComplementSolver_->apply(pNew, pRhs, result);
+ }
+ else
+ applySolverForSchurComplement_(pNew, pRhs);
+
+ // 2. Update p
+ pNew *= relaxationFactorP_;
+ p *= (1 - relaxationFactorP_);
+ p += pNew;
+
+ // 3. Solve A*uNew + B*p = f for uNew
+ auto uRhs2 = f;
+ B.mmv(p, uRhs2);
+ auto uNew = u;
+ if(useStandardSolverForA_){
+ Dune::InverseOperatorResult result1;
+ solverForA_->apply(uNew, uRhs2, result1);
+ }
+ else
+ applySolverForA_(uNew, uRhs2);
+
+ // 4. Solve (D - C * diag(A)^-1 * B) * deltaP = g - (C * uNew + D*p) for deltaP
+ auto pRhs2 = g;
+ C.mmv(uNew, pRhs2);
+ D.mmv(p, pRhs2);
+ auto deltaP = p;
+ if(useStandardSolverForSchurComplement_){
+ Dune::InverseOperatorResult result2;
+ schurComplementSolver_->apply(deltaP, pRhs2, result2);
+ }
+ else
+ applySolverForSchurComplement_(deltaP, pRhs2);
+
+ // 5. Calculate correction of u:
+ // deltaU = (-diag(A^-1) * B) * deltaP
+ auto deltaU = U(B.N());
+ B.mv(deltaP, deltaU);
+
+ applyInverseOfDiagonalOfA_(deltaU);
+
+ deltaU *= -1.0;
+
+ // 6. Update u
+ deltaU += uNew;
+ deltaU *= relaxationFactorU_;
+ u *= (1 - relaxationFactorU_);
+ u += deltaU;
+ }
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, SIMPLEType << "not supported. Use SIMPLE, SIMPLER");
+ }
+
+ /*!
+ * \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 apply
+ void applyStandardSIMPLE_()
+ {}
+
+ void initAMGSolverForA_(const Dune::ParameterTree& params)
+ {
+ using namespace Dune::Indices;
+ auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
+ amgSolverForA_ = std::make_unique<AMGForA>(linearOperator, params);
+ }
+
+ void initUMFPackForA_()
+ {
+#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 UseDirectSolverForA = false.");
+#endif
+ }
+
+ template<class Sol, class Rhs>
+ void applySolverForA_(Sol& sol, Rhs& rhs) const
+ {
+ if (useDirectSolverForA_)
+ {
+#if HAVE_UMFPACK
+ Dune::InverseOperatorResult res;
+ auto rhsTmp = rhs;
+ umfPackSolverForA_->apply(sol, rhsTmp, res);
+#endif
+ }
+ else
+ {
+ amgSolverForA_->pre(sol, rhs);
+ amgSolverForA_->apply(sol, rhs);
+ amgSolverForA_->post(sol);
+ }
+ }
+
+ void initAMGPreconditionerForA_(const Dune::ParameterTree& params)
+ {
+ using namespace Dune::Indices;
+ auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
+ amgPreconditionerForA_ = std::make_shared<AMGForA>(linearOperator, params);
+ }
+
+ void initStandardSolverForA_()
+ {
+ using namespace Dune::Indices;
+ auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
+
+ Dune::ParameterTree treeSolverForA;
+ treeSolverForA["verbose"] = params_.get<std::string>("SolverForAVerbosity", "0");
+ treeSolverForA["reduction"] = params_.get<std::string>("SolverForAResidualReduction", "1e-12");
+ treeSolverForA["maxit"] = params_.get<std::string>("SolverForAMaxIter", "1000");
+ treeSolverForA["restart"] = params_.get<std::string>("SolverForARestart", "100");
+
+ static const bool usePreconditionerForA = params_.get<bool>("UsePreconditionerForA", false);
+ if (usePreconditionerForA)
+ {
+ Dune::ParameterTree treePreconditionerA;
+ treePreconditionerA["verbose"] = params_.get<std::string>("SolverForAVerbosity", "0");
+ // preconditionerForA_ = std::make_shared<SeqJacobiMatrixFree<U, U>>(linearOperator, treePreconditionerA);
+ preconditionerForA_ = amgPreconditionerForA_; // std::make_shared<Dune::SeqJac<A, U, U>>(matrix_[_0][_0], 1, 1.0);
+ }
+ else
+ preconditionerForA_ = std::make_shared<Dune::Richardson<U, U>>();
+
+ static const auto solverForAType = params_.get<std::string>("SolverForAType", "restartedgmressolver");
+
+ if (solverForAType == "restartedgmressolver")
+ solverForA_ = std::make_unique<Dune::RestartedGMResSolver<U>>(linearOperator, preconditionerForA_, treeSolverForA);
+ else if (solverForAType == "bicgstabsolver")
+ solverForA_ = std::make_unique<Dune::BiCGSTABSolver<U>>(linearOperator, preconditionerForA_, treeSolverForA);
+ else if (solverForAType == "cgsolver")
+ solverForA_ = std::make_unique<Dune::CGSolver<U>>(linearOperator, preconditionerForA_, treeSolverForA);
+ else
+ DUNE_THROW(Dune::InvalidStateException, solverForAType << "not supported. Use restartedgmressolver, bicgstabsolver or cgsolver");
+ }
+
+ void initAMGSolverForSchur_(const Dune::ParameterTree& params)
+ {
+ schurMatrix_ = schurApproximate_();
+ auto linearOperatorSchur = std::make_shared<LinearOperatorSchur>(*schurMatrix_);
+ amgSolverForSchur_ = std::make_unique<AMGForSchur>(linearOperatorSchur, params);
+ }
+
+ void initUMFPackForSchur_()
+ {
+#if HAVE_UMFPACK
+ schurMatrix_ = schurApproximate_();
+ umfPackSolverForSchur_ = std::make_unique<Dune::UMFPack<D>>(*schurMatrix_);
+#else
+ DUNE_THROW(Dune::InvalidStateException, "UMFPack not available. Use UseDirectSolverForSchurComplement = false.");
+#endif
+ }
+
+ template<class Sol, class Rhs>
+ void applySolverForSchurComplement_(Sol& sol, Rhs& rhs) const
+ {
+ if (useDirectSolverForSchurComplement_)
+ {
+#if HAVE_UMFPACK
+ Dune::InverseOperatorResult res;
+ auto rhsTmp = rhs;
+ umfPackSolverForSchur_->apply(sol, rhsTmp, res);
+#endif
+ }
+ else
+ {
+ amgSolverForSchur_->pre(sol, rhs);
+ amgSolverForSchur_->apply(sol, rhs);
+ amgSolverForSchur_->post(sol);
+ }
+ }
+
+ void initAMGPreconditionerForSchur_(const Dune::ParameterTree& params)
+ {
+ schurMatrix_ = schurApproximate_();
+ auto linearOperatorSchur = std::make_shared<LinearOperatorSchur>(*schurMatrix_);
+ amgPreconditionerForSchur_ = std::make_shared<AMGForSchur>(linearOperatorSchur, params);
+ }
+
+ // explicit assembly of the approximate Schur complement matrix [D - C * diag(A)^-1 * B]
+ auto schurApproximate_()
+ {
+ using namespace Dune::Indices;
+ auto& A = matrix_[_0][_0];
+ auto& B = matrix_[_0][_1];
+ auto& C = matrix_[_1][_0];
+ auto& D = matrix_[_1][_1];
+
+ using AType = std::decay_t<decltype(A)>;
+ using BType = std::decay_t<decltype(B)>;
+ using DType = std::decay_t<decltype(D)>;
+
+ const std::size_t sizeA = A.N();
+
+ //get the diagonal matrix invDiagA = diag(A)^-1
+ AType invDiagA;
+ invDiagA.setBuildMode(AType::random);
+ setDiagPattern(invDiagA, sizeA);
+
+ auto row = A.begin();
+
+ for(; row != A.end(); ++row)
+ {
+ using size_type = typename AType::size_type;
+ size_type rowIdx = row.index();
+ invDiagA[rowIdx][rowIdx] = 1./(A[rowIdx][rowIdx]);
+ }
+
+ // diag(A)^-1 * B
+ BType invAB;
+ Dune::matMultMat(invAB, invDiagA, B);
+
+ // C * diag(A)^-1 * B
+ DType schurApproximate;
+ Dune::matMultMat(schurApproximate, C, invAB);
+
+ // D - C * diag(A)^-1 * B
+ schurApproximate *= -1.0;
+ schurApproximate += D;
+
+ // const std::size_t sizeS = schurApproximate.N();
+ //
+ // //get the diagonal matrix invDiagS = diag(S)^-1
+ // DType invDiagS;
+ // invDiagS.setBuildMode(DType::random);
+ // setDiagPattern(invDiagS, sizeS);
+ //
+ // auto rowS = schurApproximate.begin();
+ //
+ // for(; rowS != schurApproximate.end(); ++rowS)
+ // {
+ // using size_type = typename DType::size_type;
+ // size_type rowSIdx = rowS.index();
+ // invDiagS[rowSIdx][rowSIdx] = 1./(schurApproximate[rowSIdx][rowSIdx]);
+ // }
+
+ // return std::make_unique<DType>(invDiagS);
+
+ return std::make_unique<DType>(schurApproximate);
+ }
+
+ template <class MType>
+ void setDiagPattern(MType& diagM, std::size_t size)
+ {
+ // set the size of the sub-matrices
+ diagM.setSize(size, size);
+
+ // set occupation pattern of the coefficient matrix
+ Dune::MatrixIndexSet occupationPatternDiagM;
+ occupationPatternDiagM.resize(size, size);
+
+ // evaluate the acutal pattern
+ for (int i = 0; i < diagM.N(); ++i)
+ {
+ occupationPatternDiagM.add(i, i);
+ }
+
+ occupationPatternDiagM.exportIdx(diagM);
+ }
+
+ void initStandardSchurComplementSolver_()
+ {
+ auto linearOperator = std::make_shared<SimpleOperator>(matrix_);
+ if(useDiagonalInSchurComplement_)
+ linearOperator->setSolver([&](U& a, U& b, Dune::InverseOperatorResult result) { applyInverseOfDiagonalOfA_(a); });
+ else if(useStandardSolverForA_ && useStandardSolverForAInSchurComplement_){
+ linearOperator->setSolver([&](U& a, U& b, Dune::InverseOperatorResult result) { solverForA_->apply(a, b, result); });
+ }
+ else
+ linearOperator->setSolver([&](U& a, U& b, Dune::InverseOperatorResult result) { applySolverForA_(a, b); });
+
+ Dune::ParameterTree treeSchurComplementSolver;
+ treeSchurComplementSolver["verbose"] = params_.get<std::string>("SolverForSchurVerbosity", "0");
+ treeSchurComplementSolver["reduction"] = params_.get<std::string>("SolverForSchurResidualReduction", "1e-12");
+ treeSchurComplementSolver["maxit"] = params_.get<std::string>("SolverForSchurMaxIter", "1000");
+ treeSchurComplementSolver["restart"] = params_.get<std::string>("SolverForSchurRestart", "100");
+
+ static const bool usePreconditionerForSchurComplement = params_.get<bool>("UsePreconditionerForSchurComplement", false);
+ if (usePreconditionerForSchurComplement)
+ {
+ Dune::ParameterTree treePreconditionerSchur;
+ treePreconditionerSchur["verbose"] = params_.get<std::string>("SolverForSchurVerbosity", "0");
+ // schurComplementPreconditioner_ = std::make_shared<SeqJacobiMatrixFree<P, P>>(linearOperator, treePreconditionerSchur);
+ schurComplementPreconditioner_ = amgPreconditionerForSchur_; // std::make_shared<Dune::SeqJac<D, P, P>>(*schurMatrix_, 1, 1.0);
+ }
+ else
+ schurComplementPreconditioner_ = std::make_shared<Dune::Richardson<P, P>>();
+
+ static const auto schurComplementSolverType = params_.get<std::string>("SolverForSchurType", "restartedgmressolver");
+
+ if (schurComplementSolverType == "restartedgmressolver")
+ schurComplementSolver_ = std::make_unique<Dune::RestartedGMResSolver<P>>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver);
+ else if (schurComplementSolverType == "bicgstabsolver")
+ schurComplementSolver_ = std::make_unique<Dune::BiCGSTABSolver<P>>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver);
+ else if (schurComplementSolverType == "cgsolver")
+ schurComplementSolver_ = std::make_unique<Dune::CGSolver<P>>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver);
+ else
+ DUNE_THROW(Dune::InvalidStateException, schurComplementSolverType << "not supported. Use restartedgmressolver, bicgstabsolver or cgsolver");
+ }
+
+ void applyInverseOfDiagonalOfA_(U& x) const
+ {
+ using namespace Dune::Indices;
+ const auto& A = matrix_[_0][_0];
+
+ for (std::size_t i = 0; i < A.N(); ++i)
+ {
+ auto tmp = A[i][i];
+ tmp.invert();
+ x[i] *= tmp;
+ }
+ }
+
+ //! \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 relaxationFactorP_;
+ scalar_field_type relaxationFactorU_;
+ //! \brief The verbosity level
+ const int verbosity_;
+
+ std::unique_ptr<AMGForA> amgSolverForA_;
+ std::shared_ptr<AMGForA> amgPreconditionerForA_;
+ std::unique_ptr<Dune::IterativeSolver<U,U>> solverForA_;
+ std::shared_ptr<Dune::Preconditioner<U,U>> preconditionerForA_;
+ std::unique_ptr<AMGForSchur> amgSolverForSchur_;
+ std::shared_ptr<AMGForSchur> amgPreconditionerForSchur_;
+ std::unique_ptr<Dune::IterativeSolver<P,P>> schurComplementSolver_;
+ std::shared_ptr<Dune::Preconditioner<P,P>> schurComplementPreconditioner_;
+ std::unique_ptr<A> A_;
+ std::unique_ptr<D> schurMatrix_;
+#if HAVE_UMFPACK
+ std::unique_ptr<Dune::UMFPack<A>> umfPackSolverForA_;
+ std::unique_ptr<Dune::UMFPack<D>> umfPackSolverForSchur_;
+#endif
+ const std::string paramGroup_;
+ const bool useStandardSolverForA_;
+ const bool useDirectSolverForA_;
+ const bool useDiagonalInSchurComplement_;
+ const bool useDiagonalInUpdate_;
+ const bool useStandardSolverForSchurComplement_;
+ const bool useDirectSolverForSchurComplement_;
+ const bool useStandardSolverForAInSchurComplement_;
+
+ const Dune::ParameterTree& params_;
+};
+// DUNE_REGISTER_PRECONDITIONER("simple", defaultPreconditionerBlockLevelCreator<SeqSimple>());
+
+
+
+
+
+
+
diff --git a/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh b/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh
index 78a22d0c60f1e2e24cef2a597757176e1fec2e03..0e17bdf14fd6e67f9b6684edaf3d77b24c5f83f0 100644
--- a/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh
+++ b/dune/phasefield/utility/SEQ_Uzawa_preconditioner.hh
@@ -5,6 +5,10 @@
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>
+#if HAVE_UMFPACK
+#include <dune/istl/umfpack.hh>
+#endif
+
/*!
* \brief A preconditioner based on the Uzawa algorithm for saddle-point problems of the form
* \f$
@@ -69,8 +73,8 @@ public:
* \param mat The matrix to operate on.
* \param params Collection of paramters.
*/
- SeqUzawa(const M& op, const Dune::ParameterTree& params)
- : matrix_(op)
+ SeqUzawa(const M& mat, const Dune::ParameterTree& params)
+ : matrix_(mat)
, numIterations_(params.get<std::size_t>("NumIterations", 1))
, relaxationFactor_(params.get<scalar_field_type>("Relaxation", 1.0))
, verbosity_(params.get<int>("VerbosityLevel", 1))
@@ -191,9 +195,9 @@ private:
* 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.
+ * When facing convergence issues, you may set VerbosityLevel = 1 to see the estimate of lambdaMax.
+ * In a new simulation run, you can then set DetermineRelaxationFactor = false and set some other value
+ * for 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.
*/