Commit 21c2b3d8 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[AMGBackend] fix AMGBackend for sequential grids and HAVE_MPI == 1

So far, the AMGBackend only worked for sequential grids if HAVE_MPI
was 0. This was due to the fact that AmgTraits distinguished between
parallel and sequential by evaluating HAVE_MPI. If HAVE_MPI was 1
but the Grid was sequential, this lead to compiler errors for the
constructor calls of the Communicator, LinearOperator and
ScalarProduct.

Do this properly now by using Dune::Capabilities for distinction
between the sequential and parallel case. For Dune 2.4, the
capability isParallel has to be used, while for Dune 3.0, it is
canCommunicate. The correct AmgTraits and constructor calls are
chosen by template specialization.

This solves FS#300.
parent b9e5d439
......@@ -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_;
......
......@@ -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
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment