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, [[maybe_unused]] double reduction, [[maybe_unused]] InverseOperatorResult& res)
288  {
289  if(first_)
290  {
291  amg_.pre(x,b);
292  first_=false;
293  x_=x;
294  }
295  amg_.apply(x,b);
296  }
297 
298  void apply(X& x, X& b, InverseOperatorResult& res)
299  {
300  return apply(x,b,1e-8,res);
301  }
302 
303  virtual SolverCategory::Category category() const
304  {
305  return amg_.category();
306  }
307 
308  ~AMGInverseOperator()
309  {
310  if(!first_)
311  amg_.post(x_);
312  }
313  AMGInverseOperator(const AMGInverseOperator& other)
314  : x_(other.x_), amg_(other.amg_), first_(other.first_)
315  {
316  }
317  private:
318  X x_;
319  AMGType amg_;
320  bool first_;
321  };
322 
323 public:
325  typedef AMGInverseOperator CoarseLevelSolver;
326 
334  template<class P>
336  {
337  coarseOperator_=transferPolicy.getCoarseLevelOperator();
338  AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
339  criterion_,
340  smootherArgs_);
341 
342  return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
343 
344  }
345 
346 private:
348  std::shared_ptr<Operator> coarseOperator_;
350  SmootherArgs smootherArgs_;
352  Criterion criterion_;
353 };
354 
360 template<class FO, class CSP, class S>
362  public Preconditioner<typename FO::domain_type, typename FO::range_type>
363 {
364 public:
368  typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
373  typedef FO FineOperatorType;
377  typedef typename FineOperatorType::range_type FineRangeType;
381  typedef typename FineOperatorType::domain_type FineDomainType;
386  typedef typename CSP::Operator CoarseOperatorType;
390  typedef typename CoarseOperatorType::range_type CoarseRangeType;
394  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
398  typedef S SmootherType;
399 
415  std::shared_ptr<SmootherType> smoother,
417  CoarseOperatorType>& policy,
418  CoarseLevelSolverPolicy& coarsePolicy,
419  std::size_t preSteps=1, std::size_t postSteps=1)
420  : operator_(&op), smoother_(smoother),
421  preSteps_(preSteps), postSteps_(postSteps)
422  {
423  policy_ = policy.clone();
424  policy_->createCoarseLevelSystem(*operator_);
425  coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_);
426  }
427 
429  : operator_(other.operator_), coarseSolver_(new CoarseLevelSolver(*other.coarseSolver_)),
430  smoother_(other.smoother_), policy_(other.policy_->clone()),
431  preSteps_(other.preSteps_), postSteps_(other.postSteps_)
432  {}
433 
435  {
436  // Each instance has its own policy.
437  delete policy_;
438  delete coarseSolver_;
439  }
440 
441  void updatePreconditioner(FineOperatorType& /* op */,
442  std::shared_ptr<SmootherType> smoother,
443  CoarseLevelSolverPolicy& coarsePolicy)
444  {
445  updatePreconditioner(smoother, coarsePolicy);
446  }
447 
448  void updatePreconditioner(std::shared_ptr<SmootherType> smoother,
449  CoarseLevelSolverPolicy& coarsePolicy)
450  {
451  //assume new matrix is not reallocated the new precondition should anyway be made
452  smoother_ = smoother;
453  if (coarseSolver_) {
454  policy_->calculateCoarseEntries(*operator_);
455  coarsePolicy.setCoarseOperator(*policy_);
456  coarseSolver_->updatePreconditioner();
457  } else {
458  // we should probably not be heere
459  policy_->createCoarseLevelSystem(*operator_);
460  coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_);
461  }
462  }
463 
464  void pre(FineDomainType& x, FineRangeType& b)
465  {
466  smoother_->pre(x,b);
467  }
468 
469  void post([[maybe_unused]] FineDomainType& x)
470  {
471  }
472 
473  void apply(FineDomainType& v, const FineRangeType& d)
474  {
475  FineDomainType u(v);
476  FineRangeType rhs(d);
477  LevelContext context;
478  SequentialInformation info;
479  context.pinfo=&info;
480  context.lhs=&u;
481  context.update=&v;
482  context.smoother=smoother_;
483  context.rhs=&rhs;
484  context.matrix=operator_;
485  // Presmoothing
486  presmooth(context, preSteps_);
487  //Coarse grid correction
488  policy_->moveToCoarseLevel(*context.rhs);
489  InverseOperatorResult res;
490  coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res);
491  *context.lhs=0;
492  policy_->moveToFineLevel(*context.lhs);
493  *context.update += *context.lhs;
494  // Postsmoothing
495  postsmooth(context, postSteps_);
496 
497  }
498 // //! Category of the preconditioner (see SolverCategory::Category)
499  virtual SolverCategory::Category category() const
500  {
501  return SolverCategory::sequential;
502  }
503 
504 private:
508  struct LevelContext
509  {
511  typedef S SmootherType;
513  std::shared_ptr<SmootherType> smoother;
515  FineDomainType* lhs;
516  /*
517  * @brief The right hand side holding the current residual.
518  *
519  * This is passed to the smoother as the right hand side.
520  */
521  FineRangeType* rhs;
527  FineDomainType* update;
529  SequentialInformation* pinfo;
535  const FineOperatorType* matrix;
536  };
537  const FineOperatorType* operator_;
539  CoarseLevelSolver* coarseSolver_;
541  std::shared_ptr<S> smoother_;
543  LevelTransferPolicyCpr<FO,typename CSP::Operator>* policy_;
545  std::size_t preSteps_;
547  std::size_t postSteps_;
548 };
549 }// end namespace Amg
550 }// end namespace Dune
551 
553 #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:325
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:335
O Operator
The type of the linear operator used.
Definition: twolevelmethodcpr.hh:248
Definition: twolevelmethodcpr.hh:363
S SmootherType
The type of the fine level smoother.
Definition: twolevelmethodcpr.hh:398
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethodcpr.hh:377
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:414
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethodcpr.hh:381
CoarseOperatorType::domain_type CoarseDomainType
The type of the domain of the coarse level operator.
Definition: twolevelmethodcpr.hh:394
CSP::Operator CoarseOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:386
CSP CoarseLevelSolverPolicy
The type of the policy for constructing the coarse level solver.
Definition: twolevelmethodcpr.hh:366
CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver
The type of the coarse level solver.
Definition: twolevelmethodcpr.hh:368
CoarseOperatorType::range_type CoarseRangeType
The type of the range of the coarse level operator.
Definition: twolevelmethodcpr.hh:390
FO FineOperatorType
The linear operator of the finel level system.
Definition: twolevelmethodcpr.hh:373