dune-istl  2.7.1
twolevelmethod.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_TWOLEVELMETHOD_HH
4 #define DUNE_ISTL_TWOLEVELMETHOD_HH
5 
6 #include <tuple>
7 
9 #include"amg.hh"
10 #include"galerkin.hh"
11 #include<dune/istl/solver.hh>
12 
13 #include<dune/common/unused.hh>
14 
22 namespace Dune
23 {
24 namespace Amg
25 {
26 
36 template<class FO, class CO>
38 {
39 public:
44  typedef FO FineOperatorType;
48  typedef typename FineOperatorType::range_type FineRangeType;
52  typedef typename FineOperatorType::domain_type FineDomainType;
57  typedef CO CoarseOperatorType;
61  typedef typename CoarseOperatorType::range_type CoarseRangeType;
65  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
70  std::shared_ptr<CoarseOperatorType>& getCoarseLevelOperator()
71  {
72  return operator_;
73  }
79  {
80  return rhs_;
81  }
82 
88  {
89  return lhs_;
90  }
100  virtual void moveToCoarseLevel(const FineRangeType& fineRhs)=0;
110  virtual void moveToFineLevel(FineDomainType& fineLhs)=0;
118  virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0;
119 
121  virtual LevelTransferPolicy* clone() const =0;
122 
125 
126  protected:
132  std::shared_ptr<CoarseOperatorType> operator_;
133 };
134 
140 template<class O, class C>
142  : public LevelTransferPolicy<O,O>
143 {
145 public:
147  typedef C Criterion;
149 
151  : criterion_(crit)
152  {}
153 
154  void createCoarseLevelSystem(const O& fineOperator)
155  {
156  prolongDamp_ = criterion_.getProlongationDampingFactor();
160  Dune::Amg::EdgeProperties,Dune::IdentityMap,Dune::IdentityMap> PropertiesGraph;
161  MatrixGraph mg(fineOperator.getmat());
162  PropertiesGraph pg(mg,Dune::IdentityMap(),Dune::IdentityMap());
163  typedef NegateSet<typename ParallelInformation::OwnerSet> OverlapFlags;
164 
165  aggregatesMap_ = std::make_shared<AggregatesMap>(pg.maxVertex()+1);
166 
167  int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
168 
169  std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
170  aggregatesMap_->buildAggregates(fineOperator.getmat(), pg, criterion_, true);
171  std::cout<<"no aggregates="<<noAggregates<<" iso="<<isoAggregates<<" one="<<oneAggregates<<" skipped="<<skippedAggregates<<std::endl;
172  // misuse coarsener to renumber aggregates
174  typedef std::vector<bool>::iterator Iterator;
175  typedef Dune::IteratorPropertyMap<Iterator, Dune::IdentityMap> VisitedMap;
176  std::vector<bool> excluded(fineOperator.getmat().N(), false);
177  VisitedMap vm(excluded.begin(), Dune::IdentityMap());
178  ParallelInformation pinfo;
179  std::size_t aggregates = renumberer.coarsen(pinfo, pg, vm,
180  *aggregatesMap_, pinfo,
181  noAggregates);
182  std::vector<bool>& visited=excluded;
183 
184  typedef std::vector<bool>::iterator Iterator;
185 
186  for(Iterator iter= visited.begin(), end=visited.end();
187  iter != end; ++iter)
188  *iter=false;
189  matrix_.reset(productBuilder.build(mg, vm,
191  *aggregatesMap_,
192  aggregates,
193  OverlapFlags()));
194  productBuilder.calculate(fineOperator.getmat(), *aggregatesMap_, *matrix_, pinfo, OverlapFlags());
195  this->lhs_.resize(this->matrix_->M());
196  this->rhs_.resize(this->matrix_->N());
197  this->operator_ = std::make_shared<O>(*matrix_);
198  }
199 
200  void moveToCoarseLevel(const typename FatherType::FineRangeType& fineRhs)
201  {
203  ::restrictVector(*aggregatesMap_, this->rhs_, fineRhs, ParallelInformation());
204  this->lhs_=0;
205  }
206 
208  {
210  ::prolongateVector(*aggregatesMap_, this->lhs_, fineLhs,
211  prolongDamp_, ParallelInformation());
212  }
213 
215  {
216  return new AggregationLevelTransferPolicy(*this);
217  }
218 
219 private:
220  typename O::matrix_type::field_type prolongDamp_;
221  std::shared_ptr<AggregatesMap> aggregatesMap_;
222  Criterion criterion_;
223  std::shared_ptr<typename O::matrix_type> matrix_;
224 };
225 
232 template<class O, class S, class C>
234 {
235 public:
237  typedef O Operator;
239  typedef typename O::range_type X;
241  typedef C Criterion;
243  typedef S Smoother;
254  : smootherArgs_(args), criterion_(c)
255  {}
258  : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
259  criterion_(other.criterion_)
260  {}
261 private:
268  struct AMGInverseOperator : public InverseOperator<X,X>
269  {
270  AMGInverseOperator(const typename AMGType::Operator& op,
271  const Criterion& crit,
272  const typename AMGType::SmootherArgs& args)
273  : amg_(op, crit,args), first_(true)
274  {}
275 
276  void apply(X& x, X& b, double reduction, InverseOperatorResult& res)
277  {
278  DUNE_UNUSED_PARAMETER(reduction);
279  DUNE_UNUSED_PARAMETER(res);
280  if(first_)
281  {
282  amg_.pre(x,b);
283  first_=false;
284  x_=x;
285  }
286  amg_.apply(x,b);
287  }
288 
289  void apply(X& x, X& b, InverseOperatorResult& res)
290  {
291  return apply(x,b,1e-8,res);
292  }
293 
295  virtual SolverCategory::Category category() const
296  {
297  return amg_.category();
298  }
299 
300  ~AMGInverseOperator()
301  {
302  if(!first_)
303  amg_.post(x_);
304  }
305  AMGInverseOperator(const AMGInverseOperator& other)
306  : x_(other.x_), amg_(other.amg_), first_(other.first_)
307  {
308  }
309  private:
310  X x_;
311  AMGType amg_;
312  bool first_;
313  };
314 
315 public:
317  typedef AMGInverseOperator CoarseLevelSolver;
318 
326  template<class P>
328  {
329  coarseOperator_=transferPolicy.getCoarseLevelOperator();
330  AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
331  criterion_,
332  smootherArgs_);
333 
334  return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
335 
336  }
337 
338 private:
340  std::shared_ptr<Operator> coarseOperator_;
342  SmootherArgs smootherArgs_;
344  Criterion criterion_;
345 };
346 
352 template<class FO, class CSP, class S>
354  public Preconditioner<typename FO::domain_type, typename FO::range_type>
355 {
356 public:
360  typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
365  typedef FO FineOperatorType;
369  typedef typename FineOperatorType::range_type FineRangeType;
373  typedef typename FineOperatorType::domain_type FineDomainType;
378  typedef typename CSP::Operator CoarseOperatorType;
382  typedef typename CoarseOperatorType::range_type CoarseRangeType;
386  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
390  typedef S SmootherType;
391 
407  std::shared_ptr<SmootherType> smoother,
409  CoarseOperatorType>& policy,
410  CoarseLevelSolverPolicy& coarsePolicy,
411  std::size_t preSteps=1, std::size_t postSteps=1)
412  : operator_(&op), smoother_(smoother),
413  preSteps_(preSteps), postSteps_(postSteps)
414  {
415  policy_ = policy.clone();
416  policy_->createCoarseLevelSystem(*operator_);
417  coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_);
418  }
419 
421  : operator_(other.operator_), coarseSolver_(new CoarseLevelSolver(*other.coarseSolver_)),
422  smoother_(other.smoother_), policy_(other.policy_->clone()),
423  preSteps_(other.preSteps_), postSteps_(other.postSteps_)
424  {}
425 
427  {
428  // Each instance has its own policy.
429  delete policy_;
430  delete coarseSolver_;
431  }
432 
434  {
435  smoother_->pre(x,b);
436  }
437 
439  {
440  DUNE_UNUSED_PARAMETER(x);
441  }
442 
443  void apply(FineDomainType& v, const FineRangeType& d)
444  {
445  FineDomainType u(v);
446  FineRangeType rhs(d);
447  LevelContext context;
449  context.pinfo=&info;
450  context.lhs=&u;
451  context.update=&v;
452  context.smoother=smoother_;
453  context.rhs=&rhs;
454  context.matrix=operator_;
455  // Presmoothing
456  presmooth(context, preSteps_);
457  //Coarse grid correction
458  policy_->moveToCoarseLevel(*context.rhs);
460  coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res);
461  *context.lhs=0;
462  policy_->moveToFineLevel(*context.lhs);
463  *context.update += *context.lhs;
464  // Postsmoothing
465  postsmooth(context, postSteps_);
466 
467  }
468 
471  {
473  }
474 
475 private:
479  struct LevelContext
480  {
482  typedef S SmootherType;
484  std::shared_ptr<SmootherType> smoother;
486  FineDomainType* lhs;
487  /*
488  * @brief The right hand side holding the current residual.
489  *
490  * This is passed to the smoother as the right hand side.
491  */
492  FineRangeType* rhs;
498  FineDomainType* update;
500  SequentialInformation* pinfo;
506  const FineOperatorType* matrix;
507  };
508  const FineOperatorType* operator_;
510  CoarseLevelSolver* coarseSolver_;
512  std::shared_ptr<S> smoother_;
514  LevelTransferPolicy<FO,typename CSP::Operator>* policy_;
516  std::size_t preSteps_;
518  std::size_t postSteps_;
519 };
520 }// end namespace Amg
521 }// end namespace Dune
522 
524 #endif
Define general, extensible interface for operators. The available implementation wraps a matrix.
The AMG preconditioner.
Provides a class for building the galerkin product based on a aggregation scheme.
Define general, extensible interface for inverse operators.
G::MutableMatrix * build(G &fineGraph, V &visitedMap, const ParallelInformation &pinfo, AggregatesMap< typename G::VertexDescriptor > &aggregates, const typename G::Matrix::size_type &size, const Set &copy)
Calculates the coarse matrix via a Galerkin product.
Definition: galerkin.hh:564
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:97
Operator Operator
The matrix operator type.
Definition: amg.hh:70
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:468
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:490
void calculate(const M &fine, const AggregatesMap< V > &aggregates, M &coarse, const I &pinfo, const O &copy)
Calculate the galerkin product.
Definition: allocator.hh:7
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:556
Class representing the properties of an ede in the matrix graph.
Definition: dependency.hh:38
Class representing a node in the matrix graph.
Definition: dependency.hh:125
Definition: galerkin.hh:117
The (undirected) graph of a matrix.
Definition: graph.hh:49
Attaches properties to the edges and vertices of a graph.
Definition: graph.hh:976
VertexDescriptor maxVertex() const
Get the maximal vertex descriptor.
Definition: indicescoarsener.hh:35
Definition: pinfo.hh:26
The default class for the smoother arguments.
Definition: smoother.hh:37
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Abstract base class for transfer between levels and creation of the coarse level system.
Definition: twolevelmethod.hh:38
CO CoarseOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:57
virtual void moveToCoarseLevel(const FineRangeType &fineRhs)=0
Transfers the data to the coarse level.
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethod.hh:48
virtual void createCoarseLevelSystem(const FineOperatorType &fineOperator)=0
Algebraically creates the coarse level system.
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethod.hh:61
CoarseRangeType & getCoarseLevelRhs()
Get the coarse level right hand side.
Definition: twolevelmethod.hh:78
virtual ~LevelTransferPolicy()
Destructor.
Definition: twolevelmethod.hh:124
CoarseDomainType lhs_
The coarse level lhs.
Definition: twolevelmethod.hh:130
std::shared_ptr< CoarseOperatorType > operator_
the coarse level linear operator.
Definition: twolevelmethod.hh:132
CoarseRangeType rhs_
The coarse level rhs.
Definition: twolevelmethod.hh:128
std::shared_ptr< CoarseOperatorType > & getCoarseLevelOperator()
Get the coarse level operator.
Definition: twolevelmethod.hh:70
virtual void moveToFineLevel(FineDomainType &fineLhs)=0
Updates the fine level linear system after the correction of the coarse levels system.
FO FineOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:44
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethod.hh:65
CoarseDomainType & getCoarseLevelLhs()
Get the coarse level left hand side.
Definition: twolevelmethod.hh:87
virtual LevelTransferPolicy * clone() const =0
Clone the current object.
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethod.hh:52
A LeveTransferPolicy that used aggregation to construct the coarse level system.
Definition: twolevelmethod.hh:143
C Criterion
Definition: twolevelmethod.hh:147
AggregationLevelTransferPolicy * clone() const
Clone the current object.
Definition: twolevelmethod.hh:214
AggregationLevelTransferPolicy(const Criterion &crit)
Definition: twolevelmethod.hh:150
void moveToFineLevel(typename FatherType::FineDomainType &fineLhs)
Updates the fine level linear system after the correction of the coarse levels system.
Definition: twolevelmethod.hh:207
void moveToCoarseLevel(const typename FatherType::FineRangeType &fineRhs)
Definition: twolevelmethod.hh:200
SequentialInformation ParallelInformation
Definition: twolevelmethod.hh:148
LevelTransferPolicy< O, O > FatherType
Definition: twolevelmethod.hh:146
void createCoarseLevelSystem(const O &fineOperator)
Algebraically creates the coarse level system.
Definition: twolevelmethod.hh:154
A policy class for solving the coarse level system using one step of AMG.
Definition: twolevelmethod.hh:234
OneStepAMGCoarseSolverPolicy(const SmootherArgs &args, const Criterion &c)
Constructs the coarse solver policy.
Definition: twolevelmethod.hh:253
AMGInverseOperator CoarseLevelSolver
The type of solver constructed for the coarse level.
Definition: twolevelmethod.hh:317
OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy &other)
Copy constructor.
Definition: twolevelmethod.hh:257
O::range_type X
The type of the range and domain of the operator.
Definition: twolevelmethod.hh:239
C Criterion
The type of the crition used for the aggregation within AMG.
Definition: twolevelmethod.hh:241
Dune::Amg::SmootherTraits< S >::Arguments SmootherArgs
The type of the arguments used for constructing the smoother.
Definition: twolevelmethod.hh:245
O Operator
The type of the linear operator used.
Definition: twolevelmethod.hh:237
CoarseLevelSolver * createCoarseLevelSolver(P &transferPolicy)
Constructs a coarse level solver.
Definition: twolevelmethod.hh:327
AMG< Operator, X, Smoother > AMGType
The type of the AMG construct on the coarse level.
Definition: twolevelmethod.hh:247
S Smoother
The type of the smoother used in AMG.
Definition: twolevelmethod.hh:243
Definition: twolevelmethod.hh:355
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethod.hh:382
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethod.hh:373
TwoLevelMethod(const TwoLevelMethod &other)
Definition: twolevelmethod.hh:420
void pre(FineDomainType &x, FineRangeType &b)
Definition: twolevelmethod.hh:433
FO FineOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:365
CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver
The type of the coarse level solver.
Definition: twolevelmethod.hh:360
void apply(FineDomainType &v, const FineRangeType &d)
Definition: twolevelmethod.hh:443
CSP CoarseLevelSolverPolicy
The type of the policy for constructing the coarse level solver.
Definition: twolevelmethod.hh:358
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethod.hh:386
TwoLevelMethod(const FineOperatorType &op, std::shared_ptr< SmootherType > smoother, const LevelTransferPolicy< FineOperatorType, CoarseOperatorType > &policy, CoarseLevelSolverPolicy &coarsePolicy, std::size_t preSteps=1, std::size_t postSteps=1)
Constructs a two level method.
Definition: twolevelmethod.hh:406
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: twolevelmethod.hh:470
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethod.hh:369
~TwoLevelMethod()
Definition: twolevelmethod.hh:426
CSP::Operator CoarseOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:378
void post(FineDomainType &x)
Definition: twolevelmethod.hh:438
S SmootherType
The type of the fine level smoother.
Definition: twolevelmethod.hh:390
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Statistics about the application of an inverse operator.
Definition: solver.hh:46
Abstract base class for all solvers.
Definition: solver.hh:97
virtual void apply(X &x, X &b, InverseOperatorResult &res)=0
Apply inverse operator,.
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23