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

30 31
using namespace std;

32
namespace AMDiS { namespace Parallel {
33

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

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

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

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

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

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

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

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

91
#ifndef NDEBUG
92 93 94 95 96 97 98
    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);
99
    petsc::options_view(viewer);
100
    PetscViewerDestroy(&viewer);
101

102 103
    }
#endif
104 105 106
  }


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

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

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

119
    createMatVec(*seqMat);
120

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

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

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

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

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

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

146
    matAssembly();
147

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


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

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

166

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


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


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

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

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

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

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

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

	if (!dofMat)
222 223
	  continue;

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

328
    matAssembly();
329

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

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


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

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

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

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

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


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

400
      VecSet(getVecSolInterior(), 0.0);
401

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

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

408

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

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

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

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

	VecNormalize(nullspaceBasis, PETSC_NULL);
434 435

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

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

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

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

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

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

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

466 467

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


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

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

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

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

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


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

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

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

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

    VecDestroy(&tmp);
  }


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

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


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


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

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

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

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


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


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

595 596 597 598 599 600 601 602 603
  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
  {
    IS isrow, iscol;
    interiorMap->createIndexSet(isrow, startRow, numberOfRows);
    interiorMap->createIndexSet(iscol, startCol, numberOfCols);
    MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
    ISDestroy(&iscol);
    ISDestroy(&isrow);
  }
604

605

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

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


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


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

634

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


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

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

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

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

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

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

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

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

      int globalRowDof = rowMultiIndex.global;
677

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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