My Project
twolevelmethodcpr.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_TWOLEVELMETHODCPR_HH
4 #define DUNE_ISTL_TWOLEVELMETHODCPR_HH
5 
6 // NOTE: This file is a modified version of dune/istl/paamg/twolevelmethod.hh from
7 // dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
8 
9 #include <tuple>
10 
11 #include<dune/istl/operators.hh>
12 //#include "amg.hh"
13 //#include"galerkin.hh"
14 #include<dune/istl/paamg/amg.hh>
15 #include<dune/istl/paamg/galerkin.hh>
16 #include<dune/istl/solver.hh>
17 
18 #include<dune/common/unused.hh>
19 #include<dune/common/version.hh>
20 
28 namespace Dune
29 {
30 namespace Amg
31 {
32 
42 template<class FO, class CO>
44 {
45 public:
50  typedef FO FineOperatorType;
54  typedef typename FineOperatorType::range_type FineRangeType;
58  typedef typename FineOperatorType::domain_type FineDomainType;
63  typedef CO CoarseOperatorType;
67  typedef typename CoarseOperatorType::range_type CoarseRangeType;
71  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
76  std::shared_ptr<CoarseOperatorType>& getCoarseLevelOperator()
77  {
78  return operator_;
79  }
85  {
86  return rhs_;
87  }
88 
94  {
95  return lhs_;
96  }
106  virtual void moveToCoarseLevel(const FineRangeType& fineRhs)=0;
116  virtual void moveToFineLevel(FineDomainType& fineLhs)=0;
124  virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0;
125 
129  virtual void calculateCoarseEntries(const FineOperatorType& fineOperator) = 0;
130 
132  virtual LevelTransferPolicyCpr* clone() const =0;
133 
136 
137  protected:
143  std::shared_ptr<CoarseOperatorType> operator_;
144 };
145 
151 template<class O, class C>
153  : public LevelTransferPolicyCpr<O,O>
154 {
155  typedef Dune::Amg::AggregatesMap<typename O::matrix_type::size_type> AggregatesMap;
156 public:
158  typedef C Criterion;
159  typedef SequentialInformation ParallelInformation;
160 
161  AggregationLevelTransferPolicyCpr(const Criterion& crit)
162  : criterion_(crit)
163  {}
164 
165  void createCoarseLevelSystem(const O& fineOperator)
166  {
167  prolongDamp_ = criterion_.getProlongationDampingFactor();
168  GalerkinProduct<ParallelInformation> productBuilder;
169  typedef typename Dune::Amg::MatrixGraph<const typename O::matrix_type> MatrixGraph;
170  typedef typename Dune::Amg::PropertiesGraph<MatrixGraph,Dune::Amg::VertexProperties,
171  Dune::Amg::EdgeProperties,Dune::IdentityMap,Dune::IdentityMap> PropertiesGraph;
172  MatrixGraph mg(fineOperator.getmat());
173  PropertiesGraph pg(mg,Dune::IdentityMap(),Dune::IdentityMap());
174  typedef NegateSet<typename ParallelInformation::OwnerSet> OverlapFlags;
175 
176  aggregatesMap_.reset(new AggregatesMap(pg.maxVertex()+1));
177 
178  int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
179 
180  std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
181  aggregatesMap_->buildAggregates(fineOperator.getmat(), pg, criterion_, true);
182  std::cout<<"no aggregates="<<noAggregates<<" iso="<<isoAggregates<<" one="<<oneAggregates<<" skipped="<<skippedAggregates<<std::endl;
183  // misuse coarsener to renumber aggregates
184  Dune::Amg::IndicesCoarsener<Dune::Amg::SequentialInformation,int> renumberer;
185  typedef std::vector<bool>::iterator Iterator;
186  typedef Dune::IteratorPropertyMap<Iterator, Dune::IdentityMap> VisitedMap;
187  std::vector<bool> excluded(fineOperator.getmat().N(), false);
188  VisitedMap vm(excluded.begin(), Dune::IdentityMap());
189  ParallelInformation pinfo;
190  std::size_t aggregates = renumberer.coarsen(pinfo, pg, vm,
191  *aggregatesMap_, pinfo,
192  noAggregates);
193  std::vector<bool>& visited=excluded;
194 
195  typedef std::vector<bool>::iterator Iterator;
196 
197  for(Iterator iter= visited.begin(), end=visited.end();
198  iter != end; ++iter)
199  *iter=false;
200  matrix_.reset(productBuilder.build(mg, vm,
201  SequentialInformation(),
202  *aggregatesMap_,
203  aggregates,
204  OverlapFlags()));
205  productBuilder.calculate(fineOperator.getmat(), *aggregatesMap_, *matrix_, pinfo, OverlapFlags());
206  this->lhs_.resize(this->matrix_->M());
207  this->rhs_.resize(this->matrix_->N());
208  this->operator_.reset(new O(*matrix_));
209  }
210 
211  void moveToCoarseLevel(const typename FatherType::FineRangeType& fineRhs)
212  {
213  Transfer<std::size_t,typename FatherType::FineRangeType,ParallelInformation>
214  ::restrictVector(*aggregatesMap_, this->rhs_, fineRhs, ParallelInformation());
215  this->lhs_=0;
216  }
217 
219  {
220  Transfer<std::size_t,typename FatherType::FineRangeType,ParallelInformation>
221  ::prolongateVector(*aggregatesMap_, this->lhs_, fineLhs,
222  prolongDamp_, ParallelInformation());
223  }
224 
226  {
227  return new AggregationLevelTransferPolicyCpr(*this);
228  }
229 
230 private:
231  typename O::matrix_type::field_type prolongDamp_;
232  std::shared_ptr<AggregatesMap> aggregatesMap_;
233  Criterion criterion_;
234  std::shared_ptr<typename O::matrix_type> matrix_;
235 };
236 
243 template<class O, class S, class C>
245 {
246 public:
248  typedef O Operator;
250  typedef typename O::range_type X;
252  typedef C Criterion;
254  typedef S Smoother;
256  typedef typename Dune::Amg::SmootherTraits<S>::Arguments SmootherArgs;
258  typedef AMG<Operator,X,Smoother> AMGType;
265  : smootherArgs_(args), criterion_(c)
266  {}
269  : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
270  criterion_(other.criterion_)
271  {}
272 private:
279  struct AMGInverseOperator : public InverseOperator<X,X>
280  {
281  AMGInverseOperator(const typename AMGType::Operator& op,
282  const Criterion& crit,
283  const typename AMGType::SmootherArgs& args)
284  : amg_(op, crit,args), first_(true)
285  {}
286 
287  void apply(X& x, X& b, double reduction, InverseOperatorResult& res)
288  {
289  DUNE_UNUSED_PARAMETER(reduction);
290  DUNE_UNUSED_PARAMETER(res);
291  if(first_)
292  {
293  amg_.pre(x,b);
294  first_=false;
295  x_=x;
296  }
297  amg_.apply(x,b);
298  }
299 
300  void apply(X& x, X& b, InverseOperatorResult& res)
301  {
302  return apply(x,b,1e-8,res);
303  }
304 
305  virtual SolverCategory::Category category() const
306  {
307  return amg_.category();
308  }
309 
310  ~AMGInverseOperator()
311  {
312  if(!first_)
313  amg_.post(x_);
314  }
315  AMGInverseOperator(const AMGInverseOperator& other)
316  : x_(other.x_), amg_(other.amg_), first_(other.first_)
317  {
318  }
319  private:
320  X x_;
321  AMGType amg_;
322  bool first_;
323  };
324 
325 public:
327  typedef AMGInverseOperator CoarseLevelSolver;
328 
336  template<class P>
338  {
339  coarseOperator_=transferPolicy.getCoarseLevelOperator();
340  AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
341  criterion_,
342  smootherArgs_);
343 
344  return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
345 
346  }
347 
348 private:
350  std::shared_ptr<Operator> coarseOperator_;
352  SmootherArgs smootherArgs_;
354  Criterion criterion_;
355 };
356 
362 template<class FO, class CSP, class S>
364  public Preconditioner<typename FO::domain_type, typename FO::range_type>
365 {
366 public:
370  typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
375  typedef FO FineOperatorType;
379  typedef typename FineOperatorType::range_type FineRangeType;
383  typedef typename FineOperatorType::domain_type FineDomainType;
388  typedef typename CSP::Operator CoarseOperatorType;
392  typedef typename CoarseOperatorType::range_type CoarseRangeType;
396  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
400  typedef S SmootherType;
401 
417  std::shared_ptr<SmootherType> smoother,
419  CoarseOperatorType>& policy,
420  CoarseLevelSolverPolicy& coarsePolicy,
421  std::size_t preSteps=1, std::size_t postSteps=1)
422  : operator_(&op), smoother_(smoother),
423  preSteps_(preSteps), postSteps_(postSteps)
424  {
425  policy_ = policy.clone();
426  policy_->createCoarseLevelSystem(*operator_);
427  coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_);
428  }
429 
431  : operator_(other.operator_), coarseSolver_(new CoarseLevelSolver(*other.coarseSolver_)),
432  smoother_(other.smoother_), policy_(other.policy_->clone()),
433  preSteps_(other.preSteps_), postSteps_(other.postSteps_)
434  {}
435 
437  {
438  // Each instance has its own policy.
439  delete policy_;
440  delete coarseSolver_;
441  }
442 
443  void updatePreconditioner(FineOperatorType& /* op */,
444  std::shared_ptr<SmootherType> smoother,
445  CoarseLevelSolverPolicy& coarsePolicy)
446  {
447  updatePreconditioner(smoother, coarsePolicy);
448  }
449 
450  void updatePreconditioner(std::shared_ptr<SmootherType> smoother,
451  CoarseLevelSolverPolicy& coarsePolicy)
452  {
453  //assume new matrix is not reallocated the new precondition should anyway be made
454  smoother_ = smoother;
455  if (coarseSolver_) {
456  policy_->calculateCoarseEntries(*operator_);
457  coarsePolicy.setCoarseOperator(*policy_);
458  coarseSolver_->updatePreconditioner();
459  } else {
460  // we should probably not be heere
461  policy_->createCoarseLevelSystem(*operator_);
462  coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_);
463  }
464  }
465 
466  void pre(FineDomainType& x, FineRangeType& b)
467  {
468  smoother_->pre(x,b);
469  }
470 
471  void post(FineDomainType& x)
472  {
473  DUNE_UNUSED_PARAMETER(x);
474  }
475 
476  void apply(FineDomainType& v, const FineRangeType& d)
477  {
478  FineDomainType u(v);
479  FineRangeType rhs(d);
480  LevelContext context;
481  SequentialInformation info;
482  context.pinfo=&info;
483  context.lhs=&u;
484  context.update=&v;
485  context.smoother=smoother_;
486  context.rhs=&rhs;
487  context.matrix=operator_;
488  // Presmoothing
489  presmooth(context, preSteps_);
490  //Coarse grid correction
491  policy_->moveToCoarseLevel(*context.rhs);
492  InverseOperatorResult res;
493  coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res);
494  *context.lhs=0;
495  policy_->moveToFineLevel(*context.lhs);
496  *context.update += *context.lhs;
497  // Postsmoothing
498  postsmooth(context, postSteps_);
499 
500  }
501 // //! Category of the preconditioner (see SolverCategory::Category)
502  virtual SolverCategory::Category category() const
503  {
504  return SolverCategory::sequential;
505  }
506 
507 private:
511  struct LevelContext
512  {
514  typedef S SmootherType;
516  std::shared_ptr<SmootherType> smoother;
518  FineDomainType* lhs;
519  /*
520  * @brief The right hand side holding the current residual.
521  *
522  * This is passed to the smoother as the right hand side.
523  */
524  FineRangeType* rhs;
530  FineDomainType* update;
532  SequentialInformation* pinfo;
538  const FineOperatorType* matrix;
539  };
540  const FineOperatorType* operator_;
542  CoarseLevelSolver* coarseSolver_;
544  std::shared_ptr<S> smoother_;
546  LevelTransferPolicyCpr<FO,typename CSP::Operator>* policy_;
548  std::size_t preSteps_;
550  std::size_t postSteps_;
551 };
552 }// end namespace Amg
553 }// end namespace Dune
554 
556 #endif
A LeveTransferPolicy that used aggregation to construct the coarse level system.
Definition: twolevelmethodcpr.hh:154
void moveToFineLevel(typename FatherType::FineDomainType &fineLhs)
Updates the fine level linear system after the correction of the coarse levels system.
Definition: twolevelmethodcpr.hh:218
void createCoarseLevelSystem(const O &fineOperator)
Algebraically creates the coarse level system.
Definition: twolevelmethodcpr.hh:165
AggregationLevelTransferPolicyCpr * clone() const
Clone the current object.
Definition: twolevelmethodcpr.hh:225
Abstract base class for transfer between levels and creation of the coarse level system.
Definition: twolevelmethodcpr.hh:44
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethodcpr.hh:67
FO FineOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:50
virtual void createCoarseLevelSystem(const FineOperatorType &fineOperator)=0
Algebraically creates the coarse level system.
virtual void moveToCoarseLevel(const FineRangeType &fineRhs)=0
Transfers the data to the coarse level.
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethodcpr.hh:71
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethodcpr.hh:54
virtual ~LevelTransferPolicyCpr()
Destructor.
Definition: twolevelmethodcpr.hh:135
virtual void moveToFineLevel(FineDomainType &fineLhs)=0
Updates the fine level linear system after the correction of the coarse levels system.
CoarseDomainType & getCoarseLevelLhs()
Get the coarse level left hand side.
Definition: twolevelmethodcpr.hh:93
virtual void calculateCoarseEntries(const FineOperatorType &fineOperator)=0
???.
CoarseDomainType lhs_
The coarse level lhs.
Definition: twolevelmethodcpr.hh:141
CoarseRangeType rhs_
The coarse level rhs.
Definition: twolevelmethodcpr.hh:139
virtual LevelTransferPolicyCpr * clone() const =0
Clone the current object.
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethodcpr.hh:58
CO CoarseOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:63
std::shared_ptr< CoarseOperatorType > & getCoarseLevelOperator()
Get the coarse level operator.
Definition: twolevelmethodcpr.hh:76
CoarseRangeType & getCoarseLevelRhs()
Get the coarse level right hand side.
Definition: twolevelmethodcpr.hh:84
std::shared_ptr< CoarseOperatorType > operator_
the coarse level linear operator.
Definition: twolevelmethodcpr.hh:143
A policy class for solving the coarse level system using one step of AMG.
Definition: twolevelmethodcpr.hh:245
AMG< Operator, X, Smoother > AMGType
The type of the AMG construct on the coarse level.
Definition: twolevelmethodcpr.hh:258
OneStepAMGCoarseSolverPolicyCpr(const OneStepAMGCoarseSolverPolicyCpr &other)
Copy constructor.
Definition: twolevelmethodcpr.hh:268
AMGInverseOperator CoarseLevelSolver
The type of solver constructed for the coarse level.
Definition: twolevelmethodcpr.hh:327
S Smoother
The type of the smoother used in AMG.
Definition: twolevelmethodcpr.hh:254
C Criterion
The type of the crition used for the aggregation within AMG.
Definition: twolevelmethodcpr.hh:252
Dune::Amg::SmootherTraits< S >::Arguments SmootherArgs
The type of the arguments used for constructing the smoother.
Definition: twolevelmethodcpr.hh:256
O::range_type X
The type of the range and domain of the operator.
Definition: twolevelmethodcpr.hh:250
OneStepAMGCoarseSolverPolicyCpr(const SmootherArgs &args, const Criterion &c)
Constructs the coarse solver policy.
Definition: twolevelmethodcpr.hh:264
CoarseLevelSolver * createCoarseLevelSolver(P &transferPolicy)
Constructs a coarse level solver.
Definition: twolevelmethodcpr.hh:337
O Operator
The type of the linear operator used.
Definition: twolevelmethodcpr.hh:248
Definition: twolevelmethodcpr.hh:365
S SmootherType
The type of the fine level smoother.
Definition: twolevelmethodcpr.hh:400
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethodcpr.hh:379
TwoLevelMethodCpr(const FineOperatorType &op, std::shared_ptr< SmootherType > smoother, const LevelTransferPolicyCpr< FineOperatorType, CoarseOperatorType > &policy, CoarseLevelSolverPolicy &coarsePolicy, std::size_t preSteps=1, std::size_t postSteps=1)
Constructs a two level method.
Definition: twolevelmethodcpr.hh:416
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethodcpr.hh:383
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethodcpr.hh:396
CSP::Operator CoarseOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:388
CSP CoarseLevelSolverPolicy
The type of the policy for constructing the coarse level solver.
Definition: twolevelmethodcpr.hh:368
CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver
The type of the coarse level solver.
Definition: twolevelmethodcpr.hh:370
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethodcpr.hh:392
FO FineOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:375