PetscSolver.cc 24.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
#include <vector>
#include <set>

Thomas Witkowski's avatar
Thomas Witkowski committed
16
#include "parallel/PetscSolver.h"
17
18
#include "parallel/StdMpi.h"
#include "parallel/ParallelDebug.h"
19
20
21
#include "DOFVector.h"
#include "Debug.h"
#include "SystemVector.h"
22

23
24
25
26
#include "petscksp.h"

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
27
28
  using namespace std;

29
30
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
  {    
31
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
32
      cout << "[0]  Petsc-Iteration " << iter << ": " << rnorm << endl;
33
34
35
36

    return 0;
  }
 
37

Thomas Witkowski's avatar
Thomas Witkowski committed
38
  void PetscSolver::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
39
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
40
    FUNCNAME("PetscSolver::solve()");
41

42
43
    TEST_EXIT(meshDistributor)("Should not happen!\n");

44
    double wtime = MPI::Wtime();
45

46
    fillPetscMatrix(systemMatrix, rhs);
47
    solvePetscMatrix(*solution, adaptInfo);   
48

49
50
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  MPI::Wtime() - wtime);
51
52
53
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
  void PetscSolver::setDofMatrix(DOFMatrix* mat, int dispMult, 
				 int dispAddRow, int dispAddCol)
56
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
57
    FUNCNAME("PetscSolver::setDofMatrix()");
58
59
60
61
62
63
64
65
66
67
68
69
70

    TEST_EXIT(mat)("No DOFMatrix!\n");

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

Thomas Witkowski's avatar
Thomas Witkowski committed
71
72
    vector<int> cols;
    vector<double> values;
73
74
    cols.reserve(300);
    values.reserve(300);
75
    
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    vector<int> globalCols;
77

78

79
    // === Traverse all rows of the dof matrix and insert row wise the values ===
Thomas Witkowski's avatar
Merge    
Thomas Witkowski committed
80
    // === to the PETSc matrix.                                               ===
81
82
83
84
85

    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

      // Global index of the current row dof.
86
      int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
87
      // Test if the current row dof is a periodic dof.
88
89
      bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof);

Thomas Witkowski's avatar
Thomas Witkowski committed
90
      if (!periodicRow) {
Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
93
	// === Row DOF index is not periodic. ===

	// Calculate PETSc row index.
Thomas Witkowski's avatar
Merge    
Thomas Witkowski committed
94

Thomas Witkowski's avatar
Thomas Witkowski committed
95
	int rowIndex = globalRowDof * dispMult + dispAddRow;
96

97
98
99
	cols.clear();
	values.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
102

Thomas Witkowski's avatar
Thomas Witkowski committed
103
104
105
106
	  // Global index of the current column index.
	  int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
	  // Test if the current col dof is a periodic dof.
	  bool periodicCol = meshDistributor->isPeriodicDof(globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
	  // Calculate PETSc col index.
	  int colIndex = globalColDof * dispMult + dispAddCol;
109

Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
	  // Ignore all zero entries, expect it is a diagonal entry.
	  //	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
112
113
// 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
// 	    continue;
114

Thomas Witkowski's avatar
Thomas Witkowski committed
115
	  if (!periodicCol) {
Thomas Witkowski's avatar
Thomas Witkowski committed
116
	    // Calculate the exact position of the column index in the PETSc matrix.
117
118
	    cols.push_back(globalColDof * dispMult + dispAddCol);
	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
119
	  } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
	    // === Row index is not periodic, but column index is. ===
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
	    // Create set of all periodic associations of the column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
123
	    std::set<int> perAsc;
Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
	    std::set<int>& perColAsc = 
	      meshDistributor->getPerDofAssociations(globalColDof);
126
127
	    for (std::set<int>::iterator it = perColAsc.begin(); 
		 it != perColAsc.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
128
129
	      if (*it >= -3)
		perAsc.insert(*it);
130
    
Thomas Witkowski's avatar
Thomas Witkowski committed
131
	    // Scale value to the number of periodic associations of the column index.
132
	    double scaledValue = 
Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
137
138
139
140
141
	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));

	    
	    // === Create set of all matrix column indices due to the periodic ===
	    // === associations of the column DOF index.                       ===

	    vector<int> newCols;

	    // First, add the original matrix index.
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
	    newCols.push_back(globalColDof);
	    
Thomas Witkowski's avatar
Thomas Witkowski committed
144
	    // And add all periodic matrix indices.
145
146
	    for (std::set<int>::iterator it = perAsc.begin(); 
		 it != perAsc.end(); ++it) {
147
148
149
150
	      int nCols = static_cast<int>(newCols.size());

	      for (int i = 0; i < nCols; i++) {
 		TEST_EXIT_DBG(meshDistributor->isPeriodicDof(newCols[i], *it))
Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
 		  ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
		   *it, newCols[i]);

154
		newCols.push_back(meshDistributor->getPeriodicMapping(newCols[i], *it));
Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
	      }
	    }
157

158
	    for (unsigned int i = 0; i < newCols.size(); i++) {
159
160
	      cols.push_back(newCols[i] * dispMult + dispAddCol);
	      values.push_back(scaledValue);	      
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
	    }
	  }
163
	}
164
165
166

	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
		     &(cols[0]), &(values[0]), ADD_VALUES);	
167
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
	// === Row DOF index is periodic. ===
169

Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
173
174
175
176
177
	// Because this row is periodic, we will have to add the entries of this 
	// matrix row to multiple rows. The following maps store to each row an
	// array of column indices and values of the entries that must be added to
	// the PETSc matrix.
	map<int, vector<int> > colsMap;
	map<int, vector<double> > valsMap;

	// Traverse all column entries.
Thomas Witkowski's avatar
Thomas Witkowski committed
178
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
179
	     icursor != icend; ++icursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
182
183

	  // Global index of the current column index.
	  int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));

Thomas Witkowski's avatar
Thomas Witkowski committed
184
	  // Ignore all zero entries, expect it is a diagonal entry.
185
186
// 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
// 	    continue;
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
190
191

	  // === Add all periodic associations of both, the row and the column ===
	  // === indices to the set perAsc.                                    ===

Thomas Witkowski's avatar
Thomas Witkowski committed
192
193
	  std::set<int> perAsc;

194
195
196
197
198
199
200
201
202
	  if (meshDistributor->isPeriodicDof(globalColDof)) {
	    std::set<int>& perColAsc = 
	      meshDistributor->getPerDofAssociations(globalColDof);
	    for (std::set<int>::iterator it = perColAsc.begin(); 
		 it != perColAsc.end(); ++it)
	      if (*it >= -3)
		perAsc.insert(*it);
	  }

203
204
205
206
207
208
	  std::set<int>& perRowAsc = 
	    meshDistributor->getPerDofAssociations(globalRowDof);
	  for (std::set<int>::iterator it = perRowAsc.begin(); 
	       it != perRowAsc.end(); ++it)
	    if (*it >= -3)
	      perAsc.insert(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
	  // Scale the value with respect to the number of periodic associations.
211
	  double scaledValue = 
Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
214
215
216
217
218
219
220
221
	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));


	  // === Create all matrix entries with respect to the periodic  ===
	  // === associations of the row and column indices.             ===

	  vector<pair<int, int> > entry;
	  
	  // First, add the original entry.
	  entry.push_back(make_pair(globalRowDof, globalColDof));
Thomas Witkowski's avatar
Thomas Witkowski committed
222

Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
	  // Then, traverse the periodic associations of the row and column indices
	  // and create the corresponding entries.
Thomas Witkowski's avatar
Thomas Witkowski committed
225
226
227
	  for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
	    int nEntry = static_cast<int>(entry.size());
	    for (int i = 0; i < nEntry; i++) {
228
	      int perRowDof = 0;
229
	      if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].first))
230
		perRowDof = meshDistributor->getPeriodicMapping(entry[i].first, *it);
Thomas Witkowski's avatar
Thomas Witkowski committed
231
	      else
232
		perRowDof = entry[i].first;
Thomas Witkowski's avatar
Thomas Witkowski committed
233

234
	      int perColDof;
235
	      if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].second))
236
		perColDof = meshDistributor->getPeriodicMapping(entry[i].second, *it);
Thomas Witkowski's avatar
Thomas Witkowski committed
237
	      else
238
		perColDof = entry[i].second;	      	      
239
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241
	      entry.push_back(make_pair(perRowDof, perColDof));
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243
244
	    }
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
247
248

	  // === Translate the matrix entries to PETSc's matrix.

	  for (vector<pair<int, int> >::iterator eIt = entry.begin(); 
249
	       eIt != entry.end(); ++eIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
250
	    // Calculate row and column indices of the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
	    int rowIndex = eIt->first * dispMult + dispAddRow;
	    int colIndex = eIt->second * dispMult + dispAddCol;
Thomas Witkowski's avatar
Thomas Witkowski committed
253

254
255
	    colsMap[rowIndex].push_back(colIndex);
	    valsMap[rowIndex].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
	  }
	}
258
259


Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
262
	// === Finally, add all periodic rows to the PETSc matrix. ===

	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
263
264
265
266
267
268
269
270
	     rowIt != colsMap.end(); ++rowIt) {
	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
	    ("Should not happen!\n");

	  int rowIndex = rowIt->first;
	  MatSetValues(petscMatrix, 1, &rowIndex, rowIt->second.size(),
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      }
272
273
274
275
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
  void PetscSolver::setDofVector(Vec& petscVec, DOFVector<double>* vec, 
				 int dispMult, int dispAdd)
278
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    FUNCNAME("PetscSolver::setDofVector()");
280

281
    // Traverse all used DOFs in the dof vector.
282
283
284
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
      // Calculate global row index of the dof.
285
      DegreeOfFreedom globalRowDof = 
286
	meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex());
Thomas Witkowski's avatar
Merge    
Thomas Witkowski committed
287
      // Calculate PETSc index of the row dof.
288
      int index = globalRowDof * dispMult + dispAdd;
289

290
291
292
      if (meshDistributor->isPeriodicDof(globalRowDof)) {
	std::set<int>& perAsc = meshDistributor->getPerDofAssociations(globalRowDof);
	double value = *dofIt / (perAsc.size() + 1.0);
293
294
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);

295
	for (std::set<int>::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) {
296
	  int mappedDof = meshDistributor->getPeriodicMapping(globalRowDof, *perIt);
297
298
	  int mappedIndex = mappedDof * dispMult + dispAdd;
	  VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES);
299
300
301
302
303
304
	}
      } else {
	// The dof index is not periodic.
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
305
    }
306
307
308
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
309
  void PetscSolver::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
310
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
311
    FUNCNAME("PetscSolver::createPetscNnzStructure()");
312
313
314
315

    TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
    TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");

316
    int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
317
318
319
320
321
322
323
324
325
326
    d_nnz = new int[nRankRows];
    o_nnz = new int[nRankRows];
    for (int i = 0; i < nRankRows; i++) {
      d_nnz[i] = 0;
      o_nnz[i] = 0;
    }

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
    typedef vector<pair<int, int> > MatrixNnzEntry;
    typedef map<int, DofContainer> RankToDofContainer;
329
330

    // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
331
    // that this rank will send to. These nnz entries will be assembled on this rank,
332
333
    // but because the row DOFs are not DOFs of this rank they will be send to the
    // owner of the row DOFs.
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    map<int, MatrixNnzEntry> sendMatrixEntry;
335

336
337
    // Maps to each DOF that must be send to another rank the rank number of the
    // receiving rank.
Thomas Witkowski's avatar
Thomas Witkowski committed
338
    map<DegreeOfFreedom, int> sendDofToRank;
339

340
341

    // First, create for all ranks we send data to MatrixNnzEntry object with 0 entries.
342
343
    for (RankToDofContainer::iterator it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
344
345
      sendMatrixEntry[it->first].resize(0);

346
347
348
349
350
351
352
353
354
355
356
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	sendDofToRank[**dofIt] = it->first;
    }


    std::set<int> recvFromRank;
    for (RankToDofContainer::iterator it = meshDistributor->getSendDofs().begin();
	 it != meshDistributor->getSendDofs().end(); ++it)
      recvFromRank.insert(it->first);

357

358
359
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
360
361
362
363
 	if (!(*mat)[i][j])
	  continue;

	Matrix bmat = (*mat)[i][j]->getBaseMatrix();
364

365
366
	traits::col<Matrix>::type col(bmat);
	traits::const_value<Matrix>::type value(bmat);
367
	  
368
369
370
371
372
	typedef traits::range_generator<row, Matrix>::type cursor_type;
	typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	
	for (cursor_type cursor = begin<row>(bmat), 
	       cend = end<row>(bmat); cursor != cend; ++cursor) {
373
	  
374
	  int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
375

Thomas Witkowski's avatar
Thomas Witkowski committed
376
	  vector<int> rows;
377
	  rows.push_back(globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
378
	  vector<int> rowRank;
379
380
381
382
383
384
385
386
	  if (meshDistributor->getIsRankDof(*cursor)) {
	    rowRank.push_back(meshDistributor->getMpiRank());
	  } else {
	    // Find out who is the member of this DOF.
	    TEST_EXIT_DBG(sendDofToRank.count(*cursor))("Should not happen!\n");
	    
	    rowRank.push_back(sendDofToRank[*cursor]);
	  }
387

388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
	  // Map the local row number to the global DOF index and create from it
	  // the global PETSc row index of this DOF.
	  
	  int petscRowIdx = globalRowDof * nComponents + i;
	  
	  if (meshDistributor->getIsRankDof(*cursor)) {
	    	    
	    // === The current row DOF is a rank dof, so create the corresponding ===
	    // === nnz values directly on rank's nnz data.                        ===
	    
	    // This is the local row index of the local PETSc matrix.
	    int localPetscRowIdx = 
	      petscRowIdx - meshDistributor->getRstart() * nComponents;
	    
	    TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
	      ("Should not happen! \n Debug info: localRowIdx = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   nCompontens = %d   nRankRows = %d\n",
	       *cursor, meshDistributor->mapLocalToGlobal(*cursor), petscRowIdx, localPetscRowIdx, meshDistributor->getRstart(), nComponents, nRankRows);
	    
	    
	    // Traverse all non zero entries in this row.
	    for (icursor_type icursor = begin<nz>(cursor), 
		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
	      int petscColIdx = 
		meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
412
	      
413
	      //	      if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
414
415
416
417
418
419
420
		// The row DOF is a rank DOF, if also the column is a rank DOF, 
		// increment the d_nnz values for this row, otherwise the o_nnz value.
		if (petscColIdx >= meshDistributor->getRstart() * nComponents && 
		    petscColIdx < meshDistributor->getRstart() * nComponents + nRankRows)
		  d_nnz[localPetscRowIdx]++;
		else
		  o_nnz[localPetscRowIdx]++;
421
		//	      }    
422
423
424
425
426
427
428
429
430
431
432
433
	    }
	  } else {
	    // === The current row DOF is not a rank dof, i.e., it will be created ===
	    // === on this rank, but after this it will be send to another rank    ===
	    // === matrix. So we need to send also the corresponding nnz structure ===
	    // === of this row to the corresponding rank.                          ===
	    
	    // Send all non zero entries to the member of the row DOF.
	    int sendToRank = sendDofToRank[*cursor];
	    
	    for (icursor_type icursor = begin<nz>(cursor), 
		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
434
	      //	      if (value(*icursor) != 0.0) {
435
436
		int petscColIdx = 
		  meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
437
438
		
		sendMatrixEntry[sendToRank].
Thomas Witkowski's avatar
Thomas Witkowski committed
439
		  push_back(make_pair(petscRowIdx, petscColIdx));
440
		//	      }
441
442
443
444
	    }
	    
	  } // if (isRankDof[*cursor]) ... else ...
	} // for each row in mat[i][j]
445
446
447
448
449
      } 
    }

    // === Send and recv the nnz row structure to/from other ranks. ===

450
    StdMpi<MatrixNnzEntry> stdMpi(meshDistributor->getMpiComm(), true);
451
    stdMpi.send(sendMatrixEntry);
452
453
454
    for (std::set<int>::iterator it = recvFromRank.begin(); 
	 it != recvFromRank.end(); ++it)
      stdMpi.recv(*it);
455
    stdMpi.startCommunication();
456

457

458
459
460
    // === Evaluate the nnz structure this rank got from other ranks and add it to ===
    // === the PETSc nnz data structure.                                           ===

Thomas Witkowski's avatar
Thomas Witkowski committed
461
    for (map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
462
463
464
465
466
467
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second.size() > 0) {
	for (unsigned int i = 0; i < it->second.size(); i++) {
	  int r = it->second[i].first;
	  int c = it->second[i].second;

468
	  int localRowIdx = r - meshDistributor->getRstart() * nComponents;
469
470
471
472
473

	  TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
	    ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
	     r, localRowIdx, nRankRows, it->first);
	  
474
475
	  if (c < meshDistributor->getRstart() * nComponents || 
	      c >= meshDistributor->getRstart() * nComponents + nRankRows)
476
477
478
479
480
	    o_nnz[localRowIdx]++;
	  else
	    d_nnz[localRowIdx]++;
	}
      }
481
    }
482
483
484
485
486
487
488
489
490
491

    // The above algorithm for calculating the number of nnz per row over-
    // approximates the value, i.e., the number is always equal or larger to 
    // the real number of nnz values in the global parallel matrix. For small
    // matrices, the problem may arise, that the result is larger than the
    // number of elements in a row. This is fixed in the following.

    if (nRankRows < 100) 
      for (int i = 0; i < nRankRows; i++)
	d_nnz[i] = std::min(d_nnz[i], nRankRows);
492
493
494
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
495
  void PetscSolver::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
496
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
497
    FUNCNAME("PetscSolver::fillPetscMatrix()");
498

499
    double wtime = MPI::Wtime();
500
501
    int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
    int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents;
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516

    // === Create PETSc vector (rhs, solution and a temporary vector). ===

    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    VecSetSizes(petscSolVec, nRankRows, nOverallRows);
    VecSetType(petscSolVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
    VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
    VecSetType(petscTmpVec, VECMPI);

517
518
519
520
521
    int recvAllValues = 0;
    int sendValue = static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    meshDistributor->getMpiComm().Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

    if (!d_nnz || recvAllValues != 0) {
522
523
      if (d_nnz) {
	delete [] d_nnz;
524
	d_nnz = NULL;
525
	delete [] o_nnz;
526
	o_nnz = NULL;
527
528
      }

529
      createPetscNnzStructure(mat);
530
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();
531
    }
532

533

534
535
536
    // === Create PETSc matrix with the computed nnz data structure. ===

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
537
		    0, d_nnz, 0, o_nnz, &petscMatrix);
538
    
539
#if (DEBUG != 0)
540
    INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
541
#endif
542
543
544
545

#if (DEBUG != 0)
    int a, b;
    MatGetOwnershipRange(petscMatrix, &a, &b);
546
547
548
549
    TEST_EXIT(a == meshDistributor->getRstart() * nComponents)
      ("Wrong matrix ownership range!\n");
    TEST_EXIT(b == meshDistributor->getRstart() * nComponents + nRankRows)
      ("Wrong matrix ownership range!\n");
550
551
#endif

552

553
554
555
556
557
    // === Transfer values from DOF matrices to the PETSc matrix. === 

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j])
558
559
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

560
#if (DEBUG != 0)
561
    INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
562
#endif
563
564
565
566

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

567
#if (DEBUG != 0)
568
569
    INFO(info, 8)("Fill petsc matrix 3 needed %.5f seconds\n", 
		  TIME_USED(MPI::Wtime(), wtime));
570
#endif
571

572
573
574
575
576
577
578
579
    // === Transfer values from DOF vector to the PETSc vector. === 

    for (int i = 0; i < nComponents; i++)
      setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

580
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime);
581
582
583
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
584
  void PetscSolver::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
585
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
586
    FUNCNAME("PetscSolver::solvePetscMatrix()");
587
588
589
590
591
592
593
594
595
596

#if 0
    // Set old solution to be initiual guess for petsc solver.
    for (int i = 0; i < nComponents; i++)
      setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i);

    VecAssemblyBegin(petscSolVec);
    VecAssemblyEnd(petscSolVec);
#endif

597
    // === Init PETSc solver. ===
598
599

    KSP solver;
600
601
602
603
    PC pc;

    providePetscSolver(solver, pc);

604
605
606
    // Do not delete the solution vector, use it for the initial guess.
    //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);

607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
    
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPGetPC(solver, &pc);
    PCSetType(pc, PCFIELDSPLIT);

#if 0
    typedef map<int, DofContainer> RankToDofContainer;
    typedef map<DegreeOfFreedom, bool> DofIndexToBool;

    std::set<DegreeOfFreedom> boundaryDofsSet;
    std::set<DegreeOfFreedom> boundaryLocalDofs;
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator rankIt = sendDofs.begin(); rankIt != sendDofs.end(); ++rankIt)
      for (DofContainer::iterator dofIt = rankIt->second.begin(); dofIt != rankIt->second.end(); ++dofIt) {
	boundaryLocalDofs.insert(**dofIt);

	for (int i = 0; i < nComponents; i++)
	  boundaryDofsSet.insert(meshDistributor->mapLocalToGlobal(**dofIt) * nComponents + i);
      }

    vector<DegreeOfFreedom> boundaryDofs(boundaryDofsSet.begin(), boundaryDofsSet.end());
628

629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
    std::set<DegreeOfFreedom> otherBoundaryLocalDofs;
    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator rankIt = recvDofs.begin(); rankIt != recvDofs.end(); ++rankIt)
      for (DofContainer::iterator dofIt = rankIt->second.begin(); dofIt != rankIt->second.end(); ++dofIt)
	otherBoundaryLocalDofs.insert(**dofIt);

    std::set<DegreeOfFreedom> interiorDofsSet;
    DofIndexToBool& isRankDof = meshDistributor->getIsRankDof();
    for (DofIndexToBool::iterator dofIt = isRankDof.begin(); dofIt != isRankDof.end(); ++dofIt)
      if (dofIt->second && 
	  boundaryLocalDofs.count(dofIt->first) == 0 && 
	  otherBoundaryLocalDofs.count(dofIt->first) == 0)
	for (int i = 0; i < nComponents; i++)	    
	  interiorDofsSet.insert(meshDistributor->mapLocalToGlobal(dofIt->first) * nComponents + i);

    vector<DegreeOfFreedom> interiorDofs(interiorDofsSet.begin(), interiorDofsSet.end());

    IS interiorIs;
    ISCreateGeneral(PETSC_COMM_WORLD, interiorDofs.size(), &(interiorDofs[0]), PETSC_COPY_VALUES, &interiorIs);
    PCFieldSplitSetIS(pc, "interior", interiorIs);

    IS boundaryIs;
    ISCreateGeneral(PETSC_COMM_WORLD, boundaryDofs.size(), &(boundaryDofs[0]), PETSC_COPY_VALUES, &boundaryIs);
    PCFieldSplitSetIS(pc, "boundary", boundaryIs);
#endif

    // === Run PETSc. ===
656
657
658

    KSPSolve(solver, petscRhsVec, petscSolVec);

659

660
    // === Transfere values from PETSc's solution vectors to the dof vectors.
661

662
663
664
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

665
    int nRankDofs = meshDistributor->getNumberRankDofs();
666
    for (int i = 0; i < nComponents; i++) {
667
      DOFVector<double> &dofvec = *(vec.getDOFVector(i));
668
      for (int j = 0; j < nRankDofs; j++)
669
	dofvec[meshDistributor->mapLocalToDofIndex(j)] = 
670
	  vecPointer[j * nComponents + i]; 
671
672
673
674
675
    }

    VecRestoreArray(petscSolVec, &vecPointer);


676
    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to more ===
677
    // === than one partition.                                                 ===
678
    meshDistributor->synchVector(vec);
679
680
681
682
683
684
685


    // === Print information about solution process. ===

    int iterations = 0;
    KSPGetIterationNumber(solver, &iterations);
    MSG("  Number of iterations: %d\n", iterations);
686
687
    adaptInfo->setSolverIterations(iterations);

688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
    double norm = 0.0;
    MatMult(petscMatrix, petscSolVec, petscTmpVec);
    VecAXPY(petscTmpVec, -1.0, petscRhsVec);
    VecNorm(petscTmpVec, NORM_2, &norm);
    MSG("  Residual norm: %e\n", norm);


    // === Destroy Petsc's variables. ===

    MatDestroy(petscMatrix);
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
    KSPDestroy(solver);
  }

704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719

  void PetscSolver::providePetscSolver(KSP &solver, PC &pc)
  {
    FUNCNAME("PetscSolver::providePetscSolver()");

    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPGetPC(solver, &pc);

    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
    KSPSetFromOptions(solver);
    PCSetFromOptions(pc);
  }

720
}