seqsolverbackend.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_SEQ_SOLVER_BACKEND_HH
#define DUMUX_SEQ_SOLVER_BACKEND_HH
 
#include <type_traits>
#include <tuple>
#include <utility>
 
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
#include <dune/istl/io.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/common/typetraits/matrix.hh>
#include <dumux/common/typetraits/utility.hh>
#include <dumux/linear/solver.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/preconditioners.hh>
#include <dumux/linear/linearsolverparameters.hh>
 
namespace Dumux {

class IterativePreconditionedSolverImpl
{
public:
 
    template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
    static bool solve(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                      const std::string& modelParamGroup = "")
    {
        Preconditioner precond(A, s.precondIter(), s.relaxation());
 
        // make a linear operator from a matrix
        using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>;
        MatrixAdapter linearOperator(A);
 
        Solver solver(linearOperator, precond, s.residReduction(), s.maxIter(), s.verbosity());
 
        Vector bTmp(b);
 
        Dune::InverseOperatorResult result;
        solver.apply(x, bTmp, result);
 
        return result.converged;
    }
 
    template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
    static bool solveWithGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                               const std::string& modelParamGroup = "")
    {
        // get the restart threshold
        const int restartGMRes = getParamFromGroup<int>(modelParamGroup, "LinearSolver.GMResRestart", 10);
 
        Preconditioner precond(A, s.precondIter(), s.relaxation());
 
        // make a linear operator from a matrix
        using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>;
        MatrixAdapter linearOperator(A);
 
        Solver solver(linearOperator, precond, s.residReduction(), restartGMRes, s.maxIter(), s.verbosity());
 
        Vector bTmp(b);
 
        Dune::InverseOperatorResult result;
        solver.apply(x, bTmp, result);
 
        return result.converged;
    }
 
    template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
    static bool solveWithILU0Prec(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                                  const std::string& modelParamGroup = "")
    {
        Preconditioner precond(A, s.relaxation());
 
        using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>;
        MatrixAdapter operatorA(A);
 
        Solver solver(operatorA, precond, s.residReduction(), s.maxIter(), s.verbosity());
 
        Vector bTmp(b);
 
        Dune::InverseOperatorResult result;
        solver.apply(x, bTmp, result);
 
        return result.converged;
    }
 
    // solve with RestartedGMRes (needs restartGMRes as additional argument)
    template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
    static bool solveWithILU0PrecGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                                       const std::string& modelParamGroup = "")
    {
        // get the restart threshold
        const int restartGMRes = getParamFromGroup<int>(modelParamGroup, "LinearSolver.GMResRestart", 10);
 
        Preconditioner precond(A, s.relaxation());
 
        using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>;
        MatrixAdapter operatorA(A);
 
        Solver solver(operatorA, precond, s.residReduction(), restartGMRes, s.maxIter(), s.verbosity());
 
        Vector bTmp(b);
 
        Dune::InverseOperatorResult result;
        solver.apply(x, bTmp, result);
 
        return result.converged;
    }
 
#if DUNE_VERSION_GTE(DUNE_ISTL,2,7)
    // solve with generic parameter tree
    template<class Preconditioner, class Solver, class Matrix, class Vector>
    static bool solveWithParamTree(const Matrix& A, Vector& x, const Vector& b,
                                   const Dune::ParameterTree& params)
    {
        // make a linear operator from a matrix
        using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>;
        const auto linearOperator = std::make_shared<MatrixAdapter>(A);
 
#if DUNE_VERSION_GTE(DUNE_ISTL,2,8)
        auto precond = std::make_shared<Preconditioner>(linearOperator, params.sub("preconditioner"));
#else
        auto precond = std::make_shared<Preconditioner>(A, params.sub("preconditioner"));
#endif
        Solver solver(linearOperator, precond, params);
 
        Vector bTmp(b);
 
        Dune::InverseOperatorResult result;
        solver.apply(x, bTmp, result);
 
        return result.converged;
    }
#endif
};

template<class M>
constexpr std::size_t preconditionerBlockLevel() noexcept
{
    return isMultiTypeBlockMatrix<M>::value ? 2 : 1;
}

class ILUnBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILUn preconditioned BiCGSTAB solver";
    }
};

class SORBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SOR preconditioned BiCGSTAB solver";
    }
};

class SSORBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SSOR preconditioned BiCGSTAB solver";
    }
};

class GSBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SSOR preconditioned BiCGSTAB solver";
    }
};

class JacBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "Jac preconditioned BiCGSTAB solver";
    }
};

class ILUnCGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILUn preconditioned CG solver";
    }
};

class SORCGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SOR preconditioned CG solver";
    }
};

class SSORCGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SSOR preconditioned CG solver";
    }
};

class GSCGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "GS preconditioned CG solver";
    }
};

class JacCGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "GS preconditioned CG solver";
    }
};

class SSORRestartedGMResBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::RestartedGMResSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "SSOR preconditioned GMRes solver";
    }
};

class ILU0BiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILU0 preconditioned BiCGSTAB solver";
    }
};

class ILU0CGBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::CGSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILU0 preconditioned BiCGSTAB solver";
    }
};

class ILU0RestartedGMResBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::RestartedGMResSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solveWithILU0PrecGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILU0 preconditioned BiCGSTAB solver";
    }
};

class ILUnRestartedGMResBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
        using Solver = Dune::RestartedGMResSolver<Vector>;
 
        return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
    }
 
    std::string name() const
    {
        return "ILUn preconditioned GMRes solver";
    }
};

class ExplicitDiagonalSolver : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        Vector rhs(b);
        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
        Dune::SeqJac<Matrix, Vector, Vector, precondBlockLevel> jac(A, 1, 1.0);
        jac.pre(x, rhs);
        jac.apply(x, rhs);
        jac.post(x);
        return true;
    }
 
    std::string name() const
    {
        return "Explicit diagonal matrix solver";
    }
};
 
#if HAVE_SUPERLU
class SuperLUBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        static_assert(isBCRSMatrix<Matrix>::value, "SuperLU only works with BCRS matrices!");
        using BlockType = typename Matrix::block_type;
        static_assert(BlockType::rows == BlockType::cols, "Matrix block must be quadratic!");
        constexpr auto blockSize = BlockType::rows;
 
        Dune::SuperLU<Matrix> solver(A, this->verbosity() > 0);
 
        Vector bTmp(b);
        solver.apply(x, bTmp, result_);
 
        int size = x.size();
        for (int i = 0; i < size; i++)
        {
            for (int j = 0; j < blockSize; j++)
            {
                using std::isnan;
                using std::isinf;
                if (isnan(x[i][j]) || isinf(x[i][j]))
                {
                    result_.converged = false;
                    break;
                }
            }
        }
 
        return result_.converged;
    }
 
    std::string name() const
    {
        return "SuperLU solver";
    }
 
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
 
private:
    Dune::InverseOperatorResult result_;
};
#endif // HAVE_SUPERLU
 
#if HAVE_UMFPACK
class UMFPackBackend : public LinearSolver
{
public:
 
    UMFPackBackend(const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
    {
        ordering_ = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.UMFPackOrdering", 1);
    }

    void setOrdering(int i)
    { ordering_ = i; }

    int ordering() const
    { return ordering_; }
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        static_assert(isBCRSMatrix<Matrix>::value, "UMFPack only works with BCRS matrices!");
        using BlockType = typename Matrix::block_type;
        static_assert(BlockType::rows == BlockType::cols, "Matrix block must be quadratic!");
        constexpr auto blockSize = BlockType::rows;
 
        Dune::UMFPack<Matrix> solver;
        solver.setVerbosity(this->verbosity() > 0);
        solver.setOption(UMFPACK_ORDERING, ordering_);
        solver.setMatrix(A);
 
        Vector bTmp(b);
        solver.apply(x, bTmp, result_);
 
        int size = x.size();
        for (int i = 0; i < size; i++)
        {
            for (int j = 0; j < blockSize; j++)
            {
                using std::isnan;
                using std::isinf;
                if (isnan(x[i][j]) || isinf(x[i][j]))
                {
                    result_.converged = false;
                    break;
                }
            }
        }
 
        return result_.converged;
    }
 
    std::string name() const
    {
        return "UMFPack solver";
    }
 
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
 
private:
    Dune::InverseOperatorResult result_;
    int ordering_;
};
#endif // HAVE_UMFPACK
 

// \{

template <class LinearSolverTraits>
class UzawaBiCGSTABBackend : public LinearSolver
{
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        using Preconditioner = SeqUzawa<Matrix, Vector, Vector>;
        using Solver = Dune::BiCGSTABSolver<Vector>;
        static const auto solverParams = LinearSolverParameters<LinearSolverTraits>::createParameterTree(this->paramGroup());
        return IterativePreconditionedSolverImpl::template solveWithParamTree<Preconditioner, Solver>(A, x, b, solverParams);
    }
 
    std::string name() const
    {
        return "Uzawa preconditioned BiCGSTAB solver";
    }
};

template<class M, class X, class Y, int blockLevel = 2>
class BlockDiagILU0Preconditioner : public Dune::Preconditioner<X, Y>
{
    template<std::size_t i>
    using DiagBlockType = std::decay_t<decltype(std::declval<M>()[Dune::index_constant<i>{}][Dune::index_constant<i>{}])>;
 
    template<std::size_t i>
    using VecBlockType = std::decay_t<decltype(std::declval<X>()[Dune::index_constant<i>{}])>;
 
    template<std::size_t i>
    using BlockILU = Dune::SeqILU<DiagBlockType<i>, VecBlockType<i>, VecBlockType<i>, blockLevel-1>;
 
    using ILUTuple = typename makeFromIndexedType<std::tuple, BlockILU, std::make_index_sequence<M::N()> >::type;
 
public:
    using matrix_type = typename std::decay_t<M>;
    using domain_type = X;
    using range_type = Y;
    using field_type = typename X::field_type;

    BlockDiagILU0Preconditioner(const M& m, double w = 1.0)
    : BlockDiagILU0Preconditioner(m, w, std::make_index_sequence<M::N()>{})
    {
        static_assert(blockLevel >= 2, "Only makes sense for MultiTypeBlockMatrix!");
    }
 
    void pre (X& v, Y& d) final {}
 
    void apply (X& v, const Y& d) final
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(ilu_)), [&](const auto i)
        {
            std::get<decltype(i)::value>(ilu_).apply(v[i], d[i]);
        });
    }
 
    void post (X&) final {}

    Dune::SolverCategory::Category category() const final
    {
        return Dune::SolverCategory::sequential;
    }
 
private:
    template<std::size_t... Is>
    BlockDiagILU0Preconditioner (const M& m, double w, std::index_sequence<Is...> is)
    : ilu_(std::make_tuple(BlockILU<Is>(m[Dune::index_constant<Is>{}][Dune::index_constant<Is>{}], w)...))
    {}
 
    ILUTuple ilu_;
};
 

class BlockDiagILU0BiCGSTABSolver : public LinearSolver
{
 
public:
    using LinearSolver::LinearSolver;
 
    template<class Matrix, class Vector>
    bool solve(const Matrix& M, Vector& x, const Vector& b)
    {
        BlockDiagILU0Preconditioner<Matrix, Vector, Vector> preconditioner(M);
        Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
        Dune::BiCGSTABSolver<Vector> solver(op, preconditioner, this->residReduction(),
                                            this->maxIter(), this->verbosity());
        auto bTmp(b);
        solver.apply(x, bTmp, result_);
 
        return result_.converged;
    }
 
    const Dune::InverseOperatorResult& result() const
    {
      return result_;
    }
 
    std::string name() const
    { return "block-diagonal ILU0-preconditioned BiCGSTAB solver"; }
 
private:
    Dune::InverseOperatorResult result_;
};

class BlockDiagILU0RestartedGMResSolver : public LinearSolver
{
 
public:
    using LinearSolver::LinearSolver;
 
    template<int precondBlockLevel = 2, class Matrix, class Vector>
    bool solve(const Matrix& M, Vector& x, const Vector& b)
    {
        BlockDiagILU0Preconditioner<Matrix, Vector, Vector> preconditioner(M);
        Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
        static const int restartGMRes = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.GMResRestart");
        Dune::RestartedGMResSolver<Vector> solver(op, preconditioner, this->residReduction(), restartGMRes,
                                                  this->maxIter(), this->verbosity());
        auto bTmp(b);
        solver.apply(x, bTmp, result_);
 
        return result_.converged;
    }
 
    const Dune::InverseOperatorResult& result() const
    {
      return result_;
    }
 
    std::string name() const
    { return "block-diagonal ILU0-preconditioned restarted GMRes solver"; }
 
private:
    Dune::InverseOperatorResult result_;
};

template<class M, class X, class Y, int blockLevel = 2>
class BlockDiagAMGPreconditioner : public Dune::Preconditioner<X, Y>
{
    template<std::size_t i>
    using DiagBlockType = std::decay_t<decltype(std::declval<M>()[Dune::index_constant<i>{}][Dune::index_constant<i>{}])>;
 
    template<std::size_t i>
    using VecBlockType = std::decay_t<decltype(std::declval<X>()[Dune::index_constant<i>{}])>;
 
    template<std::size_t i>
    using Smoother = Dune::SeqSSOR<DiagBlockType<i>, VecBlockType<i>, VecBlockType<i>>;
 
    template<std::size_t i>
    using LinearOperator = Dune::MatrixAdapter<DiagBlockType<i>, VecBlockType<i>, VecBlockType<i>>;
 
    template<std::size_t i>
    using ScalarProduct = Dune::SeqScalarProduct<VecBlockType<i>>;
 
    template<std::size_t i>
    using BlockAMG = Dune::Amg::AMG<LinearOperator<i>, VecBlockType<i>, Smoother<i>>;
 
    using AMGTuple = typename makeFromIndexedType<std::tuple, BlockAMG, std::make_index_sequence<M::N()> >::type;
 
public:
    using matrix_type = typename std::decay_t<M>;
    using domain_type = X;
    using range_type = Y;
    using field_type = typename X::field_type;

    template<class LOP, class Criterion, class SmootherArgs>
    BlockDiagAMGPreconditioner(const LOP& lop, const Criterion& c, const SmootherArgs& sa)
    : BlockDiagAMGPreconditioner(lop, c, sa, std::make_index_sequence<M::N()>{})
    {
        static_assert(blockLevel >= 2, "Only makes sense for MultiTypeBlockMatrix!");
    }

    void pre (X& v, Y& d) final
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(amg_)), [&](const auto i)
        {
            std::get<decltype(i)::value>(amg_).pre(v[i], d[i]);
        });
    }

    void apply (X& v, const Y& d) final
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(amg_)), [&](const auto i)
        {
            std::get<decltype(i)::value>(amg_).apply(v[i], d[i]);
        });
    }

    void post (X& v) final
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(amg_)), [&](const auto i)
        {
            std::get<decltype(i)::value>(amg_).post(v[i]);
        });
    }

    Dune::SolverCategory::Category category() const final
    {
        return Dune::SolverCategory::sequential;
    }
 
private:
    template<class LOP, class Criterion, class SmootherArgs, std::size_t... Is>
    BlockDiagAMGPreconditioner (const LOP& lop, const Criterion& c, const SmootherArgs& sa, std::index_sequence<Is...> is)
    : amg_(std::make_tuple(BlockAMG<Is>(*std::get<Is>(lop), *std::get<Is>(c), *std::get<Is>(sa))...))
    {}
 
    AMGTuple amg_;
};

class BlockDiagAMGBiCGSTABSolver : public LinearSolver
{
    template<class M, std::size_t i>
    using DiagBlockType = std::decay_t<decltype(std::declval<M>()[Dune::index_constant<i>{}][Dune::index_constant<i>{}])>;
 
    template<class X, std::size_t i>
    using VecBlockType = std::decay_t<decltype(std::declval<X>()[Dune::index_constant<i>{}])>;
 
    template<class M, class X, std::size_t i>
    using Smoother = Dune::SeqSSOR<DiagBlockType<M, i>, VecBlockType<X, i>, VecBlockType<X, i>>;
 
    template<class M, class X, std::size_t i>
    using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother<M, X, i>>::Arguments;
 
    template<class M, std::size_t i>
    using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<DiagBlockType<M, i>, Dune::Amg::FirstDiagonal>>;
 
    template<class M, class X, std::size_t i>
    using LinearOperator = Dune::MatrixAdapter<DiagBlockType<M, i>, VecBlockType<X, i>, VecBlockType<X, i>>;
 
public:
    using LinearSolver::LinearSolver;
 
    // Solve saddle-point problem using a Schur complement based preconditioner
    template<class Matrix, class Vector>
    bool solve(const Matrix& m, Vector& x, const Vector& b)
    {
        Dune::Amg::Parameters params(15, 2000, 1.2, 1.6, Dune::Amg::atOnceAccu);
        params.setDefaultValuesIsotropic(3);
        params.setDebugLevel(this->verbosity());
 
        auto criterion = makeCriterion_<Criterion, Matrix>(params, std::make_index_sequence<Matrix::N()>{});
        auto smootherArgs = makeSmootherArgs_<SmootherArgs, Matrix, Vector>(std::make_index_sequence<Matrix::N()>{});
 
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<Matrix::N()>{}, [&](const auto i)
        {
            auto& args = std::get<decltype(i)::value>(smootherArgs);
            args->iterations = 1;
            args->relaxationFactor = 1;
        });
 
        auto linearOperator = makeLinearOperator_<LinearOperator, Matrix, Vector>(m, std::make_index_sequence<Matrix::N()>{});
 
        BlockDiagAMGPreconditioner<Matrix, Vector, Vector> preconditioner(linearOperator, criterion, smootherArgs);
 
        Dune::MatrixAdapter<Matrix, Vector, Vector> op(m);
        Dune::BiCGSTABSolver<Vector> solver(op, preconditioner, this->residReduction(),
                                            this->maxIter(), this->verbosity());
        auto bTmp(b);
        solver.apply(x, bTmp, result_);
 
        return result_.converged;
    }
 
    const Dune::InverseOperatorResult& result() const
    {
      return result_;
    }
 
    std::string name() const
    { return "block-diagonal AMG-preconditioned BiCGSTAB solver"; }
 
private:
    template<template<class M, std::size_t i> class Criterion, class Matrix, class Params, std::size_t... Is>
    auto makeCriterion_(const Params& p, std::index_sequence<Is...>)
    {
        return std::make_tuple(std::make_shared<Criterion<Matrix, Is>>(p)...);
    }
 
    template<template<class M, class X, std::size_t i> class SmootherArgs, class Matrix, class Vector, std::size_t... Is>
    auto makeSmootherArgs_(std::index_sequence<Is...>)
    {
        return std::make_tuple(std::make_shared<SmootherArgs<Matrix, Vector, Is>>()...);
    }
 
    template<template<class M, class X, std::size_t i> class LinearOperator, class Matrix, class Vector, std::size_t... Is>
    auto makeLinearOperator_(const Matrix& m, std::index_sequence<Is...>)
    {
        return std::make_tuple(std::make_shared<LinearOperator<Matrix, Vector, Is>>(m[Dune::index_constant<Is>{}][Dune::index_constant<Is>{}])...);
    }
 
    Dune::InverseOperatorResult result_;
};
 
// \}
 
} // end namespace Dumux
 
#endif