Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Christian Wehner
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__VANKA__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__VANKA__
35 :
36 : #include "common/util/smart_pointer.h"
37 : #include "lib_algebra/operator/interface/preconditioner.h"
38 :
39 : #ifdef UG_PARALLEL
40 : #include "pcl/pcl_util.h"
41 : #include "lib_algebra/parallelization/parallelization_util.h"
42 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
43 : #endif
44 :
45 : namespace ug{
46 :
47 : static const size_t MAXBLOCKSIZE = 53;
48 :
49 : template<typename Matrix_type, typename Vector_type>
50 0 : bool Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
51 : {
52 : DenseVector< VariableArray1<number> > s;
53 : DenseVector< VariableArray1<number> > localx;
54 : DenseMatrix< VariableArray2<number> > mat;
55 :
56 : size_t blockind[MAXBLOCKSIZE];
57 :
58 0 : for(size_t i=0; i < x.size(); i++)
59 : {
60 0 : x[i]=0;
61 : };
62 :
63 0 : for(size_t i=0; i < x.size(); i++)
64 : {
65 0 : if (A(i,i)!=0){
66 : /* // do usual gauss-seidel (would be needed in case of Dirichlet pressure condition)
67 : typename Vector_type::value_type def = b[i];
68 :
69 : for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
70 : // s -= it.value() * x[it.index()];
71 : MatMultAdd(def, 1.0, def, -1.0, it.value(), x[it.index()]);
72 : // x[i] = s/A(i,i)
73 : InverseMatMult(x[i], relax, A(i,i), s);*/
74 0 : continue;
75 : };
76 :
77 : size_t blocksize=0;
78 :
79 0 : for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
80 0 : blockind[blocksize] = it.index();
81 0 : blocksize++;
82 : };
83 : if (blocksize>MAXBLOCKSIZE) UG_THROW("MAXBLOCKSIZE too small\n");
84 0 : mat.resize(blocksize,blocksize);
85 0 : s.resize(blocksize);
86 0 : localx.resize(blocksize);
87 0 : for (size_t j=0;j<blocksize;j++){
88 : // fill local block matrix
89 0 : for (size_t k=0;k<blocksize;k++){
90 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
91 : };
92 : // compute rhs
93 0 : typename Vector_type::value_type sj = b[blockind[j]];
94 0 : for(typename Matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
95 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
96 : };
97 : s.subassign(j,sj);
98 : };
99 : // solve block
100 0 : InverseMatMult(localx,1,mat,s);
101 0 : for (size_t j=0;j<blocksize;j++){
102 0 : x[blockind[j]] += relax*localx[j];
103 : };
104 : }
105 :
106 0 : return true;
107 : }
108 :
109 : // Diagonal Vanka block smoother:
110 : // When setting up the local block matrix the side-diagonal entries are left away, except for the pressure.
111 : // The local block matrix therefore has the form
112 : //
113 : // x 0 ... 0 x
114 : // 0 x 0 ... 0 x
115 : // 0 0 x 0 ... x
116 : // ...
117 : // 0 ... 0 x x
118 : // x x ... x x x
119 : //
120 : // The velocity off-diagonal entries are considered in the local defect vector.
121 : // The local block system can be solved in O(n) time so that a step of this smoother is computationly cheaper than the
122 : // full Vanka smoother.
123 : template<typename Matrix_type, typename Vector_type>
124 0 : bool Diag_Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
125 : {
126 : typedef typename Vector_type::value_type vector_block_type;
127 : DenseVector< VariableArray1<vector_block_type> > s;
128 : DenseVector< VariableArray1<vector_block_type> > localx;
129 : DenseMatrix< VariableArray2<number> > mat;
130 0 : s.resize(MAXBLOCKSIZE);
131 : typedef typename Matrix_type::value_type block_type;
132 :
133 : size_t blockind[MAXBLOCKSIZE];
134 :
135 0 : for(size_t i=0; i < x.size(); i++)
136 : {
137 0 : x[i]=0;
138 : };
139 :
140 0 : for(size_t i=0; i < x.size(); i++)
141 : {
142 0 : if (A(i,i)!=0){
143 : /* // do usual gauss-seidel (would be needed in case of Dirichlet pressure condition)
144 : typename Vector_type::value_type def = b[i];
145 :
146 : for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
147 : // s -= it.value() * x[it.index()];
148 : MatMultAdd(def, 1.0, def, -1.0, it.value(), x[it.index()]);
149 : // x[i] = s/A(i,i)
150 : InverseMatMult(x[i], relax, A(i,i), def);*/
151 0 : continue;
152 : };
153 :
154 : size_t blocksize=0;
155 :
156 0 : for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
157 0 : if (it.index()==i) continue;
158 0 : blockind[blocksize] = it.index();
159 0 : s[blocksize] = b[blockind[blocksize]];
160 0 : for(typename Matrix_type::const_row_iterator rowit = A.begin_row(blockind[blocksize]); rowit != A.end_row(blockind[blocksize]) ; ++rowit){
161 0 : if ((rowit.index()==blockind[blocksize])||(rowit.index()==i)) continue;
162 : // s[blocksize] -= a_ij*x_j
163 0 : MatMultAdd(s[blocksize], 1.0, s[blocksize], -1.0, rowit.value(), x[rowit.index()]);
164 : };
165 0 : blocksize++;
166 : };
167 0 : if (blocksize>MAXBLOCKSIZE) UG_THROW("MAXBLOCKSIZE too small\n");
168 : // remark: blocksize is without pressure variable, so actual blocksize is blocksize+1
169 0 : block_type a_ii = A(i,i);
170 0 : typename Vector_type::value_type s_i = b[i];
171 : // Gauss elimination on local block matrix
172 0 : for (size_t j=0;j<blocksize;j++){
173 0 : block_type a_q = A(i,blockind[j]);
174 0 : block_type a_jj = A(blockind[j],blockind[j]);
175 0 : a_q /= a_jj;
176 : // s_i -= a_ij/a_jj*s_j
177 0 : MatMultAdd(s_i, 1.0, s_i, -1.0, a_q, s[j]);
178 : // a_ii -= a_ij/a_jj*a_ji
179 0 : a_ii-=a_q*A(blockind[j],i);
180 : }
181 : // solve diagonalized system
182 : // x[i] = s_i/a_ii
183 : InverseMatMult(x[i], 1.0, a_ii, s_i);
184 0 : for (size_t j=0;j<blocksize;j++){
185 : // s_j-=a_ji*x_i
186 0 : MatMultAdd(s[j], 1.0, s[j], -1.0, A(blockind[j],i), x[i]);
187 : // x_j=1/a_jj*s_j
188 0 : InverseMatMult(x[blockind[j]], relax, A(blockind[j],blockind[j]),s[j]);
189 : }
190 : }
191 0 : return true;
192 : }
193 :
194 :
195 : /// Vanka Preconditioner
196 : template <typename TAlgebra>
197 : class Vanka : public IPreconditioner<TAlgebra>
198 : {
199 : public:
200 : /// Algebra type
201 : typedef TAlgebra algebra_type;
202 :
203 : /// Vector type
204 : typedef typename TAlgebra::vector_type vector_type;
205 :
206 : /// Matrix type
207 : typedef typename TAlgebra::matrix_type matrix_type;
208 :
209 : /// Matrix Operator type
210 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
211 :
212 : /// Base type
213 : typedef IPreconditioner<TAlgebra> base_type;
214 :
215 : protected:
216 : using base_type::set_debug;
217 : using base_type::debug_writer;
218 : using base_type::write_debug;
219 :
220 : public:
221 : /// default constructor
222 0 : Vanka() {m_relax=1;};
223 :
224 : /// Clone
225 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
226 : {
227 0 : SmartPtr<Vanka<algebra_type> > newInst(new Vanka<algebra_type>());
228 0 : newInst->set_debug(debug_writer());
229 0 : newInst->set_damp(this->damping());
230 0 : return newInst;
231 : }
232 :
233 : /// Destructor
234 0 : virtual ~Vanka()
235 0 : {};
236 :
237 0 : virtual bool supports_parallel() const {return true;}
238 :
239 : /// set relaxation parameter
240 : public:
241 0 : void set_relax(number omega){m_relax=omega;};
242 :
243 : protected:
244 : number m_relax;
245 :
246 : protected:
247 : /// Name of preconditioner
248 0 : virtual const char* name() const {return "Vanka";}
249 :
250 : /// Preprocess routine
251 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
252 : {
253 : #ifdef UG_PARALLEL
254 : if(pcl::NumProcs() > 1)
255 : {
256 : // copy original matrix
257 : MakeConsistent(*pOp, m_A);
258 : // set zero on slaves
259 : std::vector<IndexLayout::Element> vIndex;
260 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
261 : SetDirichletRow(m_A, vIndex);
262 : }
263 : #endif
264 0 : return true;
265 : }
266 :
267 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
268 : {
269 : #ifdef UG_PARALLEL
270 : if(pcl::NumProcs() > 1)
271 : {
272 : // make defect unique
273 : // todo: change that copying
274 : vector_type dhelp;
275 : dhelp.resize(d.size()); dhelp = d;
276 : dhelp.change_storage_type(PST_UNIQUE);
277 :
278 : if(!Vanka_step(m_A, c, dhelp, m_relax)) return false;
279 :
280 : c.set_storage_type(PST_UNIQUE);
281 : return true;
282 : }
283 : else
284 : #endif
285 : {
286 0 : if(!Vanka_step(*pOp, c, d, m_relax)) return false;
287 :
288 : #ifdef UG_PARALLEL
289 : c.set_storage_type(PST_UNIQUE);
290 : #endif
291 : return true;
292 : }
293 : }
294 :
295 : /// Postprocess routine
296 0 : virtual bool postprocess() {return true;}
297 :
298 : protected:
299 : #ifdef UG_PARALLEL
300 : matrix_type m_A;
301 : #endif
302 :
303 : };
304 :
305 : /// Diagvanka Preconditioner, description see above diagvanka_step function
306 : template <typename TAlgebra>
307 : class DiagVanka : public IPreconditioner<TAlgebra>
308 : {
309 : public:
310 : /// Algebra type
311 : typedef TAlgebra algebra_type;
312 :
313 : /// Vector type
314 : typedef typename TAlgebra::vector_type vector_type;
315 :
316 : /// Matrix type
317 : typedef typename TAlgebra::matrix_type matrix_type;
318 :
319 : /// Matrix Operator type
320 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
321 :
322 : /// Base type
323 : typedef IPreconditioner<TAlgebra> base_type;
324 :
325 : protected:
326 : using base_type::set_debug;
327 : using base_type::debug_writer;
328 : using base_type::write_debug;
329 :
330 : public:
331 : /// default constructor
332 0 : DiagVanka() {m_relax=1;};
333 :
334 : /// Clone
335 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
336 : {
337 0 : SmartPtr<DiagVanka<algebra_type> > newInst(new DiagVanka<algebra_type>());
338 0 : newInst->set_debug(debug_writer());
339 0 : newInst->set_damp(this->damping());
340 0 : return newInst;
341 : }
342 :
343 0 : virtual bool supports_parallel() const {return true;}
344 :
345 : /// Destructor
346 0 : virtual ~DiagVanka()
347 0 : {};
348 :
349 : /// set relaxation parameter
350 : public:
351 0 : void set_relax(number omega){m_relax=omega;};
352 :
353 : protected:
354 : number m_relax;
355 :
356 : protected:
357 : /// Name of preconditioner
358 0 : virtual const char* name() const {return "DiagVanka";}
359 :
360 : /// Preprocess routine
361 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
362 : {
363 : #ifdef UG_PARALLEL
364 : if(pcl::NumProcs() > 1)
365 : {
366 : // copy original matrix
367 : MakeConsistent(*pOp, m_A);
368 : // set zero on slaves
369 : std::vector<IndexLayout::Element> vIndex;
370 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
371 : SetDirichletRow(m_A, vIndex);
372 : }
373 : #endif
374 0 : return true;
375 : }
376 :
377 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
378 : {
379 : #ifdef UG_PARALLEL
380 : if(pcl::NumProcs() > 1)
381 : {
382 : // make defect unique
383 : // todo: change that copying
384 : vector_type dhelp;
385 : dhelp.resize(d.size()); dhelp = d;
386 : dhelp.change_storage_type(PST_UNIQUE);
387 :
388 : if(!Diag_Vanka_step(m_A, c, dhelp, m_relax)) return false;
389 :
390 : c.set_storage_type(PST_UNIQUE);
391 : return true;
392 : }
393 : else
394 : #endif
395 : {
396 :
397 0 : if(!Diag_Vanka_step(*pOp, c, d, m_relax)) return false;
398 :
399 : #ifdef UG_PARALLEL
400 : c.set_storage_type(PST_UNIQUE);
401 : #endif
402 : return true;
403 : }
404 : }
405 :
406 : /// Postprocess routine
407 0 : virtual bool postprocess() {return true;}
408 :
409 : protected:
410 : #ifdef UG_PARALLEL
411 : matrix_type m_A;
412 : #endif
413 :
414 : };
415 :
416 :
417 : } // end namespace ug
418 :
419 : #endif
|
