PetscSolverGlobalMatrix.cc 25.9 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
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.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
15
16
17
18
19
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
20
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
21
22
23
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

24
25
26
    petscData.create(interiorMap, coarseSpaceMap, 
		     subdomainLevel, mpiCommLocal, mpiCommGlobal);

Thomas Witkowski's avatar
Thomas Witkowski committed
27
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
28
      updateSubdomainData();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
29
      fillPetscMatrixWithCoarseSpace(seqMat);
30
31
32
      return;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
33
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
34
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
35
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");
36
    
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    double wtime = MPI::Wtime();
38

39
    // === If required, recompute non zero structure of the matrix. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
40

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
41
42
    if (checkMeshChange(seqMat))
      nnzInterior.create(seqMat, mpiCommGlobal, *interiorMap, 
43
44
			 &(meshDistributor->getPeriodicMap()),
			 meshDistributor->getElementObjectDb());
Thomas Witkowski's avatar
Thomas Witkowski committed
45

46
47
    // === Create PETSc vector (solution and a temporary vector). ===

48
49
    int nRankRows = interiorMap->getRankDofs();
    int nOverallRows = interiorMap->getOverallDofs();
50

51
    VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &petscSolVec);
52
53


Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
56
57
58
59
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

#if (DEBUG != 0)
    int a, b;
60
    MatGetOwnershipRange(petscData.getInteriorMat(), &a, &b);
61
62
    TEST_EXIT(a == interiorMap->getStartDofs())("Wrong matrix ownership range!\n");
    TEST_EXIT(b == interiorMap->getStartDofs() + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
63
64
65
66
67
68
      ("Wrong matrix ownership range!\n");
#endif


    // === Transfer values from DOF matrices to the PETSc matrix. === 

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
69
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
72
73
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
74
75
76
77
78

#if (DEBUG != 0)
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

79
    petscData.assembly();
80
81
82
    
    if (printMatInfo) {
      MatInfo matInfo;
83
      MatGetInfo(petscData.getInteriorMat(), MAT_GLOBAL_SUM, &matInfo);
84
85
86
87
88
89
90
      MSG("Matrix info:\n");
      MSG("  memory usage: %e MB\n", matInfo.memory / (1024.0 * 1024.0));
      MSG("  mallocs: %d\n", static_cast<int>(matInfo.mallocs));
      MSG("  nz allocated: %d\n", static_cast<int>(matInfo.nz_allocated));
      MSG("  nz used: %d\n", static_cast<int>(matInfo.nz_used));
      MSG("  nz unneeded: %d\n", static_cast<int>(matInfo.nz_unneeded));
    }
91
92
93

    // === Remove Dirichlet BC DOFs. ===

94
    //    removeDirichletBcDofs(mat);
95
96
    

97
    // === Init PETSc solver. ===
98

99
100
    KSPCreate(mpiCommGlobal, &kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
101
102
103
104
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(), 
		    petscData.getInteriorMat(), 
		    SAME_NONZERO_PATTERN); 
105
106
107
108
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(kspInterior, KSPBCGS);
    KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
    KSPSetFromOptions(kspInterior);
109

110
    initPreconditioner(pcInterior);
111

112
113
    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
114
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
118
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
119
120
121
  }


122
123
124
125
126
127
128
129
  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
  {
    Matrix<DOFMatrix*> m(1, 1);
    m[0][0] = mat;
    fillPetscMatrix(&m);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
130
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
131
132
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
133
134

    TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
135
136
    TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize())
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
137

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
138
139
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(seqMat);

140
141
142
    int nRowsRankInterior = interiorMap->getRankDofs();
    int nRowsOverallInterior = interiorMap->getOverallDofs();

143
144
145
    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (subdomainLevel == 0);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
146
147
    if (checkMeshChange(seqMat, localMatrix)) {
      nnzInterior.create(seqMat, mpiCommGlobal, *interiorMap, NULL, 
148
149
			 meshDistributor->getElementObjectDb(),
			 localMatrix);
150

Thomas Witkowski's avatar
Thomas Witkowski committed
151
      if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
152
153
	MSG("NO COARSE SPACE NNZ!\n");
	/*
154
155
156
	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap, NULL, meshDistributor->getElementObjectDb());
	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
Thomas Witkowski's avatar
Thomas Witkowski committed
157
	*/
158
      }
159
160
    }

161

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    // === Prepare traverse of sequentially created matrices. ===

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

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

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
181
    int nComponents = seqMat->getSize();
182
183
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
184
	if (!(*seqMat)[i][j])
185
186
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[i];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[j];

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
190
191
	traits::col<Matrix>::type col((*seqMat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*seqMat)[i][j]->getBaseMatrix());
192
193
	
	// Traverse all rows.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
194
195
	for (cursor_type cursor = begin<row>((*seqMat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*seqMat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
	  bool isRowCoarse = isCoarseSpace(i, feSpaces[i], *cursor);
198
199
200
201
202
203
204
205
206
207
  
	  cols.clear();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

Thomas Witkowski's avatar
Thomas Witkowski committed
208
	    bool isColCoarse = isCoarseSpace(j, feSpaces[j], col(*icursor));
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
211
	    if (isColCoarse) {
	      if (isRowCoarse) {
212
213
214
215
216
217
218
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      } else {
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
219
	      if (isRowCoarse) {
220
221
222
223
224
225
226
227
228
229
230
231
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      }
	    }
	  }  // for each nnz in row


	  // === Set matrix values. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
232
233
	  if (isRowCoarse) {
	    int rowIndex = rowCoarseSpace->getMatIndex(i, *cursor);
234
	    for (unsigned int k = 0; k < cols.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	      cols[k] = colCoarseSpace->getMatIndex(j, cols[k]);
236

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
237
	    // matCoarseCoarse
238
	    MatSetValues(petscData.getCoarseMatComp(i),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
239
			 1, &rowIndex, cols.size(),
240
241
242
243
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      if (subdomainLevel == 0) {
244
		for (unsigned int k = 0; k < colsOther.size(); k++)
245
246
247
248
249
250
		  colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
	      } else {
		for (unsigned int k = 0; k < colsOther.size(); k++)
		  colsOther[k] = 
		    interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior;
	      }
251

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
252
	      // matCoarseInt
253
	      MatSetValues(petscData.getCoarseIntMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
254
			   1, &rowIndex, colsOther.size(),
255
256
257
258
259
260
261
262
263
264
265
266
267
268
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
	    int localRowIndex = 
	      (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(i, *cursor) :
	       interiorMap->getMatIndex(i, *cursor));

	    for (unsigned int k = 0; k < cols.size(); k++) {
	      if (subdomainLevel == 0)
		cols[k] = interiorMap->getLocalMatIndex(j, cols[k]);
	      else
		cols[k] = interiorMap->getMatIndex(j, cols[k]);
	    }
	    
269
  	    MatSetValues(petscData.getInteriorMat(), 1, &localRowIndex, cols.size(),
270
271
272
273
274
275
276
277
278
  			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      int globalRowIndex = interiorMap->getMatIndex(i, *cursor);

	      if (subdomainLevel != 0)
		globalRowIndex += rStartInterior;

	      for (unsigned int k = 0; k < colsOther.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
279
		colsOther[k] = colCoarseSpace->getMatIndex(j, colsOther[k]);
280

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
281
	      // matIntCoarse
282
  	      MatSetValues(petscData.getIntCoarseMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
283
			   1, &globalRowIndex, colsOther.size(),
284
285
286
287
288
289
290
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

291
    petscData.assembly();
292

293
294
295
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
296
297
298
299
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(),
		    petscData.getInteriorMat(),
		    SAME_NONZERO_PATTERN);
300
301
302
303
304
305
306
307
308
309
310
311
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCLU);
    if (subdomainLevel == 0)
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
    else
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    KSPSetFromOptions(kspInterior);  
  }


312
313
314
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
315

Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
319
    VecCreateMPI(mpiCommGlobal, 
		 interiorMap->getRankDofs(), 
		 nGlobalOverallInterior,
		 &rhsInterior);
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
    if (coarseSpaceMap.size()) 
Thomas Witkowski's avatar
Thomas Witkowski committed
322
      VecCreateMPI(mpiCommGlobal, 
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
		   coarseSpaceMap[0]->getRankDofs(), 
		   coarseSpaceMap[0]->getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
		   &rhsCoarseSpace);
    
327
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
328
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
329
    
330

Thomas Witkowski's avatar
Thomas Witkowski committed
331
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, rhsCoarseSpace, vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
338
    } else {
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, vec->getDOFVector(i), i);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
339

340
341
    VecAssemblyBegin(rhsInterior);
    VecAssemblyEnd(rhsInterior);
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
347
      VecAssemblyBegin(rhsCoarseSpace);
      VecAssemblyEnd(rhsCoarseSpace);
    }

348
349
350

    // === Remove Dirichlet BC DOFs. ===
    //    removeDirichletBcDofs(vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
 
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
  }


355
356
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
359
360
361
362
363
364
365
366
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
      VecSet(petscSolVec, 0.0);
      
      for (int i = 0; i < nComponents; i++)
367
	setDofVector(petscSolVec, vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
370
371
372
      
      VecAssemblyBegin(petscSolVec);
      VecAssemblyEnd(petscSolVec);
    }

373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
376
377
378
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
379
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
380

381
382
383
384
385
386
387
388
389
390
      if (constNullspaceComponent.size() > 0) {
	nullspace.clear();
	SystemVector *basisVec = new SystemVector(vec);
	basisVec->set(0.0);
	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);

	nullspace.push_back(basisVec);
      } 

Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
393
394
395
396
397
      if (nullspace.size() > 0) {
	VecDuplicate(petscSolVec, &nullspaceBasis);
	for (int i = 0; i < nComponents; i++)
	  setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
398
399

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

404
	MatMult(petscData.getInteriorMat(), nullspaceBasis, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
410
	PetscReal n;
	VecNorm(petscSolVec, NORM_2, &n);
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
411

412
      MatSetNullSpace(petscData.getInteriorMat(), matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
413
      KSPSetNullSpace(kspInterior, matNullspace);
414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
      // === Remove null space, if requested. ===
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
418
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
421
422
423
424
	
	MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
425
426
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    // PETSc.
428
    solve(rhsInterior, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
429

430
431

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
432
      MatNullSpaceDestroy(&matNullspace);
433
434
435
436
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
439
    VecGetArray(petscSolVec, &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
440

441
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    for (int i = 0; i < nComponents; i++) {
443
      DOFVector<double> &dv = *(vec.getDOFVector(i));
444

445
      DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap();
446
447
448
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dv[it->first] = vecPointer[c++];
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
451
    }

    VecRestoreArray(petscSolVec, &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
454
455

    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===  
    // === more than one partition.                                       ===  
    meshDistributor->synchVector(vec);
456
457
458
  }


459
460
461
462
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

463
464
465
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
      VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp);
    else
      VecCreateMPI(mpiCommLocal,
		   interiorMap->getRankDofs(),
		   interiorMap->getOverallDofs(),
		   &tmp);

    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
    VecGetArray(rhs, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++)
      tmpValues[i] = rhsValues[i];

    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

485
486
487
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
488
    KSPSolve(kspInterior, tmp, tmp);
489
490
491
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
492
493
494
495
496
497
498
499
500
501
502

    VecGetArray(tmp, &tmpValues);
    VecGetArray(sol, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];

    VecRestoreArray(sol, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    VecDestroy(&tmp);
503
504
    t0 += MPI::Wtime() - wtime;

505
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
506
507
508
  }


509
510
511
512
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

513
514
    exitPreconditioner(pcInterior);

515
    petscData.destroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
516

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
517
    KSPDestroy(&kspInterior);
518
519

    if (petscSolVec != PETSC_NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
      VecDestroy(&petscSolVec);
521
522
      petscSolVec = PETSC_NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
523
524
525
  }


526
527
528
529
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

530
    VecDestroy(&rhsInterior);
531

Thomas Witkowski's avatar
Thomas Witkowski committed
532
    if (coarseSpaceMap.size())
533
      VecDestroy(&rhsCoarseSpace);
534
535
536
  }


537
538
539
540
541
542
543
544
545
546
  bool PetscSolverGlobalMatrix::checkMeshChange(Matrix<DOFMatrix*> *mat,
						bool localMatrix)
  {
    FUNCNAME("PetscSolverGlobalMatrix::checkMeshChange()");

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

547
    if (!nnzInterior.dnnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
548
549
      vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
      
550
      interiorMap->setComputeMatIndex(!localMatrix);
551
552
      interiorMap->update(feSpaces);

553
      nnzInterior.clear();
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
      if (coarseSpaceMap.size()) {
556
557
558
	nnzCoarse.clear();
	nnzCoarseInt.clear();
	nnzIntCoarse.clear();
559
560
561
562
563
564
565
566
567
568
569
570
      }

      updateSubdomainData();
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();

      return true;
    }

    return false;
  }


571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    vector<string> isNames;
    Parameters::get("parallel->solver->is blocks", isNames);

    int nBlocks = isNames.size();
    if (nBlocks == 0)
      return;

    for (int i = 0; i < nBlocks; i++) {
      MSG("Create for block %s\n", isNames[i].c_str());

      vector<int> blockComponents;
      Parameters::get("parallel->solver->is block " + lexical_cast<string>(i),
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

      TEST_EXIT(nComponents > 0)("No is block for block %d defined!\n", i);

      // Check if blocks are continous
      for (int j = 0; j < nComponents; j++) {
	TEST_EXIT(blockComponents[j] == blockComponents[0] + j)
	  ("Does not yet support not continous IS blocks! Block %s\n", 
	   isNames[i].c_str());
      }

      IS is;
      interiorMap->createIndexSet(is, blockComponents[0], nComponents);
      PCFieldSplitSetIS(pc, isNames[i].c_str(), is);
      ISDestroy(&is);
    }
  }


607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
622
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
623
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
624
625
626
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
627
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
628
629

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
630
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
631
632
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
633
634
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
637
638
639
640
641
642
643
644
645

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

    vector<int> cols;
    vector<double> values;
    cols.reserve(300);
    values.reserve(300);
    
    vector<int> globalCols;

646
647
648
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
649
650
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
653
654
    DofMap& rowMap = (*interiorMap)[rowFe].getMap();
    DofMap& colMap = (*interiorMap)[colFe].getMap();

    // === Traverse all rows of the DOF matrix and insert row wise the values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
657
658
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
659
      // Global index of the current row DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
660
      int globalRowDof = rowMap[*cursor].global;
661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
      // Test if the current row DOF is a periodic DOF.
663
      bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
664

Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

668
	// Get PETSc's mat row index.
669
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
670
671
672
673
674
675
676
677

	cols.clear();
	values.clear();

	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {

	  // Global index of the current column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
678
	  int globalColDof = colMap[col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
	  // Test if the current col dof is a periodic dof.
680
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
681
	  // Get PETSc's mat col index.
682
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
685
686
687
688
689

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;

	  if (!periodicCol) {
	    // Calculate the exact position of the column index in the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
694
695
696
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
697
698
699
	    perMap.fillAssociations(colFe, globalColDof, 
				    meshDistributor->getElementObjectDb(), perAsc);

Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
702
703
704
705
706
707
708
	    // Scale value to the number of periodic associations of the column index.
	    double scaledValue = 
	      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;
709
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
710
	    for (unsigned int i = 0; i < newCols.size(); i++) {
711
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
714
715
716
	      values.push_back(scaledValue);	      
	    }
	  }
	}

717
  	MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
718
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
721
722
723
724
725
726
727
728
729
730
731
732
      } else {
	// === Row DOF index is periodic. ===

	// 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.
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
	  // Global index of the current column index.
733
	  int globalColDof = (*interiorMap)[colFe][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
736
737
738
739
740
741
742

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;

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

	  std::set<int> perAsc;
743
744
745
746
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750
751
752
753
754
755
756

	  // Scale the value with respect to the number of periodic associations.
	  double scaledValue = 
	    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;
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
757
758
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
761

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

762
	  for (unsigned int i = 0; i < entry.size(); i++) {
763
764
	    int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
765

766
767
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
770
771
772
773
774
775
776
777
778
779
	  }
	}


	// === Finally, add all periodic rows to the PETSc matrix. ===

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

	  int rowIndex = rowIt->first;
780

781
	  MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
784
785
786
787
788
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
791
					     DOFVector<double>* vec, 
792
					     int nRowVec, 
793
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
796
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

797
    const FiniteElemSpace *feSpace = vec->getFeSpace();
798
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
801
802
    
    ParallelDofMapping *rowCoarseSpace = NULL;
    if (coarseSpaceMap.size())
      rowCoarseSpace = coarseSpaceMap[nRowVec];
803

Thomas Witkowski's avatar
Thomas Witkowski committed
804
805
806
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
807

808
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
811
      if (isCoarseSpace(nRowVec, feSpace, dofIt.getDOFIndex())) {
812
813
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
814
	int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
815
816
817
818
819
820
821
822
823
824
825
826
827
828
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
	// Calculate global row index of the DOF.
	DegreeOfFreedom globalRowDof = 
	  (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
	    interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior;
	else
	  index =
	    interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
829

830
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
833
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
834

Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
837
838
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
	    int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
839

Thomas Witkowski's avatar
Thomas Witkowski committed
840
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
841
842
843
844
	  }	  
	} else {	  
	  // The DOF index is not periodic.
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
845
846
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
847
848
849
850
851
      }
    }
  }


852
853
  void PetscSolverGlobalMatrix::removeDirichletBcDofs(Matrix<DOFMatrix*> *mat)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
854
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");    
855
856
857
858
859
860
861
862
  }


  void PetscSolverGlobalMatrix::removeDirichletBcDofs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
863
}