1 #include"config.h"
2 #include "anisotropic.hh"
3 #include <dune/common/timer.hh>
4 #include <dune/common/shared_ptr.hh>
5 #include <dune/common/parallel/indexset.hh>
6 #include <dune/common/parallel/communication.hh>
7 #include <dune/istl/paamg/twolevelmethod.hh>
8 #include <dune/istl/overlappingschwarz.hh>
9 #include <dune/istl/paamg/pinfo.hh>
10 #include <dune/istl/solvers.hh>
11
12 template<class M, class V>
randomize(const M & mat,V & b)13 void randomize(const M& mat, V& b)
14 {
15 V x=b;
16
17 srand((unsigned)std::clock());
18
19 typedef typename V::iterator iterator;
20 for(iterator i=x.begin(); i != x.end(); ++i)
21 *i=(rand() / (RAND_MAX + 1.0));
22
23 mat.mv(static_cast<const V&>(x), b);
24 }
25
testTwoLevelMethod()26 void testTwoLevelMethod()
27 {
28 const int BS=1;
29 int N=100;
30 typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
31 ParallelIndexSet indices;
32 typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
33 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
34 typedef Dune::FieldVector<double,BS> VectorBlock;
35 typedef Dune::BlockVector<VectorBlock> Vector;
36 typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
37 typedef Dune::CollectiveCommunication<void*> Comm;
38 Comm c;
39 int n;
40 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, c, &n, 1);
41 Vector b(mat.N()), x(mat.M());
42 x=0;
43 randomize(mat, b);
44 #ifndef USE_OVERLAPPINGSCHWARZ
45 // Smoother used for the finest level
46 typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> FSmoother;
47 FSmoother fineSmoother(mat,1,1.0);
48 #else
49 // Smoother used for the finest level. There is an additional
50 // template parameter to provide the subdomain solver.
51 typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,
52 Dune::SymmetricMultiplicativeSchwarzMode> FSmoother;
53 // Type of the subdomain vector
54 typedef FSmoother::subdomain_vector SubdomainVector;
55 // Create subdomain.
56 std::size_t stride=2;
57 SubdomainVector subdomains((((N-1)/stride)+1)*(((N-1)/stride)+1));
58
59 for(int i=0; i<N; ++i)
60 for(int j=0; j<N; ++j)
61 {
62 int index=i/stride*(((N-1)/stride)+1)+j/stride;
63 subdomains[index].insert(i*N+j);
64 }
65 //create smoother
66 FSmoother fineSmoother(mat,subdomains, 1.0, false);
67 #endif
68 typedef Dune::SeqJac<BCRSMat,Vector,Vector> CSmoother;
69 typedef Dune::Amg::SmootherTraits<CSmoother>::Arguments SmootherArgs;
70 typedef Dune::Amg::CoarsenCriterion<
71 Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> > Criterion;
72 typedef Dune::Amg::AggregationLevelTransferPolicy<Operator,Criterion> TransferPolicy;
73 typedef Dune::Amg::OneStepAMGCoarseSolverPolicy<Operator,CSmoother, Criterion> CoarsePolicy; // Policy for coarse solver creation
74 Criterion crit;
75 CoarsePolicy coarsePolicy=CoarsePolicy(SmootherArgs(), crit);
76 TransferPolicy transferPolicy(crit);
77 Operator fop(mat);
78 Dune::Amg::TwoLevelMethod<Operator,CoarsePolicy,FSmoother> preconditioner(fop,
79 Dune::stackobject_to_shared_ptr(fineSmoother),
80 transferPolicy,
81 coarsePolicy);
82 Dune::GeneralizedPCGSolver<Vector> amgCG(fop,preconditioner,.8,80,2);
83 Dune::Amg::TwoLevelMethod<Operator,CoarsePolicy,FSmoother> preconditioner1(preconditioner);
84 Dune::InverseOperatorResult res;
85 amgCG.apply(x,b,res);
86 }
87
main()88 int main()
89 {
90 testTwoLevelMethod();
91 return 0;
92 }
93