ParallelDomainBase.cc 77.3 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 new global and local DOF numbering. ===
111

112
    createLocalGlobalNumbering();
113

114
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
117

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
405
406
    MSG("Is rank DOF: %d\n", isRankDof[6717 - rstart * nComponents]);

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

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

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

421
422
423
	    // 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;
424

425
	    if (isRankDof[*cursor]) {
426
427
428
429
430
431

	      // === 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;
432
433

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

449
450
451
452
453
		  // 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]++;
454
		  else
455
		    o_nnz[localPetscRowIdx]++;
456
457
458
		}    
	      }
	    } else {
459
460
461
462
463
	      
	      // === 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.                          ===
464

465
466
	      // Find out who is the member of this DOF.
	      int sendToRank = -1;
467
468
469
470
471
	      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) {
472
473
474
475
476
477

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

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

		if (sendToRank != -1)
		  break;
	      }

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

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

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

506
    // === Send and recv the nnz row structure to/from other ranks. ===
507

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

514
515
    // === Evaluate the nnz structure this rank got from other ranks and add it to ===
    // === the PETSc nnz data structure.                                           ===
516

517
518
519
520
521
522
    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;
523

524
	  int localRowIdx = r - rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
525

526
527
528
	  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);
529
	  
530
	  if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
531
	    o_nnz[localRowIdx]++;
532
	  else
533
	    d_nnz[localRowIdx]++;
534
535
	}
      }
536
    }  
Thomas Witkowski's avatar
Thomas Witkowski committed
537

538
    // === Create PETSc matrix with the computed nnz data structure. ===
539
540

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
541
542
543
544
545
546
547
548
		    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
549

550
    // === Transfer values from DOF matrices to the PETSc matrix. === 
551
552
553

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
554
	if ((*mat)[i][j])
555
556
557
558
559
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

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

560
561
    // === Transfer values from DOF vector to the PETSc vector. === 

562
    for (int i = 0; i < nComponents; i++)
563
      setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
564

565
566
567
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

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


572
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
573
574
575
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

576
577
578
579
580
581
582
583
584
#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

585
    // === Init Petsc solver. ===
586

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

597
598
599

    // === Run Petsc. ===

600
    KSPSolve(solver, petscRhsVec, petscSolVec);
601

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

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

    VecRestoreArray(petscSolVec, &vecPointer);

614
615
616

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

619
620
621

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

622
623
624
625
626
627
628
629
630
631
    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);

632
633
634

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

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

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

663
    stdMpi.startCommunication<int>();
664
665
666
667
#endif

#if 1

668
669
670
671
672
673
674
675
676
    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++) {
677
678
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
679
680
681

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

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

    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++) {
710
	DOFVector<double> *dofvec = vec.getDOFVector(j);
711
 	for (int k = 0; k < nRecvDOFs; k++)
712
	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
713
714
715
716
717
718
719
      }

      delete [] recvBuffers[i];
    }
    
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];
720
#endif
721
722
  }

723
724
725

  void ParallelDomainBase::checkMeshChange()
  {
726
727
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

728
729
730
731
732
    // === 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);
733

734
    if (recvAllValues == 0)
735
736
      return;

737
738
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
739

740
741
742
743
744
745
746
747
748
    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;
749

750
      // === Check the boundaries and adapt mesh if necessary. ===
751

752
753
754
755
756
757
758
759
      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);
760

761
762

#if (DEBUG != 0)
763
    debug::writeMesh(feSpace, -1, "mesh");
764
765
766
767
768
769
770
771
#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();
772
773

    exit(0);
774
775
776
777
778
779
780
781
782
783
784
785
  }

  
  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) {
786
787
788
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
789
	elCode.init(boundIt->rankObj);
790
791
792
793
794
795
796
	sendCodes[it->first].push_back(elCode);
      }
    }

    StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm);
    stdMpi.prepareCommunication(true);
    stdMpi.send(sendCodes);
797
798
799
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
800
    // === Compare received mesh structure codes. ===
801
    
802
803
    bool meshFitTogether = true;

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

831
    return meshFitTogether;
832
  }
833
  
834
835
836
837
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
838
    SerUtil::serialize(out, vecSize);
839
840
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
841
      SerUtil::serialize(out, dofIndex);
842
843
844
845
846
847
848
849
850
851
    }
  }


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

    int vecSize = 0;
852
    SerUtil::deserialize(in, vecSize);
853
854
855
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
856
      SerUtil::deserialize(in, dofIndex);
857
858
859
860
861
862
863
864
865
866
867
868

      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();
869
    SerUtil::serialize(out, mapSize);
870
871
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
872
      SerUtil::serialize(out, rank);
873
874
875
876
877
878
879
880
881
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
882
    SerUtil::deserialize(in, mapSize);
883
884
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
885
      SerUtil::deserialize(in, rank);
886
887
888
889
      deserialize(in, data[rank], dofMap);      
    }
  }

890

891
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
892
893
894
895
896
897
898
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

    TraverseStack stack;
899
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
900
901
902
903
904
905
906
907
    while (elInfo) {
      Element *element = elInfo->getElement();

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

      if (partitionData && partitionData->getPartitionStatus() == IN) {
908
	if (partitionData->getLevel() == 0)
909
	  elNum = element->getIndex();
910
	
Thomas Witkowski's avatar
Thomas Witkowski committed
911
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
912
913
914
915
916
917
918
919
920
921
922
923
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

924

925
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
926
  {
927
928
    FUNCNAME("ParallelDomainBase::partitionMesh()");

929
930
931
932
933
934
935
936
937
938
939
940
    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);
  }

941

942
  void ParallelDomainBase::createInteriorBoundaryInfo()
943
  {
944
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
945

946
947
948
949
950
951
952
953
954
955
956
957
958
    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");
    }     
959
960
961

    // === 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
962

963
    TraverseStack stack;
964
965
966
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
967
968
969
970
971
    while (elInfo) {
      Element *element = elInfo->getElement();

      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
972

973
      // Check, if the element is within rank's partition.
974
      if (partitionData->getPartitionStatus() == IN) {
975
	for (int i = 0; i < nNeighbours; i++) {
976
977
978
979
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
980
981
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
982

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

985
986
987
988
989
	    // 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. ===

990
	    AtomicBoundary bound;	    	    
991
992
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
993
	    bound.rankObj.elType = elInfo->getType();
994
	    bound.rankObj.subObj = subObj;
995
	    bound.rankObj.ithObj = i;
996
997

	    bound.neighObj.el = elInfo->getNeighbour(i);
998
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
999
	    bound.neighObj.elType = -1;