3#ifndef DUNE_AMGSMOOTHER_HH
4#define DUNE_AMGSMOOTHER_HH
9#include <dune/istl/schwarz.hh>
10#include <dune/istl/novlpschwarz.hh>
11#include <dune/common/propertymap.hh>
69 template<
class X,
class Y>
76 template<
class X,
class Y,
class C,
class T>
78 :
public SmootherTraits<T>
81 template<
class C,
class T>
82 struct SmootherTraits<NonoverlappingBlockPreconditioner<C,T> >
83 :
public SmootherTraits<T>
92 typedef typename T::matrix_type Matrix;
102 void setMatrix(
const Matrix& matrix)
106 virtual void setMatrix(
const Matrix& matrix, [[maybe_unused]]
const AggregatesMap& amap)
112 const Matrix& getMatrix()
const
123 void setComm([[maybe_unused]] T1& comm)
126 const SequentialInformation& getComm()
137 const Matrix* matrix_;
140 SequentialInformation comm_;
144 struct ConstructionArgs
148 template<
class T,
class C=SequentialInformation>
149 class DefaultParallelConstructionArgs
150 :
public ConstructionArgs<T>
153 virtual ~DefaultParallelConstructionArgs()
156 void setComm(
const C& comm)
161 const C& getComm()
const
170 template<
class X,
class Y>
171 class DefaultConstructionArgs<Richardson<X,Y>>
173 typedef Richardson<X,Y> T;
175 typedef typename SmootherTraits<T>::Arguments SmootherArgs;
178 virtual ~DefaultConstructionArgs()
181 template <
class... Args>
182 void setMatrix(
const Args&...)
185 void setArgs(
const SmootherArgs& args)
191 void setComm([[maybe_unused]] T1& comm)
194 const SequentialInformation& getComm()
199 const SmootherArgs getArgs()
const
205 const SmootherArgs* args_;
206 SequentialInformation comm_;
212 struct ConstructionTraits;
217 template<
class M,
class X,
class Y,
int l>
224 return std::make_shared<SeqSSOR<M,X,Y,l>>
233 template<
class M,
class X,
class Y,
int l>
240 return std::make_shared<SeqSOR<M,X,Y,l>>
249 template<
class M,
class X,
class Y,
int l>
256 return std::make_shared<SeqJac<M,X,Y,l>>
264 template<
class X,
class Y>
269 static inline std::shared_ptr<Richardson<X,Y>>
construct(Arguments& args)
271 return std::make_shared<Richardson<X,Y>>
272 (args.getArgs().relaxationFactor);
277 template<
class M,
class X,
class Y>
278 class ConstructionArgs<
SeqILU<M,X,Y> >
282 ConstructionArgs(
int n=0)
304 template<
class M,
class X,
class Y>
307 typedef ConstructionArgs<SeqILU<M,X,Y> >
Arguments;
309 static inline std::shared_ptr<SeqILU<M,X,Y>>
construct(Arguments& args)
311 return std::make_shared<SeqILU<M,X,Y>>
312 (args.getMatrix(), args.getN(), args.getArgs().relaxationFactor);
319 template<
class M,
class X,
class Y,
class C>
322 typedef DefaultParallelConstructionArgs<M,C>
Arguments;
324 static inline std::shared_ptr<ParSSOR<M,X,Y,C>>
construct(Arguments& args)
326 return std::make_shared<ParSSOR<M,X,Y,C>>
327 (args.getMatrix(), args.getArgs().iterations,
328 args.getArgs().relaxationFactor, args.getComm());
332 template<
class X,
class Y,
class C,
class T>
335 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
337 static inline std::shared_ptr<BlockPreconditioner<X,Y,C,T>>
construct(
Arguments& args)
339 auto seqPrec = SeqConstructionTraits::construct(args);
340 return std::make_shared<BlockPreconditioner<X,Y,C,T>> (seqPrec, args.getComm());
344 template<
class C,
class T>
345 struct ConstructionTraits<NonoverlappingBlockPreconditioner<C,T> >
347 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
348 typedef ConstructionTraits<T> SeqConstructionTraits;
349 static inline std::shared_ptr<NonoverlappingBlockPreconditioner<C,T>>
construct(
Arguments& args)
351 auto seqPrec = SeqConstructionTraits::construct(args);
352 return std::make_shared<NonoverlappingBlockPreconditioner<C,T>> (seqPrec, args.getComm());
370 typedef typename Smoother::range_type Range;
371 typedef typename Smoother::domain_type Domain;
380 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
392 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
403 template<
typename LevelContext>
404 void presmooth(LevelContext& levelContext,
size_t steps)
406 for(std::size_t i=0; i < steps; ++i) {
409 ::preSmooth(*levelContext.smoother, *levelContext.lhs,
412 *levelContext.update += *levelContext.lhs;
415 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs, *levelContext.rhs);
416 levelContext.pinfo->project(*levelContext.rhs);
425 template<
typename LevelContext>
428 for(std::size_t i=0; i < steps; ++i) {
430 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs,
433 levelContext.pinfo->project(*levelContext.rhs);
435 ::postSmooth(*levelContext.smoother, *levelContext.lhs, *levelContext.rhs);
437 *levelContext.update += *levelContext.lhs;
441 template<
class M,
class X,
class Y,
int l>
442 struct SmootherApplier<
SeqSOR<M,X,Y,l> >
445 typedef typename Smoother::range_type Range;
446 typedef typename Smoother::domain_type Domain;
448 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
450 smoother.template apply<true>(v,d);
454 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
456 smoother.template apply<false>(v,d);
460 template<
class M,
class X,
class Y,
class C,
int l>
461 struct SmootherApplier<BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > >
463 typedef BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > Smoother;
464 typedef typename Smoother::range_type Range;
465 typedef typename Smoother::domain_type Domain;
467 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
469 smoother.template apply<true>(v,d);
473 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
475 smoother.template apply<false>(v,d);
479 template<
class M,
class X,
class Y,
class C,
int l>
480 struct SmootherApplier<NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > >
482 typedef NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > Smoother;
483 typedef typename Smoother::range_type Range;
484 typedef typename Smoother::domain_type Domain;
486 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
488 smoother.template apply<true>(v,d);
492 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
494 smoother.template apply<false>(v,d);
501 template<
class M,
class X,
class MO,
class MS,
class A>
502 class SeqOverlappingSchwarz;
504 struct MultiplicativeSchwarzMode;
508 template<
class M,
class X,
class MS,
class TA>
509 struct SmootherApplier<SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,
512 typedef SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,MS,TA> Smoother;
513 typedef typename Smoother::range_type Range;
514 typedef typename Smoother::domain_type Domain;
516 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
518 smoother.template apply<true>(v,d);
522 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
524 smoother.template apply<false>(v,d);
533 struct SeqOverlappingSchwarzSmootherArgs
534 :
public DefaultSmootherArgs<T>
541 SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=
vertex,
542 bool onthefly_=
false)
543 : overlap(overlap_), onthefly(onthefly_)
547 template<
class M,
class X,
class TM,
class TS,
class TA>
548 struct SmootherTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
550 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> Arguments;
553 template<
class M,
class X,
class TM,
class TS,
class TA>
554 class ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
555 :
public DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
557 typedef DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Father;
564 typedef typename Vector::value_type Subdomain;
566 virtual void setMatrix(
const M& matrix,
const AggregatesMap& amap)
568 Father::setMatrix(matrix);
570 std::vector<bool> visited(amap.noVertices(),
false);
571 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
572 VisitedMapType visitedMap(visited.begin());
574 MatrixGraph<const M> graph(matrix);
576 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
578 switch(Father::getArgs().overlap) {
581 VertexAdder visitor(subdomains, amap);
582 createSubdomains(matrix, graph, amap, visitor, visitedMap);
585 case SmootherArgs::pairwise :
587 createPairDomains(graph);
590 case SmootherArgs::aggregate :
592 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
593 createSubdomains(matrix, graph, amap, visitor, visitedMap);
598 createSubdomains(matrix, graph, amap, visitor, visitedMap);
601 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
604 void setMatrix(
const M& matrix)
606 Father::setMatrix(matrix);
609 AggregatesMap amap(matrix.N());
610 VertexDescriptor v=0;
611 for(
typename AggregatesMap::iterator iter=amap.begin();
612 iter!=amap.end(); ++iter)
615 std::vector<bool> visited(amap.noVertices(),
false);
616 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
617 VisitedMapType visitedMap(visited.begin());
619 MatrixGraph<const M> graph(matrix);
621 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
623 switch(Father::getArgs().overlap) {
626 VertexAdder visitor(subdomains, amap);
627 createSubdomains(matrix, graph, amap, visitor, visitedMap);
630 case SmootherArgs::aggregate :
632 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
639 case SmootherArgs::pairwise :
641 createPairDomains(graph);
646 createSubdomains(matrix, graph, amap, visitor, visitedMap);
651 const Vector& getSubDomains()
659 VertexAdder(Vector& subdomains_,
const AggregatesMap& aggregates_)
660 : subdomains(subdomains_),
max(-1), subdomain(-1), aggregates(aggregates_)
663 void operator()(
const T& edge)
666 subdomains[subdomain].insert(edge.target());
668 int setAggregate(
const AggregateDescriptor& aggregate_)
670 subdomain=aggregate_;
674 int noSubdomains()
const
680 AggregateDescriptor
max;
681 AggregateDescriptor subdomain;
682 const AggregatesMap& aggregates;
687 void operator()(
const T& edge)
689 int setAggregate(
const AggregateDescriptor& aggregate_)
693 int noSubdomains()
const
700 struct AggregateAdder
702 AggregateAdder(Vector& subdomains_,
const AggregatesMap& aggregates_,
703 const MatrixGraph<const M>& graph_, VM& visitedMap_)
704 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
705 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
708 void operator()(
const T& edge)
710 subdomains[subdomain].insert(edge.target());
715 assert(aggregates[edge.target()]!=aggregate);
717 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
718 graph, vlist, adder, adder,
723 int setAggregate(
const AggregateDescriptor& aggregate_)
725 adder.setAggregate(aggregate_);
726 aggregate=aggregate_;
729 int noSubdomains()
const
735 AggregateDescriptor aggregate;
738 const AggregatesMap& aggregates;
740 const MatrixGraph<const M>& graph;
744 void createPairDomains(
const MatrixGraph<const M>& graph)
748 typedef typename M::size_type size_type;
750 std::set<std::pair<size_type,size_type> > pairs;
752 for(VIter v=graph.begin(), ve=graph.end(); ve != v; ++v)
753 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
756 if(e.source()<e.target())
757 pairs.insert(std::make_pair(e.source(),e.target()));
759 pairs.insert(std::make_pair(e.target(),e.source()));
763 subdomains.resize(pairs.size());
764 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
765 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
766 typename Vector::iterator subdomain=subdomains.begin();
768 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
770 subdomain->insert(s->first);
771 subdomain->insert(s->second);
774 std::size_t minsize=10000;
775 std::size_t maxsize=0;
777 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
778 sum+=subdomains[i].size();
779 minsize=
std::min(minsize, subdomains[i].size());
780 maxsize=
std::max(maxsize, subdomains[i].size());
782 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
783 <<
" no="<<subdomains.size()<<std::endl;
786 template<
class Visitor>
787 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
788 const AggregatesMap& amap, Visitor& overlapVisitor,
789 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
796 AggregateDescriptor maxAggregate=0;
798 for(std::size_t i=0; i < amap.noVertices(); ++i)
802 maxAggregate =
std::max(maxAggregate, amap[i]);
804 subdomains.resize(maxAggregate+1+isolated);
807 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
808 subdomains[i].clear();
813 VertexAdder aggregateVisitor(subdomains, amap);
815 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
816 if(!get(visitedMap, i)) {
817 AggregateDescriptor aggregate=amap[i];
821 subdomains.push_back(Subdomain());
822 aggregate=subdomains.size()-1;
824 overlapVisitor.setAggregate(aggregate);
825 aggregateVisitor.setAggregate(aggregate);
826 subdomains[aggregate].insert(i);
828 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
829 overlapVisitor, visitedMap);
832 std::size_t minsize=10000;
833 std::size_t maxsize=0;
835 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
836 sum+=subdomains[i].size();
837 minsize=
std::min(minsize, subdomains[i].size());
838 maxsize=
std::max(maxsize, subdomains[i].size());
840 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
841 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
850 template<
class M,
class X,
class TM,
class TS,
class TA>
851 struct ConstructionTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
853 typedef ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
Arguments;
855 static inline std::shared_ptr<SeqOverlappingSchwarz<M,X,TM,TS,TA>>
construct(
Arguments& args)
857 return std::make_shared<SeqOverlappingSchwarz<M,X,TM,TS,TA>>
859 args.getSubDomains(),
860 args.getArgs().relaxationFactor,
861 args.getArgs().onthefly);
Provides classes for the Coloring process of AMG.
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
Construction Arguments for the default smoothers.
Definition: smoother.hh:91
VertexIteratorT< const MatrixGraph< Matrix > > ConstVertexIterator
The constant vertex iterator type.
Definition: graph.hh:306
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
Block parallel preconditioner.
Definition: schwarz.hh:279
Richardson preconditioner.
Definition: preconditioners.hh:711
Sequential ILU preconditioner.
Definition: preconditioners.hh:530
The sequential jacobian preconditioner.
Definition: preconditioners.hh:410
std::vector< subdomain_type, typename std::allocator_traits< TA >::template rebind_alloc< subdomain_type > > subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:795
Sequential SOR preconditioner.
Definition: preconditioners.hh:259
Sequential SSOR preconditioner.
Definition: preconditioners.hh:139
Helper classes for the construction of classes without empty constructor.
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:479
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:504
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:578
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:569
DefaultSmootherArgs()
Default constructor.
Definition: smoother.hh:54
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:42
static void postSmooth(Smoother &smoother, Domain &v, const Range &d)
apply post smoothing in forward direction
Definition: smoother.hh:392
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:50
FieldTraits< T >::real_type RelaxationFactor
The type of the relaxation factor.
Definition: smoother.hh:40
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:590
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
RelaxationFactor relaxationFactor
The relaxation factor to use.
Definition: smoother.hh:49
static void preSmooth(Smoother &smoother, Domain &v, const Range &d)
apply pre smoothing in forward direction
Definition: smoother.hh:380
int iterations
The numbe of iterations to perform.
Definition: smoother.hh:45
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
PartitionSet<... > Overlap
Type of PartitionSet for the overlap partition.
Definition: partitionset.hh:248
Dune namespace.
Definition: alignedallocator.hh:11
Define general preconditioner interface.
Traits class for generically constructing non default constructable types.
Definition: construction.hh:37
The default class for the smoother arguments.
Definition: smoother.hh:36
Helper class for applying the smoothers.
Definition: smoother.hh:368
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64