diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh
index bc79e5d4738fe32bfaa8ce0225257ef340bd4b02..b6d0ee62d488351e41d5da166d2ef95db51161ed 100644
--- a/dumux/linear/amgbackend.hh
+++ b/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
diff --git a/dumux/linear/amgparallelhelpers.hh b/dumux/linear/amgparallelhelpers.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fa9b1b20a7198d4bfce8c51d2c61c13e40d1e50b
--- /dev/null
+++ b/dumux/linear/amgparallelhelpers.hh
@@ -0,0 +1,819 @@
+// -*- 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
diff --git a/dumux/linear/amgproperties.hh b/dumux/linear/amgproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cfc3c32aba96edbc99bafa90708a7f8af2aef926
--- /dev/null
+++ b/dumux/linear/amgproperties.hh
@@ -0,0 +1,212 @@
+// -*- 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
diff --git a/test/decoupled/2p/test_impesproblem.hh b/test/decoupled/2p/test_impesproblem.hh
index 5cc6eda7116209ee749a8ba5551ef2aa905a3ac3..672cbd279c2ae68c731f8f9b5c600575eac0f24e 100644
--- a/test/decoupled/2p/test_impesproblem.hh
+++ b/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
}
/*!
diff --git a/test/geomechanics/el1p2c/el1p2cproblem.hh b/test/geomechanics/el1p2c/el1p2cproblem.hh
index 502b7c52cf2e1c441664e01c55f5dac81a03aaa8..6cbb57fc44518655edcf50f75ccd6a0a2bdd2b2e 100644
--- a/test/geomechanics/el1p2c/el1p2cproblem.hh
+++ b/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
}
diff --git a/test/implicit/1p/1ptestproblem.hh b/test/implicit/1p/1ptestproblem.hh
index 186c8da7ec76fd743653ac80fc30319521d8f610..8a1488f5a0c4e90f19c8a5bb19dcd91660fdfa81 100644
--- a/test/implicit/1p/1ptestproblem.hh
+++ b/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);
diff --git a/test/implicit/1p/test_box1pwithamg.cc b/test/implicit/1p/test_box1pwithamg.cc
index f05fa5fb06f2957732729086fc26a28935375e49..8ee1b7a01137908b55ddc6a10594ddf0ca4c0a68 100644
--- a/test/implicit/1p/test_box1pwithamg.cc
+++ b/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
diff --git a/test/implicit/1p/test_cc1pwithamg.cc b/test/implicit/1p/test_cc1pwithamg.cc
index e3c30bc09198a3e9b934e3a8c38b4242d7af3ba3..30494c06b3cffe83e1579391eb2e179c0265ece2 100644
--- a/test/implicit/1p/test_cc1pwithamg.cc
+++ b/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