5#ifndef DUNE_AMGSMOOTHER_HH
6#define DUNE_AMGSMOOTHER_HH
11#include <dune/istl/schwarz.hh>
12#include <dune/istl/novlpschwarz.hh>
13#include <dune/common/propertymap.hh>
71 template<
class X,
class Y>
78 template<
class X,
class Y,
class C,
class T>
80 :
public SmootherTraits<T>
83 template<
class C,
class T>
84 struct SmootherTraits<NonoverlappingBlockPreconditioner<C,T> >
85 :
public SmootherTraits<T>
94 typedef typename T::matrix_type Matrix;
104 void setMatrix(
const Matrix& matrix)
108 virtual void setMatrix(
const Matrix& matrix, [[maybe_unused]]
const AggregatesMap& amap)
114 const Matrix& getMatrix()
const
125 void setComm([[maybe_unused]] T1& comm)
128 const SequentialInformation& getComm()
139 const Matrix* matrix_;
142 SequentialInformation comm_;
146 struct ConstructionArgs
150 template<
class T,
class C=SequentialInformation>
151 class DefaultParallelConstructionArgs
152 :
public ConstructionArgs<T>
155 virtual ~DefaultParallelConstructionArgs()
158 void setComm(
const C& comm)
163 const C& getComm()
const
172 template<
class X,
class Y>
173 class DefaultConstructionArgs<Richardson<X,Y>>
175 typedef Richardson<X,Y> T;
177 typedef typename SmootherTraits<T>::Arguments SmootherArgs;
180 virtual ~DefaultConstructionArgs()
183 template <
class... Args>
184 void setMatrix(
const Args&...)
187 void setArgs(
const SmootherArgs& args)
193 void setComm([[maybe_unused]] T1& comm)
196 const SequentialInformation& getComm()
201 const SmootherArgs getArgs()
const
207 const SmootherArgs* args_;
208 SequentialInformation comm_;
214 struct ConstructionTraits;
219 template<
class M,
class X,
class Y,
int l>
226 return std::make_shared<SeqSSOR<M,X,Y,l>>
235 template<
class M,
class X,
class Y,
int l>
242 return std::make_shared<SeqSOR<M,X,Y,l>>
251 template<
class M,
class X,
class Y,
int l>
258 return std::make_shared<SeqJac<M,X,Y,l>>
266 template<
class X,
class Y>
271 static inline std::shared_ptr<Richardson<X,Y>>
construct(Arguments& args)
273 return std::make_shared<Richardson<X,Y>>
274 (args.getArgs().relaxationFactor);
279 template<
class M,
class X,
class Y>
280 class ConstructionArgs<
SeqILU<M,X,Y> >
284 ConstructionArgs(
int n=0)
306 template<
class M,
class X,
class Y>
309 typedef ConstructionArgs<SeqILU<M,X,Y> >
Arguments;
311 static inline std::shared_ptr<SeqILU<M,X,Y>>
construct(Arguments& args)
313 return std::make_shared<SeqILU<M,X,Y>>
314 (args.getMatrix(), args.getN(), args.getArgs().relaxationFactor);
321 template<
class M,
class X,
class Y,
class C>
324 typedef DefaultParallelConstructionArgs<M,C>
Arguments;
326 static inline std::shared_ptr<ParSSOR<M,X,Y,C>>
construct(Arguments& args)
328 return std::make_shared<ParSSOR<M,X,Y,C>>
329 (args.getMatrix(), args.getArgs().iterations,
330 args.getArgs().relaxationFactor, args.getComm());
334 template<
class X,
class Y,
class C,
class T>
337 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
339 static inline std::shared_ptr<BlockPreconditioner<X,Y,C,T>>
construct(
Arguments& args)
341 auto seqPrec = SeqConstructionTraits::construct(args);
342 return std::make_shared<BlockPreconditioner<X,Y,C,T>> (seqPrec, args.getComm());
346 template<
class C,
class T>
347 struct ConstructionTraits<NonoverlappingBlockPreconditioner<C,T> >
349 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
350 typedef ConstructionTraits<T> SeqConstructionTraits;
351 static inline std::shared_ptr<NonoverlappingBlockPreconditioner<C,T>>
construct(
Arguments& args)
353 auto seqPrec = SeqConstructionTraits::construct(args);
354 return std::make_shared<NonoverlappingBlockPreconditioner<C,T>> (seqPrec, args.getComm());
372 typedef typename Smoother::range_type Range;
373 typedef typename Smoother::domain_type Domain;
382 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
394 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
405 template<
typename LevelContext>
406 void presmooth(LevelContext& levelContext,
size_t steps)
408 for(std::size_t i=0; i < steps; ++i) {
411 ::preSmooth(*levelContext.smoother, *levelContext.lhs,
414 *levelContext.update += *levelContext.lhs;
417 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs, *levelContext.rhs);
418 levelContext.pinfo->project(*levelContext.rhs);
427 template<
typename LevelContext>
430 for(std::size_t i=0; i < steps; ++i) {
432 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs,
435 levelContext.pinfo->project(*levelContext.rhs);
437 ::postSmooth(*levelContext.smoother, *levelContext.lhs, *levelContext.rhs);
439 *levelContext.update += *levelContext.lhs;
443 template<
class M,
class X,
class Y,
int l>
444 struct SmootherApplier<
SeqSOR<M,X,Y,l> >
447 typedef typename Smoother::range_type Range;
448 typedef typename Smoother::domain_type Domain;
450 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
452 smoother.template apply<true>(v,d);
456 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
458 smoother.template apply<false>(v,d);
462 template<
class M,
class X,
class Y,
class C,
int l>
463 struct SmootherApplier<BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > >
465 typedef BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > Smoother;
466 typedef typename Smoother::range_type Range;
467 typedef typename Smoother::domain_type Domain;
469 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
471 smoother.template apply<true>(v,d);
475 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
477 smoother.template apply<false>(v,d);
481 template<
class M,
class X,
class Y,
class C,
int l>
482 struct SmootherApplier<NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > >
484 typedef NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > Smoother;
485 typedef typename Smoother::range_type Range;
486 typedef typename Smoother::domain_type Domain;
488 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
490 smoother.template apply<true>(v,d);
494 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
496 smoother.template apply<false>(v,d);
503 template<
class M,
class X,
class MO,
class MS,
class A>
504 class SeqOverlappingSchwarz;
506 struct MultiplicativeSchwarzMode;
510 template<
class M,
class X,
class MS,
class TA>
511 struct SmootherApplier<SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,
514 typedef SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,MS,TA> Smoother;
515 typedef typename Smoother::range_type Range;
516 typedef typename Smoother::domain_type Domain;
518 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
520 smoother.template apply<true>(v,d);
524 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
526 smoother.template apply<false>(v,d);
535 struct SeqOverlappingSchwarzSmootherArgs
536 :
public DefaultSmootherArgs<T>
543 SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=
vertex,
544 bool onthefly_=
false)
545 : overlap(overlap_), onthefly(onthefly_)
549 template<
class M,
class X,
class TM,
class TS,
class TA>
550 struct SmootherTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
552 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> Arguments;
555 template<
class M,
class X,
class TM,
class TS,
class TA>
556 class ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
557 :
public DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
559 typedef DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Father;
566 typedef typename Vector::value_type Subdomain;
568 virtual void setMatrix(
const M& matrix,
const AggregatesMap& amap)
570 Father::setMatrix(matrix);
572 std::vector<bool> visited(amap.noVertices(),
false);
573 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
574 VisitedMapType visitedMap(visited.begin());
576 MatrixGraph<const M> graph(matrix);
578 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
580 switch(Father::getArgs().overlap) {
583 VertexAdder visitor(subdomains, amap);
584 createSubdomains(matrix, graph, amap, visitor, visitedMap);
587 case SmootherArgs::pairwise :
589 createPairDomains(graph);
592 case SmootherArgs::aggregate :
594 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
595 createSubdomains(matrix, graph, amap, visitor, visitedMap);
600 createSubdomains(matrix, graph, amap, visitor, visitedMap);
603 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
606 void setMatrix(
const M& matrix)
608 Father::setMatrix(matrix);
611 AggregatesMap amap(matrix.N());
612 VertexDescriptor v=0;
613 for(
typename AggregatesMap::iterator iter=amap.begin();
614 iter!=amap.end(); ++iter)
617 std::vector<bool> visited(amap.noVertices(),
false);
618 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
619 VisitedMapType visitedMap(visited.begin());
621 MatrixGraph<const M> graph(matrix);
623 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
625 switch(Father::getArgs().overlap) {
628 VertexAdder visitor(subdomains, amap);
629 createSubdomains(matrix, graph, amap, visitor, visitedMap);
632 case SmootherArgs::aggregate :
634 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
641 case SmootherArgs::pairwise :
643 createPairDomains(graph);
648 createSubdomains(matrix, graph, amap, visitor, visitedMap);
653 const Vector& getSubDomains()
661 VertexAdder(Vector& subdomains_,
const AggregatesMap& aggregates_)
662 : subdomains(subdomains_),
max(-1), subdomain(-1), aggregates(aggregates_)
665 void operator()(
const T& edge)
668 subdomains[subdomain].insert(edge.target());
670 int setAggregate(
const AggregateDescriptor& aggregate_)
672 subdomain=aggregate_;
676 int noSubdomains()
const
682 AggregateDescriptor
max;
683 AggregateDescriptor subdomain;
684 const AggregatesMap& aggregates;
689 void operator()(
const T& edge)
691 int setAggregate(
const AggregateDescriptor& aggregate_)
695 int noSubdomains()
const
702 struct AggregateAdder
704 AggregateAdder(Vector& subdomains_,
const AggregatesMap& aggregates_,
705 const MatrixGraph<const M>& graph_, VM& visitedMap_)
706 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
707 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
710 void operator()(
const T& edge)
712 subdomains[subdomain].insert(edge.target());
717 assert(aggregates[edge.target()]!=aggregate);
719 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
720 graph, vlist, adder, adder,
725 int setAggregate(
const AggregateDescriptor& aggregate_)
727 adder.setAggregate(aggregate_);
728 aggregate=aggregate_;
731 int noSubdomains()
const
737 AggregateDescriptor aggregate;
740 const AggregatesMap& aggregates;
742 const MatrixGraph<const M>& graph;
746 void createPairDomains(
const MatrixGraph<const M>& graph)
750 typedef typename M::size_type size_type;
752 std::set<std::pair<size_type,size_type> > pairs;
754 for(VIter v=graph.begin(), ve=graph.end(); ve != v; ++v)
755 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
758 if(e.source()<e.target())
759 pairs.insert(std::make_pair(e.source(),e.target()));
761 pairs.insert(std::make_pair(e.target(),e.source()));
765 subdomains.resize(pairs.size());
766 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
767 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
768 typename Vector::iterator subdomain=subdomains.begin();
770 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
772 subdomain->insert(s->first);
773 subdomain->insert(s->second);
776 std::size_t minsize=10000;
777 std::size_t maxsize=0;
779 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
780 sum+=subdomains[i].size();
784 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
785 <<
" no="<<subdomains.size()<<std::endl;
788 template<
class Visitor>
789 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
790 const AggregatesMap& amap, Visitor& overlapVisitor,
791 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
798 AggregateDescriptor maxAggregate=0;
800 for(std::size_t i=0; i < amap.noVertices(); ++i)
804 maxAggregate =
std::max(maxAggregate, amap[i]);
806 subdomains.resize(maxAggregate+1+isolated);
809 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
810 subdomains[i].clear();
815 VertexAdder aggregateVisitor(subdomains, amap);
817 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
818 if(!
get(visitedMap, i)) {
819 AggregateDescriptor aggregate=amap[i];
823 subdomains.push_back(Subdomain());
824 aggregate=subdomains.size()-1;
826 overlapVisitor.setAggregate(aggregate);
827 aggregateVisitor.setAggregate(aggregate);
828 subdomains[aggregate].insert(i);
830 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
831 overlapVisitor, visitedMap);
834 std::size_t minsize=10000;
835 std::size_t maxsize=0;
837 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
838 sum+=subdomains[i].size();
842 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
843 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
852 template<
class M,
class X,
class TM,
class TS,
class TA>
853 struct ConstructionTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
855 typedef ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
Arguments;
857 static inline std::shared_ptr<SeqOverlappingSchwarz<M,X,TM,TS,TA>>
construct(
Arguments& args)
859 return std::make_shared<SeqOverlappingSchwarz<M,X,TM,TS,TA>>
861 args.getSubDomains(),
862 args.getArgs().relaxationFactor,
863 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:560
Construction Arguments for the default smoothers.
Definition: smoother.hh:93
VertexIteratorT< const MatrixGraph< Matrix > > ConstVertexIterator
The constant vertex iterator type.
Definition: graph.hh:308
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:73
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:298
Block parallel preconditioner.
Definition: schwarz.hh:278
Richardson preconditioner.
Definition: preconditioners.hh:878
Sequential ILU preconditioner.
Definition: preconditioners.hh:697
The sequential jacobian preconditioner.
Definition: preconditioners.hh:413
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:797
Sequential SOR preconditioner.
Definition: preconditioners.hh:262
Sequential SSOR preconditioner.
Definition: preconditioners.hh:142
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:218
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:471
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:406
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:580
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:571
DefaultSmootherArgs()
Default constructor.
Definition: smoother.hh:56
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
static void postSmooth(Smoother &smoother, Domain &v, const Range &d)
apply post smoothing in forward direction
Definition: smoother.hh:394
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
FieldTraits< T >::real_type RelaxationFactor
The type of the relaxation factor.
Definition: smoother.hh:42
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:592
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:428
RelaxationFactor relaxationFactor
The relaxation factor to use.
Definition: smoother.hh:51
static void preSmooth(Smoother &smoother, Domain &v, const Range &d)
apply pre smoothing in forward direction
Definition: smoother.hh:382
int iterations
The number of iterations to perform.
Definition: smoother.hh:47
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:141
PartitionSet<... > Overlap
Type of PartitionSet for the overlap partition.
Definition: partitionset.hh:249
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Define general preconditioner interface.
Traits class for generically constructing non default constructable types.
Definition: construction.hh:39
The default class for the smoother arguments.
Definition: smoother.hh:38
Helper class for applying the smoothers.
Definition: smoother.hh:370
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:66