Loading...
Searching...
No Matches
amgcpr.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_AMG_AMG_CPR_HH
#define DUNE_AMG_AMG_CPR_HH
// NOTE: This file is a modified version of dune/istl/paamg/amg.hh from
// dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/common/TimingMacros.hpp>
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/transfer.hh>
#include <dune/istl/paamg/hierarchy.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
#include <dune/istl/solvertype.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/exceptions.hh>
#include <memory>
namespace Dune
{
namespace Amg
{
#if HAVE_MPI
template<class M, class T>
void redistributeMatrixAmg(M&, M&, SequentialInformation&, SequentialInformation&, T&)
{
// noop
}
template<class M, class PI>
typename std::enable_if<!std::is_same<PI,SequentialInformation>::value,void>::type
redistributeMatrixAmg(M& mat, M& matRedist, PI& info, PI& infoRedist,
Dune::RedistributeInformation<PI>& redistInfo)
{
OPM_TIMEBLOCK(redistributeMatrixAmg);
info.buildGlobalLookup(mat.N());
redistributeMatrixEntries(mat, matRedist, info, infoRedist, redistInfo);
info.freeGlobalLookup();
}
#endif
template<class M, class X, class S, class P, class K, class A>
class KAMG;
template<class T>
class KAmgTwoGrid;
template<class M, class X, class S, class PI=SequentialInformation,
class A=std::allocator<X> >
class AMGCPR : public PreconditionerWithUpdate<X,X>
{
template<class M1, class X1, class S1, class P1, class K1, class A1>
friend class KAMG;
friend class KAmgTwoGrid<AMGCPR>;
public:
typedef M Operator;
typedef PI ParallelInformation;
typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
typedef X Domain;
typedef X Range;
typedef InverseOperator<X,X> CoarseSolver;
typedef S Smoother;
typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs, const Parameters& parms);
template<class C>
AMGCPR(const Operator& fineOperator, const C& criterion,
const SmootherArgs& smootherArgs=SmootherArgs(),
const ParallelInformation& pinfo=ParallelInformation());
AMGCPR(const AMGCPR& amg);
~AMGCPR();
void pre(Domain& x, Range& b) override;
void apply(Domain& v, const Range& d) override;
SolverCategory::Category category() const override
{
return category_;
}
bool hasPerfectUpdate() const override
{
return false;
}
void post(Domain& x) override;
template<class A1>
void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
std::size_t levels();
std::size_t maxlevels();
void recalculateHierarchy()
{
OPM_TIMEBLOCK(recalculateHierarch);
auto copyFlags = NegateSet<typename PI::OwnerSet>();
const auto& matrices = matrices_->matrices();
const auto& aggregatesMapHierarchy = matrices_->aggregatesMaps();
const auto& infoHierarchy = matrices_->parallelInformation();
const auto& redistInfoHierarchy = matrices_->redistributeInformation();
BaseGalerkinProduct productBuilder;
auto aggregatesMap = aggregatesMapHierarchy.begin();
auto info = infoHierarchy.finest();
auto redistInfo = redistInfoHierarchy.begin();
auto matrix = matrices.finest();
auto coarsestMatrix = matrices.coarsest();
using Matrix = typename M::matrix_type;
#if HAVE_MPI
if(matrix.isRedistributed()) {
redistributeMatrixAmg(const_cast<Matrix&>(matrix->getmat()),
const_cast<Matrix&>(matrix.getRedistributed().getmat()),
const_cast<PI&>(*info), const_cast<PI&>(info.getRedistributed()),
const_cast<Dune::RedistributeInformation<PI>&>(*redistInfo));
}
#endif
for(; matrix!=coarsestMatrix; ++aggregatesMap) {
const Matrix& fine = (matrix.isRedistributed() ? matrix.getRedistributed() : *matrix).getmat();
++matrix;
++info;
++redistInfo;
productBuilder.calculate(fine, *(*aggregatesMap), const_cast<Matrix&>(matrix->getmat()), *info, copyFlags);
#if HAVE_MPI
if(matrix.isRedistributed()) {
redistributeMatrixAmg(const_cast<Matrix&>(matrix->getmat()),
const_cast<Matrix&>(matrix.getRedistributed().getmat()),
const_cast<PI&>(*info), const_cast<PI&>(info.getRedistributed()),
const_cast<Dune::RedistributeInformation<PI>&>(*redistInfo));
}
#endif
}
}
template<class C>
void updateSolver(C& criterion, const Operator& /* matrix */, const PI& pinfo);
void update() override;
bool usesDirectCoarseLevelSolver() const;
private:
template<class C>
void createHierarchies(C& criterion, Operator& matrix,
const PI& pinfo)
{
//OPM_TIMEBLOCK(createHierarchies);
// create shared_ptr with empty deleter
std::shared_ptr< Operator > op( &matrix, []( Operator* ) {});
std::shared_ptr< PI > pifo( const_cast< PI* > (&pinfo), []( PI * ) {});
createHierarchies( criterion, op, pifo);
}
template<class C>
void createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
std::shared_ptr< PI > pinfo );
void setupCoarseSolver();
struct LevelContext
{
typedef Smoother SmootherType;
typename Hierarchy<Smoother,A>::Iterator smoother;
typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
typename ParallelInformationHierarchy::Iterator pinfo;
typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
typename Hierarchy<Domain,A>::Iterator lhs;
typename Hierarchy<Domain,A>::Iterator update;
typename Hierarchy<Range,A>::Iterator rhs;
std::size_t level;
};
void mgc(LevelContext& levelContext);
void additiveMgc();
void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
bool moveToCoarseLevel(LevelContext& levelContext);
void initIteratorsWithFineLevel(LevelContext& levelContext);
std::shared_ptr<OperatorHierarchy> matrices_;
SmootherArgs smootherArgs_;
std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
std::shared_ptr<CoarseSolver> solver_;
std::shared_ptr< Hierarchy<Range,A> > rhs_;
std::shared_ptr< Hierarchy<Domain,A> > lhs_;
std::shared_ptr< Hierarchy<Domain,A> > update_;
using ScalarProduct = Dune::ScalarProduct<X>;
std::shared_ptr<ScalarProduct> scalarProduct_;
std::size_t gamma_;
std::size_t preSteps_;
std::size_t postSteps_;
bool buildHierarchy_;
bool additive;
bool coarsesolverconverged;
std::shared_ptr<Smoother> coarseSmoother_;
SolverCategory::Category category_;
std::size_t verbosity_;
};
template<class M, class X, class S, class PI, class A>
inline AMGCPR<M,X,S,PI,A>::AMGCPR(const AMGCPR& amg)
: matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
smoothers_(amg.smoothers_), solver_(amg.solver_),
rhs_(), lhs_(), update_(),
scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
buildHierarchy_(amg.buildHierarchy_),
additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
coarseSmoother_(amg.coarseSmoother_),
category_(amg.category_),
verbosity_(amg.verbosity_)
{
if(amg.rhs_)
rhs_.reset( new Hierarchy<Range,A>(*amg.rhs_) );
if(amg.lhs_)
lhs_.reset( new Hierarchy<Domain,A>(*amg.lhs_) );
if(amg.update_)
update_.reset( new Hierarchy<Domain,A>(*amg.update_) );
}
template<class M, class X, class S, class PI, class A>
AMGCPR<M,X,S,PI,A>::AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs,
const Parameters& parms)
: matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
rhs_(), lhs_(), update_(), scalarProduct_(0),
gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
additive(parms.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(),
// #warning should category be retrieved from matrices?
category_(SolverCategory::category(*smoothers_->coarsest())),
verbosity_(parms.debugLevel())
{
assert(matrices_->isBuilt());
// build the necessary smoother hierarchies
matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
}
template<class M, class X, class S, class PI, class A>
template<class C>
AMGCPR<M,X,S,PI,A>::AMGCPR(const Operator& matrix,
const C& criterion,
const SmootherArgs& smootherArgs,
const PI& pinfo)
: smootherArgs_(smootherArgs),
smoothers_(new Hierarchy<Smoother,A>), solver_(),
rhs_(), lhs_(), update_(), scalarProduct_(),
gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
additive(criterion.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(),
category_(SolverCategory::category(pinfo)),
verbosity_(criterion.debugLevel())
{
if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
}
template<class M, class X, class S, class PI, class A>
AMGCPR<M,X,S,PI,A>::~AMGCPR()
{
if(buildHierarchy_) {
if(solver_)
solver_.reset();
if(coarseSmoother_)
coarseSmoother_.reset();
}
}
template<class M, class X, class S, class PI, class A>
template<class C>
void AMGCPR<M,X,S,PI,A>::updateSolver(C& /* criterion */, const Operator& /* matrix */, const PI& /* pinfo */)
{
update();
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::update()
{
OPM_TIMEBLOCK(update);
Timer watch;
smoothers_.reset(new Hierarchy<Smoother,A>);
solver_.reset();
coarseSmoother_.reset();
scalarProduct_.reset();
buildHierarchy_= true;
coarsesolverconverged = true;
smoothers_.reset(new Hierarchy<Smoother,A>);
recalculateHierarchy();
matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
setupCoarseSolver();
if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
std::cout << "Recalculating galerkin and coarse smoothers "<< matrices_->maxlevels() << " levels "
<< watch.elapsed() << " seconds." << std::endl;
}
}
template<class M, class X, class S, class PI, class A>
template<class C>
void AMGCPR<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
std::shared_ptr< PI > pinfo )
{
OPM_TIMEBLOCK(createHierarchies);
Timer watch;
matrices_.reset(new OperatorHierarchy(matrix, pinfo));
matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
// build the necessary smoother hierarchies
matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
setupCoarseSolver();
if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
<<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
{
OPM_TIMEBLOCK(setupCoarseSolver);
// test whether we should solve on the coarse level. That is the case if we
// have that level and if there was a redistribution on this level then our
// communicator has to be valid (size()>0) as the smoother might try to communicate
// in the constructor.
if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
&& ( ! matrices_->redistributeInformation().back().isSetup() ||
matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
{
// We have the carsest level. Create the coarse Solver
SmootherArgs sargs(smootherArgs_);
sargs.iterations = 1;
typename ConstructionTraits<Smoother>::Arguments cargs;
cargs.setArgs(sargs);
if(matrices_->redistributeInformation().back().isSetup()) {
// Solve on the redistributed partitioning
cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
}else{
cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
cargs.setComm(*matrices_->parallelInformation().coarsest());
}
coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
// Use superlu if we are purely sequential or with only one processor on the coarsest level.
if( SolverSelector::isDirectSolver &&
(std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
|| matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
|| (matrices_->parallelInformation().coarsest().isRedistributed()
&& matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
&& matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
)
{ // redistribute and 1 proc
if(matrices_->parallelInformation().coarsest().isRedistributed())
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
{
// We are still participating on this level
solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
}
else
solver_.reset();
}
else
{
solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
}
if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
}
else
{
if(matrices_->parallelInformation().coarsest().isRedistributed())
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
// We are still participating on this level
// we have to allocate these types using the rebound allocator
// in order to ensure that we fulfill the alignement requirements
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
// 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));
else
solver_.reset();
}else
{
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
// 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));
// // we have to allocate these types using the rebound allocator
// // in order to ensure that we fulfill the alignement requirements
// using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other;
// Alloc alloc;
// auto p = alloc.allocate(1);
// alloc.construct(p,
// const_cast<M&>(*matrices_->matrices().coarsest()),
// *scalarProduct_,
// *coarseSmoother_, 1E-2, 1000, 0);
// solver_.reset(p,[](BiCGSTABSolver<X>* p){
// Alloc alloc;
// alloc.destroy(p);
// alloc.deallocate(p,1);
// });
}
}
}
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::pre(Domain& x, Range& b)
{
OPM_TIMEBLOCK(pre);
// Detect Matrix rows where all offdiagonal entries are
// zero and set x such that A_dd*x_d=b_d
// Thus users can be more careless when setting up their linear
// systems.
typedef typename M::matrix_type Matrix;
typedef typename Matrix::ConstRowIterator RowIter;
typedef typename Matrix::ConstColIterator ColIter;
typedef typename Matrix::block_type Block;
Block zero;
zero=typename Matrix::field_type();
const Matrix& mat=matrices_->matrices().finest()->getmat();
for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
bool isDirichlet = true;
bool hasDiagonal = false;
Block diagonal;
for(ColIter col=row->begin(); col!=row->end(); ++col) {
if(row.index()==col.index()) {
diagonal = *col;
hasDiagonal = false;
}else{
if(*col!=zero)
isDirichlet = false;
}
}
if(isDirichlet && hasDiagonal)
diagonal.solve(x[row.index()], b[row.index()]);
}
if(smoothers_->levels()>0)
smoothers_->finest()->pre(x,b);
else
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
typedef std::shared_ptr< Range > RangePtr ;
typedef std::shared_ptr< Domain > DomainPtr;
// Hierarchy takes ownership of pointers
RangePtr copy( new Range(b) );
rhs_.reset( new Hierarchy<Range,A>(copy) );
DomainPtr dcopy( new Domain(x) );
lhs_.reset( new Hierarchy<Domain,A>(dcopy) );
DomainPtr dcopy2( new Domain(x) );
update_.reset( new Hierarchy<Domain,A>(dcopy2) );
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*update_);
// Preprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Range,A>::Iterator RIterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_->coarsest();
Iterator smoother = smoothers_->finest();
RIterator rhs = rhs_->finest();
DIterator lhs = lhs_->finest();
if(smoothers_->levels()>0) {
assert(lhs_->levels()==rhs_->levels());
assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
if(smoother!=coarsest)
for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
smoother->pre(*lhs,*rhs);
smoother->pre(*lhs,*rhs);
}
// The preconditioner might change x and b. So we have to
// copy the changes to the original vectors.
x = *lhs_->finest();
b = *rhs_->finest();
}
template<class M, class X, class S, class PI, class A>
std::size_t AMGCPR<M,X,S,PI,A>::levels()
{
return matrices_->levels();
}
template<class M, class X, class S, class PI, class A>
std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
{
return matrices_->maxlevels();
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::apply(Domain& v, const Range& d)
{
OPM_TIMEBLOCK(apply);
LevelContext levelContext;
if(additive) {
*(rhs_->finest())=d;
additiveMgc();
v=*lhs_->finest();
}else{
// Init all iterators for the current level
initIteratorsWithFineLevel(levelContext);
*levelContext.lhs = v;
*levelContext.rhs = d;
*levelContext.update=0;
levelContext.level=0;
mgc(levelContext);
if(postSteps_==0||matrices_->maxlevels()==1)
levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
v=*levelContext.update;
}
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
{
OPM_TIMEBLOCK(initIteratorsWithFineLevel);
levelContext.smoother = smoothers_->finest();
levelContext.matrix = matrices_->matrices().finest();
levelContext.pinfo = matrices_->parallelInformation().finest();
levelContext.redist =
matrices_->redistributeInformation().begin();
levelContext.aggregates = matrices_->aggregatesMaps().begin();
levelContext.lhs = lhs_->finest();
levelContext.update = update_->finest();
levelContext.rhs = rhs_->finest();
}
template<class M, class X, class S, class PI, class A>
bool AMGCPR<M,X,S,PI,A>
::moveToCoarseLevel(LevelContext& levelContext)
{
OPM_TIMEBLOCK(moveToCoarseLevel);
bool processNextLevel=true;
if(levelContext.redist->isSetup()) {
levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
levelContext.rhs.getRedistributed());
processNextLevel = levelContext.rhs.getRedistributed().size()>0;
if(processNextLevel) {
//restrict defect to coarse level right hand side.
typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
++levelContext.pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
static_cast<const Range&>(fineRhs.getRedistributed()),
*levelContext.pinfo);
}
}else{
//restrict defect to coarse level right hand side.
typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
++levelContext.pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*levelContext.aggregates),
*levelContext.rhs, static_cast<const Range&>(*fineRhs),
*levelContext.pinfo);
}
if(processNextLevel) {
// prepare coarse system
++levelContext.lhs;
++levelContext.update;
++levelContext.matrix;
++levelContext.level;
++levelContext.redist;
if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
// next level is not the globally coarsest one
++levelContext.smoother;
++levelContext.aggregates;
}
// prepare the update on the next level
*levelContext.update=0;
}
return processNextLevel;
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>
::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
{
OPM_TIMEBLOCK(moveToFineLevel);
if(processNextLevel) {
if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
// previous level is not the globally coarsest one
--levelContext.smoother;
--levelContext.aggregates;
}
--levelContext.redist;
--levelContext.level;
//prolongate and add the correction (update is in coarse left hand side)
--levelContext.matrix;
//typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
--levelContext.lhs;
--levelContext.pinfo;
}
if(levelContext.redist->isSetup()) {
// Need to redistribute during prolongateVector
levelContext.lhs.getRedistributed()=0;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
levelContext.lhs.getRedistributed(),
matrices_->getProlongationDampingFactor(),
*levelContext.pinfo, *levelContext.redist);
}else{
*levelContext.lhs=0;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
matrices_->getProlongationDampingFactor(),
*levelContext.pinfo);
}
if(processNextLevel) {
--levelContext.update;
--levelContext.rhs;
}
*levelContext.update += *levelContext.lhs;
}
template<class M, class X, class S, class PI, class A>
bool AMGCPR<M,X,S,PI,A>::usesDirectCoarseLevelSolver() const
{
return IsDirectSolver< CoarseSolver>::value;
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::mgc(LevelContext& levelContext){
//OPM_TIMEBLOCK(mgc);
if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
// Solve directly
InverseOperatorResult res;
res.converged=true; // If we do not compute this flag will not get updated
if(levelContext.redist->isSetup()) {
levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
if(levelContext.rhs.getRedistributed().size()>0) {
// We are still participating in the computation
levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
levelContext.rhs.getRedistributed());
solver_->apply(levelContext.update.getRedistributed(),
levelContext.rhs.getRedistributed(), res);
}
levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
}else{
levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
solver_->apply(*levelContext.update, *levelContext.rhs, res);
}
if (!res.converged)
coarsesolverconverged = false;
}else{
// presmoothing
presmooth(levelContext, preSteps_);
#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
bool processNextLevel = moveToCoarseLevel(levelContext);
if(processNextLevel) {
// next level
for(std::size_t i=0; i<gamma_; i++)
mgc(levelContext);
}
moveToFineLevel(levelContext, processNextLevel);
#else
*lhs=0;
#endif
if(levelContext.matrix == matrices_->matrices().finest()) {
coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
if(!coarsesolverconverged){
//DUNE_THROW(MathError, "Coarse solver did not converge");
}
}
// postsmoothing
postsmooth(levelContext, postSteps_);
}
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::additiveMgc(){
OPM_TIMEBLOCK(additiveMgc);
// restrict residual to all levels
typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
++pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
}
// pinfo is invalid, set to coarsest level
//pinfo = matrices_->parallelInformation().coarsest
// calculate correction for all levels
lhs = lhs_->finest();
typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
// presmoothing
*lhs=0;
smoother->apply(*lhs, *rhs);
}
// Coarse level solve
#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
InverseOperatorResult res;
pinfo->copyOwnerToAll(*rhs, *rhs);
solver_->apply(*lhs, *rhs, res);
if(!res.converged)
DUNE_THROW(MathError, "Coarse solver did not converge");
#else
*lhs=0;
#endif
// Prologate and add up corrections from all levels
--pinfo;
--aggregates;
for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
}
}
template<class M, class X, class S, class PI, class A>
void AMGCPR<M,X,S,PI,A>::post([[maybe_unused]] Domain& x)
{
OPM_TIMEBLOCK(post);
// Postprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_->coarsest();
Iterator smoother = smoothers_->finest();
DIterator lhs = lhs_->finest();
if(smoothers_->levels()>0) {
if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
smoother->post(*lhs);
if(smoother!=coarsest)
for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
smoother->post(*lhs);
smoother->post(*lhs);
}
//delete &(*lhs_->finest());
lhs_.reset();
//delete &(*update_->finest());
update_.reset();
//delete &(*rhs_->finest());
rhs_.reset();
}
template<class M, class X, class S, class PI, class A>
template<class A1>
void AMGCPR<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
{
matrices_->getCoarsestAggregatesOnFinest(cont);
}
} // end namespace Amg
} // end namespace Dune
#endif
Dune::Amg::AMGCPR
Parallel algebraic multigrid based on agglomeration.
Definition amgcpr.hh:88
Dune::Amg::KAMG
Definition amgcpr.hh:70
Dune::Amg::KAmgTwoGrid
Definition amgcpr.hh:73
Dune::PreconditionerWithUpdate
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Dune::Amg::AMGCPR::Operator
M Operator
The matrix operator type.
Definition amgcpr.hh:96
Dune::Amg::AMGCPR::LevelContext::smoother
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition amgcpr.hh:295
Dune::Amg::AMGCPR::post
void post(Domain &x) override
Definition amgcpr.hh:963
Dune::Amg::AMGCPR::AMGCPR
AMGCPR(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition amgcpr.hh:415
Dune::Amg::AMGCPR::LevelContext::update
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition amgcpr.hh:319
Dune::Amg::AMGCPR::LevelContext::rhs
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition amgcpr.hh:323
Dune::Amg::AMGCPR::OperatorHierarchy
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition amgcpr.hh:105
Dune::Amg::AMGCPR::LevelContext::pinfo
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition amgcpr.hh:303
Dune::Amg::AMGCPR::LevelContext::matrix
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition amgcpr.hh:299
Dune::Amg::AMGCPR::getCoarsestAggregateNumbers
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition amgcpr.hh:990
Dune::Amg::AMGCPR::ParallelInformationHierarchy
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition amgcpr.hh:107
Dune::Amg::AMGCPR::recalculateHierarchy
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition amgcpr.hh:199
Dune::Amg::AMGCPR::SmootherArgs
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition amgcpr.hh:123
Dune::Amg::AMGCPR::apply
void apply(Domain &v, const Range &d) override
Definition amgcpr.hh:713
Dune::Amg::AMGCPR::Range
X Range
The range type.
Definition amgcpr.hh:112
Dune::Amg::AMGCPR::pre
void pre(Domain &x, Range &b) override
Definition amgcpr.hh:619
Dune::Amg::AMGCPR::Domain
X Domain
The domain type.
Definition amgcpr.hh:110
Dune::Amg::AMGCPR::update
void update() override
Update the coarse solver and the hierarchies.
Definition amgcpr.hh:475
Dune::Amg::AMGCPR::Smoother
S Smoother
The type of the smoother.
Definition amgcpr.hh:120
Dune::Amg::AMGCPR::updateSolver
void updateSolver(C &criterion, const Operator &, const PI &pinfo)
Update the coarse solver and the hierarchies.
Definition amgcpr.hh:469
Dune::Amg::AMGCPR::LevelContext::aggregates
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition amgcpr.hh:311
Dune::Amg::AMGCPR::ParallelInformation
PI ParallelInformation
The type of the parallel information.
Definition amgcpr.hh:103
Dune::Amg::AMGCPR::CoarseSolver
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition amgcpr.hh:114
Dune::Amg::AMGCPR::LevelContext::level
std::size_t level
The level index.
Definition amgcpr.hh:327
Dune::Amg::AMGCPR::LevelContext::lhs
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition amgcpr.hh:315
Dune::Amg::AMGCPR::usesDirectCoarseLevelSolver
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition amgcpr.hh:852
Dune::Amg::AMGCPR::LevelContext::redist
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition amgcpr.hh:307
Dune::Amg::AMGCPR::category
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category).
Definition amgcpr.hh:167