PetscSolverGlobalMatrix.cc 25.9 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14 15 16 17 18 19
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
20
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
21 22 23
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

24 25 26
    petscData.create(interiorMap, coarseSpaceMap, 
		     subdomainLevel, mpiCommLocal, mpiCommGlobal);

Thomas Witkowski's avatar
Thomas Witkowski committed
27
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
28
      updateSubdomainData();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
29
      fillPetscMatrixWithCoarseSpace(seqMat);
30 31 32
      return;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
33
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
34
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
35
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");
36
    
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    double wtime = MPI::Wtime();
38

39
    // === If required, recompute non zero structure of the matrix. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
40

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
41 42
    if (checkMeshChange(seqMat))
      nnzInterior.create(seqMat, mpiCommGlobal, *interiorMap, 
43 44
			 &(meshDistributor->getPeriodicMap()),
			 meshDistributor->getElementObjectDb());
Thomas Witkowski's avatar
Thomas Witkowski committed
45

46 47
    // === Create PETSc vector (solution and a temporary vector). ===

48 49
    int nRankRows = interiorMap->getRankDofs();
    int nOverallRows = interiorMap->getOverallDofs();
50

51
    VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &petscSolVec);
52 53


Thomas Witkowski's avatar
Thomas Witkowski committed
54 55 56 57 58 59
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

#if (DEBUG != 0)
    int a, b;
60
    MatGetOwnershipRange(petscData.getInteriorMat(), &a, &b);
61 62
    TEST_EXIT(a == interiorMap->getStartDofs())("Wrong matrix ownership range!\n");
    TEST_EXIT(b == interiorMap->getStartDofs() + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
63 64 65 66 67 68
      ("Wrong matrix ownership range!\n");
#endif


    // === Transfer values from DOF matrices to the PETSc matrix. === 

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
69
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
70 71
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
72 73
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
74 75 76 77 78

#if (DEBUG != 0)
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

79
    petscData.assembly();
80 81 82
    
    if (printMatInfo) {
      MatInfo matInfo;
83
      MatGetInfo(petscData.getInteriorMat(), MAT_GLOBAL_SUM, &matInfo);
84 85 86 87 88 89 90
      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));
    }
91 92 93

    // === Remove Dirichlet BC DOFs. ===

94
    //    removeDirichletBcDofs(mat);
95 96
    

97
    // === Init PETSc solver. ===
98

99 100
    KSPCreate(mpiCommGlobal, &kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
101 102 103 104
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(), 
		    petscData.getInteriorMat(), 
		    SAME_NONZERO_PATTERN); 
105 106 107 108
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(kspInterior, KSPBCGS);
    KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
    KSPSetFromOptions(kspInterior);
109

110
    initPreconditioner(pcInterior);
111

112 113
    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
114
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116 117 118
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
119 120 121
  }


122 123 124 125 126 127 128 129
  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
  {
    Matrix<DOFMatrix*> m(1, 1);
    m[0][0] = mat;
    fillPetscMatrix(&m);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
130
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
131 132
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
133 134

    TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
135 136
    TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize())
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
137

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
138 139
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(seqMat);

140 141 142
    int nRowsRankInterior = interiorMap->getRankDofs();
    int nRowsOverallInterior = interiorMap->getOverallDofs();

143 144 145
    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (subdomainLevel == 0);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
146 147
    if (checkMeshChange(seqMat, localMatrix)) {
      nnzInterior.create(seqMat, mpiCommGlobal, *interiorMap, NULL, 
148 149
			 meshDistributor->getElementObjectDb(),
			 localMatrix);
150

Thomas Witkowski's avatar
Thomas Witkowski committed
151
      if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
152 153
	MSG("NO COARSE SPACE NNZ!\n");
	/*
154 155 156
	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap, NULL, meshDistributor->getElementObjectDb());
	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
Thomas Witkowski's avatar
Thomas Witkowski committed
157
	*/
158
      }
159 160
    }

161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
    // === 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);

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
181
    int nComponents = seqMat->getSize();
182 183
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
184
	if (!(*seqMat)[i][j])
185 186
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
187 188 189
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[i];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[j];

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
190 191
	traits::col<Matrix>::type col((*seqMat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*seqMat)[i][j]->getBaseMatrix());
192 193
	
	// Traverse all rows.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
194 195
	for (cursor_type cursor = begin<row>((*seqMat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*seqMat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
	  bool isRowCoarse = isCoarseSpace(i, feSpaces[i], *cursor);
198 199 200 201 202 203 204 205 206 207
  
	  cols.clear();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
208
	    bool isColCoarse = isCoarseSpace(j, feSpaces[j], col(*icursor));
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210 211
	    if (isColCoarse) {
	      if (isRowCoarse) {
212 213 214 215 216 217 218
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      } else {
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
219
	      if (isRowCoarse) {
220 221 222 223 224 225 226 227 228 229 230 231
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      }
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
232 233
	  if (isRowCoarse) {
	    int rowIndex = rowCoarseSpace->getMatIndex(i, *cursor);
234
	    for (unsigned int k = 0; k < cols.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	      cols[k] = colCoarseSpace->getMatIndex(j, cols[k]);
236

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
237
	    // matCoarseCoarse
238
	    MatSetValues(petscData.getCoarseMatComp(i),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
239
			 1, &rowIndex, cols.size(),
240 241 242 243
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      if (subdomainLevel == 0) {
244
		for (unsigned int k = 0; k < colsOther.size(); k++)
245 246 247 248 249 250
		  colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
	      } else {
		for (unsigned int k = 0; k < colsOther.size(); k++)
		  colsOther[k] = 
		    interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior;
	      }
251

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
252
	      // matCoarseInt
253
	      MatSetValues(petscData.getCoarseIntMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
254
			   1, &rowIndex, colsOther.size(),
255 256 257 258 259 260 261 262 263 264 265 266 267 268
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
	    int localRowIndex = 
	      (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(i, *cursor) :
	       interiorMap->getMatIndex(i, *cursor));

	    for (unsigned int k = 0; k < cols.size(); k++) {
	      if (subdomainLevel == 0)
		cols[k] = interiorMap->getLocalMatIndex(j, cols[k]);
	      else
		cols[k] = interiorMap->getMatIndex(j, cols[k]);
	    }
	    
269
  	    MatSetValues(petscData.getInteriorMat(), 1, &localRowIndex, cols.size(),
270 271 272 273 274 275 276 277 278
  			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      int globalRowIndex = interiorMap->getMatIndex(i, *cursor);

	      if (subdomainLevel != 0)
		globalRowIndex += rStartInterior;

	      for (unsigned int k = 0; k < colsOther.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
279
		colsOther[k] = colCoarseSpace->getMatIndex(j, colsOther[k]);
280

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
281
	      // matIntCoarse
282
  	      MatSetValues(petscData.getIntCoarseMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
283
			   1, &globalRowIndex, colsOther.size(),
284 285 286 287 288 289 290
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

291
    petscData.assembly();
292

293 294 295
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
296 297 298 299
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(),
		    petscData.getInteriorMat(),
		    SAME_NONZERO_PATTERN);
300 301 302 303 304 305 306 307 308 309 310 311
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCLU);
    if (subdomainLevel == 0)
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
    else
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    KSPSetFromOptions(kspInterior);  
  }


312 313 314
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
315

Thomas Witkowski's avatar
Thomas Witkowski committed
316 317 318 319
    VecCreateMPI(mpiCommGlobal, 
		 interiorMap->getRankDofs(), 
		 nGlobalOverallInterior,
		 &rhsInterior);
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
    if (coarseSpaceMap.size()) 
Thomas Witkowski's avatar
Thomas Witkowski committed
322
      VecCreateMPI(mpiCommGlobal, 
Thomas Witkowski's avatar
Thomas Witkowski committed
323 324
		   coarseSpaceMap[0]->getRankDofs(), 
		   coarseSpaceMap[0]->getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
325 326
		   &rhsCoarseSpace);
    
327
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
328
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
329
    
330

Thomas Witkowski's avatar
Thomas Witkowski committed
331
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
333 334
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, rhsCoarseSpace, vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
335 336 337 338
    } else {
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, vec->getDOFVector(i), i);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
339

340 341
    VecAssemblyBegin(rhsInterior);
    VecAssemblyEnd(rhsInterior);
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
344 345 346 347
      VecAssemblyBegin(rhsCoarseSpace);
      VecAssemblyEnd(rhsCoarseSpace);
    }

348 349 350

    // === Remove Dirichlet BC DOFs. ===
    //    removeDirichletBcDofs(vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
 
Thomas Witkowski's avatar
Thomas Witkowski committed
352 353 354
  }


355 356
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
357 358 359 360 361 362 363 364 365 366
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
      VecSet(petscSolVec, 0.0);
      
      for (int i = 0; i < nComponents; i++)
367
	setDofVector(petscSolVec, vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
368 369 370 371 372
      
      VecAssemblyBegin(petscSolVec);
      VecAssemblyEnd(petscSolVec);
    }

373

Thomas Witkowski's avatar
Thomas Witkowski committed
374 375
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
376 377 378
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
379
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
380

381 382 383 384 385 386 387 388 389 390
      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);
      } 

Thomas Witkowski's avatar
Thomas Witkowski committed
391 392 393 394 395 396 397
      if (nullspace.size() > 0) {
	VecDuplicate(petscSolVec, &nullspaceBasis);
	for (int i = 0; i < nComponents; i++)
	  setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
398 399

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
400 401 402 403
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

404
	MatMult(petscData.getInteriorMat(), nullspaceBasis, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
405 406 407 408 409 410
	PetscReal n;
	VecNorm(petscSolVec, NORM_2, &n);
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
411

412
      MatSetNullSpace(petscData.getInteriorMat(), matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
413
      KSPSetNullSpace(kspInterior, matNullspace);
414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
      // === Remove null space, if requested. ===
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
418
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
419 420 421 422 423 424
	
	MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
425 426
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    // PETSc.
428
    solve(rhsInterior, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
429

430 431

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
432
      MatNullSpaceDestroy(&matNullspace);
433 434 435 436
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
437 438
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
439
    VecGetArray(petscSolVec, &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
440

441
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    for (int i = 0; i < nComponents; i++) {
443
      DOFVector<double> &dv = *(vec.getDOFVector(i));
444

445
      DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap();
446 447 448
      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
449 450 451
    }

    VecRestoreArray(petscSolVec, &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
452 453 454 455

    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===  
    // === more than one partition.                                       ===  
    meshDistributor->synchVector(vec);
456 457 458
  }


459 460 461 462
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

463 464 465
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
      VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp);
    else
      VecCreateMPI(mpiCommLocal,
		   interiorMap->getRankDofs(),
		   interiorMap->getOverallDofs(),
		   &tmp);

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

485 486 487
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
488
    KSPSolve(kspInterior, tmp, tmp);
489 490 491
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
492 493 494 495 496 497 498 499 500 501 502

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

    for (int i = 0; i < interiorMap->getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];

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

    VecDestroy(&tmp);
503 504
    t0 += MPI::Wtime() - wtime;

505
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
506 507 508
  }


509 510 511 512
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

513 514
    exitPreconditioner(pcInterior);

515
    petscData.destroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
516

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
517
    KSPDestroy(&kspInterior);
518 519

    if (petscSolVec != PETSC_NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
      VecDestroy(&petscSolVec);
521 522
      petscSolVec = PETSC_NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
523 524 525
  }


526 527 528 529
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

530
    VecDestroy(&rhsInterior);
531

Thomas Witkowski's avatar
Thomas Witkowski committed
532
    if (coarseSpaceMap.size())
533
      VecDestroy(&rhsCoarseSpace);
534 535 536
  }


537 538 539 540 541 542 543 544 545 546
  bool PetscSolverGlobalMatrix::checkMeshChange(Matrix<DOFMatrix*> *mat,
						bool localMatrix)
  {
    FUNCNAME("PetscSolverGlobalMatrix::checkMeshChange()");

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

547
    if (!nnzInterior.dnnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
548 549
      vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
      
550
      interiorMap->setComputeMatIndex(!localMatrix);
551 552
      interiorMap->update(feSpaces);

553
      nnzInterior.clear();
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
      if (coarseSpaceMap.size()) {
556 557 558
	nnzCoarse.clear();
	nnzCoarseInt.clear();
	nnzIntCoarse.clear();
559 560 561 562 563 564 565 566 567 568 569 570
      }

      updateSubdomainData();
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();

      return true;
    }

    return false;
  }


571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606
  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;
      Parameters::get("parallel->solver->is block " + lexical_cast<string>(i),
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

      TEST_EXIT(nComponents > 0)("No is block for block %d defined!\n", i);

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

      IS is;
      interiorMap->createIndexSet(is, blockComponents[0], nComponents);
      PCFieldSplitSetIS(pc, isNames[i].c_str(), is);
      ISDestroy(&is);
    }
  }


607 608 609 610 611 612 613 614 615 616 617 618 619 620 621
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
622
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
623
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
624 625 626
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
627
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
628 629

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
630
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
631 632
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
633 634
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
635 636 637 638 639 640 641 642 643 644 645

    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);
    
    vector<int> globalCols;

646 647 648
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
649 650
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
651 652 653 654
    DofMap& rowMap = (*interiorMap)[rowFe].getMap();
    DofMap& colMap = (*interiorMap)[colFe].getMap();

    // === Traverse all rows of the DOF matrix and insert row wise the values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
655 656
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
657 658
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
659
      // Global index of the current row DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
660
      int globalRowDof = rowMap[*cursor].global;
661

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

Thomas Witkowski's avatar
Thomas Witkowski committed
665 666 667
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

668
	// Get PETSc's mat row index.
669
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
670 671 672 673 674 675 676 677

	cols.clear();
	values.clear();

	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {

	  // Global index of the current column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
678
	  int globalColDof = colMap[col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
	  // Test if the current col dof is a periodic dof.
680
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
681
	  // Get PETSc's mat col index.
682
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
683 684 685 686 687 688 689

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;

	  if (!periodicCol) {
	    // Calculate the exact position of the column index in the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
690 691
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
692 693 694 695 696
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
697 698 699
	    perMap.fillAssociations(colFe, globalColDof, 
				    meshDistributor->getElementObjectDb(), perAsc);

Thomas Witkowski's avatar
Thomas Witkowski committed
700 701 702 703 704 705 706 707 708
	    // Scale value to the number of periodic associations of the column index.
	    double scaledValue = 
	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));

	    
	    // === Create set of all matrix column indices due to the periodic ===
	    // === associations of the column DOF index.                       ===

	    vector<int> newCols;
709
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
710
	    for (unsigned int i = 0; i < newCols.size(); i++) {
711
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
712 713 714 715 716
	      values.push_back(scaledValue);	      
	    }
	  }
	}

717
  	MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
718
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
719 720 721 722 723 724 725 726 727 728 729 730 731 732
      } else {
	// === Row DOF index is periodic. ===

	// Because this row is periodic, we will have to add the entries of this 
	// matrix row to multiple rows. The following maps store to each row an
	// array of column indices and values of the entries that must be added to
	// the PETSc matrix.
	map<int, vector<int> > colsMap;
	map<int, vector<double> > valsMap;

	// Traverse all column entries.
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
	  // Global index of the current column index.
733
	  int globalColDof = (*interiorMap)[colFe][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
734 735 736 737 738 739 740 741 742

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;

	  // === Add all periodic associations of both, the row and the column ===
	  // === indices to the set perAsc.                                    ===

	  std::set<int> perAsc;
743 744 745 746
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
747 748 749 750 751 752 753 754 755 756

	  // Scale the value with respect to the number of periodic associations.
	  double scaledValue = 
	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));


	  // === Create all matrix entries with respect to the periodic  ===
	  // === associations of the row and column indices.             ===

	  vector<pair<int, int> > entry;
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
757 758
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
759 760 761

	  // === Translate the matrix entries to PETSc's matrix.

762
	  for (unsigned int i = 0; i < entry.size(); i++) {
763 764
	    int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
765

766 767
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
768 769 770 771 772 773 774 775 776 777 778 779
	  }
	}


	// === Finally, add all periodic rows to the PETSc matrix. ===

	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
	     rowIt != colsMap.end(); ++rowIt) {
	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
	    ("Should not happen!\n");

	  int rowIndex = rowIt->first;
780

781
	  MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
782 783 784 785 786 787 788
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
789 790
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
791
					     DOFVector<double>* vec, 
792
					     int nRowVec, 
793
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
794 795 796
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

797
    const FiniteElemSpace *feSpace = vec->getFeSpace();
798
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
799 800 801 802
    
    ParallelDofMapping *rowCoarseSpace = NULL;
    if (coarseSpaceMap.size())
      rowCoarseSpace = coarseSpaceMap[nRowVec];
803

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

808
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
809 810
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
811
      if (isCoarseSpace(nRowVec, feSpace, dofIt.getDOFIndex())) {
812 813
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
814
	int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
815 816 817 818 819 820 821 822 823 824 825 826 827 828
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
	// Calculate global row index of the DOF.
	DegreeOfFreedom globalRowDof = 
	  (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
	    interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior;
	else
	  index =
	    interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
829

830
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
831 832 833
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
834

Thomas Witkowski's avatar
Thomas Witkowski committed
835 836 837 838
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
	    int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
839

Thomas Witkowski's avatar
Thomas Witkowski committed
840
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
841 842 843 844
	  }	  
	} else {	  
	  // The DOF index is not periodic.
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
845 846
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
847 848 849 850 851
      }
    }
  }


852 853
  void PetscSolverGlobalMatrix::removeDirichletBcDofs(Matrix<DOFMatrix*> *mat)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
854
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");    
855 856 857 858 859 860 861 862
  }


  void PetscSolverGlobalMatrix::removeDirichletBcDofs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
863
}