1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #include <config.h>
4
5 #include <dune/common/parallel/mpihelper.hh>
6 #include <dune/common/parallel/mpicommunication.hh>
7 #include <dune/istl/paamg/matrixhierarchy.hh>
8 #include <dune/istl/paamg/smoother.hh>
9 #include <dune/istl/preconditioners.hh>
10 #include <dune/istl/owneroverlapcopy.hh>
11 #include <dune/istl/schwarz.hh>
12 #include "anisotropic.hh"
13
14 template<int blockSize>
testHierarchy(int N)15 void testHierarchy(int N)
16 {
17 typedef int LocalId;
18 typedef int GlobalId;
19 typedef Dune::OwnerOverlapCopyCommunication<LocalId,GlobalId> Communication;
20 typedef Communication::ParallelIndexSet ParallelIndexSet;
21 typedef Dune::FieldMatrix<double,blockSize,blockSize> MatrixBlock;
22 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
23 typedef Dune::FieldVector<double,blockSize> VectorBlock;
24 typedef Dune::BlockVector<VectorBlock> Vector;
25
26 int n;
27 Communication pinfo(Dune::MPIHelper::getCommunicator());
28 ParallelIndexSet& indices = pinfo.indexSet();
29
30 typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
31 RemoteIndices& remoteIndices = pinfo.remoteIndices();
32
33 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, Dune::MPIHelper::getCollectiveCommunication(), &n);
34 Vector b(indices.size());
35
36 remoteIndices.rebuild<false>();
37
38 typedef Dune::NegateSet<Communication::OwnerSet> OverlapFlags;
39 typedef Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication> Operator;
40 typedef Dune::Amg::MatrixHierarchy<Operator,Communication> Hierarchy;
41 typedef Dune::Amg::Hierarchy<Vector> VHierarchy;
42
43 Operator op(mat, pinfo);
44 Hierarchy hierarchy(
45 Dune::stackobject_to_shared_ptr(op),
46 Dune::stackobject_to_shared_ptr(pinfo));
47 VHierarchy vh(
48 Dune::stackobject_to_shared_ptr(b));
49
50 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
51 Criterion;
52
53 Criterion criterion(100,4);
54 Dune::Timer timer;
55
56 hierarchy.template build<OverlapFlags>(criterion);
57 hierarchy.coarsenVector(vh);
58
59 std::cout<<"Building hierarchy took "<<timer.elapsed()<<std::endl;
60
61 std::cout<<"=== Vector hierarchy has "<<vh.levels()<<" levels! ==="<<std::endl;
62 timer.reset();
63
64 hierarchy.recalculateGalerkin(OverlapFlags());
65
66 std::cout<<"Recalculation took "<<timer.elapsed()<<std::endl;
67 std::vector<std::size_t> data;
68
69 hierarchy.getCoarsestAggregatesOnFinest(data);
70 }
71
72
main(int argc,char ** argv)73 int main(int argc, char** argv)
74 {
75 Dune::MPIHelper::instance(argc, argv);
76
77 constexpr int blockSize = 1;
78 int N=10;
79
80 if(argc>1)
81 N = atoi(argv[1]);
82
83 testHierarchy<blockSize>(N);
84 }
85