dune-istl  2.8.0
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 
20 namespace Dune
21 {
22 namespace Amg
23 {
24 
34 template<class FO, class CO>
36 {
37 public:
42  typedef FO FineOperatorType;
46  typedef typename FineOperatorType::range_type FineRangeType;
50  typedef typename FineOperatorType::domain_type FineDomainType;
55  typedef CO CoarseOperatorType;
59  typedef typename CoarseOperatorType::range_type CoarseRangeType;
63  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
68  std::shared_ptr<CoarseOperatorType>& getCoarseLevelOperator()
69  {
70  return operator_;
71  }
77  {
78  return rhs_;
79  }
80 
86  {
87  return lhs_;
88  }
98  virtual void moveToCoarseLevel(const FineRangeType& fineRhs)=0;
108  virtual void moveToFineLevel(FineDomainType& fineLhs)=0;
116  virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0;
117 
119  virtual LevelTransferPolicy* clone() const =0;
120 
123 
124  protected:
130  std::shared_ptr<CoarseOperatorType> operator_;
131 };
132 
138 template<class O, class C>
140  : public LevelTransferPolicy<O,O>
141 {
143 public:
145  typedef C Criterion;
147 
149  : criterion_(crit)
150  {}
151 
152  void createCoarseLevelSystem(const O& fineOperator)
153  {
154  prolongDamp_ = criterion_.getProlongationDampingFactor();
158  Dune::Amg::EdgeProperties,Dune::IdentityMap,Dune::IdentityMap> PropertiesGraph;
159  MatrixGraph mg(fineOperator.getmat());
160  PropertiesGraph pg(mg,Dune::IdentityMap(),Dune::IdentityMap());
161  typedef NegateSet<typename ParallelInformation::OwnerSet> OverlapFlags;
162 
163  aggregatesMap_ = std::make_shared<AggregatesMap>(pg.maxVertex()+1);
164 
165  int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
166 
167  std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
168  aggregatesMap_->buildAggregates(fineOperator.getmat(), pg, criterion_, true);
169  std::cout<<"no aggregates="<<noAggregates<<" iso="<<isoAggregates<<" one="<<oneAggregates<<" skipped="<<skippedAggregates<<std::endl;
170  // misuse coarsener to renumber aggregates
172  typedef std::vector<bool>::iterator Iterator;
173  typedef Dune::IteratorPropertyMap<Iterator, Dune::IdentityMap> VisitedMap;
174  std::vector<bool> excluded(fineOperator.getmat().N(), false);
175  VisitedMap vm(excluded.begin(), Dune::IdentityMap());
176  ParallelInformation pinfo;
177  std::size_t aggregates = renumberer.coarsen(pinfo, pg, vm,
178  *aggregatesMap_, pinfo,
179  noAggregates);
180  std::vector<bool>& visited=excluded;
181 
182  typedef std::vector<bool>::iterator Iterator;
183 
184  for(Iterator iter= visited.begin(), end=visited.end();
185  iter != end; ++iter)
186  *iter=false;
187  matrix_.reset(productBuilder.build(mg, vm,
189  *aggregatesMap_,
190  aggregates,
191  OverlapFlags()));
192  productBuilder.calculate(fineOperator.getmat(), *aggregatesMap_, *matrix_, pinfo, OverlapFlags());
193  this->lhs_.resize(this->matrix_->M());
194  this->rhs_.resize(this->matrix_->N());
195  this->operator_ = std::make_shared<O>(*matrix_);
196  }
197 
198  void moveToCoarseLevel(const typename FatherType::FineRangeType& fineRhs)
199  {
201  ::restrictVector(*aggregatesMap_, this->rhs_, fineRhs, ParallelInformation());
202  this->lhs_=0;
203  }
204 
206  {
208  ::prolongateVector(*aggregatesMap_, this->lhs_, fineLhs,
209  prolongDamp_, ParallelInformation());
210  }
211 
213  {
214  return new AggregationLevelTransferPolicy(*this);
215  }
216 
217 private:
218  typename O::matrix_type::field_type prolongDamp_;
219  std::shared_ptr<AggregatesMap> aggregatesMap_;
220  Criterion criterion_;
221  std::shared_ptr<typename O::matrix_type> matrix_;
222 };
223 
230 template<class O, class S, class C>
232 {
233 public:
235  typedef O Operator;
237  typedef typename O::range_type X;
239  typedef C Criterion;
241  typedef S Smoother;
252  : smootherArgs_(args), criterion_(c)
253  {}
256  : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
257  criterion_(other.criterion_)
258  {}
259 private:
266  struct AMGInverseOperator : public InverseOperator<X,X>
267  {
268  AMGInverseOperator(const typename AMGType::Operator& op,
269  const Criterion& crit,
270  const typename AMGType::SmootherArgs& args)
271  : amg_(op, crit,args), first_(true)
272  {}
273 
274  void apply(X& x, X& b, [[maybe_unused]] double reduction, [[maybe_unused]] InverseOperatorResult& res)
275  {
276  if(first_)
277  {
278  amg_.pre(x,b);
279  first_=false;
280  x_=x;
281  }
282  amg_.apply(x,b);
283  }
284 
285  void apply(X& x, X& b, InverseOperatorResult& res)
286  {
287  return apply(x,b,1e-8,res);
288  }
289 
291  virtual SolverCategory::Category category() const
292  {
293  return amg_.category();
294  }
295 
296  ~AMGInverseOperator()
297  {
298  if(!first_)
299  amg_.post(x_);
300  }
301  AMGInverseOperator(const AMGInverseOperator& other)
302  : x_(other.x_), amg_(other.amg_), first_(other.first_)
303  {
304  }
305  private:
306  X x_;
307  AMGType amg_;
308  bool first_;
309  };
310 
311 public:
313  typedef AMGInverseOperator CoarseLevelSolver;
314 
322  template<class P>
324  {
325  coarseOperator_=transferPolicy.getCoarseLevelOperator();
326  AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
327  criterion_,
328  smootherArgs_);
329 
330  return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
331 
332  }
333 
334 private:
336  std::shared_ptr<Operator> coarseOperator_;
338  SmootherArgs smootherArgs_;
340  Criterion criterion_;
341 };
342 
348 template<class FO, class CSP, class S>
350  public Preconditioner<typename FO::domain_type, typename FO::range_type>
351 {
352 public:
356  typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
361  typedef FO FineOperatorType;
365  typedef typename FineOperatorType::range_type FineRangeType;
369  typedef typename FineOperatorType::domain_type FineDomainType;
374  typedef typename CSP::Operator CoarseOperatorType;
378  typedef typename CoarseOperatorType::range_type CoarseRangeType;
382  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
386  typedef S SmootherType;
387 
403  std::shared_ptr<SmootherType> smoother,
405  CoarseOperatorType>& policy,
406  CoarseLevelSolverPolicy& coarsePolicy,
407  std::size_t preSteps=1, std::size_t postSteps=1)
408  : operator_(&op), smoother_(smoother),
409  preSteps_(preSteps), postSteps_(postSteps)
410  {
411  policy_ = policy.clone();
412  policy_->createCoarseLevelSystem(*operator_);
413  coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_);
414  }
415 
417  : operator_(other.operator_), coarseSolver_(new CoarseLevelSolver(*other.coarseSolver_)),
418  smoother_(other.smoother_), policy_(other.policy_->clone()),
419  preSteps_(other.preSteps_), postSteps_(other.postSteps_)
420  {}
421 
423  {
424  // Each instance has its own policy.
425  delete policy_;
426  delete coarseSolver_;
427  }
428 
430  {
431  smoother_->pre(x,b);
432  }
433 
434  void post([[maybe_unused]] FineDomainType& x)
435  {}
436 
437  void apply(FineDomainType& v, const FineRangeType& d)
438  {
439  FineDomainType u(v);
440  FineRangeType rhs(d);
441  LevelContext context;
443  context.pinfo=&info;
444  context.lhs=&u;
445  context.update=&v;
446  context.smoother=smoother_;
447  context.rhs=&rhs;
448  context.matrix=operator_;
449  // Presmoothing
450  presmooth(context, preSteps_);
451  //Coarse grid correction
452  policy_->moveToCoarseLevel(*context.rhs);
454  coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res);
455  *context.lhs=0;
456  policy_->moveToFineLevel(*context.lhs);
457  *context.update += *context.lhs;
458  // Postsmoothing
459  postsmooth(context, postSteps_);
460 
461  }
462 
465  {
467  }
468 
469 private:
473  struct LevelContext
474  {
476  typedef S SmootherType;
478  std::shared_ptr<SmootherType> smoother;
480  FineDomainType* lhs;
481  /*
482  * @brief The right hand side holding the current residual.
483  *
484  * This is passed to the smoother as the right hand side.
485  */
486  FineRangeType* rhs;
492  FineDomainType* update;
494  SequentialInformation* pinfo;
500  const FineOperatorType* matrix;
501  };
502  const FineOperatorType* operator_;
504  CoarseLevelSolver* coarseSolver_;
506  std::shared_ptr<S> smoother_;
508  LevelTransferPolicy<FO,typename CSP::Operator>* policy_;
510  std::size_t preSteps_;
512  std::size_t postSteps_;
513 };
514 }// end namespace Amg
515 }// end namespace Dune
516 
518 #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:561
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:98
Operator Operator
The matrix operator type.
Definition: amg.hh:71
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
void calculate(const M &fine, const AggregatesMap< V > &aggregates, M &coarse, const I &pinfo, const O &copy)
Calculate the galerkin product.
Definition: allocator.hh:9
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
Class representing the properties of an ede in the matrix graph.
Definition: dependency.hh:37
Class representing a node in the matrix graph.
Definition: dependency.hh:124
Definition: galerkin.hh:116
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:34
Definition: pinfo.hh:26
The default class for the smoother arguments.
Definition: smoother.hh:36
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:36
CO CoarseOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:55
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:46
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:59
CoarseRangeType & getCoarseLevelRhs()
Get the coarse level right hand side.
Definition: twolevelmethod.hh:76
virtual ~LevelTransferPolicy()
Destructor.
Definition: twolevelmethod.hh:122
CoarseDomainType lhs_
The coarse level lhs.
Definition: twolevelmethod.hh:128
std::shared_ptr< CoarseOperatorType > operator_
the coarse level linear operator.
Definition: twolevelmethod.hh:130
CoarseRangeType rhs_
The coarse level rhs.
Definition: twolevelmethod.hh:126
std::shared_ptr< CoarseOperatorType > & getCoarseLevelOperator()
Get the coarse level operator.
Definition: twolevelmethod.hh:68
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:42
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethod.hh:63
CoarseDomainType & getCoarseLevelLhs()
Get the coarse level left hand side.
Definition: twolevelmethod.hh:85
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:50
A LeveTransferPolicy that used aggregation to construct the coarse level system.
Definition: twolevelmethod.hh:141
C Criterion
Definition: twolevelmethod.hh:145
AggregationLevelTransferPolicy * clone() const
Clone the current object.
Definition: twolevelmethod.hh:212
AggregationLevelTransferPolicy(const Criterion &crit)
Definition: twolevelmethod.hh:148
void moveToFineLevel(typename FatherType::FineDomainType &fineLhs)
Updates the fine level linear system after the correction of the coarse levels system.
Definition: twolevelmethod.hh:205
void moveToCoarseLevel(const typename FatherType::FineRangeType &fineRhs)
Definition: twolevelmethod.hh:198
SequentialInformation ParallelInformation
Definition: twolevelmethod.hh:146
LevelTransferPolicy< O, O > FatherType
Definition: twolevelmethod.hh:144
void createCoarseLevelSystem(const O &fineOperator)
Algebraically creates the coarse level system.
Definition: twolevelmethod.hh:152
A policy class for solving the coarse level system using one step of AMG.
Definition: twolevelmethod.hh:232
OneStepAMGCoarseSolverPolicy(const SmootherArgs &args, const Criterion &c)
Constructs the coarse solver policy.
Definition: twolevelmethod.hh:251
AMGInverseOperator CoarseLevelSolver
The type of solver constructed for the coarse level.
Definition: twolevelmethod.hh:313
OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy &other)
Copy constructor.
Definition: twolevelmethod.hh:255
O::range_type X
The type of the range and domain of the operator.
Definition: twolevelmethod.hh:237
C Criterion
The type of the crition used for the aggregation within AMG.
Definition: twolevelmethod.hh:239
Dune::Amg::SmootherTraits< S >::Arguments SmootherArgs
The type of the arguments used for constructing the smoother.
Definition: twolevelmethod.hh:243
O Operator
The type of the linear operator used.
Definition: twolevelmethod.hh:235
CoarseLevelSolver * createCoarseLevelSolver(P &transferPolicy)
Constructs a coarse level solver.
Definition: twolevelmethod.hh:323
AMG< Operator, X, Smoother > AMGType
The type of the AMG construct on the coarse level.
Definition: twolevelmethod.hh:245
S Smoother
The type of the smoother used in AMG.
Definition: twolevelmethod.hh:241
Definition: twolevelmethod.hh:351
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethod.hh:378
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethod.hh:369
TwoLevelMethod(const TwoLevelMethod &other)
Definition: twolevelmethod.hh:416
void pre(FineDomainType &x, FineRangeType &b)
Definition: twolevelmethod.hh:429
FO FineOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:361
CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver
The type of the coarse level solver.
Definition: twolevelmethod.hh:356
void apply(FineDomainType &v, const FineRangeType &d)
Definition: twolevelmethod.hh:437
CSP CoarseLevelSolverPolicy
The type of the policy for constructing the coarse level solver.
Definition: twolevelmethod.hh:354
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethod.hh:382
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:402
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: twolevelmethod.hh:464
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethod.hh:365
~TwoLevelMethod()
Definition: twolevelmethod.hh:422
CSP::Operator CoarseOperatorType
The linear operator of the finel level system. Has to be derived from AssembledLinearOperator.
Definition: twolevelmethod.hh:374
void post([[maybe_unused]] FineDomainType &x)
Definition: twolevelmethod.hh:434
S SmootherType
The type of the fine level smoother.
Definition: twolevelmethod.hh:386
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