PetscSolverGlobalMatrix.cc 29.2 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 25 26
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
27
#include "solver/PetscTypes.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28

29 30
using namespace std;

31
namespace AMDiS { namespace Parallel {
32

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

39
    bool matSolverPackage = false;
40 41
    if (setOptions) {
      PetscParameters params;
42

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

49
      std::string kspSolver = params.solverMap[solverName];
50

51 52
      if (params.matSolverPackage.find(kspSolver) != params.matSolverPackage.end()) {
	// direct solvers
53 54 55
	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());
56 57 58
	setMaxIterations(1);
	zeroStartVector = true;
	matSolverPackage = true;
59
      } else if (params.emptyParam.find(kspSolver) == params.emptyParam.end() && solverName != "petsc") {
60
	// other solvers
61
	petsc_options_insert_string(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
62
      }
63

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

76 77 78
      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());
79

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

90
#ifndef NDEBUG
91 92 93 94 95 96 97
    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);
98 99 100
#if (PETSC_VERSION_MINOR >= 7)
    PetscOptionsView(PETSC_NULL, viewer);
#else
101
    PetscOptionsView(viewer);
102
#endif
103
    PetscViewerDestroy(&viewer);
104

105 106
    }
#endif
107 108 109
  }


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

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

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

122
    createMatVec(*seqMat);
123

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

129 130
    // === Create PETSc vector (solution and a temporary vector). ===

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

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

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

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

149
    matAssembly();
150

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


163
    // === Init PETSc solver and preconditioner objects. ===
164

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

169

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


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


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

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

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

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

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

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

	if (!dofMat)
225 226
	  continue;

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

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

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

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

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

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

249 250
	  cols.clear();
	  colsOther.clear();
251
	  values.clear();
252 253 254
	  valuesOther.clear();

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

331
    matAssembly();
332

333 334
    // === Create solver for the non primal (thus local) variables. ===

335
    KSPCreate(domainComm, &kspInterior);
336 337 338 339 340
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior());
#else
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
#endif
341 342 343
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
344 345
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
346
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
347 348
    } else {
      PCSetType(pcInterior, PCLU);
349
      if (localMatrix)
350 351 352 353
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
354
    KSPSetFromOptions(kspInterior);
355 356 357
  }


358 359
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
360
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
361

362
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
363
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
364 365

    // === Transfer values from DOF vector to the PETSc vector. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
366
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
367
      for (int i = 0; i < vec->getSize(); i++)
368
	setDofVector(getVecRhsInterior(),
Thomas Witkowski's avatar
Thomas Witkowski committed
369
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
370 371
    } else {
      for (int i = 0; i < vec->getSize(); i++)
372
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
373 374
    }

375
    vecRhsAssembly();
376 377 378 379 380 381 382 383 384 385 386 387

    // === 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
388 389 390
  }


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

407
      VecSet(getVecSolInterior(), 0.0);
408

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

412
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
413 414
    }

415

Thomas Witkowski's avatar
Thomas Witkowski committed
416 417
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
418
    if (nullspace.size() > 0 ||
419 420
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
421
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
422

423 424 425 426 427 428 429 430
      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);
431
      }
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
      if (nullspace.size() > 0) {
434
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
435
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437 438
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
439 440

	VecNormalize(nullspaceBasis, PETSC_NULL);
441 442

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

445
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
446
	PetscReal n;
447
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
448 449
	MSG("NORM IS: %e\n", n);
      } else {
450
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
451
      }
452

Thomas Witkowski's avatar
Thomas Witkowski committed
453 454 455
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
      // === Remove null space, if requested. ===
458

Thomas Witkowski's avatar
Thomas Witkowski committed
459
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
460
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
461 462

	MSG("Remove nullspace from rhs vector.\n");
463 464 465
#if (PETSC_VERSION_MINOR >= 5)
	MatNullSpaceRemove(matNullspace, getVecRhsInterior());
#else
466
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
467
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
468 469 470 471
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
472 473
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    // PETSc.
475
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
476

477 478

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
479
      MatNullSpaceDestroy(&matNullspace);
480 481 482 483
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
484 485
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
486
    VecGetArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
487

488
    int c = 0;
489 490
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
491

492
      DofMap& d = (*interiorMap)[component].getMap();
493 494 495
      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
496 497
    }

498
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
499

500 501
    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
    // === more than one partition.                                       ===
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    meshDistributor->synchVector(vec);
503 504 505
  }


506 507 508
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    Vec tmp;
509
    if (domainComm.Get_size() == 1)
510
      createLocalVec(*interiorMap, tmp);
511
    else
512
      createVec(*interiorMap, tmp);
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528

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

529
    for (int i = 0; i < interiorMap->getRankDofs(); i++)
530 531 532 533 534 535 536 537 538
      rhsValues[i] = tmpValues[i];

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

    VecDestroy(&tmp);
  }


539 540
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
541
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
542

543 544
    exitPreconditioner(pcInterior);
    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
545 546 547
  }


548 549
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
550
    vecDestroy();
551 552 553
  }


554 555 556 557 558 559 560 561 562 563 564 565 566 567 568
  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;
569
      Parameters::get("parallel->solver->is block " + boost::lexical_cast<string>(i),
570 571 572
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

573
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
574 575 576 577

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

582
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
583 584 585 586
    }
  }


587 588
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc,
						 const char* splitName,
589 590 591 592
						 vector<int> &components)
  {
    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
593
    PCFieldSplitSetIS(pc, splitName, is);
594 595
    ISDestroy(&is);
  }
596 597


598 599 600 601 602 603 604
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
605

606 607 608 609 610 611 612 613 614
  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);
    MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
    ISDestroy(&iscol);
    ISDestroy(&isrow);
  }
615

616

617
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
618
  {
619
    KSPCreate(domainComm, &ksp);
620 621 622 623 624
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
#endif
625 626 627
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
628
    MatSetOptionsPrefix(getMatInterior(), kspPrefix.c_str());
629 630 631 632 633
    KSPSetFromOptions(ksp);

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


637
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
638 639
  {
    KSPDestroy(&ksp);
640
  }
641 642


643 644 645 646 647 648
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

649

650
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
651
  { }
652 653


654
  /// called by \ref fillPetscMatrix
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
655
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
656
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
657 658 659
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
660
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
661 662

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
663
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
664 665
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
666 667
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
668 669 670 671 672 673 674 675

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

677 678
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
679
    std::map<int, int> associatedRows;
Thomas Witkowski's avatar
Thomas Witkowski committed
680
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
681

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
682 683
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
684

685
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()),
686
         cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
687
      MultiIndex rowMultiIndex;
688
      if ((*interiorMap)[rowComp].find(cursor.value(), rowMultiIndex) == false)
689
        continue;
690 691

      int globalRowDof = rowMultiIndex.global;
692

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

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

700 701
      int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
      int assRowIndex = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
702

703 704 705 706
      if (periodicRow) {
        std::set<int> perAsc;
        perMap.fillAssociations(rowFe, globalRowDof,
                          meshDistributor->getElementObjectDb(), perAsc);
707

708 709
        vector<int> newRows;
        perMap.mapDof(rowFe, globalRowDof, perAsc, newRows);
Thomas Witkowski's avatar
Thomas Witkowski committed
710

711 712 713
        assRowIndex = interiorMap->getMatIndex(rowComp,
                                *std::max_element(newRows.begin(),newRows.end()));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
714

715 716
      cols.clear();
      values.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
717

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

721 722 723 724
        // 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
725

726 727 728
        int globalColDof = colMultiIndex.global;
        // Get PETSc's mat col index.
        int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
729

730 731 732
        cols.push_back(colIndex);
        values.push_back(value(*icursor));
          }
Thomas Witkowski's avatar
Thomas Witkowski committed
733

734 735 736 737 738 739
          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
740

741 742 743
        associatedRows.insert(std::make_pair(rowIndex, assRowIndex));
          }
        }
Thomas Witkowski's avatar
Thomas Witkowski committed
744

745 746
        int globalPeriodic = associatedRows.size();
        Parallel::mpi::globalAdd(globalPeriodic);
Thomas Witkowski's avatar
Thomas Witkowski committed
747

748
        if (globalPeriodic) {
Thomas Witkowski's avatar
Thomas Witkowski committed
749

750 751
          MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
          MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
752

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

755 756 757 758
          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
759

760
        MatSetValues(getMatInterior(), 1, &c[0], 2, c, v, INSERT_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
761
      }
762 763 764

      MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
      MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
765 766 767 768
    }
  }


769
  /// called by \ref fillPetscRhs
770
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
771
					     Vec vecCoarse,
772 773
					     DOFVector<double>* vec,
					     int rowComp,
774
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
775
  {
776
    FUNCNAME_DBG("PetscSolverGlobalMatrix::setDofVector()");
Thomas Witkowski's avatar
Thomas Witkowski committed
777

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

781 782
    std::set<int> associatedRows;

783
    ParallelDofMapping *rowCoarseSpace =
784
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
785 786

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
792
      DegreeOfFreedom dof = dofIt.getDOFIndex();
793

Thomas Witkowski's avatar
Thomas Witkowski committed
794
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
795
        continue;
796

797 798
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
799
      // Dirichlet rows can be set only be the owner ranks.</