PetscSolverGlobalMatrix.cc 28.7 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
20

21
#include <mpi.h>
22
// #include "DirichletBC.h"
23
#include "DOFVector.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
#include "parallel/PetscSolverGlobalMatrix.h"
25
#include "parallel/PetscHelper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
27
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
28
#include "solver/PetscTypes.h"
29
#include "utility/PetscWrapper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
30

31
32
using namespace std;

33
namespace AMDiS { namespace Parallel {
34

35
  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name, bool setOptions)
36
37
38
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
39
  { FUNCNAME_DBG("PetscSolverGlobalMatrix()");
40

41
    bool matSolverPackage = false;
42
43
    if (setOptions) {
      PetscParameters params;
44

45
46
47
      // set the solver
      std::string solverName = "petsc";
      Parameters::get(name, solverName);
48
      if (solverName == "petsc")
49
	Parameters::get(name + "->ksp_type", solverName);
50

51
      std::string kspSolver = params.solverMap[solverName];
52

53
54
      if (params.matSolverPackage.find(kspSolver) != params.matSolverPackage.end()) {
	// direct solvers
55
56
57
	petsc::options_insert_string(("-" + kspPrefix + "ksp_type preonly").c_str());
	petsc::options_insert_string(("-" + kspPrefix + "pc_type lu").c_str());
	petsc::options_insert_string(("-" + kspPrefix + "pc_factor_mat_solver_package " + kspSolver).c_str());
58
59
60
	setMaxIterations(1);
	zeroStartVector = true;
	matSolverPackage = true;
61
      } else if (params.emptyParam.find(kspSolver) == params.emptyParam.end() && solverName != "petsc") {
62
	// other solvers
63
	petsc::options_insert_string(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
64
      }
65

66
67
68
69
70
      // set the preconditioner
      string precon = "";
      Parameters::get(name + "->pc_type", precon);
      if (!precon.size())
	Parameters::get(name + "->left precon", precon);
71
72
      if (!precon.size())
	Parameters::get(name + "->right precon", precon);
73
74
      if (!matSolverPackage && params.emptyParam.find(precon) == params.emptyParam.end()) {
	precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
75
	petsc::options_insert_string(("-" + kspPrefix + "pc_type " + precon).c_str());
76
      }
77

78
79
80
      petsc::options_insert_string(("-" + kspPrefix + "ksp_max_it " + boost::lexical_cast<std::string>(getMaxIterations())).c_str());
      petsc::options_insert_string(("-" + kspPrefix + "ksp_rtol " + boost::lexical_cast<std::string>(getRelative())).c_str());
      petsc::options_insert_string(("-" + kspPrefix + "ksp_atol " + boost::lexical_cast<std::string>(getTolerance())).c_str());
81

82
      if (getInfo() >= 20)
83
	petsc::options_insert_string(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
84
      else if (getInfo() >= 10)
85
	petsc::options_insert_string(("-" + kspPrefix + "ksp_monitor").c_str());
86
87
    }
    if (!matSolverPackage) {
88
      Parameters::get(name + "->use zero start vector", zeroStartVector);
89
    }
90
    Parameters::get("parallel->print matrix info", printMatInfo);
91

92
#ifndef NDEBUG
93
94
95
96
97
98
99
    bool printOptionsInfo = false;
    Parameters::get("parallel->debug->print options info", printOptionsInfo);
    if (printOptionsInfo) {
    MSG("PetscOptionsView:\n");
    PetscViewer viewer;
    PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
    PetscViewerSetType(viewer, PETSCVIEWERASCII);
100
    petsc::options_view(viewer);
101
    PetscViewerDestroy(&viewer);
102

103
104
    }
#endif
105
106
107
  }


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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
112
113
114
    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");
115

116
#ifndef NDEBUG
117
    Timer t;
118
#endif
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
119

120
    createMatVec(*seqMat);
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
123
      fillPetscMatrixWithCoarseSpace(seqMat);
124
125
      return;
    }
126

127
128
    // === Create PETSc vector (solution and a temporary vector). ===

129
#ifndef NDEBUG
130
131
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
#endif

134
    // === Transfer values from DOF matrices to the PETSc matrix. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
135

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
136
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
139
140
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
141

142
#ifndef NDEBUG
143
144
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
145
146
#endif

147
    matAssembly();
148

149
150
    if (printMatInfo) {
      MatInfo matInfo;
151
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
152
153
154
155
156
157
158
      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));
    }
159
160


161
    // === Init PETSc solver and preconditioner objects. ===
162

163
164
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
Praetorius, Simon's avatar
Praetorius, Simon committed
165
    initPreconditioner(*seqMat, mat[0][0]);
166

167

168
#ifndef NDEBUG
169
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", t.elapsed());
Thomas Witkowski's avatar
Thomas Witkowski committed
170
#endif
171
172


173
174
175
176
177
178
179
180
181
182
183
    // === 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);
    }
184
185
186
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
187
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
188
  {
189
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
190

191
    TEST_EXIT_DBG(interiorMap)("No interiorMap! Should not happen!\n");
192
    TEST_EXIT_DBG(coarseSpaceMap.size() == (size_t)seqMat->getSize())
193
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
194

195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
    // === 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);

211
    bool localMatrix =
212
213
      (meshDistributor->getMeshLevelData().getMpiComm(meshLevel) == MPI::COMM_SELF);

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
217
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
220
221
222
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
223
224
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
225
226
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
227

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

Thomas Witkowski's avatar
Thomas Witkowski committed
230
231
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
232

233
	// Traverse all rows.
234
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()),
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
236

237
	  bool isRowCoarse = isCoarseSpace(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
238

239
	  // For the case, this is a dirichlet row we have to check whether the
Thomas Witkowski's avatar
Thomas Witkowski committed
240
	  // rank is also owner of this row DOF.
241
242
243
	  if (dirichletRows.count(cursor.value())) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(cursor.value())) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(cursor.value())))
244
	      continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	  }
246

247
248
	  cols.clear();
	  colsOther.clear();
249
	  values.clear();
250
251
252
	  valuesOther.clear();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
256
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
257

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
258
	    if (isColCoarse == false)
259
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
260
261
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
264
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
265
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
268
269
270
271
272
273
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
274
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
277

278
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
279
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
280
281
	    		 1, &rowIndex, cols.size(),
	    		 &(cols[0]), &(values[0]), ADD_VALUES);
282
283

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
284
	      for (unsigned int i = 0; i < colsOther.size(); i++)
285
286
287
		colsOther[i] =
		  interiorMap->getMatIndex(colComponent, colsOther[i]) +
		  rStartInterior;
288

289
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent),
290
291
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
292
293
	    }
	  } else {
294
	    if ((*interiorMap)[rowComponent].isSet(cursor.value()) == false)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
295
296
	      continue;

297
298
	    int localRowIndex =
	      (localMatrix ?
299
300
	       interiorMap->getLocalMatIndex(rowComponent, cursor.value()) :
	       interiorMap->getMatIndex(rowComponent, cursor.value()));
301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
	    for (unsigned int i = 0; i < cols.size(); i++) {
303
	      if (localMatrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
304
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
305
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
306
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
307
	    }
308

309
310
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
311

312
	    if (colsOther.size()) {
313
	      int globalRowIndex =
314
		interiorMap->getMatIndex(rowComponent, cursor.value()) + rStartInterior;
315

Thomas Witkowski's avatar
Thomas Witkowski committed
316
	      for (unsigned int i = 0; i < colsOther.size(); i++)
317
		colsOther[i] =
Thomas Witkowski's avatar
Thomas Witkowski committed
318
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
319

320
	      MatSetValues(getMatInteriorCoarseByComponent(colComponent),
321
322
			   1, &globalRowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
323
324
	    }
	  }
325
	}
326
327
328
      }
    }

329
    matAssembly();
330

331
332
    // === Create solver for the non primal (thus local) variables. ===

333
    KSPCreate(domainComm, &kspInterior);
334
    petsc::ksp_set_operators(kspInterior, getMatInterior(), getMatInterior());
335
336
337
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
338
339
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
340
      petsc::pc_factor_set_mat_solver_package(pcInterior, MATSOLVERMUMPS);
341
342
    } else {
      PCSetType(pcInterior, PCLU);
343
      if (localMatrix)
344
	petsc::pc_factor_set_mat_solver_package(pcInterior, MATSOLVERUMFPACK);
345
      else
346
	petsc::pc_factor_set_mat_solver_package(pcInterior, MATSOLVERMUMPS);
347
    }
348
    KSPSetFromOptions(kspInterior);
349
350
351
  }


352
353
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
354
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
355

356
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
357
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
358
359

    // === Transfer values from DOF vector to the PETSc vector. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
361
      for (int i = 0; i < vec->getSize(); i++)
362
	setDofVector(getVecRhsInterior(),
Thomas Witkowski's avatar
Thomas Witkowski committed
363
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
    } else {
      for (int i = 0; i < vec->getSize(); i++)
366
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
    }

369
    vecRhsAssembly();
370
371
372
373
374
375
376
377
378
379
380
381

    // === 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
382
383
384
  }


385
386
387
388
389
  /// 1.) set startsolution
  /// 2.) create null-space
  /// 3.) solve Ax=b
  /// 4.) destroy null-space
  /// 5.) transfer solution back to DOFVector
390
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec,
391
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
397
398
  {
    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
399
400
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

401
      VecSet(getVecSolInterior(), 0.0);
402

Thomas Witkowski's avatar
Thomas Witkowski committed
403
      for (int i = 0; i < nComponents; i++)
404
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
405

406
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
407
408
    }

409

Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
Reuther, Sebastian's avatar
Reuther, Sebastian committed
412
    SystemVector *basisVec = NULL;
413
    if (nullspace.size() > 0 ||
414
415
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
416
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
417

418
419
      if (constNullspaceComponent.size() > 0) {
	nullspace.clear();
420
        basisVec = new SystemVector(vec);
421
422
423
424
425
	basisVec->set(0.0);
	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);

	nullspace.push_back(basisVec);
426
      }
427

Thomas Witkowski's avatar
Thomas Witkowski committed
428
      if (nullspace.size() > 0) {
429
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
430
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
434
435

	VecNormalize(nullspaceBasis, PETSC_NULL);
436
437

	MatNullSpaceCreate(domainComm, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE),
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
			   1, &nullspaceBasis, &matNullspace);

440
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	PetscReal n;
442
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
	MSG("NORM IS: %e\n", n);
      } else {
445
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
446
      }
447

Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
450
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
451

Thomas Witkowski's avatar
Thomas Witkowski committed
452
      // === Remove null space, if requested. ===
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
455
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
456
457

	MSG("Remove nullspace from rhs vector.\n");
458
	petsc::mat_nullspace_remove(matNullspace, getVecRhsInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
461
462
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
463
464
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
465
    // PETSc.
466
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
467

468
    // delete PETSc stuff
469
    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
      MatNullSpaceDestroy(&matNullspace);
471
472
      VecDestroy(&nullspaceBasis);
    }
473

474
    // delete allocated memory
Reuther, Sebastian's avatar
Reuther, Sebastian committed
475
    delete basisVec;
476

Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
479
    VecGetArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
480

481
    int c = 0;
482
483
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
484

485
      DofMap& d = (*interiorMap)[component].getMap();
486
487
488
      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
489
490
    }

491
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
492

493
494
    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
    // === more than one partition.                                       ===
Thomas Witkowski's avatar
Thomas Witkowski committed
495
    meshDistributor->synchVector(vec);
496
497
498
  }


499
500
501
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    Vec tmp;
502
    if (domainComm.Get_size() == 1)
503
      createLocalVec(*interiorMap, tmp);
504
    else
505
      createVec(*interiorMap, tmp);
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521

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

    KSPSolve(kspInterior, tmp, tmp);

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

522
    for (int i = 0; i < interiorMap->getRankDofs(); i++)
523
524
525
526
527
528
529
530
531
      rhsValues[i] = tmpValues[i];

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

    VecDestroy(&tmp);
  }


532
533
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
534
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
535

536
537
    exitPreconditioner(pcInterior);
    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
540
  }


541
542
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
543
    vecDestroy();
544
545
546
  }


547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
  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;
562
      Parameters::get("parallel->solver->is block " + boost::lexical_cast<string>(i),
563
564
565
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

566
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
567
568
569
570

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

575
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
576
577
578
579
    }
  }


580
581
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc,
						 const char* splitName,
582
583
584
585
						 vector<int> &components)
  {
    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
586
    PCFieldSplitSetIS(pc, splitName, is);
587
588
    ISDestroy(&is);
  }
589
590


591
592
593
594
595
596
597
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
598

599
600
601
602
603
  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
  {
    IS isrow, iscol;
    interiorMap->createIndexSet(isrow, startRow, numberOfRows);
    interiorMap->createIndexSet(iscol, startCol, numberOfCols);
604
    petsc::mat_get_sub_matrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
605
606
607
    ISDestroy(&iscol);
    ISDestroy(&isrow);
  }
608

609

610
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
611
  {
612
    KSPCreate(domainComm, &ksp);
613
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
614
615
616
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
617
    MatSetOptionsPrefix(getMatInterior(), kspPrefix.c_str());
618
619
620
621
622
    KSPSetFromOptions(ksp);

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


626
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
627
628
  {
    KSPDestroy(&ksp);
629
  }
630
631


632
633
634
635
636
637
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

638

639
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
640
  { }
641
642


643
  /// called by \ref fillPetscMatrix
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
644
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
645
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
646
647
648
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
649
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
652
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
655
656
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
660
661
662
663
664

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

666
667
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
668
    std::map<int, int> associatedRows;
Thomas Witkowski's avatar
Thomas Witkowski committed
669
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
670

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
671
672
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
673

674
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()),
675
         cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
676
      MultiIndex rowMultiIndex;
677
      if ((*interiorMap)[rowComp].find(cursor.value(), rowMultiIndex) == false)
678
        continue;
679
680

      int globalRowDof = rowMultiIndex.global;
681

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

Thomas Witkowski's avatar
Thomas Witkowski committed
685
      // Dirichlet rows can be set only be the owner ranks.
686
      if (dirichletRows.count(cursor.value()) && !((*interiorMap)[rowComp].isRankDof(cursor.value())))
687
        continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
688

689
690
      int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
      int assRowIndex = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
691

692
693
694
695
      if (periodicRow) {
        std::set<int> perAsc;
        perMap.fillAssociations(rowFe, globalRowDof,
                          meshDistributor->getElementObjectDb(), perAsc);
696

697
698
        vector<int> newRows;
        perMap.mapDof(rowFe, globalRowDof, perAsc, newRows);
Thomas Witkowski's avatar
Thomas Witkowski committed
699

700
701
702
        assRowIndex = interiorMap->getMatIndex(rowComp,
                                *std::max_element(newRows.begin(),newRows.end()));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
703

704
705
      cols.clear();
      values.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
706

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

710
711
712
713
        // Global index of the current column index.
        MultiIndex colMultiIndex;
        if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
          continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
714

715
716
717
        int globalColDof = colMultiIndex.global;
        // Get PETSc's mat col index.
        int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
718

719
720
721
        cols.push_back(colIndex);
        values.push_back(value(*icursor));
          }
Thomas Witkowski's avatar
Thomas Witkowski committed
722

723
724
725
726
727
728
          if (!periodicRow || assRowIndex == rowIndex)
        MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(),
                &(cols[0]), &(values[0]), ADD_VALUES);
          else {
        MatSetValues(getMatInterior(), 1, &assRowIndex, cols.size(),
                &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
729

730
731
732
        associatedRows.insert(std::make_pair(rowIndex, assRowIndex));
          }
        }
Thomas Witkowski's avatar
Thomas Witkowski committed
733

734
735
        int globalPeriodic = associatedRows.size();
        Parallel::mpi::globalAdd(globalPeriodic);
Thomas Witkowski's avatar
Thomas Witkowski committed
736

737
        if (globalPeriodic) {
Thomas Witkowski's avatar
Thomas Witkowski committed
738

739
740
          MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
          MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
741

742
          int c[2]; double v[2] = {1.0, -1.0};
Thomas Witkowski's avatar
Thomas Witkowski committed
743

744
745
746
747
          std::map<int, int>::iterator mapIt = associatedRows.begin();
          for (;mapIt != associatedRows.end(); mapIt++) {
        c[0] = mapIt->first;
        c[1] = mapIt->second;
Thomas Witkowski's avatar
Thomas Witkowski committed
748

749
        MatSetValues(getMatInterior(), 1, &c[0], 2, c, v, INSERT_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
750
      }
751
752
753

      MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
      MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
756
757
    }
  }


758
  /// called by \ref fillPetscRhs
759
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
760
					     Vec vecCoarse,
761
762
					     DOFVector<double>* vec,
					     int rowComp,
763
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
764
  {
765
    FUNCNAME_DBG("PetscSolverGlobalMatrix::setDofVector()");
Thomas Witkowski's avatar
Thomas Witkowski committed
766

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

770
771
    std::set<int> associatedRows;

772
    ParallelDofMapping *rowCoarseSpace =
773
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
774
775

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
781
      DegreeOfFreedom dof = dofIt.getDOFIndex();
782

Thomas Witkowski's avatar
Thomas Witkowski committed
783
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
784
        continue;
785

786
787
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
788
      // Dirichlet rows can be set only be the owner ranks.
789
      if (dirichletValues.count(dof)) {
790
791
792
        if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
            (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
          continue;
793
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
794

795
      if (isCoarseDof) {
796
        TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("vecCoarse not set! Should not happen!\n");
797

798
799
        int index = rowCoarseSpace->getMatIndex(rowComp, dof);
        VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
800
      } else {
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
        if ((*interiorMap)[rowComp].isSet(dof) == false)
          continue;

        // Calculate global row index of the DOF.
        DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;

        // Get PETSc's mat index of the row DOF.
        int index = 0, mappedIndex = -1;
        if (interiorMap->isMatIndexFromGlobal())
          index = interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
        else
          index =
            interiorMap->getMatIndex(rowComp, dof) + rStartInterior;

        bool periodicRow = perMap.isPeriodic(feSpace, globalRowDof);
        if (periodicRow) {
          std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
          vector<int> newRows;
          perMap.mapDof(feSpace, globalRowDof, perAsc, newRows);
          mappedIndex = interiorMap->getMatIndex(rowComp,
                              *std::max_element(newRows.begin(),newRows.end()));
        }

        if (!periodicRow || mappedIndex == index)
          VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
        else {
          VecSetValue(vecInterior, mappedIndex, *dofIt, ADD_VALUES);
          associatedRows.insert(index);
        }
      }
    }
Praetorius, Simon's avatar