PetscSolverGlobalMatrix.cc 26.8 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
7
8
9
10
11
//
// 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.

12
13
#include <mpi.h>
#include "DirichletBC.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 {

20
21
22
23
24
25
26
27
28
29
30

  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
  {
    Parameters::get("parallel->use zero start vector", zeroStartVector);
    Parameters::get("parallel->print matrix info", printMatInfo);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
31
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

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

41
    createMatVec(*seqMat);
42

Thomas Witkowski's avatar
Thomas Witkowski committed
43
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
44
      fillPetscMatrixWithCoarseSpace(seqMat);
45
46
47
      return;
    }

48
49
    // === Create PETSc vector (solution and a temporary vector). ===

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
56
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
59
60
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
63
64
65

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

66
    matAssembly();
67

68
69
    if (printMatInfo) {
      MatInfo matInfo;
70
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
71
72
73
74
75
76
77
      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));
    }
78
79


80
    // === Init PETSc solver and preconditioner objects. ===
81

82
83
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
84
    initPreconditioner(pcInterior);
85

86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
88
89
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
90
91


92
93
94
95
96
97
98
99
100
101
102
    // === For debugging allow to write the matrix to a file. ===

    bool dbgWriteMatrix = false;
    Parameters::get("parallel->debug->write matrix", dbgWriteMatrix);
    if (dbgWriteMatrix) {
      PetscViewer matView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.mat",
			    FILE_MODE_WRITE, &matView);
      MatView(getMatInterior(), matView);
      PetscViewerDestroy(&matView);
    }
103
104
105
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
106
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
107
108
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
109
110

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

114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    // === 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
133
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
137
138
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
139
140
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
143

Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
	std::set<DegreeOfFreedom> &dirichletRows = dofMat->getDirichletRows();

Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
148
149
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
156
157
158
159

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
	  if (dirichletRows.count(*cursor)) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor)))
160
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
161
	  }
162
163
164
165
166
167
168
169
170
171
  
	  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
172
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
173

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
174
	    if (isColCoarse == false)
175
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
176
177
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
181
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
184
185
186
187
188
189
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
190
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
195
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
196
			 1, &rowIndex, cols.size(),
197
198
199
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
202
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
203
		  rStartInterior;	      
204

Thomas Witkowski's avatar
Thomas Witkowski committed
205
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
206
207
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
208
209
	    }
	  } else {
210
	    if ((*interiorMap)[rowComponent].isSet(*cursor) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
211
212
	      continue;

213
	    int localRowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
	      (subdomainLevel == 0 ? 
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218
	    for (unsigned int i = 0; i < cols.size(); i++) {
219
	      if (subdomainLevel == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
220
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
221
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
222
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
223
224
	    }
	    
225
226
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
227
	    
228
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
231
	      int globalRowIndex = 
		interiorMap->getMatIndex(rowComponent, *cursor) + rStartInterior;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
232
233
234
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
  	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
237
238
	       		   1, &globalRowIndex, colsOther.size(),
  	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
239
240
241
242
243
244
	    }
	  }
	} 
      }
    }

245
    matAssembly();
246

247

248
249
250
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
251
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
252
		    SAME_NONZERO_PATTERN);
253
254
255
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
256
257
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
258
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
259
260
261
262
263
264
265
    } else {
      PCSetType(pcInterior, PCLU);
      if (subdomainLevel == 0)
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
266
    KSPSetFromOptions(kspInterior);
267
268
269
  }


270
271
272
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
273

274
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
275
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    
Thomas Witkowski's avatar
Thomas Witkowski committed
277
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
279
      for (int i = 0; i < vec->getSize(); i++)
280
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
281
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
    } else {
      for (int i = 0; i < vec->getSize(); i++)
284
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
    }

287
    vecRhsAssembly();
288
289
290
291
292
293
294
295
296
297
298
299

    // === For debugging allow to write the rhs vector to a file. ===

    bool dbgWriteRhs = false;
    Parameters::get("parallel->debug->write rhs", dbgWriteRhs);
    if (dbgWriteRhs) {
      PetscViewer vecView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.vec",
			    FILE_MODE_WRITE, &vecView);
      VecView(getVecRhsInterior(), vecView);
      PetscViewerDestroy(&vecView);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302
  }


303
304
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
309
310
311
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
312
313
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

314
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
      
      for (int i = 0; i < nComponents; i++)
317
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
318

319
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
    }

322

Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
325
326
327
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
328
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
329

330
331
332
333
334
335
336
337
338
339
      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
340
      if (nullspace.size() > 0) {
341
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
342
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
346
347

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
350
351
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

352
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	PetscReal n;
354
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
357
358
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364
      // === Remove null space, if requested. ===
365

Thomas Witkowski's avatar
Thomas Witkowski committed
366
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
367
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
368
369

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
370
	
371
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
375
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
376
377
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    // PETSc.
379
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
380

381
382

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
383
      MatNullSpaceDestroy(&matNullspace);
384
385
386
387
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
390
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
391

392
    int c = 0;
393
394
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
395

396
      DofMap& d = (*interiorMap)[component].getMap();
397
398
399
      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
400
401
    }

402
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
405
406

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


410
411
412
413
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

414
415
416
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

417
418
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
419
      interiorMap->createLocalVec(tmp);
420
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
421
      interiorMap->createVec(tmp);
422
423
424
425
426
427
428
429
430
431
432

    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);

433
434
435
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
436
    KSPSolve(kspInterior, tmp, tmp);
437
438
439
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
440
441
442
443
444
445
446
447
448
449
450

    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);
451
452
    t0 += MPI::Wtime() - wtime;

453
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
454
455
456
  }


457
458
459
460
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

461
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
462

463
464
465
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
468
  }


469
470
471
472
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

473
    vecDestroy();
474
475
476
  }


477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
  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());

496
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
497
498
499
500
501
502
503
504

      // 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());
      }

505
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
506
507
508
509
    }
  }


510
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
511
						 const char* splitName, 
512
513
514
515
516
517
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
518
    PCFieldSplitSetIS(pc, splitName, is);
519
520
521
522
    ISDestroy(&is);
  }


523
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

    KSPCreate(mpiCommGlobal, &ksp);
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
    KSPSetFromOptions(ksp);

    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }  


  void PetscSolverGlobalMatrix::exitSolver(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
  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
564
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
565
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
566
567
568
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
569
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
570
571

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
572
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
575
576
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
577
578
579
580
581
582
583
584
585
586
587

    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;

588
589
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
592
593
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
594
595

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
598
599
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
600
      // Global index of the current row DOF.
601
      MultiIndex rowMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
602
      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
603
604
605
	continue;

      int globalRowDof = rowMultiIndex.global;
606

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

Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
612
613
      // Dirichlet rows can be set only be the owner ranks.
      if (dirichletRows.count(*cursor) && !((*interiorMap)[rowComp].isRankDof(*cursor)))
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
614
615
616
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

617
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
618
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
619
620
621
622
623
624
625
626

	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.
627
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
628
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
629
630
631
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
632
	  // Test if the current col dof is a periodic dof.
633
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
634
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
635
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
636
637
638
639
640
641
642

	  // 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
643
644
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
647
648
649
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
655
656
657
658
659
660
661
	    // 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;
662
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
663
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
664
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
668
669
	      values.push_back(scaledValue);	      
	    }
	  }
	}

670
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
671
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
674
675
676
677
678
679
680
681
682
683
684
685
      } 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.
Thomas Witkowski's avatar
Thomas Witkowski committed
686
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
689
690
691
692
693
694
695

	  // 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;
696
697
698
699
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
702
703
704
705
706
707
708
709

	  // 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
710
711
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
714

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

715
	  for (unsigned int i = 0; i < entry.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
	    int rowIdx = interiorMap->getMatIndex(rowComp, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(colComp, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
718

719
720
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
723
724
725
726
727
728
729
730
731
732
	  }
	}


	// === 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;
733

734
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
737
738
739
740
741
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
744
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
745
					     int rowComp, 
746
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

750
    const FiniteElemSpace *feSpace = vec->getFeSpace();
751
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
752

753
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
756
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

    map<DegreeOfFreedom, double> &dirichletValues = vec->getDirichletValues();
757

Thomas Witkowski's avatar
Thomas Witkowski committed
758
759
760
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
763
764
765
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
      
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
766

767
768
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
769
      // Dirichlet rows can be set only be the owner ranks.
770
771
772
773
774
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
775

776
      if (isCoarseDof) {
777
778
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
779
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
780
781
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
782
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
783
784
	  continue;

785
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
786
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
787
788
789
790
791
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
793
794
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
795
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
796

797
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
798
799
800
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
801

Thomas Witkowski's avatar
Thomas Witkowski committed
802
803
804
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
805
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
806

Thomas Witkowski's avatar
Thomas Witkowski committed
807
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
808
809
810
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
813
814
815
816
      }
    }
  }

817
818
819
820
821
822
823
824
825
826

  PetscSolver* PetscSolverGlobalMatrix::createSubSolver(int component,
							string kspPrefix)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createSubSolver()");

    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
827
    subSolver->setKspPrefix(kspPrefix);
828
    subSolver->setMeshDistributor(meshDistributor);
829
830
831
832
833
834
835
836
837
    subSolver->init(fe, fe);

    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();

    return subSolver;
  }

838
839
840
841

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    Vec nullSpaceBasis;
    VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
    
    SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
    basisVec.set(0.0);
    basisVec.getDOFVector(constFeSpace)->set(1.0);
    
    setDofVector(nullSpaceBasis, basisVec, true);
    VecAssemblyBegin(nullSpaceBasis);
    VecAssemblyEnd(nullSpaceBasis);
    VecNormalize(nullSpaceBasis, PETSC_NULL);
    
    if (test) {
      Vec tmp;
      MatGetVecs(getMatInterior(), &tmp, PETSC_NULL);
      MatMult(getMatInterior(), nullSpaceBasis, tmp);
      PetscReal n;
      VecNorm(tmp, NORM_2, &n);
      MSG("NORM IS: %e\n", n);
      VecDestroy(&tmp);
865
    }
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
866
867
868
869
870
871
872
873
    
    
    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }

874

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
875
876
877
878
879
880
881
882
883
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
884

Thomas Witkowski's avatar
Thomas Witkowski committed
885
}