Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

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