ParallelDomainBase.cc 87.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 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

Thomas Witkowski's avatar
Thomas Witkowski committed
506
    StdMpi<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.                                                 ===
Thomas Witkowski's avatar
Thomas Witkowski committed
614
    synchVector(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
    KSPDestroy(solver);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640

  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
641
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
642
    StdMpi<std::vector<double> > stdMpi(mpiComm);
643
644

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
645
	 sendIt != sendDofs.end(); ++sendIt) {
646
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
647
648
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
649
      
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
652
653
654
655

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

Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
658
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
659

Thomas Witkowski's avatar
Thomas Witkowski committed
660
    stdMpi.startCommunication<double>(MPI_DOUBLE);
661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
663
664
665
666
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      for (unsigned int i = 0; i < recvIt->second.size(); i++)
	vec[*(recvIt->second)[i]] = stdMpi.getRecvData(recvIt->first)[i];
  }
667
668


Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
  void ParallelDomainBase::synchVector(SystemVector &vec)
  {
    StdMpi<std::vector<double> > stdMpi(mpiComm);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
685
      stdMpi.send(sendIt->first, dofs);
686
687
688
    }

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents);
691

Thomas Witkowski's avatar
Thomas Witkowski committed
692
    stdMpi.startCommunication<double>(MPI_DOUBLE);
693
694

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
695
696
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
697
698

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
 	for (int j = 0; j < nRecvDofs; j++)
	  (*dofvec)[*(recvIt->second)[j]] = 
	    stdMpi.getRecvData(recvIt->first)[counter++];
704
705
706
707
      }
    }
  }

708
709
710

  void ParallelDomainBase::checkMeshChange()
  {
711
712
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

713
714
715
716
717
    // === 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);
718

719
    if (recvAllValues == 0)
720
721
      return;

722
723
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
724

725
726
727
728
729
730
731
732
733
    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;
734

735
      // === Check the boundaries and adapt mesh if necessary. ===
736

737
738
739
740
741
742
743
744
      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);
745

746
747

#if (DEBUG != 0)
748
    debug::writeMesh(feSpace, -1, "mesh");
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
#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) {
769
770
771
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
772
	elCode.init(boundIt->rankObj);
773
774
775
776
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
777
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
778
    stdMpi.send(sendCodes);
779
780
781
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
782
    // === Compare received mesh structure codes. ===
783
    
784
785
    bool meshFitTogether = true;

786
787
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
788
789
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
790
      
791
792
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
793
794
795
796
	
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
797
798
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
799
800
801
802
803
804
805
806
	  
	  bool b = fitElementToMeshCode(refineManager, 
					recvCodes[i], 
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
					boundIt->rankObj.elType);  
	  if (b)
	    meshFitTogether = false;
807
	}
808
	
809
	i++;
810
811
812
      }
    }

813
    return meshFitTogether;
814
  }
815
  
816
817
818
819
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
820
    SerUtil::serialize(out, vecSize);
821
822
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
823
      SerUtil::serialize(out, dofIndex);
824
825
826
827
828
829
830
831
832
833
    }
  }


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

    int vecSize = 0;
834
    SerUtil::deserialize(in, vecSize);
835
836
837
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
838
      SerUtil::deserialize(in, dofIndex);
839
840
841
842
843
844
845
846
847
848
849
850

      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();
851
    SerUtil::serialize(out, mapSize);
852
853
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
854
      SerUtil::serialize(out, rank);
855
856
857
858
859
860
861
862
863
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
864
    SerUtil::deserialize(in, mapSize);
865
866
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
867
      SerUtil::deserialize(in, rank);
868
869
870
871
      deserialize(in, data[rank], dofMap);      
    }
  }

872

873
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
874
875
876
877
878
879
880
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

    TraverseStack stack;
881
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
882
883
884
885
886
887
888
889
    while (elInfo) {
      Element *element = elInfo->getElement();

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

      if (partitionData && partitionData->getPartitionStatus() == IN) {
890
	if (partitionData->getLevel() == 0)
891
	  elNum = element->getIndex();
892
	
Thomas Witkowski's avatar
Thomas Witkowski committed
893
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
894
895
896
897
898
899
900
901
902
903
904
905
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

906

907
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
908
  {
909
910
    FUNCNAME("ParallelDomainBase::partitionMesh()");

911
912
913
914
915
916
917
918
919
920
921
922
    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);
  }

923

924
  void ParallelDomainBase::createInteriorBoundaryInfo()
925
  {
926
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
927

Thomas Witkowski's avatar
Thomas Witkowski committed
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
    // === First, create the interior boundaries based on macro element's  ===
    // === neighbour informations.                                         ===

    createMacroElementInteriorBoundaryInfo();

    // === Second, search the whole mesh for interior boundaries that consists of ===
    // === substructures of the macro elements.                                   ===

    createSubstructureInteriorBoundaryInfo();


    // === Once we have this information, we must care about the order of the atomic ===
    // === bounds in the three boundary handling object. Eventually all the bound-   ===
    // === aries have to be in the same order on both ranks that share the bounday.  ===

    StdMpi<std::vector<AtomicBoundary> > stdMpi(mpiComm);
    stdMpi.send(myIntBoundary.boundary);
    stdMpi.recv(otherIntBoundary.boundary);
    stdMpi.startCommunication<int>(MPI_INT);
   
    // === The information about all neighbouring boundaries has been received. So ===
    // === the rank tests if its own atomic boundaries are in the same order. If   ===
    // === not, the atomic boundaries are swaped to the correct order.             ===

    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {

      // === We have received from rank "rankIt->first" the ordered list of element ===
      // === indices. We now have to sort the corresponding list in this rank to    ===
      // === get the same order.                                                    ===
      
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	// If the expected object is not at place, search for it.

	BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;

	if ((rankIt->second)[j].neighObj != recvedBound) {
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
 	    if ((rankIt->second)[k].neighObj == recvedBound)
	      break;

	  // The element must always be found, because the list is just in another order.
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	    ("Should never happen!\n");

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }
    }

    // === Do the same for the periodic boundaries. ===

    if (periodicBoundary.boundary.size() > 0) {
      stdMpi.clear();

      InteriorBoundary sendBounds, recvBounds;     
      for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	   rankIt != periodicBoundary.boundary.end(); ++rankIt) {

	TEST_EXIT_DBG(rankIt->first != mpiRank)
	  ("It is no possible to have an interior boundary within a rank partition!\n");

	if (rankIt->first < mpiRank)
	  sendBounds.boundary[rankIt->first] = rankIt->second;
	else
	  recvBounds.boundary[rankIt->first] = rankIt->second;
      }

      stdMpi.send(sendBounds.boundary);
      stdMpi.recv(recvBounds.boundary);
      stdMpi.startCommunication<int>(MPI_INT);

      for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	   rankIt != periodicBoundary.boundary.end(); ++rankIt) {
	if (rankIt->first <= mpiRank)
	  continue;
	  
	for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	  
	  BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;
	  
	  if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) {    
	    int k = j + 1;	    
	    for (; k < static_cast<int>(rankIt->second.size()); k++)
	      if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound)
		break;
	    
	    // The element must always be found, because the list is just in 
	    // another order.
	    TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	      ("Should never happen!\n");
	    
	    // Swap the current with the found element.
	    AtomicBoundary tmpBound = (rankIt->second)[k];
	    (rankIt->second)[k] = (rankIt->second)[j];
	    (rankIt->second)[j] = tmpBound;	
	  }  	  
	}
      }     
    } // periodicBoundary.boundary.size() > 0
  }


  void ParallelDomainBase::createMacroElementInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createMacroElementInteriorBoundaryInfo()");

1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
    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");
    }