Supported AMG as a linear solver without dune-pdelab

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/strip-pdelab@12914 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/linear/amgbackend.hh

@@ -24,132 +24,18 @@
 #ifndef DUMUX_AMGBACKEND_HH
 #define DUMUX_AMGBACKEND_HH
 
-#if HAVE_DUNE_PDELAB
-
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-#include <dune/pdelab/backend/novlpistlsolverbackend.hh>
-#include <dune/pdelab/backend/ovlpistlsolverbackend.hh>
-#include <dune/pdelab/backend/seqistlsolverbackend.hh>
-#include <dune/pdelab/backend/istlvectorbackend.hh>
-
-#include <dune/pdelab/finiteelementmap/q1fem.hh>
+#include <dune/common/parallel/indexset.hh>
+#include <dune/common/parallel/mpicollectivecommunication.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/solvers.hh>
+
+#include <dumux/linear/amgproperties.hh>
+#include <dumux/linear/amgparallelhelpers.hh>
 #include <dumux/linear/p0fem.hh>
 
-#include <dumux/implicit/box/boxproperties.hh>
-#include <dumux/implicit/cellcentered/ccproperties.hh>
-#include <dumux/decoupled/common/pressureproperties.hh>
-#include "linearsolverproperties.hh"
-
 namespace Dumux {
 
-// forward declaration for the property definitions
-template <class TypeTag> class AMGBackend;
-
-namespace Properties
-{
-//! the PDELab finite element map used for the gridfunctionspace
-NEW_PROP_TAG(AMGLocalFemMap);
-
-//! box: use the (multi-)linear local FEM space associated with cubes by default
-SET_PROP(BoxModel, AMGLocalFemMap)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    enum{dim = GridView::dimension};
-public:
-    typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type;
-};
-
-//! cell-centered: use the element-wise constant local FEM space by default
-SET_PROP(CCModel, AMGLocalFemMap)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    enum{dim = GridView::dimension};
-public:
-    typedef Dumux::P0LocalFiniteElementMap<Scalar,Scalar,dim>  type;
-};
-
-//! decoupled models: use the element-wise constant local FEM space by default
-SET_PROP(DecoupledModel, AMGLocalFemMap)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    enum{dim = GridView::dimension};
-public:
-    typedef Dumux::P0LocalFiniteElementMap<Scalar,Scalar,dim>  type;
-};
-
-//! the type of the employed PDELab backend
-NEW_PROP_TAG(AMGPDELabBackend);
-
-//! box: use the non-overlapping PDELab AMG backend
-SET_PROP(BoxModel, AMGPDELabBackend)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_AMG_SSOR<GridOperator> type;
-};
-
-//! cell-centered: use the overlapping PDELab AMG backend
-SET_PROP(CCModel, AMGPDELabBackend)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<GridOperator> type;
-};
-
-//! decoupled model: use the overlapping PDELab AMG backend
-SET_PROP(DecoupledModel, AMGPDELabBackend)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<GridOperator> type;
-};
-
-//! box: reset the type of solution vector to be PDELab conforming
-SET_PROP(BoxModel, SolutionVector)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef typename GridOperator::Traits::Domain type;
-};
-
-//! cell-centered: reset the type of solution vector to be PDELab conforming
-SET_PROP(CCModel, SolutionVector)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef typename GridOperator::Traits::Domain type;
-};
-
-
-//! decoupled model: reset the type of solution vector to be PDELab conforming
-SET_PROP(DecoupledModel, PressureSolutionVector)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef typename GridOperator::Traits::Domain type;
-};
-
-//! decoupled model: reset the type of solution vector to be PDELab conforming
-SET_PROP(DecoupledModel, PressureRHSVector)
-{
-    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
-public:
-    typedef typename GridOperator::Traits::Domain type;
-};
-
-//! set a property JacobianMatrix also for the decoupled models
-SET_PROP(DecoupledModel, JacobianMatrix)
-{
-public:
-    typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) type;
-};
-
-
-}
-
 /*!
  * \brief Scale the linear system by the inverse of 
  * its (block-)diagonal entries.
@@ -188,36 +74,20 @@ class AMGBackend
 {
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, AMGLocalFemMap) LocalFemMap;
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    enum { dim = GridView::dimension };
-    typedef typename Dune::PDELab::NoConstraints Constraints;
     typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
     enum { numEq = JacobianMatrix::block_type::rows};
-    typedef Dune::PDELab::GridFunctionSpace<GridView, 
-                                            LocalFemMap, 
-                                            Constraints, 
-                                            Dune::PDELab::ISTLVectorBackend<numEq> 
-                                           > ScalarGridFunctionSpace;
-    typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, 
-                                                 numEq, 
-                                                 Dune::PDELab::GridFunctionSpaceBlockwiseMapper
-                                                > GridFunctionSpace;
-    typedef typename GridFunctionSpace::template ConstraintsContainer<Scalar>::Type ConstraintsTrafo;
-    typedef int LocalOperator;
 
-public:
-    typedef typename Dune::PDELab::GridOperator<GridFunctionSpace,
-                                       GridFunctionSpace,
-                                       LocalOperator,
-                                       Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-                                       Scalar, Scalar, Scalar,
-                                       ConstraintsTrafo,
-                                       ConstraintsTrafo,
-                                       true
-                                      > GridOperator;
-    typedef typename GET_PROP_TYPE(TypeTag, AMGPDELabBackend) PDELabBackend;
+    typedef typename GET_PROP(TypeTag, AMGPDELabBackend) PDELabBackend;
+    typedef typename PDELabBackend::LinearOperator LinearOperator;
+    typedef typename PDELabBackend::VType VType;
+    typedef typename PDELabBackend::Comm Comm;
+    typedef typename PDELabBackend::Smoother Smoother;
+    typedef Dune::Amg::AMG<typename PDELabBackend::LinearOperator, VType,
+                           Smoother,Comm> AMGType;
+    typedef typename PDELabBackend::LinearOperator::matrix_type BCRSMat;
 
+public:    
     /*!
      * \brief Construct the backend.
      * 
@@ -226,13 +96,6 @@ public:
     AMGBackend(const Problem& problem)
     : problem_(problem)
     {
-        fem_ = Dune::make_shared<LocalFemMap>();
-        constraints_ = Dune::make_shared<Constraints>();
-        scalarGridFunctionSpace_ = Dune::make_shared<ScalarGridFunctionSpace>(problem.gridView(), *fem_, *constraints_);
-        gridFunctionSpace_ = Dune::make_shared<GridFunctionSpace>(*scalarGridFunctionSpace_);
-        int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
-        int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
-        imp_ = Dune::make_shared<PDELabBackend>(*gridFunctionSpace_, maxIt, verbosity);
     }
 
     /*!
@@ -245,17 +108,49 @@ public:
     template<class Matrix, class Vector>
     bool solve(Matrix& A, Vector& x, Vector& b)
     {
+        int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
+        int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
         static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
-        imp_->apply(A, x, b, residReduction);
-
-        result_.converged  = imp_->result().converged;
-        result_.iterations = imp_->result().iterations;
-        result_.elapsed    = imp_->result().elapsed;
-        result_.reduction  = imp_->result().reduction;
-        result_.conv_rate  = imp_->result().conv_rate;
+        
+#if HAVE_MPI        
+        typename PDELabBackend::Comm comm(problem_.model().gridView().comm(),
+                                 PDELabBackend::isNonOverlapping?
+                                 Dune::SolverCategory::nonoverlapping:
+                                 Dune::SolverCategory::overlapping);        
+
+        if(PDELabBackend::isNonOverlapping)
+        {
+            // extend the matrix pattern such that it is usable for AMG        
+
+        }
+        typename PDELabBackend::LinearOperator fop(A, comm);
+        typename PDELabBackend::ScalarProduct sp(comm);
+        int rank = comm_.communicator().rank();
+#else
+        typename PDELabBackend::Comm  comm;
+        typename PDELabBackend::LinearOperator fop(A);
+        typename PDELabBackend::ScalarProduct sp;
+        int rank=0;
+#endif
+        
+        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments
+            SmootherArgs;
+        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,
+                                                                          Dune::Amg::FirstDiagonal> >
+            Criterion;
+        Criterion criterion(15, 2000);
+        criterion.setDefaultValuesIsotropic(2);
+        SmootherArgs smootherArgs;
+        smootherArgs.iterations = 1;
+        smootherArgs.relaxationFactor = 1;
+        
+        AMGType amg(fop, criterion, smootherArgs, 1, 1, 1, false, comm);
+        Dune::BiCGSTABSolver<typename PDELabBackend::VType> solver(fop, sp, amg, residReduction, maxIt, 
+                                                                 rank==0?verbosity: 0);
+        solver.apply(x, b, result_);
 
         return result_.converged;
-    }
+        }
 
     /*!
      * \brief The result containing the convergence history.
@@ -267,11 +162,6 @@ public:
     
 private:
     const Problem& problem_;
-    Dune::shared_ptr<LocalFemMap> fem_;
-    Dune::shared_ptr<Constraints> constraints_;
-    Dune::shared_ptr<ScalarGridFunctionSpace> scalarGridFunctionSpace_;
-    Dune::shared_ptr<GridFunctionSpace> gridFunctionSpace_;
-    Dune::shared_ptr<PDELabBackend> imp_;
     Dune::InverseOperatorResult result_;
 };
 
@@ -282,8 +172,15 @@ template <class TypeTag>
 class SeqAMGBackend
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename AMGBackend<TypeTag>::GridOperator GridOperator;
-    typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GridOperator> PDELabBackend;
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
+    typedef Dune::Amg::SequentialInformation Comm;
+    typedef Dune::AssembledLinearOperator<MType,VType, VType> LinearOperator;
+    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+    typedef Dune::Amg::AMG<LinearOperator,VType,Smoother> AMGType;
 public:
 
     /*!
@@ -307,19 +204,22 @@ public:
     {
         int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
         int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
-        imp_ = Dune::make_shared<PDELabBackend>(maxIt, verbosity);
-
         static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
-        imp_->apply(A, x, b, residReduction);
-
-        result_.converged  = imp_->result().converged;
-        result_.iterations = imp_->result().iterations;
-        result_.elapsed    = imp_->result().elapsed;
-        result_.reduction  = imp_->result().reduction;
-        result_.conv_rate  = imp_->result().conv_rate;
-
-        imp_.template reset<PDELabBackend>(0);
-        
+        LinearOperator fop(A);
+        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<typename LinearOperator::matrix_type,
+                                                                          Dune::Amg::FirstDiagonal> >
+            Criterion;
+        Criterion criterion(15,2000);
+        criterion.setDefaultValuesIsotropic(2);
+        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+        SmootherArgs smootherArgs;
+        smootherArgs.iterations = 1;
+        smootherArgs.relaxationFactor = 1;
+
+        AMGType amg(fop, criterion, smootherArgs, 1, 1, 1);
+        Dune::BiCGSTABSolver<VType> solver(fop, amg, residReduction, maxIt, 
+                                          verbosity);
+        solver.apply(x, b, result_);        
         return result_.converged;
     }
 
@@ -333,7 +233,6 @@ public:
 
 private:
     const Problem& problem_;
-    Dune::shared_ptr<PDELabBackend> imp_;
     Dune::InverseOperatorResult result_;
 };
 
@@ -347,8 +246,15 @@ template <class TypeTag>
 class ScaledSeqAMGBackend
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename AMGBackend<TypeTag>::GridOperator GridOperator;
-    typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GridOperator> PDELabBackend;
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
+    typedef Dune::Amg::SequentialInformation Comm;
+    typedef Dune::MatrixAdapter<MType,VType, VType> LinearOperator;
+    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+    typedef Dune::Amg::AMG<LinearOperator,VType,Smoother> AMGType;
 public:
 
     /*!
@@ -374,18 +280,20 @@ public:
 
         int maxIt = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations);
         int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
-        imp_ = Dune::make_shared<PDELabBackend>(maxIt, verbosity);
-        
         static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
-        imp_->apply(A, x, b, residReduction);
-
-        result_.converged  = imp_->result().converged;
-        result_.iterations = imp_->result().iterations;
-        result_.elapsed    = imp_->result().elapsed;
-        result_.reduction  = imp_->result().reduction;
-        result_.conv_rate  = imp_->result().conv_rate;
-
-        imp_.template reset<PDELabBackend>(0);
+        LinearOperator fop(A);
+        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MType,
+                                                                          Dune::Amg::FirstDiagonal> >
+            Criterion;
+        Criterion criterion(15,2000);
+        criterion.setDefaultValuesIsotropic(2);
+        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+        SmootherArgs smootherArgs;
+        smootherArgs.iterations = 1;
+        smootherArgs.relaxationFactor = 1;
+        AMGType amg(fop, criterion, smootherArgs, 1, 1, 1);
+        Dune::BiCGSTABSolver<VType> solver(fop, amg, residReduction, maxIt, 
+                                          verbosity);
         
         return result_.converged;
     }
@@ -400,11 +308,9 @@ public:
 
 private:
     const Problem& problem_;
-    Dune::shared_ptr<PDELabBackend> imp_;
     Dune::InverseOperatorResult result_;
 };
 
 } // namespace Dumux
 
-#endif // HAVE_DUNE_PDELAB
 #endif // DUMUX_AMGBACKEND_HH

dumux/linear/amgparallelhelpers.hh

+// -*- 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 2 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_AMGPARALLELHELPERS_HH
+#define DUMUX_AMGPARALLELHELPERS_HH
+
+#include <dumux/implicit/box/boxproperties.hh>
+#include <dumux/implicit/cellcentered/ccproperties.hh>
+#include <dumux/decoupled/common/pressureproperties.hh>
+
+
+namespace Dumux
+{
+
+//========================================================
+// A parallel helper class providing a nonoverlapping
+// decomposition of all degrees of freedom
+//========================================================
+
+// operator that resets result to zero at constrained DOFS
+template<class TypeTag>
+class ParallelISTLHelper
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, AMGLocalFemMap) LocalFemMap;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? GridView::dimension : 0 };
+
+    class BaseGatherScatter
+    {
+    public:
+        BaseGatherScatter(const Problem& problem)
+        : problem_(problem)
+        {}
+
+        bool contains(int dim, int codim)
+        {
+            return dofCodim==codim;
+        }
+
+        bool fixedSize(int dim, int codim)
+        {
+            return true;
+
+        }
+        template<class EntityType>
+        int map(const EntityType& e)
+        {
+            return problem_.model().dofMapper().map(e);
+        }
+
+    private:
+        Problem& problem_;
+    };
+
+
+    /**
+     * @brief Writes 1<<24 to each data item (of the container) that is gathered or scattered
+     * and is neither interior nor border.
+     *
+     * Can be used to mark ghost cells.
+     */
+    class GhostGatherScatter
+        : public BaseGatherScatter, Dune::CommDataHandleIF<GhostGatherScatter,std::size_t>
+    {
+    public:
+        typedef std::size_t DataType;
+
+        GhostGatherScatter(std::vector<std::size_t>& ranks,
+                           const Problem& problem)
+            : BaseGatherScatter(problem), ranks_(ranks)
+        {}
+
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType& e)
+        {
+            std::size_t& data= ranks_[this->map(e)];
+            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
+                data = (1<<24);
+            buff.write(data);
+        }
+
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType& e)
+        {
+            std::size_t x;
+            std::size_t& data = ranks_[this->map(e)];
+            buff.read(x);
+            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
+                data= (1<<24);
+            else
+                data = x;
+        }
+    private:
+        std::vector<std::size_t>& ranks_;
+    };
+
+    /**
+     * @brief GatherScatter handle that sets 1<<24 for data items neither associated to
+     * the interior or border and take the minimum when scattering.
+     *
+     * Used to compute an owner rank for each unknown.
+     */
+    class InteriorBorderGatherScatter
+        : public BaseGatherScatter, Dune::CommDataHandleIF<InteriorBorderGatherScatter,std::size_t>
+    {
+    public:
+        typedef std::size_t DataType;
+
+        InteriorBorderGatherScatter(std::vector<std::size_t>& ranks,
+                                    const Problem& problem)
+            : BaseGatherScatter(problem), ranks_(ranks)
+        {}
+
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType& e)
+        {
+
+            std::size_t& data = ranks_[this->map(e)];
+            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
+                data = (1<<24);
+            buff.write(data);
+        }
+
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType& e)
+        {
+            std::size_t x;
+            std::size_t& data = ranks_[this->map(e)];
+            buff.read(x);
+            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
+                data = x;
+            else
+                data = std::min(data,x);
+        }
+    private:
+        std::vector<std::size_t>& ranks_;
+    };
+
+    /**
+     * @brief GatherScatter handle for finding out about neighbouring processor ranks.
+     *
+     */
+    struct NeighbourGatherScatter
+        : public BaseGatherScatter, Dune::CommDataHandleIF<NeighbourGatherScatter,int>
+    {
+        typedef int DataType;
+
+        NeighbourGatherScatter(int rank_, std::set<int>& neighbours_)
+            : myrank(rank_), neighbours(neighbours_)
+        {}
+
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType &e)
+        {
+            buff.write(myrank);
+        }
+
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType &e)
+        {
+            int x;
+            buff.read(x);
+            neighbours.insert(x);
+        }
+        int myrank;
+        std::set<int>& neighbours;
+    };
+
+
+    /**
+     * @brief GatherScatter handle for finding out about neighbouring processor ranks.
+     *
+     */
+    struct SharedGatherScatter
+        : public BaseGatherScatter, Dune::CommDataHandleIF<SharedGatherScatter,int>
+    {
+        typedef int DataType;
+
+        SharedGatherScatter(std::vector<int>& shared,
+                            const Problem& problem)
+            : BaseGatherScatter(problem), shared_(shared)
+        {}
+
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, EntityType& e)
+        {
+            bool data=true;
+            buff.write(data);
+        }
+
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType &e)
+        {
+            bool x;
+            buff.read(x);
+            bool& data= shared_[this->map(e)];
+            data = data || x;
+        }
+    private:
+        std::vector<int>& shared_;
+
+    };
+
+    /**
+     * @brief GatherScatter handle for finding out about neighbouring processor ranks.
+     *
+     */
+    template<typename GI>
+    struct GlobalIndexGatherScatter
+        : public BaseGatherScatter, Dune::CommDataHandleIF<GlobalIndexGatherScatter<GI>, GI>
+    {
+        typedef GI DataType;
+        GlobalIndexGatherScatter(std::vector<GI>& gindices,
+                                 const Problem& problem)
+            : BaseGatherScatter(problem), gindices_(gindices_)
+        {}
+
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType& e)
+        {
+            buff.write(gindices_[this->map(e)]);
+        }
+
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType& e)
+        {
+            DataType x;
+            buff.read(x);
+            gindices_[this->map(e)] = std::min(gindices_[this->map(e)], x);
+        }
+    private:
+        std::vector<GI>& gindices_;
+    };
+
+public:
+
+    ParallelISTLHelper (const Problem& problem, int verbose=1)
+        : problem_(problem), owner_(problem.dofMapper().size(),
+                                    problem.gridView().comm().rank()),
+          isGhost_(problem.dofMapper().size(),0.0), verbose_(verbose)
+    {
+        // find out about ghosts
+        GhostGatherScatter ggs(owner_,problem);
+
+        if (problem.gridView().comm().size()>1)
+            problem.gridView().communicate(ggs,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
+
+        isGhost_ = owner_;
+
+        // partition interior/border
+        InteriorBorderGatherScatter dh(owner_, problem_);
+
+        if (problem.gridView().comm().size()>1)
+            problem.gridView().communicate(dh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
+
+        // convert vector into mask vector
+        for(auto v=owner_.begin(), vend=owner_.end(); v!=vend;++v)
+            if(*v=problem.gridView().comm().rank())
+                *v=1.0;
+            else
+                *v=0.0;
+    }
+
+    // keep only DOFs assigned to this processor
+    template<typename W>
+    void mask (W& w) const
+    {
+        auto v1=w.begin();
+
+        for(auto v2=owner_.begin(), vend=owner_.end(); v2!=vend;++v1,++v2)
+            v1*=v2;
+    }
+
+    // access to mask vector
+    double mask (std::size_t i) const
+    {
+        return owner_[i];
+    }
+
+    // access to ghost vector
+    double ghost (std::size_t i) const
+    {
+        return isGhost_[i];
+    }
+
+#if HAVE_MPI
+
+    /**
+     * @brief Creates a matrix suitable for parallel AMG and the parallel information
+     *
+     * It is silently assumed that the unknows are associated with vertices.
+     *
+     * @tparam MatrixType The type of the ISTL matrix used.
+     * @tparam Comm The type of the OwnerOverlapCopyCommunication
+     * @param m The local matrix.
+     * @param c The parallel information object providing index set, interfaces and
+     * communicators.
+     */
+    template<typename MatrixType, typename Comm>
+    void createIndexSetAndProjectForAMG(MatrixType& m, Comm& c);
+#endif
+private:
+    const Problem& problem_;
+    std::vector<std::size_t> owner_; // vector to identify unique decomposition
+    std::vector<std::size_t> isGhost_; //vector to identify ghost dofs
+    int verbose_; //verbosity
+};
+
+/**
+ * @brief Helper class for adding up matrix entries on border.
+ * @tparam GridOperator The grid operator to work on.
+ * @tparam MatrixType The MatrixType.
+ */
+template<class TypeTag>
+class EntityExchanger
+{
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, AMGLocalFemMap) LocalFemMap;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > Matrix;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    enum {dim = GridView::dimension};
+    typedef typename GridView::Traits::Grid Grid;
+    typedef typename Matrix::block_type BlockType;
+    typedef typename GridView::template Codim<LocalFemMap::dof_codim>::Iterator EntityIterator;
+    typedef typename Grid::Traits::GlobalIdSet IDS;
+    typedef typename IDS::IdType IdType;
+    typedef typename Matrix::RowIterator RowIterator;
+    typedef typename Matrix::ColIterator ColIterator;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dim : 0 };
+
+public:
+    /*! \brief Constructor. Sets up the local to global relations.
+
+      \param[in] gridView The grid view to operate on.
+    */
+    EntityExchanger(const Problem& problem)
+    : problem_(problem)
+    {
+        gid2Index_.clear();
+        index2GID_.clear();
+
+        const GridView& gridView = problem.model().gridView();
+
+        EntityIterator entityEndIt = gridView.template end<dim>();
+        for (EntityIterator entityIt = gridView.template begin<dim>(); entityIt != entityEndIt; ++entityIt)
+        {
+            if (entityIt->partitionType() == Dune::BorderEntity)
+            {
+                int localIdx = problem_.model().dofMapper().map(*entityIt);
+                IdType globalIdx = gridView.grid().globalIdSet().id(*entityIt);
+
+                std::pair<IdType,int> g2iPair(globalIdx, localIdx);
+                gid2Index_.insert(g2iPair);
+
+                std::pair<int,IdType> i2gPair(localIdx, globalIdx);
+                index2GID_.insert(i2gPair);
+
+            }
+        }
+    }
+
+    /**
+     * @brief A DataHandle class to exchange matrix sparsity patterns.
+     *
+     *  We look at a 2D example with a nonoverlapping grid,
+     *  two processes and no ghosts with Q1 discretization.
+     *  Process 0 has the left part of the domain
+     *  with three cells and eight vertices (1-8),
+     *  Process 1 the right part with three cells
+     *  and eight vertices (2,4,7-12).
+     *  <pre>
+     *  1 _ 2        2 _ 9 _ 10
+     *  |   |        |   |   |
+     *  3 _ 4 _ 7    4 _ 7 _ 11
+     *  |   |   |        |   |
+     *  5 _ 6 _ 8        8 _ 12
+     *  </pre>
+     *  If we look at vertex 7 and the corresponding entries in the matrix for P0,
+     *  there will be entries for (7,4) and (7,8), but not for (7,2).
+     *  The MatPatternExchange class will find these entries and returns a vector "sparsity",
+     *  that contains all missing connections.
+     */
+    class MatPatternExchange
+        : public Dune::CommDataHandleIF<MatPatternExchange,IdType> {
+        typedef typename Matrix::RowIterator RowIterator;
+        typedef typename Matrix::ColIterator ColIterator;
+    public:
+        //! Export type of data for message buffer
+        typedef IdType DataType;
+
+        /** @brief Returns true if data for given valid codim should be communicated
+         */
+        bool contains (int dim, int codim) const
+        {
+            return (codim==dofCodim);
+        }
+
+        /** @brief Returns true if size of data per entity of given dim and codim is a constant
+         */
+        bool fixedsize (int dim, int codim) const
+        {
+            return false;
+        }
+
+        /** @brief How many objects of type DataType have to be sent for a given entity
+         */
+        template<class EntityType>
+        size_t size (EntityType& e) const
+        {
+            int i = problem_.model().dofMapper().map(e);
+            int n = 0;
+            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+            {
+                typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
+                if (it != index2GID_.end())
+                    n++;
+            }
+
+            return n;
+        }
+
+        /** @brief Pack data from user to message buffer
+         */
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType& e) const
+        {
+            int i = problem_.model().dofMapper().map(e);
+            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+            {
+                typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
+                if (it != index2GID_.end())
+                    buff.write(it->second);
+            }
+
+        }
+
+        /** @brief Unpack data from message buffer to user
+         */
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
+        {
+            int i = problem_.model().dofMapper().map(e);
+            for (size_t k = 0; k < n; k++)
+            {
+                IdType id;
+                buff.read(id);
+                // only add entries corresponding to border entities
+                typename std::map<IdType,int>::const_iterator it = gid2Index_.find(id);
+                if (it != gid2Index_.end()
+                    && sparsity_[i].find(it->second) == sparsity_[i].end()
+                    && helper_.ghost(it->second) != 1<<24)
+                    sparsity_[i].insert(it->second);
+            }
+        }
+
+        /**
+         * @brief Get the communicated sparsity pattern
+         * @return the vector with the sparsity pattern
+         */
+        std::vector<std::set<int> >& sparsity ()
+        {
+            return sparsity_;
+        }
+
+        /** @brief Constructor
+            @param[in] gridView Grid view.
+            @param[in] g2i Global to local index map.
+            @param[in] i2g Local to global index map.
+            @param[in] A Matrix to operate on.
+            @param[in] helper parallel istl helper.
+        */
+        MatPatternExchange (const Problem& problem,
+                            const std::map<IdType,int>& g2i,
+                            const std::map<int,IdType>& i2g, Matrix& A,
+                            const ParallelISTLHelper<TypeTag>& helper)
+            : problem_(problem), gid2Index_(g2i), index2GID_(i2g),
+              sparsity_(A.N()), A_(A), helper_(helper)
+        {}
+
+    private:
+        const Problem& problem_;
+        const std::map<IdType,int>& gid2Index_;
+        const std::map<int,IdType>& index2GID_;
+        std::vector<std::set<int> > sparsity_;
+        Matrix& A_;
+        const ParallelISTLHelper<TypeTag>& helper_;
+    };
+
+    //! Local matrix blocks associated with the global id set
+    struct MatEntry
+    {
+        IdType first;
+        BlockType second;
+        MatEntry (const IdType& f, const BlockType& s) : first(f),second(s) {}
+        MatEntry () {}
+    };
+
+    //! A DataHandle class to exchange matrix entries
+    class MatEntryExchange
+        : public Dune::CommDataHandleIF<MatEntryExchange,MatEntry> {
+        typedef typename Matrix::RowIterator RowIterator;
+        typedef typename Matrix::ColIterator ColIterator;
+    public:
+        //! Export type of data for message buffer
+        typedef MatEntry DataType;
+
+        /** @brief Returns true if data for given valid codim should be communicated
+         */
+        bool contains (int dim, int codim) const
+        {
+            return (codim==dofCodim);
+        }
+
+        /** @brief Returns true if size of data per entity of given dim and codim is a constant
+         */
+        bool fixedsize (int dim, int codim) const
+        {
+            return false;
+        }
+
+        /** @brief How many objects of type DataType have to be sent for a given entity
+         */
+        template<class EntityType>
+        size_t size (EntityType& e) const
+        {
+            int i = problem_.model().dofMapper().map(e);
+            int n = 0;
+            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+            {
+                typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
+                if (it != index2GID_.end())
+                    n++;
+            }
+
+            return n;
+        }
+
+        /** @brief Pack data from user to message buffer
+         */
+        template<class MessageBuffer, class EntityType>
+        void gather (MessageBuffer& buff, const EntityType& e) const
+        {
+            int i = problem_.model().dofMapper().map(e);
+            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+            {
+                typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
+                if (it != index2GID_.end())
+                    buff.write(MatEntry(it->second,*j));
+            }
+
+        }
+
+        /** @brief Unpack data from message buffer to user
+         */
+        template<class MessageBuffer, class EntityType>
+        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
+        {
+            int i = problem_.model().dofMapper().map(e);
+            for (size_t k = 0; k < n; k++)
+            {
+                MatEntry m;
+                buff.read(m);
+                // only add entries corresponding to border entities
+                typename std::map<IdType,int>::const_iterator it = gid2Index_.find(m.first);
+                if (it != gid2Index_.end())
+                    if (A_[i].find(it->second) != A_[i].end())
+                        A_[i][it->second] += m.second;
+            }
+        }
+
+        /** @brief Constructor
+            @param[in] gridView Grid view.
+            @param[in] g2i Global to local index map.
+            @param[in] i2g Local to global index map.
+            @param[in] A Matrix to operate on.
+        */
+        MatEntryExchange (const Problem& problem, const std::map<IdType,int>& g2i,
+                          const std::map<int,IdType>& i2g,
+                          Matrix& A)
+            : problem_(problem), gid2Index_(g2i), index2GID_(i2g), A_(A)
+        {}
+
+    private:
+        const Problem& problem_;
+        const std::map<IdType,int>& gid2Index_;
+        const std::map<int,IdType>& index2GID_;
+        Matrix& A_;
+    };
+
+    /** @brief communicates values for the sparsity pattern of the new matrix.
+        @param A Matrix to operate on.
+        @param helper ParallelelISTLHelper.
+    */
+    void getextendedmatrix (Matrix& A,const ParallelISTLHelper<TypeTag>& helper)
+    {
+        if (problem_.model().gridView().comm().size() > 1) {
+            Matrix tmp(A);
+            std::size_t nnz=0;
+            // get entries from other processes
+            MatPatternExchange datahandle(problem_, gid2Index_, index2GID_, A, helper);
+            problem_.model().gridView().communicate(datahandle,
+                                                    Dune::InteriorBorder_InteriorBorder_Interface,
+                                                    Dune::ForwardCommunication);
+            std::vector<std::set<int> >& sparsity = datahandle.sparsity();
+            // add own entries, count number of nonzeros
+            for (RowIterator i = A.begin(); i != A.end(); ++i){
+                for (ColIterator j = A[i.index()].begin(); j != A[i.index()].end(); ++j){
+                    if (sparsity[i.index()].find(j.index()) == sparsity[i.index()].end())
+                        sparsity[i.index()].insert(j.index());
+                }
+                nnz += sparsity[i.index()].size();
+            }
+            A.setSize(tmp.N(), tmp.N(), nnz);
+            A.setBuildMode(Matrix::row_wise);
+            typename Matrix::CreateIterator citer = A.createbegin();
+            typedef typename std::vector<std::set<int> >::const_iterator Iter;
+            for (Iter i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer){
+                typedef typename std::set<int>::const_iterator SIter;
+                for (SIter si = i->begin(), send = i->end(); si!=send; ++si)
+                    citer.insert(*si);
+            }
+            // set matrix old values
+            A = 0;
+            for (RowIterator i = tmp.begin(); i != tmp.end(); ++i)
+                for (ColIterator j = tmp[i.index()].begin(); j != tmp[i.index()].end(); ++j){
+                    A[i.index()][j.index()] = tmp[i.index()][j.index()];
+                }
+        }
+    }
+
+    /** @brief Sums up the entries corresponding to border vertices.
+        @param A Matrix to operate on.
+    */
+    void sumEntries (Matrix& A)
+    {
+        if (problem_.model().gridView().comm().size() > 1)
+        {
+            MatEntryExchange datahandle(problem_.model().gridView(), gid2Index_, index2GID_, A);
+            problem_.model().gridView().communicate(datahandle,
+                                                    Dune::InteriorBorder_InteriorBorder_Interface,
+                                                    Dune::ForwardCommunication);
+        }
+    }
+
+private:
+    const Problem& problem_;
+    std::map<IdType,int> gid2Index_;
+    std::map<int,IdType> index2GID_;
+};
+
+#if HAVE_MPI
+
+void getextendedmatrix (M& A) const
+{
+    const GridView& gridView = problem_.model().gridView();
+    if (gridView.comm().size() > 1) {
+        Matrix tmp(A);
+        std::size_t nnz=0;
+        // get entries from other processes
+        MatPatternExchange datahandle(problem_, gid2Index_, index2GID_, A, *this);
+        gridView.communicate(datahandle,
+                             Dune::InteriorBorder_InteriorBorder_Interface,
+                             Dune::ForwardCommunication);
+        std::vector<std::set<int> >& sparsity = datahandle.sparsity();
+        // add own entries, count number of nonzeros
+        for (RowIterator i = A.begin(); i != A.end(); ++i){
+            for (ColIterator j = A[i.index()].begin(); j != A[i.index()].end(); ++j){
+                if (sparsity[i.index()].find(j.index()) == sparsity[i.index()].end())
+                    sparsity[i.index()].insert(j.index());
+            }
+            nnz += sparsity[i.index()].size();
+        }
+        A.setSize(tmp.N(), tmp.N(), nnz);
+        A.setBuildMode(Matrix::row_wise);
+        typename Matrix::CreateIterator citer = A.createbegin();
+        typedef typename std::vector<std::set<int> >::const_iterator Iter;
+        for (Iter i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer){
+            typedef typename std::set<int>::const_iterator SIter;
+            for (SIter si = i->begin(), send = i->end(); si!=send; ++si)
+                citer.insert(*si);
+        }
+        // set matrix old values
+        A = 0;
+        for (RowIterator i = tmp.begin(); i != tmp.end(); ++i)
+            for (ColIterator j = tmp[i.index()].begin(); j != tmp[i.index()].end(); ++j){
+                A[i.index()][j.index()] = tmp[i.index()][j.index()];
+            }
+    }
+}
+
+template<typename M, typename C>
+void ParallelISTLHelper<GFS,skipBlocksizeCheck>::createIndexSetAndProjectForAMG(M& m, C& c)
+{
+    const GridView& gridview = problem_.model().gridView();
+
+    // First find out which dofs we share with other processors
+    std::vector<int> sharedDofs(problem_.model().dofMapper().size(), false);
+
+    SharedGatherScatter sgs(sharedDofs, problem_);
+
+    if (gridview.comm().size()>1)
+        gridview.communicate(sgs,Dune::All_All_Interface,
+                             Dune::ForwardCommunication);
+
+    // Count shared dofs that we own
+    typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
+    GlobalIndex count=0;
+
+    for(auto v=sharedDofs.begin(), vend=sharedDofs.end(); v != vend; ++v)
+        if(*v)
+            ++count;
+
+    dverb<<gridview.comm().rank()<<": shared count is "<< count.touint()
+         <<std::endl;
+
+    dverb<<gridview.comm().rank()<<": shared block count is "
+         << count.touint()<<std::endl;
+
+
+    std::vector<GlobalIndex> counts(gridview.comm().size());
+    gridview.comm().allgather(&count, 1, &(counts[0]));
+
+    // Compute start index start_p = \sum_{i=0}^{i<p} counts_i
+    GlobalIndex start=0;
+    for(int i=0; i<gridview.comm().rank(); ++i)
+        start=start+counts[i];
+    //std::cout<<gv.comm().rank()<<": start index = "<<start.touint()<<std::endl;
+
+
+    typedef std::vector<GlobalIndex> GIVector;
+    GIVector scalarIndices(problem_.model().dofMapper().size(),
+                           std::numeric_limits<GlobalIndex>::max());
+
+    auto shared=sharedDofs.begin();
+    auto index=scalarIndices.begin();
+
+    for(auto i=vi.begin(), ien=vi.end(); i!=vi; ++i, ++shared)
+        if(*i==1.0 && *shared){
+            *index=start;
+        }
+
+    // publish global indices for the shared DOFS to other processors.
+    typedef GlobalIndexGatherScatter<GlobalIndex> GIGS;
+    GIGS gigs(scalarIndices, problem_);
+    if (gridview.comm().size()>1)
+        gridview.communicate(gigs,Dune::All_All_Interface,
+                             Dune::ForwardCommunication);
+
+
+    // Setup the index set
+    c.indexSet().beginResize();
+    index=scalarIndices.begin();
+    auto ghost=isGhost.begin();
+
+    for(auto i=v.begin(), iend=vend(); i!=iend; ++i)
+    {
+        Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
+        if(*index!=std::numeric_limits<GlobalIndex>::max()){
+            // global index exist in index set
+            if(*i>0){
+                // This dof is managed by us.
+                attr = Dune::OwnerOverlapCopyAttributeSet::owner;
+            }
+            else if ( *g==(1<<24) && ( c.getSolverCategory() ==
+                                       static_cast<int>(SolverCategory::nonoverlapping)) ){
+                //use attribute overlap for ghosts in novlp grids
+                attr = Dune::OwnerOverlapCopyAttributeSet::overlap;
+            }
+            else {
+                attr = Dune::OwnerOverlapCopyAttributeSet::copy;
+            }
+            c.indexSet().add(index, typename C::ParallelIndexSet::LocalIndex(i, attr)));
+    }
+    c.indexSet().endResize();
+    //std::cout<<gv.comm().rank()<<": index set size = "<<c.indexSet().size()<<std::endl;
+    //std::cout<<gv.comm().rank()<<": "<<c.indexSet()<<std::endl;
+
+    // Compute neighbours using communication
+    typedef NeighbourGatherScatter NeighbourGS;
+    std::set<int> neighbours;
+    NeighbourGatherScatter ngs(gridview.comm().rank(),
+                               neighbours);
+
+    if (gridview.comm().size()>1)
+        gridview.communicate(ngs,Dune::All_All_Interface,
+                             Dune::ForwardCommunication);
+    c.remoteIndices().setNeighbours(neighbours);
+    //std::cout<<gv.comm().rank()<<": no neighbours="<<neighbours.size()<<std::endl;
+
+    c.remoteIndices().template rebuild<false>();
+    //std::cout<<c.remoteIndices()<<std::endl;
+
+}
+#endif
+}
+#endif

dumux/linear/amgproperties.hh

+// -*- 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 2 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 DUMUXAMGPROPERTIES_HH
+#define DUMUXAMGPROPERTIES_HH
+
+
+#include <dune/istl/schwarz.hh>
+#include <dune/istl/novlpschwarz.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/preconditioners.hh>
+
+#include <dumux/implicit/box/boxproperties.hh>
+#include <dumux/implicit/cellcentered/ccproperties.hh>
+#include <dumux/decoupled/common/decoupledproperties.hh>
+#include <dumux/decoupled/common/pressureproperties.hh>
+#include "linearsolverproperties.hh"
+
+namespace Dumux
+{
+
+// forward declaration for the property definitions
+template <class TypeTag> class AMGBackend;
+
+namespace Properties
+{
+//! the PDELab finite element map used for the gridfunctionspace
+NEW_PROP_TAG(AMGLocalFemMap);
+
+//! box: use the (multi-)linear local FEM space associated with cubes by default
+SET_PROP(BoxModel, AMGLocalFemMap)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    enum{dim = GridView::dimension};
+ public:
+    enum{
+        //! \brief The codimension that the degrees of freedom are attached to.
+        dof_codimension = GridView::dimension
+    };
+};
+
+//! cell-centered: use the element-wise constant local FEM space by default
+SET_PROP(CCModel, AMGLocalFemMap)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    enum{dim = GridView::dimension};
+ public:
+    enum{
+        //! \brief The codimension that the degrees of freedom are attached to.
+        dof_codimension = 0
+    };
+};
+
+//! decoupled models: use the element-wise constant local FEM space by default
+SET_PROP(DecoupledModel, AMGLocalFemMap)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    enum{dim = GridView::dimension};
+ public:
+    enum{
+        //! \brief The codimension that the degrees of freedom are attached to.
+        dof_codimension = 0
+    };
+};
+/*
+//! set a property JacobianMatrix also for the decoupled models
+SET_PROP(DecoupledModel, JacobianMatrix)
+{
+public:
+typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) type;
+};
+*/
+ //! the type of the employed PDELab backend
+NEW_PROP_TAG(AMGPDELabBackend);
+
+//! box: use the non-overlapping AMG backend
+SET_PROP(BoxModel, AMGPDELabBackend)
+{
+ public:
+    //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_AMG_SSOR<GridOperator> type;
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum {
+        numEq = JacobianMatrix::block_type::rows,
+        isNonOverlapping = true
+    };
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
+#if HAVE_MPI
+    typedef Dune::OwnerOverlapCopyCommunication<bigunsignedint<96>,int> Comm;
+    typedef Dune::NonoverlappingSchwarzOperator<MType,VType, VType> LinearOperator;
+    typedef Dune::NonoverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
+    typedef Dune::NonoverlappingBlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
+#else
+    typedef Dune::Amg::SequentialInformation Comm;
+    typedef Dune::MatrixAdapter<MType,VType,VType> LinearOperator;
+    typedef Dune::SeqScalarProduct<VType> ScalarProduct;
+    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+#endif
+};
+
+//! cell-centered: use the overlapping PDELab AMG backend
+SET_PROP(CCModel, AMGPDELabBackend)
+{
+ public:
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum {
+        numEq = JacobianMatrix::block_type::rows,
+        isNonOverlapping = true
+    };
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
+#if HAVE_MPI
+    typedef Dune::OwnerOverlapCopyCommunication<bigunsignedint<96>,int> Comm;
+    typedef Dune::OverlappingSchwarzOperator<MType,VType, VType> LinearOperator;
+    typedef Dune::OverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
+    typedef Dune::BlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
+#else
+    typedef Dune::Amg::SequentialInformation Comm;
+    typedef Dune::MatrixAdapter<MType,VType,VType> LinearOperator;
+    typedef Dune::SeqScalarProduct<VType> ScalarProduct;
+    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+#endif
+};
+
+//! decoupled model: use the overlapping PDELab AMG backend
+SET_PROP(DecoupledModel, AMGPDELabBackend)
+{
+    typedef typename Dumux::AMGBackend<TypeTag>::GridOperator GridOperator;
+ public:
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum {
+        numEq = JacobianMatrix::block_type::rows,
+        isNonOverlapping = false
+    };
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> > MType;
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > VType;
+#if HAVE_MPI
+    typedef Dune::OwnerOverlapCopyCommunication<bigunsignedint<96>,int> Comm;
+    typedef Dune::OverlappingSchwarzOperator<MType,VType, VType> LinearOperator;
+    typedef Dune::OverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
+    typedef Dune::BlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
+#else
+    typedef Dune::Amg::SequentialInformation Comm;
+    typedef Dune::MatrixAdapter<MType,VType,VType> LinearOperator;
+    typedef Dune::SeqScalarProduct<VType> ScalarProduct;
+    typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+#endif
+};
+
+//! box: reset the type of solution vector to be PDELab conforming
+SET_PROP(BoxModel, SolutionVector)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > type;
+};
+
+//! cell-centered: reset the type of solution vector to be PDELab conforming
+SET_PROP(CCModel, SolutionVector)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > type;
+};
+
+
+//! decoupled model: reset the type of solution vector to be PDELab conforming
+SET_PROP(DecoupledModel, PressureSolutionVector)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > type;
+};
+
+//! decoupled model: reset the type of solution vector to be PDELab conforming
+SET_PROP(DecoupledModel, PressureRHSVector)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    enum { numEq = JacobianMatrix::block_type::rows};
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+    typedef Dune::BlockVector<Dune::FieldVector<Scalar,numEq> > type;
+};
+
+} // namespace Properties
+} // namespace Dumux
+#endif

test/decoupled/2p/test_impesproblem.hh

@@ -49,14 +49,7 @@
 
 #include<dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh>
 
-#if HAVE_DUNE_PDELAB
-
-// Check if DUNE-PDELab has been patched for our needs.
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 #include <dumux/linear/amgbackend.hh>
-#endif // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#endif // HAVE_DUNE_PDELAB
 
 namespace Dumux
 {
@@ -134,7 +127,6 @@ SET_TYPE_PROP(IMPESTestProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<Typ
 
 SET_SCALAR_PROP(IMPESTestProblem, ImpetCFLFactor, 0.95);
 
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 // set up an additional problem where the AMG backend is used
 NEW_TYPE_TAG(IMPESTestProblemWithAMG, INHERITS_FROM(IMPESTestProblem));
 // use the AMG backend for the corresponding test
@@ -143,7 +135,6 @@ SET_TYPE_PROP(IMPESTestProblemWithAMG, LinearSolver, Dumux::AMGBackend<TypeTag>)
 SET_TYPE_PROP(IMPESTestProblemWithAMG, Grid, Dune::YaspGrid<2>);
 // set the GridCreator property
 SET_TYPE_PROP(IMPESTestProblemWithAMG, GridCreator, Dumux::DgfGridCreator<TypeTag>);
-#endif
 }
 
 /*!

test/geomechanics/el1p2c/el1p2cproblem.hh

@@ -86,14 +86,8 @@ namespace Dumux
     SET_BOOL_PROP(El1P2CProblem, ImplicitWithStabilization, true);
 
 // Check if DUNE-PDELab is available and has been patched for our needs.
-#if HAVE_DUNE_PDELAB &&  DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
     SET_TYPE_PROP(El1P2CProblem, LinearSolver, Dumux::AMGBackend<TypeTag> );
   // #endif // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-#elif HAVE_SUPERLU
-    SET_TYPE_PROP(El1P2CProblem, LinearSolver, Dumux::SuperLUBackend<TypeTag> );
-#else
-#warning Neither DUNE_PDELAB nor SuperLU detected, you are not able to run this test!
-#endif // HAVE_DUNE_PDELAB
 
 }
 

test/implicit/1p/1ptestproblem.hh

@@ -40,14 +40,7 @@
 #include <dumux/material/components/simpleh2o.hh>
 #include <dumux/material/fluidsystems/liquidphase.hh>
 
-#if HAVE_DUNE_PDELAB
-
-// Check if DUNE-PDELab has been patched for our needs.
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 #include <dumux/linear/amgbackend.hh>
-#endif // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#endif // HAVE_DUNE_PDELAB
 
 #include "1ptestspatialparams.hh"
 
@@ -93,13 +86,11 @@ SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 0);
 SET_INT_PROP(OnePTestProblem, LinearSolverPreconditionerIterations, 1);
 SET_SCALAR_PROP(OnePTestProblem, LinearSolverPreconditionerRelaxation, 1.0);
 
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 NEW_TYPE_TAG(OnePTestBoxProblemWithAMG, INHERITS_FROM(OnePTestBoxProblem));
 NEW_TYPE_TAG(OnePTestCCProblemWithAMG, INHERITS_FROM(OnePTestCCProblem));
 // Solver settings for the tests using AMG
 SET_TYPE_PROP(OnePTestBoxProblemWithAMG, LinearSolver, Dumux::AMGBackend<TypeTag> );
 SET_TYPE_PROP(OnePTestCCProblemWithAMG, LinearSolver, Dumux::AMGBackend<TypeTag> );
-#endif
 
 // Enable gravity
 SET_BOOL_PROP(OnePTestProblem, ProblemEnableGravity, true);

test/implicit/1p/test_box1pwithamg.cc

@@ -23,10 +23,8 @@
  */
 #include "config.h"
 
-#if HAVE_DUNE_PDELAB
 
 // Check if DUNE-PDELab has been patched for our needs.
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 
 #include "1ptestproblem.hh"
 #include <dumux/common/start.hh>
@@ -68,30 +66,3 @@ int main(int argc, char** argv)
     typedef TTAG(OnePTestBoxProblemWithAMG) ProblemTypeTag;
     return Dumux::start<ProblemTypeTag>(argc, argv, usage);
 }
-#else // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#warning You need to have  a patched dune-pdelab to run this test, see ../../../patches/README for details.
-
-#include <iostream>
-
-int main()
-{
-    std::cerr << "You need to have a patched dune-pdelab to run this test, "
-                 "see ../../../patches/README for details." << std::endl;
-    return 77;
-}
-
-#endif // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#else // HAVE_DUNE_PDELAB
-
-#warning You need to have dune-pdelab installed and patched to run this test.
-
-#include <iostream>
-
-int main()
-{
-    std::cerr << "You need to have dune-pdelab installed and patched to run this test.\n";
-    return 77;
-}
-#endif // HAVE_DUNE_PDELAB

test/implicit/1p/test_cc1pwithamg.cc

@@ -23,10 +23,6 @@
  */
 #include "config.h"
 
-#if HAVE_DUNE_PDELAB
-
-// Check if DUNE-PDELab has been patched for our needs.
-#ifdef DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
 
 #include "1ptestproblem.hh"
 #include <dumux/common/start.hh>
@@ -68,30 +64,4 @@ int main(int argc, char** argv)
     typedef TTAG(OnePTestCCProblemWithAMG) ProblemTypeTag;
     return Dumux::start<ProblemTypeTag>(argc, argv, usage);
 }
-#else // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#warning You need to have  a patched dune-pdelab to run this test, see ../../../patches/README for details.
-
-#include <iostream>
-
-int main()
-{
-    std::cerr << "You need to have a patched dune-pdelab to run this test, "
-                 "see ../../../patches/README for details." << std::endl;
-    return 77;
-}
 
-#endif // DUNE_PDELAB_IS_PATCHED_FOR_DUMUX
-
-#else // HAVE_DUNE_PDELAB
-
-#warning You need to have dune-pdelab installed and patched to run this test.
-
-#include <iostream>
-
-int main()
-{
-    std::cerr << "You need to have dune-pdelab installed and patched to run this test.\n";
-    return 77;
-}
-#endif // HAVE_DUNE_PDELAB