Support DUNE 2.5 with CPR reusing hierarchy

This commit is contained in:
Markus Blatt 2019-04-10 14:40:11 +02:00
parent 29e893dc8e
commit 81b7af1d8a
2 changed files with 47 additions and 16 deletions

View File

@ -149,7 +149,7 @@ namespace Opm
constexpr int category = Dune::SolverCategory::overlapping; constexpr int category = Dune::SolverCategory::overlapping;
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser; typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer; typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
SPPointer sp(ScalarProductChooser::construct(info)); SPPointer sp(ScalarProductChooser::construct(*comm_));
sp_ = std::move(sp); sp_ = std::move(sp);
#endif #endif
@ -227,6 +227,9 @@ namespace Opm
// Since DUNE 2.2 we also need to pass the smoother args instead of steps directly // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
using AmgType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value, using AmgType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
BlackoilAmgType, ParallelBlackoilAmgType>::type; BlackoilAmgType, ParallelBlackoilAmgType>::type;
using SpType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
Dune::SeqScalarProduct<Vector>,
Dune::OverlappingSchwarzScalarProduct<Vector, Comm> >::type;
using OperatorType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value, using OperatorType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
MatrixAdapter, ParallelMatrixAdapter>::type; MatrixAdapter, ParallelMatrixAdapter>::type;
typedef typename AmgType::Smoother Smoother; typedef typename AmgType::Smoother Smoother;
@ -273,7 +276,7 @@ namespace Opm
verbosity_linsolve = this->parameters_.linear_solver_verbosity_; verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
} }
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, *sp_, *amg_, linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, reinterpret_cast<SpType&>(*sp_), reinterpret_cast<AmgType&>(*amg_),
this->parameters_.linear_solver_reduction_, this->parameters_.linear_solver_reduction_,
this->parameters_.linear_solver_maxiter_, this->parameters_.linear_solver_maxiter_,
verbosity_linsolve)); verbosity_linsolve));
@ -295,11 +298,6 @@ namespace Opm
protected: protected:
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
typedef std::shared_ptr< Dune::ScalarProduct<Vector> > SPPointer;
#else
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
#endif
///! \brief The dune-istl operator (either serial or parallel ///! \brief The dune-istl operator (either serial or parallel
std::unique_ptr< Dune::LinearOperator<Vector, Vector> > opA_; std::unique_ptr< Dune::LinearOperator<Vector, Vector> > opA_;
///! \brief Serial well matrix adapter ///! \brief Serial well matrix adapter
@ -309,11 +307,12 @@ namespace Opm
///! \brief The preconditoner to use (either serial or parallel CPR with AMG) ///! \brief The preconditoner to use (either serial or parallel CPR with AMG)
std::unique_ptr< Preconditioner > amg_; std::unique_ptr< Preconditioner > amg_;
using SPPointer = std::shared_ptr< Dune::ScalarProduct<Vector> >;
SPPointer sp_; SPPointer sp_;
std::shared_ptr< Dune::BiCGSTABSolver<Vector> > linsolve_; std::shared_ptr< Dune::BiCGSTABSolver<Vector> > linsolve_;
const void* oldMat; const void* oldMat;
using POrComm = Dune::OwnerOverlapCopyCommunication<int,int>; using POrComm = Dune::OwnerOverlapCopyCommunication<int,int>;
std::shared_ptr<POrComm> comm_; std::shared_ptr<POrComm> comm_;
}; // end ISTLSolver }; // end ISTLSolver
} // namespace Opm } // namespace Opm

View File

@ -137,11 +137,19 @@ namespace Dune
/** \copydoc Preconditioner::apply */ /** \copydoc Preconditioner::apply */
void apply(Domain& v, const Range& d); void apply(Domain& v, const Range& d);
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
//! Category of the preconditioner (see SolverCategory::Category) //! Category of the preconditioner (see SolverCategory::Category)
virtual SolverCategory::Category category() const virtual SolverCategory::Category category() const
{ {
return category_; return category_;
} }
#else
enum {
//! \brief The category the preconditioner is part of.
category = std::is_same<PI,Dune::Amg::SequentialInformation>::value?
Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
};
#endif
/** \copydoc Preconditioner::post */ /** \copydoc Preconditioner::post */
void post(Domain& x); void post(Domain& x);
@ -299,8 +307,10 @@ namespace Dune
bool additive; bool additive;
bool coarsesolverconverged; bool coarsesolverconverged;
std::shared_ptr<Smoother> coarseSmoother_; std::shared_ptr<Smoother> coarseSmoother_;
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
/** @brief The solver category. */ /** @brief The solver category. */
SolverCategory::Category category_; SolverCategory::Category category_;
#endif
/** @brief The verbosity level. */ /** @brief The verbosity level. */
std::size_t verbosity_; std::size_t verbosity_;
}; };
@ -315,7 +325,9 @@ namespace Dune
buildHierarchy_(amg.buildHierarchy_), buildHierarchy_(amg.buildHierarchy_),
additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged), additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
coarseSmoother_(amg.coarseSmoother_), coarseSmoother_(amg.coarseSmoother_),
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
category_(amg.category_), category_(amg.category_),
#endif
verbosity_(amg.verbosity_) verbosity_(amg.verbosity_)
{ {
if(amg.rhs_) if(amg.rhs_)
@ -337,8 +349,10 @@ namespace Dune
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false), postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
additive(parms.getAdditive()), coarsesolverconverged(true), additive(parms.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(), coarseSmoother_(),
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
// #warning should category be retrieved from matrices? // #warning should category be retrieved from matrices?
category_(SolverCategory::category(*smoothers_->coarsest())), category_(SolverCategory::category(*smoothers_->coarsest())),
#endif
verbosity_(parms.debugLevel()) verbosity_(parms.debugLevel())
{ {
assert(matrices_->isBuilt()); assert(matrices_->isBuilt());
@ -360,14 +374,15 @@ namespace Dune
postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true), postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
additive(criterion.getAdditive()), coarsesolverconverged(true), additive(criterion.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(), coarseSmoother_(),
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
category_(SolverCategory::category(pinfo)), category_(SolverCategory::category(pinfo)),
#endif
verbosity_(criterion.debugLevel()) verbosity_(criterion.debugLevel())
{ {
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
if(SolverCategory::category(matrix) != SolverCategory::category(pinfo)) if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!"); DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
// TODO: reestablish compile time checks. #endif
//static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
// "Matrix and Solver must match in terms of category!");
createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo); createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
} }
@ -458,7 +473,16 @@ namespace Dune
} }
coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs)); coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category()); scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
#else
typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
ScalarProductChooser;
// the scalar product.
scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
#endif
typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector; typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
@ -498,15 +522,23 @@ namespace Dune
// we have to allocate these types using the rebound allocator // we have to allocate these types using the rebound allocator
// in order to ensure that we fulfill the alignement requirements // in order to ensure that we fulfill the alignement requirements
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()), solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
*scalarProduct_, // Cast needed for Dune <=2.5
reinterpret_cast<typename
std::conditional<std::is_same<PI,SequentialInformation>::value,
Dune::SeqScalarProduct<X>,
Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
*coarseSmoother_, 1E-2, 1000, 0)); *coarseSmoother_, 1E-2, 1000, 0));
else else
solver_.reset(); solver_.reset();
}else }else
{ {
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()), solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
*scalarProduct_, // Cast needed for Dune <=2.5
*coarseSmoother_, 1E-2, 1000, 0)); reinterpret_cast<typename
std::conditional<std::is_same<PI,SequentialInformation>::value,
Dune::SeqScalarProduct<X>,
Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
*coarseSmoother_, 1E-2, 1000, 0));
// // we have to allocate these types using the rebound allocator // // we have to allocate these types using the rebound allocator
// // in order to ensure that we fulfill the alignement requirements // // in order to ensure that we fulfill the alignement requirements
// using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other; // using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other;