diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh
index 672987a3dd4ab7f352198900eaf6bfbdb2da915a..104da6a78e260178fa0c819555f552cf2a61cda9 100644
--- a/dumux/linear/amgbackend.hh
+++ b/dumux/linear/amgbackend.hh
@@ -74,6 +74,7 @@ template <class TypeTag>
class AMGBackend
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP(TypeTag, AmgTraits) AmgTraits;
@@ -111,41 +112,18 @@ public:
int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity);
static const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction);
-#if HAVE_MPI
- Dune::SolverCategory::Category category = AmgTraits::isNonOverlapping?
- Dune::SolverCategory::nonoverlapping : Dune::SolverCategory::overlapping;
-
- if(AmgTraits::isNonOverlapping && firstCall_)
- {
- phelper_.initGhostsAndOwners();
- }
-
- typename AmgTraits::Comm comm(problem_.gridView().comm(), category);
-
- if(AmgTraits::isNonOverlapping)
- {
- // extend the matrix pattern such that it is usable for AMG
- EntityExchanger<TypeTag> exchanger(problem_);
- exchanger.getExtendedMatrix(A, phelper_);
- exchanger.sumEntries(A);
- }
- phelper_.createIndexSetAndProjectForAMG(A, comm);
-
- typename AmgTraits::LinearOperator fop(A, comm);
- typename AmgTraits::ScalarProduct sp(comm);
- int rank = comm.communicator().rank();
-
- // Make rhs consistent
- if(AmgTraits::isNonOverlapping)
- {
- phelper_.makeNonOverlappingConsistent(b);
- }
+ int rank = 0;
+ std::shared_ptr<typename AmgTraits::Comm> comm;
+ std::shared_ptr<typename AmgTraits::LinearOperator> fop;
+ std::shared_ptr<typename AmgTraits::ScalarProduct> sp;
+#if DUNE_VERSION_NEWER(DUNE_GRID, 3, 0)
+ static const int dofCodim = AmgTraits::dofCodim;
+ static const bool isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v;
#else
- typename AmgTraits::Comm comm;
- typename AmgTraits::LinearOperator fop(A);
- typename AmgTraits::ScalarProduct sp;
- int rank=0;
+ static const bool isParallel = Dune::Capabilities::isParallel<Grid>::v;
#endif
+ prepareLinearAlgebra_<Matrix, Vector, isParallel>(A, b, rank, comm, fop, sp);
+
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments
SmootherArgs;
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,
@@ -160,15 +138,14 @@ public:
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = 1;
- AMGType amg(fop, criterion, smootherArgs, comm);
- Dune::BiCGSTABSolver<typename AmgTraits::VType> solver(fop, sp, amg, residReduction, maxIt,
- rank==0?verbosity: 0);
-
+ AMGType amg(*fop, criterion, smootherArgs, *comm);
+ Dune::BiCGSTABSolver<typename AmgTraits::VType> solver(*fop, *sp, amg, residReduction, maxIt,
+ rank == 0 ? verbosity : 0);
solver.apply(x, b, result_);
firstCall_ = false;
return result_.converged;
- }
+ }
/*!
* \brief The result containing the convergence history.
@@ -182,7 +159,35 @@ public:
{
return problem_;
}
+
private:
+
+ /*!
+ * \brief Prepare the linear algebra member variables.
+ *
+ * At compile time, correct constructor calls have to be chosen,
+ * depending on whether the setting is parallel or sequential.
+ * Since several template parameters are present, this cannot be solved
+ * by a full function template specialization. Instead, class template
+ * specialization has to be used.
+ * This adapts example 4 from http://www.gotw.ca/publications/mill17.htm.
+
+ * The function is called from the solve function. The call is
+ * forwarded to the corresponding function of a class template.
+ *
+ * \tparam isParallel decides if the setting is parallel or sequential
+ */
+ template<class Matrix, class Vector, bool isParallel>
+ void prepareLinearAlgebra_(Matrix& A, Vector& b, int& rank,
+ std::shared_ptr<typename AmgTraits::Comm>& comm,
+ std::shared_ptr<typename AmgTraits::LinearOperator>& fop,
+ std::shared_ptr<typename AmgTraits::ScalarProduct>& sp)
+ {
+ LinearAlgebraPreparator<TypeTag, isParallel>
+ ::prepareLinearAlgebra(A, b, rank, comm, fop, sp,
+ problem_, phelper_, firstCall_);
+ }
+
const Problem& problem_;
ParallelISTLHelper<TypeTag> phelper_;
Dune::InverseOperatorResult result_;
diff --git a/dumux/linear/amgparallelhelpers.hh b/dumux/linear/amgparallelhelpers.hh
index 8416a8decec7242370118a3aff089728fca410bd..ee7b1139397b43cc4e44a8ebeffc00e0141d407e 100644
--- a/dumux/linear/amgparallelhelpers.hh
+++ b/dumux/linear/amgparallelhelpers.hh
@@ -975,6 +975,102 @@ void ParallelISTLHelper<TypeTag>::createIndexSetAndProjectForAMG(M& m, C& c)
}
#endif // HAVE_MPI
+
+/*!
+ * \brief Prepare the linear algebra member variables.
+ *
+ * At compile time, correct constructor calls have to be chosen,
+ * depending on whether the setting is parallel or sequential.
+ * Since several template parameters are present, this cannot be solved
+ * by a full function template specialization. Instead, class template
+ * specialization has to be used.
+ * This adapts example 4 from http://www.gotw.ca/publications/mill17.htm.
+ *
+ * This class template implements the function for the sequential case.
+ *
+ * \tparam isParallel decides if the setting is parallel or sequential
+ */
+template<class TypeTag, bool isParallel>
+struct LinearAlgebraPreparator
+{
+ using AmgTraits = typename GET_PROP(TypeTag, AmgTraits);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using ParallelHelper = ParallelISTLHelper<TypeTag>;
+ using Comm = typename AmgTraits::Comm;
+ using LinearOperator = typename AmgTraits::LinearOperator;
+ using ScalarProduct = typename AmgTraits::ScalarProduct;
+
+ template<class Matrix, class Vector>
+ static void prepareLinearAlgebra(Matrix& A, Vector& b,
+ int& rank,
+ std::shared_ptr<Comm>& comm,
+ std::shared_ptr<LinearOperator>& fop,
+ std::shared_ptr<ScalarProduct>& sp,
+ const Problem& problem,
+ ParallelHelper& pHelper,
+ const bool firstCall)
+ {
+ comm = std::make_shared<Comm>();
+ fop = std::make_shared<LinearOperator>(A);
+ sp = std::make_shared<ScalarProduct>();
+ }
+};
+
+/*!
+ * \brief Specialization for the parallel case.
+ */
+template<class TypeTag>
+struct LinearAlgebraPreparator<TypeTag, true>
+{
+ using AmgTraits = typename GET_PROP(TypeTag, AmgTraits);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using ParallelHelper = ParallelISTLHelper<TypeTag>;
+ using Comm = typename AmgTraits::Comm;
+ using LinearOperator = typename AmgTraits::LinearOperator;
+ using ScalarProduct = typename AmgTraits::ScalarProduct;
+
+ template<class Matrix, class Vector>
+ static void prepareLinearAlgebra(Matrix& A, Vector& b,
+ int& rank,
+ std::shared_ptr<Comm>& comm,
+ std::shared_ptr<LinearOperator>& fop,
+ std::shared_ptr<ScalarProduct>& sp,
+ const Problem& problem,
+ ParallelHelper& pHelper,
+ const bool firstCall)
+ {
+ Dune::SolverCategory::Category category = AmgTraits::isNonOverlapping ?
+ Dune::SolverCategory::nonoverlapping : Dune::SolverCategory::overlapping;
+
+ if(AmgTraits::isNonOverlapping && firstCall)
+ {
+ pHelper.initGhostsAndOwners();
+ }
+
+ comm = std::make_shared<Comm>(problem.gridView().comm(), category);
+
+ if(AmgTraits::isNonOverlapping)
+ {
+ // extend the matrix pattern such that it is usable for AMG
+ EntityExchanger<TypeTag> exchanger(problem);
+ exchanger.getExtendedMatrix(A, pHelper);
+ exchanger.sumEntries(A);
+ }
+ pHelper.createIndexSetAndProjectForAMG(A, *comm);
+
+ fop = std::make_shared<LinearOperator>(A, *comm);
+ sp = std::make_shared<ScalarProduct>(*comm);
+ rank = comm->communicator().rank();
+
+ // Make rhs consistent
+ if(AmgTraits::isNonOverlapping)
+ {
+ pHelper.makeNonOverlappingConsistent(b);
+ }
+ }
+};
+
+
} // end namespace Dumux
#endif // DUMUX_AMGPARALLELHELPERS_HH
diff --git a/dumux/linear/amgproperties.hh b/dumux/linear/amgproperties.hh
index 9723d3d4f0f9cd6ce44341b214f23523aa0e60c8..a0526b4aefa863610e6c4b0f2f4a2cddfd652d6d 100644
--- a/dumux/linear/amgproperties.hh
+++ b/dumux/linear/amgproperties.hh
@@ -32,6 +32,7 @@
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/preconditioners.hh>
+#include <dune/grid/common/capabilities.hh>
#include <dumux/implicit/box/boxproperties.hh>
#include <dumux/implicit/cellcentered/ccproperties.hh>
@@ -50,31 +51,70 @@ namespace Properties
//! The type traits required for using the AMG backend
NEW_PROP_TAG(AmgTraits);
+template <class MType, class VType, bool isParallel>
+class NonoverlappingSolverTraits
+{
+public:
+ typedef Dune::Amg::SequentialInformation Comm;
+ typedef Dune::MatrixAdapter<MType,VType,VType> LinearOperator;
+ typedef Dune::SeqScalarProduct<VType> ScalarProduct;
+ typedef Dune::SeqSSOR<MType,VType, VType> Smoother;
+};
+
+template <class MType, class VType>
+class NonoverlappingSolverTraits<MType, VType, true>
+{
+public:
+ typedef Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int> Comm;
+ typedef Dune::NonoverlappingSchwarzOperator<MType,VType, VType,Comm> LinearOperator;
+ typedef Dune::NonoverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
+ typedef Dune::NonoverlappingBlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
+};
+
//! Box: use the non-overlapping AMG
SET_PROP(BoxModel, AmgTraits)
{
public:
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
enum {
numEq = JacobianMatrix::block_type::rows,
- dofCodim = GridView::dimension,
- isNonOverlapping = true
+ dofCodim = Grid::dimension,
+ isNonOverlapping = true,
+#if DUNE_VERSION_NEWER(DUNE_GRID, 3, 0)
+ isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
+#else
+ isParallel = Dune::Capabilities::isParallel<Grid>::v
+#endif
};
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<Dune::bigunsignedint<96>,int> Comm;
- typedef Dune::NonoverlappingSchwarzOperator<MType,VType, VType,Comm> LinearOperator;
- typedef Dune::NonoverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
- typedef Dune::NonoverlappingBlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
-#else
+ typedef NonoverlappingSolverTraits<MType, VType, isParallel> SolverTraits;
+ typedef typename SolverTraits::Comm Comm;
+ typedef typename SolverTraits::LinearOperator LinearOperator;
+ typedef typename SolverTraits::ScalarProduct ScalarProduct;
+ typedef typename SolverTraits::Smoother Smoother;
+};
+
+template <class MType, class VType, bool isParallel>
+class OverlappingSolverTraits
+{
+public:
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
+};
+
+template <class MType, class VType>
+class OverlappingSolverTraits<MType, VType, true>
+{
+public:
+ typedef Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int> Comm;
+ typedef Dune::OverlappingSchwarzOperator<MType,VType, VType,Comm> LinearOperator;
+ typedef Dune::OverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
+ typedef Dune::BlockPreconditioner<VType,VType,Comm,Dune::SeqSSOR<MType,VType, VType> > Smoother;
};
//! Cell-centered: use the overlapping AMG
@@ -82,25 +122,25 @@ SET_PROP(CCModel, AmgTraits)
{
public:
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
enum {
numEq = JacobianMatrix::block_type::rows,
dofCodim = 0,
- isNonOverlapping = false
+ isNonOverlapping = false,
+#if DUNE_VERSION_NEWER(DUNE_GRID, 3, 0)
+ isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
+#else
+ isParallel = Dune::Capabilities::isParallel<Grid>::v
+#endif
};
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<Dune::bigunsignedint<96>,int> Comm;
- typedef Dune::OverlappingSchwarzOperator<MType,VType, VType,Comm> LinearOperator;
- typedef Dune::OverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
- typedef Dune::BlockPreconditioner<VType,VType,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
+ typedef OverlappingSolverTraits<MType, VType, isParallel> SolverTraits;
+ typedef typename SolverTraits::Comm Comm;
+ typedef typename SolverTraits::LinearOperator LinearOperator;
+ typedef typename SolverTraits::ScalarProduct ScalarProduct;
+ typedef typename SolverTraits::Smoother Smoother;
};
//! Decoupled model: use the overlapping AMG
@@ -108,25 +148,25 @@ SET_PROP(DecoupledModel, AmgTraits)
{
public:
typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
enum {
numEq = JacobianMatrix::block_type::rows,
dofCodim = 0,
- isNonOverlapping = false
+ isNonOverlapping = false,
+#if DUNE_VERSION_NEWER(DUNE_GRID, 3, 0)
+ isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
+#else
+ isParallel = Dune::Capabilities::isParallel<Grid>::v
+#endif
};
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<Dune::bigunsignedint<96>,int> Comm;
- typedef Dune::OverlappingSchwarzOperator<MType,VType, VType,Comm> LinearOperator;
- typedef Dune::OverlappingSchwarzScalarProduct<VType,Comm> ScalarProduct;
- typedef Dune::BlockPreconditioner<VType,VType,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
+ typedef OverlappingSolverTraits<MType, VType, isParallel> SolverTraits;
+ typedef typename SolverTraits::Comm Comm;
+ typedef typename SolverTraits::LinearOperator LinearOperator;
+ typedef typename SolverTraits::ScalarProduct ScalarProduct;
+ typedef typename SolverTraits::Smoother Smoother;
};
} // namespace Properties