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 28.5 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
    petsc::options_view(viewer);
99
    PetscViewerDestroy(&viewer);
100

101 102
    }
#endif
103 104 105
  }


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

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

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

118
    createMatVec(*seqMat);
119

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

125 126
    // === Create PETSc vector (solution and a temporary vector). ===

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

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

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

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

145
    matAssembly();
146

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


159
    // === Init PETSc solver and preconditioner objects. ===
160

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

165

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


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


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

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

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

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

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

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

	if (!dofMat)
221 222
	  continue;

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

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

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

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

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

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

245 246
	  cols.clear();
	  colsOther.clear();
247
	  values.clear();
248 249 250
	  valuesOther.clear();

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

327
    matAssembly();
328

329 330
    // === Create solver for the non primal (thus local) variables. ===

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


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

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

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

367
    vecRhsAssembly();
368 369 370 371 372 373 374 375 376 377 378 379

    // === 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
380 381 382
  }


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

399
      VecSet(getVecSolInterior(), 0.0);
400

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

404
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
405 406
    }

407

Thomas Witkowski's avatar
Thomas Witkowski committed
408 409
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
410
    if (nullspace.size() > 0 ||
411 412
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
413
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
414

415 416 417 418 419 420 421 422
      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);
423
      }
424

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

Thomas Witkowski's avatar
Thomas Witkowski committed
429 430
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
431 432

	VecNormalize(nullspaceBasis, PETSC_NULL);
433 434

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449
      // === Remove null space, if requested. ===
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
452
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
453 454

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

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    // PETSc.
463
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
464

465 466

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
467
      MatNullSpaceDestroy(&matNullspace);
468 469 470 471
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
472 473
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
474
    VecGetArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
475

476
    int c = 0;
477 478
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
479

480
      DofMap& d = (*interiorMap)[component].getMap();
481 482 483
      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
484 485
    }

486
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
487

488 489
    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
    // === more than one partition.                                       ===
Thomas Witkowski's avatar
Thomas Witkowski committed
490
    meshDistributor->synchVector(vec);
491 492 493
  }


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

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

517
    for (int i = 0; i < interiorMap->getRankDofs(); i++)
518 519 520 521 522 523 524 525 526
      rhsValues[i] = tmpValues[i];

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

    VecDestroy(&tmp);
  }


527 528
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
529
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
530

531 532
    exitPreconditioner(pcInterior);
    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
533 534 535
  }


536 537
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
538
    vecDestroy();
539 540 541
  }


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

561
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
562 563 564 565

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

570
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
571 572 573 574
    }
  }


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


586 587 588 589 590 591 592
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
593

594 595 596 597 598 599 600 601 602
  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);
  }
603

604

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

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


621
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
622 623
  {
    KSPDestroy(&ksp);
624
  }
625 626


627 628 629 630 631 632
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

633

634
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
635
  { }
636 637


638
  /// called by \ref fillPetscMatrix
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
639
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
640
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
641 642 643
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
644
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
645 646

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
647
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
648 649
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
650 651
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
652 653 654 655 656 657 658 659

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

661 662
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
663
    std::map<int, int> associatedRows;
Thomas Witkowski's avatar
Thomas Witkowski committed
664
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
665

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
666 667
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
668

669
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()),
670
         cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
671
      MultiIndex rowMultiIndex;
672
      if ((*interiorMap)[rowComp].find(cursor.value(), rowMultiIndex) == false)
673
        continue;
674 675

      int globalRowDof = rowMultiIndex.global;
676

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

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

684 685
      int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
      int assRowIndex = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
686

687 688 689 690
      if (periodicRow) {
        std::set<int> perAsc;
        perMap.fillAssociations(rowFe, globalRowDof,
                          meshDistributor->getElementObjectDb(), perAsc);
691

692 693
        vector<int> newRows;
        perMap.mapDof(rowFe, globalRowDof, perAsc, newRows);
Thomas Witkowski's avatar
Thomas Witkowski committed
694

695 696 697
        assRowIndex = interiorMap->getMatIndex(rowComp,
                                *std::max_element(newRows.begin(),newRows.end()));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
698

699 700
      cols.clear();
      values.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
701

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

705 706 707 708
        // 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
709

710 711 712
        int globalColDof = colMultiIndex.global;
        // Get PETSc's mat col index.
        int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
713

714 715 716
        cols.push_back(colIndex);
        values.push_back(value(*icursor));
          }
Thomas Witkowski's avatar
Thomas Witkowski committed
717

718 719 720 721 722 723
          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
724

725 726 727
        associatedRows.insert(std::make_pair(rowIndex, assRowIndex));
          }
        }
Thomas Witkowski's avatar
Thomas Witkowski committed
728

729 730
        int globalPeriodic = associatedRows.size();
        Parallel::mpi::globalAdd(globalPeriodic);
Thomas Witkowski's avatar
Thomas Witkowski committed
731

732
        if (globalPeriodic) {
Thomas Witkowski's avatar
Thomas Witkowski committed
733

734 735
          MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
          MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
736

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

739 740 741 742
          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
743

744
        MatSetValues(getMatInterior(), 1, &c[0], 2, c, v, INSERT_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
745
      }
746 747 748

      MatAssemblyBegin(getMatInterior(), MAT_FLUSH_ASSEMBLY);
      MatAssemblyEnd(getMatInterior(), MAT_FLUSH_ASSEMBLY);