ParallelDomainBase.cc 84.2 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
#include <algorithm>

3
#include "ParallelDomainBase.h"
4
5
6
7
8
9
10
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
11
12
#include "DOFMatrix.h"
#include "DOFVector.h"
13
#include "SystemVector.h"
14
#include "VtkWriter.h"
15
#include "ElementDofIterator.h"
16
17
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
18
#include "ElementFileWriter.h"
19
#include "VertexVector.h"
20
#include "StdMpi.h"
21
#include "MeshStructure.h"
22
#include "Debug.h"
23
24

#include "petscksp.h"
25
26
27

namespace AMDiS {

28
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
29
  {    
30
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
31
      std::cout << "  Petsc-Iteration " << iter << ": " << rnorm << std::endl;
32
33
34
35

    return 0;
  }

36
37
38
39
40
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

41
  ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF,
42
43
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
44
					 RefinementManager *refinementManager)
45
46
    : iterationIF(iIF),
      timeIF(tIF),
47
      name(iIF->getName()),
48
49
      feSpace(fe),
      mesh(fe->getMesh()),
50
      refineManager(refinementManager),
51
      initialPartitionMesh(true),
52
      nRankDofs(0),
53
      rstart(0),
54
      nComponents(1),
55
56
      deserialized(false),
      lastMeshChangeIndex(0)
57
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
61
62
63
64
    FUNCNAME("ParallelDomainBase::ParalleDomainBase()");

    TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
      ("Only meshes with one DOFAdmin are supported!\n");
    TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
      ("Wrong pre dof number for DOFAdmin!\n");

65
66
67
68
69
70
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

71
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
72
  {
73
74
75
76
77
78
79
    FUNCNAME("ParallelDomainBase::initParallelization()");

    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

    // If the problem has been already read from a file, we do not need to do anything.
    if (deserialized)
80
81
      return;

82
83
84
85
86
    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

87
88
    MSG("START PARTITIONING!\n");

89
90
91
92
93
94
95
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
    setElemWeights(adaptInfo);
    // and now partition the mesh
    partitionMesh(adaptInfo);   

Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
98
    dbgCreateElementMap(elMap);
99
100
101
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
102

103
104
105
106
107
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
108
#endif
109

110
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
111

112
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
113

114
115
116
117
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
118
119
120
121
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

122
123
124
125
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

126
127
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
128
#if (DEBUG != 0)
129
    MSG("AMDiS runs in debug mode, so make some test ...\n");
130
131
    dbgTestElementMap(elMap);
    dbgTestInteriorBoundary();
132
    dbgTestCommonDofs(true);
133
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
134

135
136
    debug::writeMesh(feSpace, -1, "macromesh");   
#endif
137

138
139
140
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
141

142
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
143

Thomas Witkowski's avatar
Thomas Witkowski committed
144
    int globalRefinement = 0;
145
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
Thomas Witkowski's avatar
Thomas Witkowski committed
146

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (globalRefinement > 0) {
148
      refineManager->globalRefine(mesh, globalRefinement);
149

150
      updateLocalGlobalNumbering();
151

152
      // === Update periodic mapping, if there are periodic boundaries. ===
153

154
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
155
    }
156
157
  }

158
159

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
160
  {}
161

162
163
164
  
  void ParallelDomainBase::updateDofAdmins()
  {
165
166
167
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
168
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
169
170
171

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
172
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
173
	admin.enlargeDOFLists();
174
175
176
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
177
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
178
179
 	admin.setDOFFree(j, false);

180
181
182
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
183
184
185
    }
  }

186

187
188
189
190
191
192
193
194
195
196
  void ParallelDomainBase::testForMacroMesh()
  {
    FUNCNAME("ParallelDomainBase::testForMacroMesh()");

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
197
	("Mesh is already refined! This does not work with parallelization!\n");
198
199
200
201
202
203
204
205
206
207
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

    TEST_EXIT(nMacroElements >= mpiSize)
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

208

209
210
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
211
  {
212
213
214
215
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");

216
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
217
218
219
220
221
222
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

223
    typedef traits::range_generator<row, Matrix>::type cursor_type;
224
225
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

226
227
228
229
230
    std::vector<int> cols;
    std::vector<double> values;
    cols.reserve(300);
    values.reserve(300);

231
232
233
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the petsc matrix.                                               ===

234
235
236
237
238
239
    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

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

240
      // Global index of the current row dof.
241
      int globalRowDof = mapLocalGlobalDofs[*cursor];
242
      // Test if the current row dof is a periodic dof.
243
      bool periodicRow = (periodicDof.count(globalRowDof) > 0);
244

245
246
247

      // === Traverse all non zero entries of the row and produce vector cols ===
      // === with the column indices of all row entries and vector values     ===
248
      // === with the corresponding values.                                   ===
249

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

	// Set only non null values.
254
	if (value(*icursor) != 0.0) {
255
	  // Global index of the current column index.
256
	  int globalColDof = mapLocalGlobalDofs[col(*icursor)];
257
	  // Calculate the exact position of the column index in the petsc matrix.
258
259
	  int colIndex = globalColDof * dispMult + dispAddCol;

260
261
	  // If the current row is not periodic, but the current dof index is periodic,
	  // we have to duplicate the value to the other corresponding periodic columns.
262
 	  if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
263
264
265
266
267
	    // The value is assign to n matrix entries, therefore, every entry 
	    // has only 1/n value of the original entry.
	    double scalFactor = 1.0 / (periodicDof[globalColDof].size() + 1.0);

	    // Insert original entry.
268
 	    cols.push_back(colIndex);
269
 	    values.push_back(value(*icursor) * scalFactor);
270

271
272
273
	    // Insert the periodic entries.
 	    for (std::set<DegreeOfFreedom>::iterator it = 
		   periodicDof[globalColDof].begin();
274
275
 		 it != periodicDof[globalColDof].end(); ++it) {
 	      cols.push_back(*it * dispMult + dispAddCol);
276
277
 	      values.push_back(value(*icursor) * scalFactor);
	    }
278
 	  } else {
279
	    // Neigher row nor column dof index is periodic, simple add entry.
280
281
282
	    cols.push_back(colIndex);
	    values.push_back(value(*icursor));
	  }
283
	}
284
285
      }

286
287
288
289
290

      // === Up to now we have assembled on row. Now, the row must be send to the ===
      // === corresponding rows to the petsc matrix.                              ===

      // Calculate petsc row index.
291
      int rowIndex = globalRowDof * dispMult + dispAddRow;
292
      
293
      if (periodicRow) {
294
295
296
	// The row dof is periodic, so send dof to all the corresponding rows.

	double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0);
297
	
298
	int diagIndex = -1;
299
	for (int i = 0; i < static_cast<int>(values.size()); i++) {
300
301
302
303
304
	  // Change only the non diagonal values in the col. For the diagonal test
	  // we compare the global dof indices of the dof matrix (not of the petsc
	  // matrix!).
	  if ((cols[i] - dispAddCol) / dispMult != globalRowDof)
	    values[i] *= scalFactor;
305
306
307
	  else
	    diagIndex = i;
	}
308
309
310
311
312
313
314
	
	// Send the main row to the petsc matrix.
	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
		     &(cols[0]), &(values[0]), ADD_VALUES);	
 
	// Set diagonal element to zero, i.e., the diagonal element of the current
	// row is not send to the periodic row indices.
315
316
317
	if (diagIndex != -1)
	  values[diagIndex] = 0.0;

318
	// Send the row to all periodic row indices.
319
320
321
	for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
	     it != periodicDof[globalRowDof].end(); ++it) {
	  int perRowIndex = *it * dispMult + dispAddRow;
322
323
	  MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(), 
		       &(cols[0]), &(values[0]), ADD_VALUES);
324
325
326
	}

      } else {
327
328
329
	// The row dof is not periodic, simply send the row to the petsc matrix.
	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
		     &(cols[0]), &(values[0]), ADD_VALUES);
330
      }    
331
    }
332
  }
333

334

335
  void ParallelDomainBase::setDofVector(Vec& petscVec, DOFVector<double>* vec, 
336
337
					int dispMult, int dispAdd)
  {
338
    // Traverse all used dofs in the dof vector.
339
340
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
341
      // Calculate global row index of the dof.
342
      int globalRow = mapLocalGlobalDofs[dofIt.getDOFIndex()];
343
      // Calculate petsc index of the row dof.
344
345
346
      int index = globalRow * dispMult + dispAdd;

      if (periodicDof.count(globalRow) > 0) {
347
348
349
	// The dof index is periodic, so devide the value to all dof entries.

	double value = *dofIt / (periodicDof[globalRow].size() + 1.0);
350
351
352
353
354
355
356
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);

	for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRow].begin();
	     it != periodicDof[globalRow].end(); ++it) {
	  index = *it * dispMult + dispAdd;
	  VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
	}
357

358
      } else {
359
	// The dof index is not periodic.
360
361
362
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
363
    }    
364
365
  }

366

367
368
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
369
370
371
372
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

373
374
    // === Create PETSc vector (rhs, solution and a temporary vector). ===

375
376
377
378
379
380
381
382
383
384
385
386
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    VecSetSizes(petscSolVec, nRankRows, nOverallRows);
    VecSetType(petscSolVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
    VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
    VecSetType(petscTmpVec, VECMPI);

387
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
388
    namespace traits = mtl::traits;
389
    typedef DOFMatrix::base_matrix_type Matrix;
390
    typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
391
392

    int d_nnz[nRankRows];
393
394
    int o_nnz[nRankRows];
    for (int i = 0; i < nRankRows; i++) {
395
      d_nnz[i] = 0;
396
397
      o_nnz[i] = 0;
    }
398

399
400
401
402
403
404
    // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
    // that this rank will send to. This nnz entries will be assembled on this rank,
    // but because the row DOFs are not DOFs of this rank they will be send to the
    // owner of the row DOFs.
    std::map<int, MatrixNnzEntry> sendMatrixEntry;

405
406
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
407
408
409
410
411
412
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
413
	  typedef traits::range_generator<row, Matrix>::type cursor_type;
414
415
	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	  
416
417
418
	  for (cursor_type cursor = begin<row>(bmat), 
		 cend = end<row>(bmat); cursor != cend; ++cursor) {

419
420
421
	    // Map the local row number to the global DOF index and create from it
	    // the global PETSc row index of this DOF.
	    int petscRowIdx = mapLocalGlobalDofs[*cursor] * nComponents + i;
422

423
	    if (isRankDof[*cursor]) {
424
425
426
427
428
429

	      // === The current row DOF is a rank dof, so create the corresponding ===
	      // === nnz values directly on rank's nnz data.                        ===

	      // This is the local row index of the local PETSc matrix.
	      int localPetscRowIdx = petscRowIdx - rstart * nComponents;
430
431

#if (DEBUG != 0)    
432
	      if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
433
		std::cout << "ERROR in rank: " << mpiRank << std::endl;
434
435
		std::cout << "  Wrong r = " << localPetscRowIdx << " " << *cursor 
			  << " " << mapLocalGlobalDofs[*cursor] << " " 
436
437
438
439
			  << nRankRows << std::endl;
		ERROR_EXIT("Should not happen!\n");
	      }
#endif
440
	      
441
	      // Traverse all non zero entries in this row.
442
443
444
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
445
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
446

447
448
449
450
451
		  // The row DOF is a rank DOF, if also the column is a rank DOF, 
		  // increment the d_nnz values for this row, otherwise the o_nnz value.
		  if (petscColIdx >= rstart * nComponents && 
		      petscColIdx < rstart * nComponents + nRankRows)
		    d_nnz[localPetscRowIdx]++;
452
		  else
453
		    o_nnz[localPetscRowIdx]++;
454
455
456
		}    
	      }
	    } else {
457
458
459
460
461
	      
	      // === The current row DOF is not a rank dof, i.e., it will be created ===
	      // === on this rank, but after this it will be send to another rank    ===
	      // === matrix. So we need to send also the corresponding nnz structure ===
	      // === of this row to the corresponding rank.                          ===
462

463
464
	      // Find out who is the member of this DOF.
	      int sendToRank = -1;
465
466
467
468
469
	      for (RankToDofContainer::iterator it = recvDofs.begin();
		   it != recvDofs.end(); ++it) {
		for (DofContainer::iterator dofIt = it->second.begin();
		     dofIt != it->second.end(); ++dofIt) {
		  if (**dofIt == *cursor) {
470
471
472
473
474
475

		    if (petscRowIdx == 6717) {
		      debug::writeDofMesh(mpiRank, *cursor, feSpace);
		      MSG("SEND DOF TO: %d/%d\n", it->first, *cursor);
		    }

476
477
478
479
480
481
482
483
484
485
486
		    sendToRank = it->first;
		    break;
		  }
		}

		if (sendToRank != -1)
		  break;
	      }

	      TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");

487
	      // Send all non zero entries to the member of the row DOF.
488
489
490
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
491
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
492
		  
493
494
		  sendMatrixEntry[sendToRank].
		    push_back(std::make_pair(petscRowIdx, petscColIdx));
495
496
497
		}
	      }

498
499
500
501
	    } // if (isRankDof[*cursor]) ... else ...
	  } // for each row in mat[i][j]
	} // if mat[i][j]
      } 
502
503
    }

504
    // === Send and recv the nnz row structure to/from other ranks. ===
505

506
    StdMpi<MatrixNnzEntry, MatrixNnzEntry> stdMpi(mpiComm, true);
507
508
509
    stdMpi.send(sendMatrixEntry);
    stdMpi.recv(sendDofs);
    stdMpi.startCommunication<int>(MPI_INT);
Thomas Witkowski's avatar
Thomas Witkowski committed
510

511
512
    // === Evaluate the nnz structure this rank got from other ranks and add it to ===
    // === the PETSc nnz data structure.                                           ===
513

514
515
516
517
518
519
    for (std::map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second.size() > 0) {
	for (unsigned int i = 0; i < it->second.size(); i++) {
	  int r = it->second[i].first;
	  int c = it->second[i].second;
520

521
	  int localRowIdx = r - rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
522

523
524
525
	  TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
	    ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
	     r, localRowIdx, nRankRows, it->first);
526
	  
527
	  if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
528
	    o_nnz[localRowIdx]++;
529
	  else
530
	    d_nnz[localRowIdx]++;
531
532
	}
      }
533
    }  
Thomas Witkowski's avatar
Thomas Witkowski committed
534

535
    // === Create PETSc matrix with the computed nnz data structure. ===
536
537

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
538
539
540
541
542
543
544
545
		    0, d_nnz, 0, o_nnz, &petscMatrix);

#if (DEBUG != 0)
    int a, b;
    MatGetOwnershipRange(petscMatrix, &a, &b);
    TEST_EXIT(a == rstart * nComponents)("Wrong matrix ownership range!\n");
    TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n");
#endif
546

547
    // === Transfer values from DOF matrices to the PETSc matrix. === 
548
549
550

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
551
	if ((*mat)[i][j])
552
553
554
555
556
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

557
558
    // === Transfer values from DOF vector to the PETSc vector. === 

559
    for (int i = 0; i < nComponents; i++)
560
      setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
561

562
563
564
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

565
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
566
567
568
  }


569
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
570
571
572
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

573
574
575
576
577
578
579
580
581
#if 0
    // Set old solution to be initiual guess for petsc solver.
    for (int i = 0; i < nComponents; i++)
      setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i);

    VecAssemblyBegin(petscSolVec);
    VecAssemblyEnd(petscSolVec);
#endif

582
    // === Init Petsc solver. ===
583

584
    KSP solver;
585
586
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
587
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
588
589
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
590
    KSPSetFromOptions(solver);
591
592
    // Do not delete the solution vector, use it for the initial guess.
    //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
593

594
595
596

    // === Run Petsc. ===

597
    KSPSolve(solver, petscRhsVec, petscSolVec);
598

599
    // === Transfere values from Petsc's solution vectors to the dof vectors.
600
601
602
603
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
604
      DOFVector<double> *dofvec = vec.getDOFVector(i);
605
      for (int j = 0; j < nRankDofs; j++)
606
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];      
607
608
609
610
    }

    VecRestoreArray(petscSolVec, &vecPointer);

611
612
613

    // === Synchronize dofs at common dofs, i.e., dofs that correspond to more ===
    // === than one partition.                                                 ===
614
    synchVectors(vec);
615

616
617
618

    // === Print information about solution process. ===

619
620
621
622
623
624
625
626
627
628
    int iterations = 0;
    KSPGetIterationNumber(solver, &iterations);
    MSG("  Number of iterations: %d\n", iterations);
    
    double norm = 0.0;
    MatMult(petscMatrix, petscSolVec, petscTmpVec);
    VecAXPY(petscTmpVec, -1.0, petscRhsVec);
    VecNorm(petscTmpVec, NORM_2, &norm);
    MSG("  Residual norm: %e\n", norm);

629
630
631

    // === Destroy Petsc's variables. ===

632
    MatDestroy(petscMatrix);
633
634
635
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
636
637
638
639
640
    KSPDestroy(solver);
  }
  
  void ParallelDomainBase::synchVectors(SystemVector &vec)
  {
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
#if 0
    StdMpi<std::vector<double>, std::vector<double> > stdMpi(mpiComm);

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
      std::vector<double> dofs;
      int nSendDOFs = sendIt->second.size();
      dofs.reserve(nComponents * sendIt->second.size());
      
      for (int j = 0; j < nComponents; j++) {
	DOFVector<double> *dofvec = vec.getDOFVector(j);
	for (int k = 0; k < nSendDOFs; k++)
	  dofs.push_back((*dofvec)[*((sendIt->second)[k])]);
      }

      stdMpi.send(sendIt->first, dofs);
    }

659
    stdMpi.startCommunication<int>();
660
661
662
663
#endif

#if 1

664
665
666
667
668
669
670
671
672
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
    
    int i = 0;
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
673
674
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
675
676
677

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
678
	DOFVector<double> *dofvec = vec.getDOFVector(j);
679
680
681
682
	for (int k = 0; k < nSendDOFs; k++)
	  sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])];
      }

683
684
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
    }

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size() * nComponents;
      recvBuffers[i] = new double[nRecvDOFs];

      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
    }

    MPI::Request::Waitall(requestCounter, request);

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size();

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
706
	DOFVector<double> *dofvec = vec.getDOFVector(j);
707
 	for (int k = 0; k < nRecvDOFs; k++)
708
	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
709
710
711
712
713
714
715
      }

      delete [] recvBuffers[i];
    }
    
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];
716
#endif
717
718
  }

719
720
721

  void ParallelDomainBase::checkMeshChange()
  {
722
723
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

724
725
726
727
728
    // === If mesh has not been changed on all ranks, return. ===

    int recvAllValues = 0;
    int sendValue = static_cast<int>(mesh->getChangeIndex() != lastMeshChangeIndex);
    mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
729

730
    if (recvAllValues == 0)
731
732
      return;

733
734
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
735

736
737
738
739
740
741
742
743
744
    clock_t first = clock();
    
    do {
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
      RankToBoundMap allBound = myIntBoundary.boundary;
      for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); 
	   it != otherIntBoundary.boundary.end(); ++it)
	allBound[it->first] = it->second;
745

746
      // === Check the boundaries and adapt mesh if necessary. ===
747

748
749
750
751
752
753
754
755
      bool meshChanged = checkAndAdaptBoundary(allBound);

      // === Check on all ranks if at least one rank's mesh has changed. ===

      int sendValue = static_cast<int>(!meshChanged);
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
    } while (recvAllValues != 0);
756

757
758

#if (DEBUG != 0)
759
    debug::writeMesh(feSpace, -1, "mesh");
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
#endif

    INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", 
		  TIME_USED(first, clock()));

    // === Because the mesh has been changed, update the DOF numbering and mappings. ===
   
    updateLocalGlobalNumbering();
  }

  
  bool ParallelDomainBase::checkAndAdaptBoundary(RankToBoundMap &allBound)
  {
    FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
780
781
782
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
783
	elCode.init(boundIt->rankObj);
784
785
786
787
	sendCodes[it->first].push_back(elCode);
      }
    }

788
    StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm, true);
789
    stdMpi.send(sendCodes);
790
791
792
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
793
    // === Compare received mesh structure codes. ===
794
    
795
796
    bool meshFitTogether = true;

797
798
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
799
800
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
801
      
802
803
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
804
805
806
807
	
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
808
809
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
810
811
812
813
814
815
816
817
	  
	  bool b = fitElementToMeshCode(refineManager, 
					recvCodes[i], 
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
					boundIt->rankObj.elType);  
	  if (b)
	    meshFitTogether = false;
818
	}
819
	
820
	i++;
821
822
823
      }
    }

824
    return meshFitTogether;
825
  }
826
  
827
828
829
830
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
831
    SerUtil::serialize(out, vecSize);
832
833
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
834
      SerUtil::serialize(out, dofIndex);
835
836
837
838
839
840
841
842
843
844
    }
  }


  void ParallelDomainBase::deserialize(std::istream &in, DofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    FUNCNAME("ParallelDomainBase::deserialize()");

    int vecSize = 0;
845
    SerUtil::deserialize(in, vecSize);
846
847
848
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
849
      SerUtil::deserialize(in, dofIndex);
850
851
852
853
854
855
856
857
858
859
860
861

      TEST_EXIT_DBG(dofMap.count(dofIndex) != 0)
	("Dof index could not be deserialized correctly!\n");

      data[i] = dofMap[dofIndex];
    }
  }


  void ParallelDomainBase::serialize(std::ostream &out, RankToDofContainer &data)
  {
    int mapSize = data.size();
862
    SerUtil::serialize(out, mapSize);
863
864
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
865
      SerUtil::serialize(out, rank);
866
867
868
869
870
871
872
873
874
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
875
    SerUtil::deserialize(in, mapSize);
876
877
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
878
      SerUtil::deserialize(in, rank);
879
880
881
882
      deserialize(in, data[rank], dofMap);      
    }
  }

883

884
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
885
886
887
888
889
890
891
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

    TraverseStack stack;
892
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
893
894
895
896
897
898
899
900
    while (elInfo) {
      Element *element = elInfo->getElement();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

      if (partitionData && partitionData->getPartitionStatus() == IN) {
901
	if (partitionData->getLevel() == 0)
902
	  elNum = element->getIndex();
903
	
Thomas Witkowski's avatar
Thomas Witkowski committed
904
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
905
906
907
908
909
910
911
912
913
914
915
916
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

917

918
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
919
  {
920
921
    FUNCNAME("ParallelDomainBase::partitionMesh()");

922
923
924
925
926
927
928
929
930
931
932
933
    if (initialPartitionMesh) {
      initialPartitionMesh = false;
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
      partitioner->partition(&elemWeights, INITIAL);
    } else {
      oldPartitionVec = partitionVec;
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
    }    

    partitioner->fillCoarsePartitionVec(&partitionVec);
  }

934

935
  void ParallelDomainBase::createInteriorBoundaryInfo()
936
  {
937
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
938

939
940
941
942
943
944
945
946
947
948
949
950
951
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
    GeoIndex subObj = CENTER;
    switch (dim) {
    case 2:
      subObj = EDGE;
      break;
    case 3:
      subObj = FACE;
      break;
    default:
      ERROR_EXIT("What is this?\n");
    }     
952
953
954

    // === First, traverse the mesh and search for all elements that have an  ===
    // === boundary with an element within another partition.                 ===
Thomas Witkowski's avatar
Thomas Witkowski committed
955

956
    TraverseStack stack;
957
958
959
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
960
961
962
963
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
964

965
      // Check, if the element is within rank's partition.
966
      if (partitionData->getPartitionStatus() == IN) {
967
	for (int i = 0; i < nNeighbours; i++) {
968
969
970
971
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
972
973
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
974

975
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
Thomas Witkowski's avatar
Thomas Witkowski committed
976

977
978
979
980
981
	    // We have found an element that is in rank's partition, but has a 
	    // neighbour outside of the rank's partition.

	    // === Create information about the boundary between the two elements. ===

982
	    AtomicBoundary bound;	    	    
983
984
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
985
	    bound.rankObj.elType = elInfo->getType();
986
	    bound.rankObj.subObj = subObj;
987
	    bound.rankObj.ithObj = i;
988
989

	    bound.neighObj.el = elInfo->getNeighbour(i);
990
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
991
	    bound.neighObj.elType = -1;
992
	    bound.neighObj.subObj = subObj;
993
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
994
	    
995
996