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
28namespace Dune
29{
30namespace Amg
31{
32
42template<class FO, class CO>
44{
45public:
50 typedef FO FineOperatorType;
54 typedef typename FineOperatorType::range_type FineRangeType;
58 typedef typename FineOperatorType::domain_type FineDomainType;
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
151template<class O, class C>
153 : public LevelTransferPolicyCpr<O,O>
154{
155 typedef Dune::Amg::AggregatesMap<typename O::matrix_type::size_type> AggregatesMap;
156public:
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
230private:
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
243template<class O, class S, class C>
245{
246public:
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 {}
272private:
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
323public:
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
346private:
348 std::shared_ptr<Operator> coarseOperator_;
350 SmootherArgs smootherArgs_;
352 Criterion criterion_;
353};
354
360template<class FO, class CSP, class S>
362 public Preconditioner<typename FO::domain_type, typename FO::range_type>
363{
364public:
368 typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
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
434 ~TwoLevelMethodCpr()
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
504private:
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
CoarseDomainType & getCoarseLevelLhs()
Get the coarse level left hand side.
Definition: twolevelmethodcpr.hh:93
virtual LevelTransferPolicyCpr * clone() const =0
Clone the current object.
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
CoarseRangeType & getCoarseLevelRhs()
Get the coarse level right hand side.
Definition: twolevelmethodcpr.hh:84
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.
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
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
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
CoarseLevelSolver * createCoarseLevelSolver(P &transferPolicy)
Constructs a coarse level solver.
Definition: twolevelmethodcpr.hh:335
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
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