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 #ifdef TEST_AGGLO
5 #define UNKNOWNS 10
6 #endif
7 #include "anisotropic.hh"
8 #include <dune/common/timer.hh>
9 #include <dune/common/parallel/indexset.hh>
10 #include <dune/common/parallel/mpicommunication.hh>
11 #include <dune/istl/paamg/amg.hh>
12 #include <dune/istl/paamg/pinfo.hh>
13 #include <dune/istl/schwarz.hh>
14 #include <dune/istl/owneroverlapcopy.hh>
15 #include <string>
16
17 template<class T, class C>
18 class DoubleStepPreconditioner
19 : public Dune::Preconditioner<typename T::domain_type, typename T::range_type>
20 {
21 public:
22 typedef typename T::domain_type X;
23 typedef typename T::range_type Y;
24
25 enum {category = T::category};
26
DoubleStepPreconditioner(T & preconditioner_,C & comm)27 DoubleStepPreconditioner(T& preconditioner_, C& comm)
28 : preconditioner(&preconditioner_), comm_(comm)
29 {}
30
pre(X & x,Y & b)31 virtual void pre (X& x, Y& b)
32 {
33 preconditioner->pre(x,b);
34 }
35
apply(X & v,const Y & d)36 virtual void apply(X& v, const Y& d)
37 {
38 preconditioner->apply(v,d);
39 comm_.copyOwnerToAll(v,v);
40 }
41
post(X & x)42 virtual void post (X& x)
43 {
44 preconditioner->post(x);
45 }
46 private:
47 T* preconditioner;
48 C& comm_;
49 };
50
51
52 class MPIError {
53 public:
54 /** @brief Constructor. */
MPIError(std::string s,int e)55 MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
56 /** @brief The error string. */
57 std::string errorstring;
58 /** @brief The mpi error code. */
59 int errorcode;
60 };
61
MPI_err_handler(MPI_Comm * comm,int * err_code,...)62 void MPI_err_handler([[maybe_unused]] MPI_Comm *comm, int *err_code, ...)
63 {
64 char *err_string=new char[MPI_MAX_ERROR_STRING];
65 int err_length;
66 MPI_Error_string(*err_code, err_string, &err_length);
67 std::string s(err_string, err_length);
68 std::cerr << "An MPI Error occurred:"<<std::endl<<s<<std::endl;
69 delete[] err_string;
70 throw MPIError(s, *err_code);
71 }
72
73 template<int BS>
testAmg(int N,int coarsenTarget)74 void testAmg(int N, int coarsenTarget)
75 {
76 std::cout<<"==================================================="<<std::endl;
77 std::cout<<"BS="<<BS<<" N="<<N<<" coarsenTarget="<<coarsenTarget<<std::endl;
78
79 int procs, rank;
80 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
81 MPI_Comm_size(MPI_COMM_WORLD, &procs);
82
83 typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
84 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
85 typedef Dune::FieldVector<double,BS> VectorBlock;
86 typedef Dune::BlockVector<VectorBlock> Vector;
87 typedef int GlobalId;
88 typedef Dune::OwnerOverlapCopyCommunication<GlobalId> Communication;
89 typedef Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication> Operator;
90 int n;
91
92 N/=BS;
93
94 Communication comm(MPI_COMM_WORLD);
95
96 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, comm.indexSet(), comm.communicator(), &n, 1);
97
98 const BCRSMat& cmat = mat;
99
100 comm.remoteIndices().template rebuild<false>();
101
102 Vector b(cmat.N()), x(cmat.M());
103
104 b=0;
105 x=100;
106
107 setBoundary(x, b, N, comm.indexSet());
108
109 Vector b1=b, x1=x;
110
111 if(N<=6) {
112 std::ostringstream name;
113 name<<rank<<": row";
114
115 Dune::printmatrix(std::cout, cmat, "A", name.str().c_str());
116 Dune::printvector(std::cout, x, "x", name.str().c_str());
117 //Dune::printvector(std::cout, b, "b", name.str().c_str());
118 //Dune::printvector(std::cout, b1, "b1", "row");
119 //Dune::printvector(std::cout, x1, "x1", "row");
120 }
121
122 Dune::Timer watch;
123
124 watch.reset();
125 Operator fop(cmat, comm);
126
127 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
128 Criterion;
129 typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
130 //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
131 //typedef Dune::SeqILU0<BCRSMat,Vector,Vector> Smoother;
132 //typedef Dune::SeqILUn<BCRSMat,Vector,Vector> Smoother;
133 typedef Dune::BlockPreconditioner<Vector,Vector,Communication,Smoother> ParSmoother;
134 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
135
136 Dune::OverlappingSchwarzScalarProduct<Vector,Communication> sp(comm);
137
138 Dune::InverseOperatorResult r, r1;
139
140 double buildtime;
141
142 SmootherArgs smootherArgs;
143
144 smootherArgs.iterations = 1;
145
146
147 Criterion criterion(15,coarsenTarget);
148 criterion.setDefaultValuesIsotropic(2);
149
150
151 typedef Dune::Amg::AMG<Operator,Vector,ParSmoother,Communication> AMG;
152
153 AMG amg(fop, criterion, smootherArgs, comm);
154
155 buildtime = watch.elapsed();
156
157 if(rank==0)
158 std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
159
160 Dune::CGSolver<Vector> amgCG(fop, sp, amg, 10e-8, 300, (rank==0) ? 2 : 0);
161 watch.reset();
162
163 amgCG.apply(x,b,r);
164 amg.recalculateHierarchy();
165
166
167 MPI_Barrier(MPI_COMM_WORLD);
168 double solvetime = watch.elapsed();
169
170 b=0;
171 x=100;
172
173 setBoundary(x, b, N, comm.indexSet());
174
175 Dune::CGSolver<Vector> amgCG1(fop, sp, amg, 10e-8, 300, (rank==0) ? 2 : 0);
176 amgCG1.apply(x,b,r);
177
178 if(!r.converged && rank==0)
179 std::cerr<<" AMG Cg solver did not converge!"<<std::endl;
180
181 if(rank==0) {
182 std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
183
184 std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
185 std::cout<<"AMG building together with slving took "<<buildtime+solvetime<<std::endl;
186 }
187
188 }
189
190 template<int BSStart, int BSEnd, int BSStep=1>
191 struct AMGTester
192 {
testAMGTester193 static void test(int N, int coarsenTarget)
194 {
195 testAmg<BSStart>(N, coarsenTarget);
196 const int next = (BSStart+BSStep>BSEnd) ? BSEnd : BSStart+BSStep;
197 AMGTester<next,BSEnd,BSStep>::test(N, coarsenTarget);
198 }
199 }
200 ;
201
202 template<int BSStart,int BSStep>
203 struct AMGTester<BSStart,BSStart,BSStep>
204 {
testAMGTester205 static void test(int N, int coarsenTarget)
206 {
207 testAmg<BSStart>(N, coarsenTarget);
208 }
209 };
210
211
main(int argc,char ** argv)212 int main(int argc, char** argv)
213 {
214 MPI_Init(&argc, &argv);
215 MPI_Errhandler handler;
216 MPI_Comm_create_errhandler(MPI_err_handler, &handler);
217 MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler);
218
219 int N=100;
220
221 int coarsenTarget=200;
222
223 if(argc>1)
224 N = atoi(argv[1]);
225
226 if(argc>2)
227 coarsenTarget = atoi(argv[2]);
228
229 #ifdef TEST_AGGLO
230 N=UNKNOWNS;
231 #endif
232 AMGTester<1,1>::test(N, coarsenTarget);
233 //AMGTester<1,5>::test(N, coarsenTarget);
234 // AMGTester<10,10>::test(N, coarsenTarget);
235
236 MPI_Finalize();
237 }
238